Dr Vermillard Patrick

Sabine

Introduction 1.1 Context

Rainfall is a primary input for numerous applications in hydrometeorology such as flood warning, water resource management, weather forecasting and emergency planning systems.

In urban areas, these applications operate at small scales, therefore rainfall information with a good accuracy at high spatial and temporal resolutions is required. So far, traditional sensors such as rain gauge and weather radar are used for measuring rainfall. Rain gauges are yet the most widely used device for providing an accurate point measurement at a ground level while weather radars are capable of providing a detailed rainfall information over a large area [START_REF] Berne | Radar for hydrology: Unfulfilled promise or unrecognized potential?[END_REF]. However, a lack of such sensors is a major problem in many regions of the globe. Nowadays, most cities worldwide are well-equipped with commercial microwave links that are mainly operated by cellular communication companies. Very recent studies conducted by [START_REF] Messer | Environmental monitoring by wireless communication networks[END_REF]Leijnse et al., 2007b) suggest that such links could become a promising sensor for measuring space-time rainfall at fine-scale. The measurement principle is based on the phenomenon that a transmitted signal along a microwave link is attenuated by rainfall [START_REF] Atlas | Path-and area-integrated rainfall measurement by microwave attenuation in the 1-3 cm band[END_REF]. In telecommunication engineering, this phenomenon is well-understood, and a relationship between rainfall intensity and the microwave signal attenuation has been es-tablished [START_REF] Olsen | The aR b relation in the calculation of rain attenuation[END_REF][START_REF] Crane | Prediction of attenuation by rain[END_REF]. Using such relationship, averaged rainfall intensity along a single microwave link path can be obtained from the signal attenuation data. The density of commercial microwave links is usually high in cities; therefore, this fact raises a new question about the feasibility of retrieving rainfall maps over urban areas using signal attenuation data from cellular communication companies.

The use of commercial microwave links for rainfall measurement could bring two benefits: i) an alternative to measure rainfall in locations where no rainfall information is available, and ii) a complement to traditional sensors, in particular weather radar. In spite of scientific progress on this field, conclusions drawn from the literature review (chapter 3 of this thesis) recognise three ongoing issues that are still subject to further investigations: (i) Understanding and quantifying possible measurement errors sources that might occur along a single or multiple links;

(ii) Combining a single or multiple links with traditional sensors in order to improve the accuracy of measured rainfall;

(iii) The development of spatial rainfall retrieval algorithms that can be used to convert signal attenuation data from an existing network of commercial microwave links into rainfall maps.

Among those issues, the third category is the central focus of the present study. We believe that this research study is one of the first attempts so far to thoroughly examine the capability of commercial microwave links employed by cellular communication networks for rainfall retrieval at the urban scale. It is true that cellular communication networks are designed to provide a better communication quality. A geometry of the network topology is complex and arbitrary, their location and operating frequencies are inhomogeneous. Such pattern of the network geometry is almost unchanged and it requires a special processing algorithm that is able to retrieve rainfall map using the recorded signal attenuation data from cellular network.

The importance and originality of this study is in the development of rainfall retrieval algorithms. Up to now, a limited number of studies have investigated the relative importance of rainfall retrieval algorithms that can be applied to real existing complex geometry of network topology. We develop, examine as well as compare two novel retrieval algorithms, namely inverse and tomographic (introduced in chapter 5 and chapter 6 of this study, respectively) that can be used in the application of cellular communication networks for rainfall monitoring over urban areas.

PURPOSE AND SCOPE OF THE PRESENT STUDY

Purpose and scope of the present study

The main purpose of this study is to assess the feasibility of rainfall retrieval over urban areas using signal attenuation data from commercial microwave links. We adopted a simulation framework applied to a realistic case study. The simulation framework allows to understand different aspects of the main purpose of this study such as dealing with errors and uncertainties in the measurement of signal attenuation data which is the key challenge in experimental frameworks that use real microwave signal data from cellular communication companies. This thesis is based on a network of commercial microwave links operated by cellular communication companies in Nantes city and uses a set of rainfall images provided by Météo France.

In this simulation framework, we address two major challenges to achieve the main purpose of this study:

-To generate rain attenuation data along microwave links. In this step, we generate path-integrated attenuation data in realistic conditions. This helps to substitute the generated rain attenuation data for real one that can be obtained by telecommunication operators.

-To retrieve rainfall map using the generated rain attenuation data. The aim of this step is to retrieve rainfall maps based on the generated rain attenuation data.

The reader should bear in mind that the present study is based on a simulation framework which is applied to a realistic case study. This study does not use real signal attenuation data recorded at microwave antenna stations of cellular companies because commercial data were not available. Therefore, a full discussion of investigating error sources along microwave links lies beyond the scope of this study although we did take into account this aspect based on the reviewed literature. In addition, the main interest in this research is limited only to rainfall measurement at the city scale.

What is microwave link?

A microwave link is the main core part of wireless communication systems such as satellites, cellular phone networks, broadcasting, military radars. The microwave link uses an electromagnetic signal in the frequency range that carries information for short (few meters) or long (several km) distances between two fixed locations [START_REF] Freeman | Radio system design for telecommunication[END_REF]. Operating frequencies are .

roughly between 5 and 50 GHz, selected by telecommunication engineers. Typically, the microwave link consists of transmitter, receiver antennas and transmission line Figure 1.1. Aside from many types (horn, slot, dipole or dielectric), parabolic antennas are very common in the microwave communication system. The microwave antennas are usually mounted on the top of buildings, towers so that the signal can be safely transmitted in a clear path.

Microwave signals can be attenuated due to various factors such as wave reflection on the ground, effects of buildings and trees and other environmental issues. Apart from its installation in a clear path, the operating frequency and the length of transmission line along these antennas are often chosen on the basis of climate conditions for a geographic location [START_REF] Freeman | Radio system design for telecommunication[END_REF].

The reason is that the transmitted signal power can also be attenuated by different atmospheric variables such gas, water vapour, hail, fog, snow and rain. Among them, rain is the most influential factor, especially at frequencies above 10 GHz [START_REF] Oguchi | Electromagnetic wave propagation and scattering in rain and other hydrometeors[END_REF][START_REF] Freeman | Radio system design for telecommunication[END_REF]; Gurung and Zhao, 2011).

Thesis Outline

The structure of this thesis consists of 7 chapters:

Chapter 2 provides the physical background on rain induced attenuation phenomenon that takes place during the signal transmission along microwave links. In this context, this chapter highlights the key theoretical concepts regarding the context of this research which will serve as a fundamental basis throughout entire thesis. These concepts are related to electromagnetic signal scattering and micro-physical structure of rainfall. Further, a formulation of specific rain attenuation as a function of signal frequency and rain intensity is also discussed.

Chapter 3 contextualises the present research by providing a detailed review on the state of the art knowledge in microwave link based rainfall measurement. This chapter is subdivided into three major sections. The first section discusses recent advances in the application of microwave links to rainfall measurement as well as their integration with weather radar and rain gauges. The second section reviews ongoing issues regarding measurement error sources along microwave links. The objective of the third section is to give an overall summary of existing retrieval algorithms can be used to retrieve rainfall based on signal attenuation data from microwave links.

Chapter 4 describes the details of our case study used to generate signal attenuation data due to rainfall. This chapter consists of three major sections. The first two sections present data sets that include a description of microwave links of cellular communication network and rainfall images obtained from weather radar. The last section focuses on a methodology for generating rain attenuation data along the microwave links. Error sources regarding environmental and hardware-equipments of microwave antenna stations are simulated in order to mimic the real nature of the signal attenuation data. The attenuation data generated in this chapter will be used as real signal attenuation data that can be obtained from cellular company operators.

Chapters 5 presents the first retrieval algorithm that will be applied to convert signal attenuation data into rainfall map. This chapter utilises microwave signal data generated in Chapter 4 in order to retrieve rainfall map using cellular network. First, the retrieval algorithm and its application conditions are presented. The presented algorithm here employs the principle of inverse problems in which a definition of apriori knowledge is of utmost importance. Second, the sensitivity analysis test to various parameters of the algorithm is performed including conditions such as the presence of different magnitudes of measurement and model errors. Retrieval performance of the algorithm depending on the network density is analysed. Capabilities as well as limitations of the proposed algorithm in capturing spatial variability of rainfall are also discussed.

Chapter 6 presents the second retrieval algorithm which uses the principle of tomography and in our case, discrete tomography. The objective is the same as stated in Chapter 5, that is to

CHAPTER 1. INTRODUCTION
reconstruct rainfall map on the basis of microwave rain attenuation data generated in Chapter 4. The present chapter is grouped by three major sections. In the first section, a core principle of tomography and a mathematical background of the tomographic algorithm used in this chapter are defined. In the context of tomography problems, arbitrary geometry of network topology, inhomogeneity of link frequencies and lengths are the key challenges to be addressed in the application conditions of the algorithm. Therefore, in the second section, specific procedures to adjust the parameters of the algorithm are presented. The last major section deals with retrieval tests carried out over different regions of the network system. In addition to this, retrieval performance of this algorithm is compared with the one to be presented in Chapter 5.

Chapter 7 summarises the main findings of this project and discusses the implications of the findings of this study to future research.

2 Physical Background

Introduction

This chapter presents a physical basis of rain induced attenuation phenomenon that takes place during the signal transmission along the microwave links. Since rain attenuation is an intrinsic characteristic of microwave link signal it is important to understand the relation between micro-physics of rainfall and attenuation of electromagnetic signal in the atmosphere. Rainfall is not the only source that attenuates the signal strength, atmospheric gases also have an impact on the microwave signal. However, attenuation caused by rain is the most influential factor among them, especially above 10 GHz frequency band employed by microwave links [START_REF] Oguchi | Electromagnetic wave propagation and scattering in rain and other hydrometeors[END_REF][START_REF] Freeman | Radio system design for telecommunication[END_REF][START_REF] Gurung | Attenuation of Microwave signal and its Impacts on Communication system[END_REF]. The magnitude of the attenuation increases from lower to higher rain intensity which indicates that there is a strong relation between these variables. The issue here is how to quantify that attenuation caused by rainfall. Therefore, the main objective of the present chapter is to give a theoretical background that explains the relation between rainfall intensity and microwave signal attenuation.

This chapter is organised as follows. First, section 2.2 gives a detailed description about the electromagnetic scattering phenomenon with an emphasis on rainfall. Second, section 2.3 focuses on the description of rainfall micro-physics, in particular rain drop size distribution characteristics. Third, section 2.6 presents the attenuation model that is used to quantify a CHAPTER 2. PHYSICAL BACKGROUND specific rain attenuation. Finally, section 2.7 summarizes the chapter.

Electromagnetic signal scattering

Electromagnetic signal scattering is explained as the redirection of electromagnetic waves that takes place when they encounter an obstacle or an inhomogeneity, scattering particle [START_REF] Hahn | Light scattering theory[END_REF]; Bohren and Huffman, 2008) as depicted in Figure 2.1. In the context of our study, this obstacle is rain drop which causes the part of the microwave signal to be scattered and partially absorbed as well. Since this takes place among all drops which are interacting with the incident signal, it results in the signal attenuation in the propagation direction of the incident signal. [START_REF] Hahn | Light scattering theory[END_REF] There are different approaches such as T-matrix [START_REF] Waterman | Matrix formulation of electromagnetic scattering[END_REF][START_REF] Mishchenko | T-matrix computations of light scattering by nonspherical particles: a review[END_REF][START_REF] Mishchenko | Calculation of the amplitude matrix for a nonspherical particle in a fixed orientation[END_REF], Mie method [START_REF] Mie | Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen[END_REF] to compute the electromagnetic signal scattering and absorption by rain drops. For example, T-matrix is the complex method which takes into account the fact that the shape of the rain drop is oblate. We will explain the rain drop characteristics later in section 2.3. The T-matrix method is considered to be computationally expensive. If the drop size is almost equal or smaller than the electromagnetic wavelength, 'Mie' scattering can be applied that encompasses the general spherical scattering solution (absorbing or non-absorbing) without a particular bound on a particle size above 10 GHz [START_REF] Mie | Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen[END_REF]. This is the scattering of electromagnetic radiation by primarily spherical particles whose diameters are comparable to the wavelength of the incident radiation. Since the minimum value at the operating frequencies of commercial microwave links is 5 GHz or greater, Mie scattering becomes important for calculating the electromagnetic signal attenuation due to rain. The ap-plication of both approaches may depend on the climate, temperature as well as the shape of the rain drop. In the next paragraph, we will provide some details about the Mie scattering theory that will help to understand the physical background behind the phenomenon of electromagnetic signal attenuation due to rainfall.

Mie scattering parameters and formulations

Scattering amplitude function is obtained from the solution of the boundary value problem at the surface of a raindrop. Considering two scattering amplitudes, i.e. S 1 and S 2 , the intensity and the state of the polarization of the scattered radiation in any direction by angle θ are given by the following equations: and Hankel function of the first kind evaluated at size parameter α (or sometimes referred as x and complex refractive index, m × α, [START_REF] Hulst | Light scattering by small particles[END_REF][START_REF] Bohren | Absorption and scattering by a sphere[END_REF]). The notation α is the size parameter of rain drop and is expressed as follows:

S 1 (θ) = ∞ n=1 2n + 1 n(n + 1) [a n π n (cos θ) + b n τ n (cos θ)] (2.1) S 2 (θ) = ∞ n=1 (2n + 1) n(n + 1) [b n π n (cos θ) + a n τ n (cos θ)] (2.2)
x = α = kD/2 = ( π λ )D = circumf erence of sphere wavelength (2.3)
Angular coefficients τ n and π n are the functions of cos θ. They are defined in terms of Legendre polynomials and their derivatives are expressed as followins:

π n (cos θ) = 1 sin θ P 1 n (cos θ) and τ n (cos θ) = d dθ P 1 n (cos θ) (2.4)
Where, P 1 n are associated Legendre functions. It is assumed that the complex forward scattering angle direction (θ) is the same as the incident electromagnetic wave direction; therefore, it is equal to zero. Using this assumption we obtain S 1 (0 o ) = S 2 (0 o ) at θ=0 o . As a result, the scattering amplitudes in the equations (2.1, 2.2) take the following form:

S(θ = 0 o ) = 1 2 N max n=1 (2n + 1)(a n + b n ) (2.5)
Where, N max -maximum size parameter; N max ≈ α + α 1/3 + 2; α -size parameter (or x).

Mie coefficients a n and b n are defined by the following expression:

a n = ψ n (α)ψ n (mα) -mψ n (mα)ψ n (α) ζ(α)ψ n (mα) -mψ n (mα)ζ n (α) (2.6) b n = mψ n (α)ψ n (mα) -ψ n (mα)ψ n (α) mζ(α)ψ n (mα) -mψ n (mα)ζ n (α) (2.7)
Where, m is the complex refractive index of the spherical rain drop and computed using Debye formula [START_REF] Ray | Broadband complex refractive indices of ice and water[END_REF] which depends on wavelength λ and temperature t.

The complex refractive index is given by m = m 1 -m 2 i, where m 1 and m 2 are the real and imaginary parts of the index of refraction, respectively. The notations (ψ and ζ) are Ricatti-Bessel functions, defined in terms of the half-integer-order Bessel function of the first kind

(J n+1/2 (z)): Ψ n (z) = ( πz 2 ) 1/2 J n+1/2 (z) (2.8) ζ n (z) = ( πz 2 ) 1/2 H n+1/2 (z) (2.9) = Ψ n (z) + iX n (z) (2.10)
Where, H n+1/2 (z) is the half-integer-order Hankel function of the second kind, where the parameter X n is defined in terms of the half-integer-order Bessel function of the second kind,

Y n+1/2 (z) X n (z) = -( πz 2 ) 1/2 Y n+1/2 (z) (2.11)
Where, z = α or mα.

Using the expression (2.5), we are able to compute the forward scattering amplitude as a function of rain drop radius. As can be seen in Figure 2.4 that the imaginary as well as real part of the scattering functions are increasing with respect to an increase of the drop radius size. The imaginary part is directly proportional to the signal scattering. This clearly represents the electromagnetic signal scattering dependence on a drop size.

Mie Efficiency factor and Cross Sections

The Mie efficiency factors are derived from the scattering amplitudes. The scattering efficiency Q sca follows from the integration of the scattered power over all directions (see the arrows of scattering wave in Figure 2.1). The extinction efficiency Q ext follows from the extinction theorem [START_REF] Hulst | Light scattering by small particles[END_REF][START_REF] Hahn | Light scattering theory[END_REF], also called forward-scattering theorem, leading to extinction efficiency:

Q ext = 2 x 2 N max n=1 (2n + 1)Re(a n + b n ) (2.12) = 4π k 2 Re{S(0)} (2.13) Q sca = 2 x 2 N max n=1 (2n + 1)Re(|a n | 2 + |b n | 2 ) (2.14)
Where, Q ext and Q sca stands for extinction and scattering efficiencies. The notation Re{S(0)} stands for real part of the complex forward scattering function.

The extinction cross section, denoted as σ ext , is defined by the following expression:

σ ext = λ 2 π Re{S(0)} (2.15)
The expression (2.15) is considered as a function of drop size diameter, wavelength and complex refractive index of water content. Interdependence between efficiency and cross section factors are as follows:

Q ext = Q sca + Q abs (2.16)
and

σ ext = σ sca + σ abs (2.17)
The formulas (2.16) and (2.17) define the amount of extinct wave coming out of the drop after the scattering and absorption effects.

Mie scattering computation

The algorithm and numerical methods for computing Mie scattering can be found in various literatures [START_REF] Shah | Numerical Methods for Mie Theory of Scattering by a Sphere[END_REF][START_REF] Wiscombe | Improved Mie scattering algorithms[END_REF][START_REF] Hajny | Computation OF Raindrops Scattering[END_REF][START_REF] Du | Mie-scattering calculation[END_REF][START_REF] Gogoi | Mie scattering computation of spherical particles with very large size parameters using an improved program with variable speed and accuracy[END_REF].

For the simplicity, we use the following steps to perform Mie scattering computation:

Step 1. Initial values for the calculation as input parameters: rain drop temperature (t = 0 o C), frequencies at 18, 23, 38 GHz;

Step 2. Compute refractive index of water m at a given temperature from Debye formula;

Step 3. Compute a n and b n for n = 1...N max , where N max ≈ α + α 1/3 + 2, from size parameter x and index of refraction m (uses recursion relations for the spherical Bessel functions), the equations (2.6, 2.7);

Step 4. Compute Extinction Q ext , scattering Q sca efficiencies given in (2.12, 2.14).

Then, we compute extinction, scattering and absorption cross sections (2.16, 2.17);

Step 5. Optionally, compute S 1 (θ) and S 2 (θ) at desired scattering angles from a n and b n , π n and τ n from recursion.

The algorithm could have been quite simple if the Bessel and Legendre polynomials were replaced by simple complex goniometric functions. The infinite series can be limited to the n being about 10 (or even less) of a sufficient accuracy. There are available program codes to perform the algorithm in several programming languages such as Java, C, Fortran, Matlab, Octave, Pascal [START_REF] Shah | Numerical Methods for Mie Theory of Scattering by a Sphere[END_REF][START_REF] Mätzler | MATLAB functions for Mie scattering and absorption, version 2[END_REF][START_REF] Mishchenko | Gustav Mie and the evolving discipline of electromagnetic scattering by particles[END_REF][START_REF] Wriedt | Mie theory 1908, on the mobile phone 2008[END_REF][START_REF] Gogoi | Mie scattering computation of spherical particles with very large size parameters using an improved program with variable speed and accuracy[END_REF]. All of those programs differ from each other in terms of speed and complexity.

We use the codes written in Python 2.7 considered to be a high-level language and user-friendly that has standard scientific library 2 to compute functions such as Bessel and Hankel.

Comparison of Mie parameters: Scattering, Absorption, Extinction

To test the impact of drop size on the electromagnetic signal amount we make a comparison between Mie parameters defined in the previous section. This helps understand the frequency dependence on the amount of scattering, absorption and extinction as a function of drop size.

To carry out this test, we use 18, 23 and 38 GHz which are the frequencies of the microwave links chosen in our study. The conditions for this comparison are as follows:

1. Frequency: 18, 23, 38 GHz; the extinction is gradually increasing, especially at 38 GHz. This results in a sharper increase compared to lower frequency, i.e. 18 GHz. Also, we can see that rain drop radius smaller than 0.2 cm represents a small change in the cross sections at all used frequencies.

CHAPTER 2. PHYSICAL BACKGROUND

The extinction cross section rises sharply when the drop radius is greater than 0.2 cm at all frequencies whereas the absorption cross section increases less slowly.

Rain drop characteristics 2.3.1 Shape and size of the rain drop

Rain drop, sometimes called as a water droplet, is formed by condensation of water vapour in a cloud, that is heavy enough to fall from the cloud and large enough to reach the surface of land or sea before evaporating in the unsaturated air beneath the cloud. The shape of a drop is more spherical and horizontally oblate due to the force of air in vertical direction as shown in

Drop size distribution

One of the main factors of the specific attenuation is considered as Drop Size Distribution (DSD). Its analytical formulations are mainly used for describing the rain drop size concentration in the atmosphere. The main reason why there is a need for computing DSD is that it helps to classify precipitation types. There are various types of distrometers such as electromechanical, optical, video to measure DSD [START_REF] Best | Drop-size distribution in cloud and fog[END_REF][START_REF] Mason | Drop-size distributions from various types of rain[END_REF][START_REF] Fišer | First results of DSD measurement by videodistrometer in the Czech Republic in 1998-1999[END_REF][START_REF] Tokay | Measurements of drop size distribution in the southwestern Amazon basin[END_REF][START_REF] Brandes | Drop size distribution retrieval with polarimetric radar: Model and application[END_REF]. The formulation of drop size distribution and its modelling were investigated in [START_REF] Marshall | The distribution of raindrops with size[END_REF][START_REF] Best | Drop-size distribution in cloud and fog[END_REF][START_REF] Joss | The variation of rain drop-size distributions at Locarno, paper presented at[END_REF][START_REF] Waldvogel | The N 0 jump of raindrop spectra[END_REF][START_REF] Torres | A general formulation for raindrop size distribution[END_REF][START_REF] Tokay | Measurements of drop size distribution in the southwestern Amazon basin[END_REF][START_REF] Bringi | Raindrop size distribution in different climatic regimes from disdrometer and dual-polarized radar analysis[END_REF][START_REF] Tapiador | A probabilistic view on raindrop size distribution modeling: a physical interpretation of rain microphysics[END_REF].

The

Gamma distribution

Due to its suitability, gamma distribution is often used to characterize drop size distribution. Based on the real recorded data by distrometer, [START_REF] Mallet | Estimation of gamma raindrop size distribution parameters: Statistical fluctuations and estimation errors[END_REF] confirm that more than 90 % of the drop size distributions follow gamma distribution. There are different types of gamma distribution functions such as gamma and modified gamma. [START_REF] Ulbrich | Natural variations in the analytical form of the raindrop size distribution[END_REF] suggested that gamma function using 3 parameters is able to describe most DSD models and each parameter can simply be computed from estimated moments. The general expression for that is as follows:

N (D) = N 0 D µ exp[-ΛD] (2.18)
Where, Λ -slope of drops; µ -shape of drops; N 0 -number of drops;

Gamma DSD is the general form of exponential form. If rain drop sizes are mainly average, the parameters µ, N 0 and Λ show large values. Theses parameters differ from one DSD type to another, see Table 2.1 and Figure 2.7.

In addition, there are some other models for the drop size spectra such as lognormal, mod- ified gamma, exponential [START_REF] Joss | The variation of rain drop-size distributions at Locarno, paper presented at[END_REF][START_REF] Willis | Functional fits to some observed drop size distributions and parameterization of rain[END_REF]. Overall, [START_REF] Torres | A general formulation for raindrop size distribution[END_REF] investigated a general formulation for describing such models using scaling law parameters in a wide variety of rainy conditions.

Marshal Palmer distribution

The classical study conducted by Marshall-Palmer (MP) is very well-known and often used for computing DSD [START_REF] Marshall | The distribution of raindrops with size[END_REF]. The study provided an exponential formulation for N(D) according to the analysis of two datasets measured with the filter paper method.

In fact, MP DSD is a special case of Gamma DSD where the parameter µ is equal to zero. MP model is very suitable for widespread rain type in continental temperate climate. Therefore, MP size distribution is widely used. However, it shows some inadequacies in expressing other observed spectra and overestimates the number of both the smallest drops [START_REF] Waldvogel | The N 0 jump of raindrop spectra[END_REF]. 2.2).

MP DSD model is expressed by the exponential form of DSD as follows:

N (D, R) = N 0 exp[-Λ(R)D] (2.19)
Where, D -rain drop diameter, mm; N 0 = 8000, mm -4 ; Λ is supposed to have a power-law dependence on the rain rate R: The main reason for the use of these two models i.e. Gamma and MP, is their suitability in temperate climate. Above mentioned studies confirm that Gamma DSD is mainly used for the evaluation of DSD function and provides more adaptable result to express the distribution of the drop sizes in many precipitation types.

Λ(R) = 4.1R -0.
R = 6 * π * 10 -4 ∞ 0 D 3 v(D)N (D)dD (2.20)
Where, v(D) is the fall velocity of raindrops (often called terminal velocity) expressed as a function of drop size in m.sec -1 . As an example, the terminal velocity function can be derived from [START_REF] Atlas | Path-and area-integrated rainfall measurement by microwave attenuation in the 1-3 cm band[END_REF]:

v(D) = cD γ (2.21)
Where, the coefficients c and γ are 3.78 and 0.67, respectively.

Besides the fact that the expression (2.21) has been shown to be a good fit for a wide range of rain drop sizes it also makes computation less complicated. More comprehensive and complex relation for describing terminal velocity can be found in [START_REF] Beard | Terminal velocity and shape of cloud and precipitation drops aloft[END_REF]).

Scaling law and self-consistency relationship

In hydrometeorological studies, most DSD models used by now such as the exponential distribution [START_REF] Marshall | The distribution of raindrops with size[END_REF], Weibull distribution [START_REF] Best | Drop-size distribution in cloud and fog[END_REF], gamma distribution [START_REF] Ulbrich | Natural variations in the analytical form of the raindrop size distribution[END_REF][START_REF] Willis | Functional fits to some observed drop size distributions and parameterization of rain[END_REF] are particular cases of a general expression (2.22) where the number of drops per unit of air volume in the size range D to D + dD, N (D, Ψ), depends on D and on the reference variable Ψ:

N (D, Ψ) = Ψ α Ψ * g( D R β Ψ ) (2.22) With x = D R β Ψ we can get g(x) = k * exp(-λx).
In this general expression Ψ can be any integral rainfall variable although R has generally been used [START_REF] Torres | A general formulation for raindrop size distribution[END_REF]. For given Ψ, α Ψ and β Ψ are constants (they do not have any functional dependence on Ψ and g(x) is a function that is independent of the value of Ψ and that will be called the general distribution function.

In fact, the expression (2.22) is a scaling law. In our context, the scaling law is theoretical approach valid over a wide range of scales (i.e. over a wide range of intensities or liquid water content). The described two DSD models in subsection 2.3.2 can be expressed by (2.22).

An important requirement of sets of power law relationships between rainfall related variables is that they should be consistent. Consistency means that power law relationships between variables should satisfy the definitions of these variables in terms of the parameters of the raindrop size distribution. This so-called self-consistency requirement has been considered explicitly by [START_REF] Bennett | The relationship between N0 and λ for Marshall-Palmer type raindrop-size distributions[END_REF]. The resulting constraints on the coefficients of power law relationships between rainfall related variables were treated recently in much more general fashion by [START_REF] Torres | A general formulation for raindrop size distribution[END_REF], as a part of their general formulation for the raindrop size distribution. Substituting (2.22), (2.21) into the definition of R in terms of the raindrop size distribution equation (2.20) leads to the self -consistency constraints:

6π * 10 -4 c ∞ 0 x 3+γ g(x)dx = 1 (2.23) a + (4 + γ)β = 1 (2.24)
Hence, g(x) must satisfy an integral equation (which reduces its degrees of freedom by one), and there is only one free scaling exponent. These self-consistency constraints guarantee that substitution of the parametrization for the raindrop size distribution (2.22) into the defining expression for rainfall intensity (2.20), R, leads to R = R.

Specific rain attenuation

The formulation of specific attenuation due to rainfall, denoted as k in dB per km, is described as follows:

k = log 10 e * Dmax D min σ ext (D)N (D, R)dD (2.25)
Where, σ ext (D) -extinction cross section, mm 2 , defined in the expression (2.15) of section 2.2; The frequency range in x axis is plotted against the specific rain attenuation in y axis.

N (D,
These models parameters are under the conditions that temperature is 0 o C and rain drop shape is spherical. It is not surprising to see that the specific attenuation in both cases monotonously increases with the correspondence of increase in frequency values. This clearly indicates that the signal attenuation is directly proportional to frequency. Similarly, higher rainfall intensities are also the cause for higher rain attenuation since a rain drop size increases at higher rain rates (e.g., 50 mm.hour -1 ) and is able to absorb and scatter significant amount of signal.

Specific attenuation as a function of rain intensity

Here, the specific attenuation is computed as a function of the rain intensity. However, it is more sensitive to the frequency value. As the rain becomes more intense (higher rate), the specific attenuation based on Gamma DSD rises faster than that of MP DSD at 38

GHz. At frequencies of 18 and 23 GHz, the same trend represents the opposite scenario. The explanation is that MP DSD is suitable only for widespread rain type which mostly consists of lower intensity values. 2.2 and the one demonstrated by [START_REF] Zhang | Project IEEE 802.16 Broadband Wireless Access Working Group[END_REF], respectively.

Summary

In this chapter, the theoretical background of rain attenuation along the microwave links has been discussed. In particular, electromagnetic signal scattering, absorption dependence on drop sizes taking into account different rain intensities have been presented. Next, micro physical structure of rain drops and different models which describe the distribution of rain drop sizes have been discussed. Moreover, we have shown the rain attenuation model to be a function of two variables: (i) frequency and (ii) rainfall rate. Further, we compared the specific rain attenuation obtained on the basis of different drop size distribution models. Overall, the present chapter will be considered as a fundamental basis for the electromagnetic signal interaction with rain drops in the atmosphere.

State of the art: Rainfall Measurement

Using Microwave Links

Introduction

Microwave links use electromagnetic radio waves which operate roughly between 5 and 50 GHz. At these frequencies, the signal travelling through the atmosphere gets attenuated by different forms of hydrometeors such as hail, fog, wind, snow and rain. As it has been discussed in chapter 2 rain causes the most significant attenuation, particularly at frequencies bands above 10 GHz used in wireless communication systems [START_REF] Freeman | Radio system design for telecommunication[END_REF]. We have seen that this phenomenon occurs because the wavelength of microwave signal is comparable with a diameter of a rain drop at these frequencies, see section 2.2. Thus, the higher rain intensity is directly responsible for signal attenuation increase. Once such a relationship has been established, the average rain intensity along the link can be estimated [START_REF] Olsen | The aR b relation in the calculation of rain attenuation[END_REF].

Today, microwave links are very high in density and already cover large parts of most urban areas. Therefore, the introduction of the microwave links in cities can be a promising and potentially valuable way of rainfall monitoring. Rainfall estimation using commercial microwave links has become an active research subject. In this context, the objective of this chapter is to provide state-of-the-art knowledge in recent advances, limitations and advantages of the mi-crowave link, and the algorithms and methods developed for rainfall monitoring.

This chapter is organised as follows. In section 3.2, applications and recent advances in microwave links approach are presented. Then, in section 3.3 we discuss the issues regarding existing challenges which are related to error sources and uncertainties in applications of the microwave links. In section 3.4, our focus is to give an overall summary of existing algorithms used to retrieve rainfall based on attenuation data from microwave links. Finally, our conclusions regarding the state-of-the-art knowledge in microwave link based rainfall measurement are presented in section 3.5. In the context of the applicability of the microwave links for rainfall measurement, a number of studies using different methodologies in both theoretical and experimental framework were carried out. These studies can be classified into two groups (i) research and (ii) operational microwave link based investigations.

Using research dedicated links, [START_REF] Holt | Measurement of rainfall by dual-wavelength microwave attenuation[END_REF][START_REF] Holt | Comparison of the use of dualfrequency and single-frequency attenuation for the measurement of path-averaged rainfall along a microwave link[END_REF] studied the capabilities of dual-frequency links for rainfall measurement. The proposed idea was to transform a difference of the signal attenuation measured by the dual-frequency link into a path-averaged rain intensity. Here, a path-averaged rain intensity is referred to rain rate which is assumed to be an average value of rain rate along the microwave link path. Similarly, (Rahimi et al., 2003) con-

RECENT ADVANCES IN MICROWAVE LINK BASED RAINFALL MEASUREMENT41

ducted an experimental study in the north-west of England. To show how the dual-frequency link can be used to measure path-averaged rainfall rate, the authors proposed step-by-step procedures for extracting rain attenuation data from measured signal level before using the actual signal itself for estimating rainfall intensity. The general performance of the approach was tested over a long hour period in two cases (112 and 52 events). The results show a good agreement between the link estimates and the one obtained from rain gauges as well as C-band weather radar both of which were assumed to be the 'ground-truth'. It is important to note that the link estimate is different from the rain gauge based one which is a point measurement at a ground level. The first 'ground-truth' used in this study came from a network of 22 tipping-bucket rain gauges. In order to be able to compare with the link-based estimates, true path-averaged rain intensity was obtained using the weighted average rain rate records of nearest rain gauges along the link. Similarly, comparisons with weather radar, which is considered to be the second ground truth, have been made at the resolution of 2 × 2 km 2 by weighting the link length with radar grid (Cartesian). Later, the suitability of the dual-frequency link was also tested in urban rainfall measurement [START_REF] Rahimi | Dual-frequency linksa complement to gauges and radar for the measurement of rain[END_REF][START_REF] Upton | Microwave links: The future for urban rainfall measurement?[END_REF]. [START_REF] Minda | High temporal resolution path-average rain gauge with 50-GHz band microwave[END_REF] employed horizontally polarized short link (820 meters) at 50 GHz for measuring the path averaged rainfall intensity. It was shown that the link approach can be used as 'path-averaged-raingauge' sensor. Time series data sets recorded by a rain gauge and a disdrometer were used to validate the system performance. According to [START_REF] Fenicia | Microwave links for rainfall estimation in an urban environment: Insights from an experimental setup in Luxembourg-City[END_REF] the dual-frequency-based estimate does not seem to be accurate as compared to one estimated by a single-frequency link.

They came to such conclusion based on the experimental set-up, located in Luxembourg city.

The data sets comprise two horizontally polarized dual-frequency links with different lengths, and 13 rain gauges closely placed along those two links. The data collected during 1.5 year by those rain gauges was used to test the performance of the system in space and time. Here, it is worth mentioning that the characteristics of those link pairs chosen in this study are based on the recommendations by [START_REF] Rahimi | Use of dual-frequency microwave links for measuring path-averaged rainfall[END_REF] where the relationship between rain intensity and attenuation difference by the link is expected to be linear and not influenced by drop size distribution variability. Besides, the authors give emphasis on the fact that the performance of the link approach is limited due to the uncertainty in rain gauge data. In addition to those investigations, (Leijnse et al., 2007b) proposed that a single-frequency link at 27 GHz is also suitable for rainfall estimation. In fact, the 'suitability' range of the link depends on the linearity of the relationship between microwave signal attenuation and rain intensity [START_REF] Atlas | Path-and area-integrated rainfall measurement by microwave attenuation in the 1-3 cm band[END_REF][START_REF] Olsen | The aR b relation in the calculation of rain attenuation[END_REF]. The authors confirm the linearity of such relationship at 27 GHz using large data set during long period (more than a year) was used to establish that relationship under the assumption of spherical rain drop in the Netherlands climate condition. Here, R is considered to be a function of drop size and velocity and k is the specific attenuation expressed as a function of the extinction cross section of the rain drop at a given frequency [START_REF] Hahn | Light scattering theory[END_REF]. The detailed description about that relationship can be found in [START_REF] Atlas | Path-and area-integrated rainfall measurement by microwave attenuation in the 1-3 cm band[END_REF][START_REF] Olsen | The aR b relation in the calculation of rain attenuation[END_REF]. Reference rainfall intensities were obtained by a rain gauge located nearby the link and by converting reflectivity maps recorded by two C-band weather radars with a spatial (2.5 × 2.5 km 2 ) and temporal (5 min) resolutions. These two instruments were separately used to validate the link approach. It should be noted that the second reference data set, which came from weather radar, is the weighted average of the radar pixels along the link. Authors used 8 rain events to test the performance of the method. It was found that the estimates based on the single-frequency approach have a closer agreement with those measured by nearby rain gauge if rainfall variability along the link is weak. On the other hand, when the rainfall variability increases, the link based measurement was found to be closer to the one estimated by weather radar. This could indicate that merging the microwave links with traditional measurement techniques may better capture the spatio-temporal variability of rainfall. Studies regarding this aspect will be further discussed in subsection 3. Microwave links are not only capable of measuring rainfall, but also applicable to measure other forms of precipitation. For example, [START_REF] Ostrometzky | Accumulated Mixed Precipitation Estimation Using Measurements from Multiple Microwave Links[END_REF] applied the rain attenuation data recorded by multiple commercial microwave links for measuring the precipitation amount regardless of its type, i.e. whether it is rain, sleet or their mixture. Even, their application in fog monitoring can be found in (David et al., 2015).

Combining microwave links with weather radar and rain gauges

Ideally, using microwave links in combination with a rain gauge or a weather radar can improve the reliability of rainfall estimates. This is certainly beneficial to locations where the traditional techniques often suffer from drawbacks. In fact, this is one of the motivations for exploring the applicability of microwave links for rainfall measurement. There has been a few investigations focused on this aspect with different objectives and methodologies which we will discuss below.

Usually weather radar based rainfall estimation is adjusted on the basis of rain gauge(s).

However, it is a fact that such adjustment is a subject to various sources of uncertainties. These are related to malfunctioning of rain gauges caused by wind and temperature effects, and due to interpolation techniques used to distribute point measurement values over a surface area. The error caused by the latter factor becomes severe in very high rain intensity with strong variability. Besides, finding a convenient location to install a rain gauge is not easy in urban areas.

In this context, using 30 km long vertically polarized dual-frequency microwave link, [START_REF] Rahimi | Attenuation calibration of an X-band weather radar using a microwave link[END_REF]) calibrated X-band weather radar located in Essen (Germany). In their experiment, the signal attenuation measured along the link was applied to correct radar reflectivity data that has been degraded because of attenuation of radar beam. The corrected radar reflectivity was then converted to rainfall rates on the basis of Z-R Marshall Palmer relation [START_REF] Marshall | The distribution of raindrops with size[END_REF]) to be able to compare with ground-truth measurement. The authors obtained the ground-truth rain rate from a network of 5 rain gauges closely placed along the link path.

Based on the comparison between the path averaged rain rate computed along the link and the ground-truth, the approach was found to be effective, especially, in convective rain events to calibrate the radar reflectivity data. However, the authors state that the capability of the link for radar calibration is limited since it corrects the radar sector that are only located in the link path.

In addition to this finding, [START_REF] Cummings | Using microwave links to adjust the radar rainfall field[END_REF]

Measurement error sources and uncertainties

The encouraging results from different applications of microwave links must not hide the presence of several error sources whose influences were investigated in many articles. Signal attenuation measured along the microwave link is not only due to rainfall, but also to other factors which are essentially related to:

-Instrumental impairments such as signal fluctuations by antenna wetting, baseline effects (gas, water vapour, temperature, wind) and signal quantization (round-off errors) in hardware equipment of both transmitter and receiver; -A measurement model, i.e. k-R relation whose parameters depend on the signal frequency and polarization, the variability of drop size distribution as well as temperature.

Here, we discuss the issues regarding the impact of both factors on the accuracy of rainfall measurement by microwave links.

Error sources due to instrumental impairments

Microwave signal fluctuations are not only due to rainfall, because other instrumental impairments may also cause a signal variation. These impairments include the effects of baseline (i.e. gas, water vapour, temperature), antenna wetting, and quantization of signal. Therefore, a key issue here is the possibility of separating rain induced attenuation from these factors.

Various studies addressed this issue on the basis of both research dedicated (Holt et Next error source is signal quantization which arises from the hardware equipment of the microwave antenna. The power resolution in this equipment depends on the accuracy of the measurement device. In general, the signal power is quantized from 0.1 dB upto several dBs, but very extreme cases, namely 0.1 and 1 dB are commonly used in practice. The higher signal power resolution is the more accurate the measurement device will be. According to [START_REF] Leijnse | Microwave link rainfall estimation: Effects of link length and frequency, temporal sampling, power resolution, and wet antenna attenuation[END_REF] the degradation of the power resolution at 1 dB causes severe errors and uncertainties for short links and at low frequencies. [START_REF] Zinevich | Prediction of rainfall intensity measurement errors using commercial microwave communication links[END_REF]) demonstrated that the quantization error for shorter links is larger than for the longer links. It is due to the fact that the path integrated attenuation measured over longer links is larger compared to along shorter links.

However, it depends on the spatial variability of rainfall. The main difficulty in identifying the quantization error can be due to low rain intensities (Leijnse et al., 2007c). This becomes severe if the quantization error at 1 dB is considered in light rain event. For example, it has been shown by (Leijnse et al., 2007c) that the effect of quantization error at 1 dB is considerably large in low intensity and long lasting rain event. However, the error in rainfall estimation along the link with a certain quantization interval does not exceed half of that interval value (Widrow and Kollár, 2008). In general, the k-R relation applied at the higher frequencies results in larger attenuation compared to the one applied at lower frequencies. The reason is that microwave links with high frequencies are more sensitive to rain rate, see 2.9 in chapter 2. In addition, the impact of signal polarization (either horizontal or vertical) becomes important if the shape of rain drops is non-spherical. However, this effect seems to be negligible for small rain drops.

Error sources due to a measurement model

Various studies recognise the critical role played by the influence of drop size distribution on the k-R model. The findings, conclusions and arguments regarding this aspect slightly differ from each other as follows:

(Atlas and [START_REF] Atlas | Path-and area-integrated rainfall measurement by microwave attenuation in the 1-3 cm band[END_REF] (Berne and Uijlenhoet, 2007) studied the impact of the spatial variability of drop size distribution including the influence of link length and frequency on the estimated path averaged rain rate along the link. In their experiment, the path averaged rainfall rates were obtained using the microwave links at frequencies from 5 to 50 GHz with link lengths from 0.5 km to 30 km.

The rainfall data is based on the simulated profiles of drop size distribution which were recorded during an intense Mediterranean rain events over 2 years period. It has been found that the link length does not have a significant impact on the k-R relation parameters used to estimate rainfall rate, but the frequency does. In contrast to [START_REF] Atlas | Path-and area-integrated rainfall measurement by microwave attenuation in the 1-3 cm band[END_REF], the authors reported that the exponent of the k-R relation, i.e. b, is almost equal to 1 at 30 GHz. Further, the authors state that such linearity can be achieved at 37 GHz as well. GHz to 18 GHz as the link length grows. However, these conclusions are under the assumption that the effects of the sampling errors on the k-R relation coefficients are negligible.

Overall, the conclusions drawn from these investigations are in agreement in terms of quantifying the influence of each of those error sources that will be introduced in the attenuation measurement in chapter 4. However, it should be noted that the true nature of these error sources

is not known whether it can be due to the measurement model (subsection 3.3.2) or due to the instrumental impairments (subsection 3.3.1). In the scope of our study, quantifying the error sources in connection with rainfall induced attenuation measurement along the microwave link can be made on the basis of these findings. Therefore, our conclusion regarding the influence of these error source is generalized rather than relying on a single conclusion of a particular finding. The nature of the introduced error source will be simulated in chapter 4. That is to say, the influence of error sources on the attenuation measurement grouped by instrumental (subsection 3.3.1) and the measurement model, i.e. k-R relation (subsection 3.3.2) is formulated depending on the link length and frequency ranges.

Rainfall mapping using multiple links

An important and interesting issue is the mapping of rainfall fields using path integrated attenuation data measured from a network of commercial microwave links. As confirmed in the work previously discussed in section 3.2, it is possible to convert the signal attenuation into path averaged rainfall rate along the link using k-R relation. process is not straightforward and requires a specific method to derive rainfall maps based on path integrated attenuation measurements along the microwave links.

The application of such specific methods to a network of commercial microwave links can be challenging due to many reasons. It is a fact that their characteristics, location, topological structure as well as frequencies bands are inhomogeneous and designed to maintain the communication quality. Studies attempted to map rainfall fields based on the real microwave attenuation data differ from each other in terms of processing techniques and algorithms.

The idea of rainfall mapping based on the rain attenuation data from multiple links has been first investigated by [START_REF] Giuli | Tomographic reconstruction of rainfall fields through microwave attenuation measurements[END_REF]. The authors applied a tomographic inversion algorithm [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]) which uses a linear combination of basis functions to reconstruct 2D specific attenuation map from one dimensional (1D) rain attenuation vector. Then, the retrieved specific attenuation map was directly converted into rainfall map on the basis of the k-R relation. However, the proposed retrieval system was limited due to the assumptions that (i) the k-R relation is linear (ii) the geometry of the microwave link topology is homogeneous.

On the one hand, these assumptions are the main conditions for the applicability of the tomography framework. On the other hand, they are not valid if real commercial microwave links are considered. An important challenge in the context of tomography is that the reconstruction problem is ill-conditioned if the number of attenuation measurements is much smaller than the estimation area. However, the proposed framework gave a deep insight into the introduction of the tomography principle in rainfall mapping by microwave links. Later, this approach was improved by [START_REF] Giuli | Microwave tomographic inversion technique based on stochastic approach for rainfall fields monitoring[END_REF]. In the same line of these findings, [START_REF] Cuccoli | Tomography applied to radiobase network for real time estimation of the rainfall rate fields[END_REF] proposed to reconstruct rainfall fields from path integrated attenuation measurements on existing communication radio links that operate at 18, 23 and 38 GHz by means of tomographic processing method that they refer as "Combined Deterministic-Stochastic Retrieval Technique".

The specific attenuation fields obtained by the algorithm were converted into rainfall fields on the basis of the k-R relation after the reconstruction process. The authors obtained interesting results when testing this method by a simulation approach applied to the real case study: rainfall fields observed by C-band radar and actual link network of the city of Florence (Italy). proposed a stochastic interpolation algorithm, which is based on the Inverse Distance Weighting function [START_REF] Shepard | A two-dimensional interpolation function for irregularlyspaced data[END_REF] to map rainfall over a irregular space grid. It was shown that the algorithm is capable of estimating the spatial rainfall over an arbitrary geometry network in the presence of large quantization noise, i.e. 1 dB. The authors state that the rainfall information by a rain gauge can be fed into the proposed algorithm to improve the retrieval accuracy.

At a very different scale, [START_REF] Overeem | Country-wide rainfall maps from cellular communication networks[END_REF] demonstrated the capability of a huge network of 2400 links to retrieve the space-time dynamics of rainfall over the 35500 km 2 area of the Netherlands by using Ordinary Kriging interpolation method [START_REF] Cressie | The origins of kriging[END_REF]). This method is suitable for interpolation of highly irregular-spaced rainfall estimate points and takes into account the local variability of the rainfall. The performance of the algorithm highly depends on the spatial variances of the rainfall which is also called semivariogram. Since it was reported in 1991, rainfall mapping by multiple microwave links has been attracting a considerable interest. However, a limited number of studies have been devoted to the mapping of rainfall fields from the path integrated attenuation on these links and the feasibility of the rainfall mapping has not been fully established yet. In particular, on the suitability of the tomographic approach, the relative importance of the issues (i.e. adequate methods for tomographic reconstruction and adjustment procedures and the influence of the nonlinearity of the k-R relation on the tomographic system) has been subject to considerable discussion.

So far, there has been very little research directly investigating the tomographic approach and yet no attempts have been made on transform-based tomography. In this context, this study aims to contribute to this growing area of research by exploring the capabilities of commercial microwave links for rainfall mapping at the urban scale. In particular, this research examines two completely different mapping techniques which will be presented in chapter 5 and chapter 6.

Summary

In conclusion, this chapter has attempted to provide a detailed review on the state of the art relating to microwave link based rainfall measurement. All of the studies reviewed here

Case Study

Introduction

This chapter presents the methodology of the current study which includes a study area, devices, tools and assumptions used to obtain rain attenuation data along microwave links. For this purpose, we adopted a simulation approach corresponding to a realistic case scenario (i.e to mimic the signal attenuation data measured along the microwave links) to study the feasibility of rainfall mapping. The case study is carried out using microwave links employed by cellular network and rainfall maps observed by weather radar at high spatio-temporal resolutions. The generated rain attenuation data in this chapter is considered as 'artificial data' which will then be used to map rainfall (chapter 5 and chapter 6). The fundamental purpose for doing this is to substitute the generated rain attenuation data for signal levels measured by a real microwave link network.

Usually, signal attenuation data can be obtained by the difference of signal levels recorded between transmitter and receiver stations. However, a major concern in the use of real signal data is to separate rain induced attenuation from other effects which are often a subject to special instrumental set-up procedures as discussed in subsection 3.3.1. Even after having performed the cleaning procedures (e.g. identifying the dry and wet periods), the most important issue is to compare the reconstructed rain map with true rain field which, in fact, does not exist even though weather radar or rain gauges are able to provide high precision rainfall information as a 'ground-truth'.

We adopted a simulation approach that does not suffer from above-mentioned challenges and it gives the following advantages:

1. The main advantage is the existence of true rain field (i.e. this is based on weather radar in our study) which enables to validate reconstructed rain field. So, we know the true rainfall. In this case, the bias related to the 'ground-truth' rain field does not exist;

2. It is convenient for verifying the feasibility of rainfall mapping using the geometry of a real network topology without any data processing and cleaning steps which are challenging in the case of real data;

3. Applicability and adaptability in realistic conditions by changing the scenario and objective of the study: (i) data collection processes can be manipulated to understand the influence of the quality of data on the mapped rainfall fields, (ii) If needed, the real signal data recorded by microwave links can be replaced with the generated rain attenuation data for the same rainfall mapping process.

However, the main difficulty associated with the simulation approach is the fact that it might not reflect a nature of real data entirely. In addition, the limitation is clearly the gap between so called 'artificial data' generated by applying a measurement attenuation model, i.e k-R relation with simulated various measurement errors, and real data measured by antenna stations of microwave links. The main issue is to decrease that gap as much as possible. Therefore, the objective of this chapter is to include such challenges in the rain attenuation data generation process based on microwave links.

This chapter is organised as follows. First, the presentation of the study area and the spatial distribution of microwave links are given in section 4.2. Then, rainfall fields captured by weather radar and characterisation of each rain event type are given in section 4.3. Next, the rain attenuation data generation process is formulated in section 4.4, and the generated attenuation data and its statistical analysis are discussed in section 4.5. Finally, section 4.6 summarizes the chapter.

Presentation of microwave links network 4.2.1 Study area and microwave links

The case study is carried out in Nantes located in Loire Atlantique department in western France (see For the convenience, one can look at the zoomed-in part of the figure in which the selected number of microwave links at 18, 23 and 38 GHz are depicted in blue, red and green lines, respectively.

The lengths of the microwave links range from 0.3 to 16.8 km with 4.2 km on average. It can be seen in Figure 4.1 that the density of the links decreases from the city center towards the countryside. It is mainly because the cellular network system is designed for increasing the communication quality in more populated areas, i.e. in the city center. Therefore, shorter links are more concentrated in the urban areas whereas the longer ones are mainly located in the suburbs.

To have a better understanding about the length distribution one can see Overall statistical description of the microwave link lengths can be found in

Classification of the microwave links density

It is important to define the density level of the microwave links network because it allows to classify dense and sparse regions in the study area. The reason for this classification is to understand the influence of the network topology on the accuracy of the rainfall retrieval.

Below, we will demonstrate how the network density level is computed. The result of this computation will be called pixel density map which can be used to validate the network system performance over different regions in rainfall mapping process.

A basic assumption is that the higher density of the links should give the higher weights over the surface. To obtain the pixel density map the following steps are performed:

1. Divide the study area into n number of pixels with the same size. In our case, the pixel size is chosen to be 2 km × 2 km resolution, see the grid size of 2. Compute the intersection length of each link with each pixel: L j = {l j1 , l j2 , l j3 , ...l ji } 3. Define the weights of the each crossed part of the link:

W ji = l j1 L j , l j2 L j , ... , l ji L j (4.1)
4. Compute the pixel density:

D i = m j=1 W ji (4.2)
where, j -link index, i -pixel index, m -number of links, W -weight of the link along its crossed part, D -pixel density, l -length of the crossed part of the link, km, L -Link length, km.

The output of the calculation procedure presented above gives the pixel density, i. higher than its median. We consider the density level to be low if the condition D ≤ D median is fulfilled, in our case, D median = 0.5. Similarly, the density level is moderate if D is between 0.5 and 1. Beyond the moderate density level, it is considered to be sufficiently high. The condition of D ≥ 1 was chosen for the higher density level because we assume that it is acceptable if about 30 % of the D values is greater than 1. A more reasonable choice for defining these thresholds can be 1st, 2nd and 3rd quartiles of D. Our choice is very close to these thresholds.

For example, upper threshold in 3rd quartiles for the given D values is 1.4 while our choice is 1. The radar data was cleaned from the effects of ground clutter, isolated pixels and partial beam blocking. The detailed description about the scanning strategy of this radar and reflectivity conversion into rain rate can be found in [START_REF] Emmanuel | Temporal and spatial variability of rainfall at the urban hydrological scale[END_REF]. The weather radar rainfall events are used to generate attenuation measurements along the links at a resolution of 0.25 km × 0.25 km and considered as reference rainfall fields.

Rainfall events

We selected 7 rain events, regrouping 207 rainfall maps in total, which represent 4 types of rainfall, namely light rain, shower, organised and unorganised storm, characterized by different spatial and temporal variabilities. A description about duration, period and the number of rainfall maps in each event is given in Table 4.4. Table 4.5 gives a summary statistics of rainfall data set obtained for light rain, shower, organised and unorganised storm events, separately. This statistics was calculated excluding the part of the rain maps which consist of zero values. It is worth mentioning that the statistics for each column was computed by flattening the values of rainfall maps in each event. Here, the rows count, mean, std, max stand for the number of maps, average, standard deviation, maximum value of rain rate for each event, respectively. The rows q1, q2 and q3 give the 1st, 2nd, 3rd quartiles of the rainfall data. We can see that the variability of rainfall is quite strong in organised, unorganised storm and shower events with a standard deviation being equal to 6.8, 10.8 and 8.2 mm.h -1 , respectively. However, 75% of all the data in total do not exceed 3 mm.h -1 . Overall, the strongest variability can be observed in unorganised storm event while the weakest variability belongs to light rain.

Computing path integrated rain attenuation along the link

We compute a path integrated rain attenuation along the link using the principle of ray intersection problem presented by [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]. This computation is based on the assumptions defined in the previous subsubsection 4.4.1.1. In the context of our study, a ray is understood as a straight line that connects a transmitter and receiver of the microwave link. This ray is intersected by certain pixels R i which divide the link L j , into different lengths l ji , see the zoomed-in part of the Figure 4.6. Note that the pixel index is the same as the crossed part of the link index. Then, the path integrated attenuation A j due to rain rate R i along a microwave link L j can be computed based on the k -R relation as follows:

A j = a j n(j) i=1 l ji R i b j + ε j (4.3)
Where, n(j) -number of pixels crossed by the jth link; R i -rain rate in the ith pixel, mm.hour The first type of error, which is due to the measurement model uncertainty can be denoted as ε j is supposed to be unbiased, complied with zero-mean (µ = 0) Gaussian distribution whose standard deviation σ m j is equal to α percent of the computed rain attenuation A j :

ε j = µ + σ m j N (0, 1) (4.4)
Where, µ -zero mean of the Gaussian distribution; σ m j = αA j ; α -the percentage variable;

N (0, 1) -Standard normal distribution function.

The measurement error caused by the k-R relation has been generated using Equation 4.4.

Introducing variable standard deviation along each link, i.e the value of σ m j , is justifiable and realistic because error source affecting the attenuation measurement at one antenna station is totally independent from the one occurring at another station.

Obviously, two questions can arise from the assumptions introduced about the nature of the error source in Equation 4.4: (i) why it is Gaussian and (ii) the core reason of variable α percentage of the standard deviation σ m j . The assumption of Gaussian distribution is selfevident because a real nature of possible error sources caused by measurement model, as a whole, is unknown. Therefore, a reasonable assumption is to assume that this error is 'Gaussian' type. Regarding the 2nd assumption, i.e. α percentage, our justification is explained as follows:

The value of α is chosen to be 5% because 95% of our confidence about the accuracy of the k-R relation is 2 * α which gives 10% average error. It has been confirmed in many studies (section 3.3 in chapter 3) that the influence of drop size distribution on the k-R relation is limited, showing 10% average error, if the linearity of this relation is dominant, i.e. the exponent of the formula 2.26 (chapter 2) is close to 1. In fact, this is clearly applicable to our case because the exponents of the k-R relations at frequencies (18, 23 and 38 GHz) are almost 1. Larger value of α indicates that our confidence about the accuracy of the rain attenuation measurement is lower.

However, we will also test the higher percentages of α to show that the attenuation measurement might also be influenced by larger amount of error source. Thus, the percentages of α, i.e. 5%

and 20%, will be considered to be two extremes cases of all the possible error sources in rain attenuation data generation process.

The second type of error is caused by signal quantization and it is introduced after computing the rain attenuation A j in Equation 4.3. The quantization error occurs when the signal is truncated at a certain precision. Typical precision values are 0.1 and 1 dB at which the actual signal is truncated. The variance of the signal quantization error is uniform and equal to σ q 2 = ∆ 2 12 (Widrow and Kollár, 2008). We quantized A i at ∆ = 0.1 dB to reflect the signal quantization impact on the rain attenuation measured along the microwave link. It has been well-understood in many studies (see, subsection 3.3.1 in chapter 3) that the quantization error source is a major trouble for shorter links in lower intensity rain events. The larger quantization value is an indication to the lower accuracy and precision of the hardware equipment of the antenna station.

It is important to note that the path integrated rain attenuation A i can get negative values if randomly generated error ε j exceeds the actual rain attenuation value. Therefore, such values are set to be 0.001 dB to avoid negative attenuation. Thus, the signal attenuation caused by rain was generated in the presence of those two common error sources.

Rain attenuation generation protocol

In summary, to generate the rain attenuation data over the given microwave links network in Figure 4.1, the following steps are performed:

Step 1. Discretize the study area with 0.25 km × 0.25 km grid which is the spatial resolution of rainfall map;

Step 2. Compute the intersection of microwave links with the discretized area, i.e l ji ;

Step 3. Extract a rainfall map of weather radar (given in subsection 4.3.2) and super-impose it over the discretized study area. In this case, the pixel index of the rainfall map is expected to be the same as the index of the discretized area;

Step 4. Compute the path integrated attenuation along each link (Equation 4.3) with the presence of k-R model error at α = 5% and 20%. These two cases are independent from each other;

Step 5. Quantize the measured rain attenuation at 0.1 dB resolution.

The generated rain attenuation data will be used in rainfall mapping process chapter 5 and chapter 6.

Generated rain attenuation data

We generated the rain attenuation data using the attenuation measurement protocol described in subsection 4.4.3. These 5 steps are repeated for 207 rainfall maps which are grouped by light rain, shower, organised and unorganised storm events. shows a 'close-to-reality' type of the measurement error sources. The maximum error on average in each event can be observed for the link lengths between 5 to 17 km as expected in the literature, see, section 3.3 in chapter 3. Another interesting pattern of the generated error is its increasing standard deviation value in higher variability of rainfall. Evidentially, the maximum error along all the links in light rain event does not exceed 0.5 dB while it is nearly 2 dB in unorganised storm event. For example, this is clearly visible if we compare the values of box plots between each event, see the difference between light rain and shower.

However, one unexpected finding was the extent to which the nature of the error introduced in this experiment does not mimic properly if the link frequency variability is considered (see, A possible explanation may be due to the simplified assumption (i.e. random and unknown) about the true nature of the error in terms of frequencies in the generated attenuation.

Simulating the error in the presence of such behaviour still helps understand the influence of error magnitudes on the attenuation measurement process. However, this can be considered as a limitation of the error nature introduced in our experiment. To compare the 5% error magnitude as discussed previously, we present the generated attenuation data in the presence of α = 20% error magnitude with same quantization at 0.1 dB. Figure 4.10 shows the standard deviation of the quantified errors at higher magnitude(α = 20%). It can be seen that in all cases the average standard deviation is 4 times bigger than the average error shown in Figure 4.9. However, we argue that introducing larger error magnitude at 20% seems to be unrealistic since the real nature of the error itself is not known. Therefore, we consider that the parameter α at 5% is acceptable to mimic the influence of the possible error sources on average. Overall, the average trend of error generated for this experiment is not far from reality.

Summary

In this chapter, the main objective was to describe the case study that uses the data sets (weather radar rainfall maps, microwave links) and a methodology for generating rain attenuation data along the microwave links. In particular, we presented the procedure of how the signal attenuation due to rainfall is generated along the microwave link in a realistic case. The generated signal attenuation data obtained from this case study will be substituted for real signal attenuation data recorded by microwave antenna stations of the cellular network. For this purpose, we used cellular network and hundreds of weather radar rainfall images in order to prepare the signal attenuation data. Besides, we also showed how the density of microwave links network can be classified into different levels. This will be used to validate the dense and sparse regions of the microwave link network separately.

Retrieval Model 1. Inverse Algorithm

Introduction

This chapter presents the first retrieval model that uses the signal attenuation data generated along multiple microwave links to retrieve rainfall fields. As far as rainfall retrieval based on microwave links is concerned, the signal attenuation data generated in chapter 4 is considered to be data, and rain rates along the link which triggered that attenuation are called parameters. In this regard, the retrieval of rainfall fields and rain attenuation measurement are inverse and direct problems, respectively. The reason is that the path integrated attenuation is partly representative of the spatial distribution of rainfall along the microwave link.

The retrieval model to be presented here employs the principle of inverse problem to find those unknown parameters based on the data (Menke, 2012). Inverse problems are often subject to non-uniqueness meaning that the solution we are seeking can be more than one. In order to single out one solution among other possible ones a priori knowledge is applied. The inverse algorithm to be presented in this chapter uses such a priori knowledge.

This chapter is organized as follows. The formulation of the rain retrieval model on the basis of inverse problem is presented in section 5.2. The algorithm used in this model is detailed in section 5.3 and the application conditions before applying the algorithm are demonstrated in section 5.4. The application of the inverse algorithm is presented in section 5.5. The statis-tical metrics used to evaluate the performance of the algorithm are given in section 5.6. The sensitivity analysis to the application conditions are given in section 5.7.

Formulation of the problem

Main assumptions

The microwave link network forming a two dimensional discrete space has been shown in Figure 4.6. The rainfall vector R which groups the rain rates in all crossed pixels of the study area are the unknowns to be determined based on the attenuation data vector A, see more details of these notations in (4.

3). The physical model which relates the path integrated attenuation

A j along the jth microwave link to the rain R i in ith pixel is given in the formula (4.3), in chapter 4, using which A vector has been generated over the microwave links network based on the available weather radar rainfall maps. The retrieval problem is formulated under the following assumptions:

-The physical model are known (see, eq. 4.

3), and this model is affected by different sources of errors.

-The nature and shape of those error sources are random, and assumed to comply with zero mean Gaussian probability distribution.

-The retrieval problem is non-linear since the operating frequencies of the microwave link network are not constant [START_REF] Das | Rain attenuation modeling in the 10-100 GHz frequency using drop size distributions for different climatic zones in tropical India[END_REF]23 and 38 GHz).

Based on such assumptions, the problem is to retrieve those pixels crossed by microwave links as shown in Figure 4.6. In order to include uncrossed pixels by the links we need to introduce additional equations to the formula (4.3) which are not available. Therefore, a reasonable approach can be to formulate the problem based on only the sampled pixels. Then, interpolation techniques such as Inverse Distance Weighting [START_REF] Shepard | A two-dimensional interpolation function for irregularlyspaced data[END_REF] or Kriging [START_REF] Creutin | Objective analyses and mapping techniques for rainfall fields: an objective comparison[END_REF], which takes into account the spatial structure of rainfall, could be applied to fill in the unsampled pixels. However, application of interpolation techniques in such cases is beyond the scope of this current study. Instead, we introduce a new approach, alternative to standard interpolation techniques, for filling in the unsampled pixels based on the nearby the sampled ones.

By doing so, we increase the rainfall retrieval area. We will give more detailed description regarding this aspect in section 5.5.

Determining the state of the problem

The length of data A and the number of unknown parameters to be estimated R are denoted as L A and L R , respectively. Then, the state of the problem can be described as follows:

1. L R < L A : The problem is over-determined. For example, we can see a part of the network located in the city center where the richness of information in terms of spatial density is very high. More than one microwave link carry rainfall information for the same pixels in many locations, see the city center in Figure 4.1. 3. L R = L A : The problem is even-determined. We could have seen even-determinedness in a denser part of the network if the pixel is equally crossed by many links. However, it is nearly impossible to find such type of the state of the problem in the network area separately.

L R > L

Globally, the state of the problem depends on the pixel resolution and the errors affecting the measured data since the assumption of constant rainfall is valid for each pixel area. The problem type would be certainly over-determined if the study area was described at a lower resolution and with more links intersecting all the pixels. However, we do not see such case stated whether under or over-determined type; therefore, it depends on the network density and the pixel resolution. In this study, the emphasis is given on more under-determined context.

As far as under-determined problems are concerned, the definition of apriori knowledge is of great importance. Therefore, the solution is usually complemented by the apriori information. We demonstrate such a solution vector using the nonlinear inverse algorithm proposed by [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF]). This will be discussed in the next section.

Inverse algorithm

This section describes the retrieval algorithm applied to identify the rainfall vector R demonstrating the best compromise with respect to the maximum likelihood between the observed values of attenuation data and staying close to an a priori guess to be defined. To do so, the identification is performed by an iterative non-linear algorithm [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF] which consists of minimizing the following expression:

φ(A, R) = (R -R 0 ) T • C R -1 • (R -R 0 ) + (A(R) -A 0 ) T • C A 0 -1 • (A(R) -A 0 ) (5.1)
where, A(R) -Model which relates the rainfall intensity to attenuation; φ(A, R) -Likelihood function; A 0 -Total attenuation vector, i.e. path integrated attenuation data generated along the microwave links; C A 0 -Uncertainty covariance matrix of vector grouping path integrated attenuation along the microwave links; R 0 -A priori rainfall vector; C R 0 -Uncertainty covariance matrix of a priori rainfall vector.

The exponents of '-1' and 'T' stand for inverse and transpose, respectively. The covariance matrices (denoted as C) indicate the respective confidence allocated to the measured attenuation data (i.e. generated in chapter 4), and to the values of R vector in the form of a priori knowledge.

Their specification allows establishing a balance between attenuation data and a priori rain rate information. We will discuss the details of these parameters, i.e. R 0 , A 0 , C R 0 , C A 0 in section 5.4. The statistical distributions of both (R -R 0 ) and (A(R) -A 0 ) are Gaussian and unbiased.

The solution vector that minimizes (5.1) is demonstrated by [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF] as follows:

R = R 0 -C R 0 • G T • (C A 0 + G • C R 0 • G T ) -1 • [(A(R) -A 0 ) -G • (R -R 0 )] (5.2)
where, G is the matrix of partial derivatives (first order) of the rain attenuation model. On the left and the right handside of (5.2), we can see the same vector R which indicates that the first guess about the solution is needed. Iterative form of (5.2) can be expressed as:

R t+1 = R 0 -C R 0 • G t T • (C A 0 + G t • C R 0 • G t T ) -1 • [(A(R t ) -A 0 ) -G t • (R t -R 0 )] (5.3)
where t -iteration number. The values of the matrix G t are defined as follows:

G t (j, i) = ∂A j (R ti ) ∂R ti (5.4)
where, A j (R ti ) is the attenuation model applied along the jth link to obtain attenuation value at R ti rain rate of the ith pixel.

Diagram in Figure 5.1 illustrates the minimization procedure of (5.1) in which the terms

Yes

Set inputs:

1. ϕ 0 -initial likelihood 2. R 0 -a priori rain vector R t and φ t are estimated. The iterative procedure starts at an arbitrary point. The first guess is needed to initialize the algorithm at 1st iteration, i.e. t=1. It is important to mention that the first guess vector R t at t=1 is different from apriori rainfall vector R 0 . However, this first guess can be accepted as R 0 . Therefore, a reasonable choice for starting the iteration is to assume that the first guess R 1 is equal to R 0 [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF].

The iteration process of the inverse algorithm is shown in the diagram Figure 5.1. The algorithm is stable if the problem is not too non-linear and does not display many local minima in the minimization process. The a priori knowledge heavily influences the solution vector R t+1 if the stated problem tends to be under-determined. In case, the apriori information is sufficient enough to describe the system, the iteration starts and immediately stops at t = 1 point. Because, the relative difference between estimated likelihoods at 0th and 1st iteration (i.e. φ 0 and φ 1 , respectively) becomes negligible. However, this 'negligible' interval depends on a user choice. We stopped the iteration if the relative difference is equal to or less than 1%, see the notation ε in Figure 5.1. If the difference between R t+1 and R 0 is reasonably small, it means that covariance operators smoothed the solution enough. Detailed description about the stability and convergence of the algorithm can be found in [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF].

The conditions of the a priori knowledge and the parametrization of (5.3) will be explained in section 5.4.

Application conditions of the inverse algorithm

In this section, application conditions of the algorithm are presented. More specifically, we discuss the details of the apriori information used to initialize the algorithm (5.3). The apriori represents an initial knowledge and our confidence about data and parameter. It is a driving part of the inverse algorithm; therefore it is important to parametrize it and associated uncertainties by which the algorithm can be influenced.

We consider that uncertainties in generated attenuation data are totally independent of uncertainties in rainfall rate vector. Thus, the apriori knowledge can be classified by two groups:

-The apriori on parameters. This represents the apriori rainfall vector, R 0 , and its associated error structure which is described in the form of covariance matrix C R 0 . It will be detailed in subsection 5.4.1.

-The apriori on data. This denotes the attenuation data, A 0 , generated along the microwave links, and our confidence about this data expressed in the form of covariance matrix C A 0 . This will be described in subsection 5.4.2.

The a priori rainfall and associated covariance matrix

The apriori knowledge on algorithm parameters is described by the rainfall vector, R 0 in mm.hour -1 and by our confidence, which is expressed as an associated error covariance matrix C R 0 . A careful choice of R 0 heavily influences the algorithm outcome if the state of the problem is under-determined. The apriori rainfall vector R 0 can be obtained by either a result of a previous inverse problem run or using a theoretical model by which a relation between attenuation data and rainfall rate vector can be established. The choice of apriori knowledge should not be complicated, so a simple approach adopted in this study is to apply k ∼ R relation along the link in order to compute R 0 from the generated data A 0 itself. However, applying such approach can only give one scalar value of average rainfall rate along the link, but not for each value of R 0 vector. If the variability of rainfall is strong along the link this approach is clearly not stable due to the constant rainfall assumption. Therefore, we initialize the algorithm with different apriori choices.

We consider two choices of the apriori rainfall:

1. Global a priori: The a priori rainfall vector R 0 values are the same everywhere and equal to the mean value of rainfall rates estimated from all the links. In this case, the initialization of a priori rainfall is called 'global apriori'. It is computed in the following order:

(a) Compute the rain rate using the attenuation value A 0j generated along jth link:

R j = ( A 0j a * L j ) 1/b
(5.5)

The formula (5.5) applies to all the links. The term R j is the rain rate (assumed to be average) along the jth link with a length of L km. Then, the initialization type is called 'local apriori'. The only condition is that k should remain constant for each value in the pixels. Its computation is performed as follows:

(a) Obtain rain rate R j along the jth link using the (5.5);

(b) Group k number of closest links for each pixel: The Euclidean distance is computed between the ith pixel center (see the black dot labelled as R 0i ) and microwave links.

This distance is depicted as the dashed-line in Figure 5.2. Only k number of links, which are the closest to the ith pixel, are chosen: {R i , ... , R ik }.

(c) Compute the mean value of rain rates obtained from those grouped links for the ith pixel:

R 0i = k j=1 1 k R ij (5.8)
Next term to be defined is C R 0 which is the error covariance matrix of the apriori rainfall.

Each element of C R 0 represents our confidence denoted as σ R 0 . In the microwave link area, the error in one pixel (denoted as x) is spatially correlated to another pixel (denoted as y) in terms of distance. This spatial structure of the a priori rainfall is assumed to be exponential:

C R 0 (x, y) = σ R 0x × σ R 0y × exp (- 1 2 × D 2 xy ∆ 2 ) (5.9)
where, σ R 0x and σ R 0y are the standard deviations of the apriori rainfall error in xth and yth pixels, respectively; D xy -Euclidean distance between xth and yth pixels, km; ∆ -decorrelation distance, km.

It is important to note that the calculation of σ R 0 for global apriori differs from the local apriori. In the global apriori, σ R 0 is the same in each pixel and computed as follows:

σ R 0i = m j=1 (R j -R) 2 m
(5.10)

The terms R j and R are given in the formulas (5.5) and (5.6), respectively.

If the local apriori is considered, σ R 0 is obtained as follows:

σ R 0i = k j=1 (R ij -R 0i ) 2 k (5.11)
The terms R ij and R 0i are obtained from the formula (5.8). It should be noted that the index i in the standard deviation can be replaced by either x or y to define the corresponding σ R 0x or σ R 0y in (5.9), respectively.

The attenuation data and associated covariance matrix

The generated data, A 0 , is considered to be equivalent to real signal attenuation recorded from the cellular network, expressed in decibel [dB]. In our study, the result of the generated data represents the path-integrated attenuation along the link, and its standard deviation can be denoted as σ A . The a priori on the attenuation data is defined by two terms of the algorithm (5.3): (i) The total attenuation vector A 0 ; (ii) covariance matrix C A 0 . It is summarized in the form of covariance matrix:

C A 0 =      σ 2 A 1 • • • 0 . . . . . . . . . 0 • • • σ 2 A j     
(5.12)

The reason why covariance has zero values except diagonal term is due to the fact that the microwave link in one location is separate equipment and not affected by the uncertainties occurred in another link. The diagonal terms of C A 0 represent a variance of the generated data.

Generally, this variance consists of the uncertainties sin the generated attenuation data σ 2 A 0 and in the model σ 2 model :

σ 2 A = σ 2 A 0 + σ 2 model (5.13) Where, σ 2 A 0 is the variance of A 0 ; σ 2 model = (β × A 0 ) 2 , β -a magnitude of the error, [%].
The main reason for adding the σ 2 model is that the k-R relation used as a model (i.e. A(R) in (5.2)) could be itself erroneous. [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF] explains that this model error concerns the state model of the algorithm and is fully independent of the measurement error.

We consider that true value and/or percentage of the error σ 2 data is unknown. For example, the cellular network may not work correctly when the signal attenuation is recorded in a real time. Therefore, it is subject to understanding the influence of the model error on the retrieval accuracy of the algorithm.

Testing different percentages of β is useful for checking the stability and consistency of the algorithm. This test will be detailed in section 5.7.

Application of the inverse algorithm

The density of microwave links usually decreases from the city center to the suburbs. This can be clearly seen from the spatial structure of the network in Figure 5.2 and associated pixel density map in Figure 4.4. Such geometry of the network is useful at least for 3 reasons: (i) increasing the estimation area based on previously obtained solution without using interpolation, (ii) addressing the under-determined status of the problem (subsection 5.2.2), (iii) refining the resolution according to the network density.

To initialize the algorithm, we start the rain retrieval at a lower resolution, which should be set in advance, and to finish the retrieval process at a desired resolution. This is implemented using Grid Nesting Procedure that consists of following stages:

Stage 1. The inverse algorithm (5.3) identifies the solution vector R at a given resolution. In this stage, the apriori rainfall (R 0 and C R 0 ) is computed.

Stage 2. The inverse algorithm continues the identification of R. The retrieval resolution in this stage is two times higher than previous one which means that a number of pixels increase four times. The apriori rainfall R 0 is reinitialized by R obtained in stage 1. As a result, the corresponding standard deviation in C R 0 is updated by the variance of the reinitialized R 0 .

Stage 3. The algorithm terminates the retrieval at a desired resolution which can be set by the user. In this stage, the apriori rainfall R 0 is reinitialized by the solution vector R obtained in stage 2 and the corresponding standard deviation in C R 0 is also updated by the variance of R 0 obtained in this stage.

The reason for calling this procedure as 'grid nesting' is that the algorithm applied in each stage takes the apriori knowledge from previous stages except the first one. The only condition for implementing the Grid Nesting Procedure is that the retrieval resolutions at each stage must be known. In our study, the retrieval resolutions in stage 1, 2 and 3 are 2 × 2 km 2 , 1 × 1 km 2 and 0.5 × 0.5 km 2 , respectively. In each stage, R 0 and C R 0 are updated in the algorithm. When the resolution changes from lower to higher, the pixels are refined by the new retrieved vector over stages. Initializing the algorithm in this way is a key feature of our retrieval model.

To better understand the functioning of the Grid Nesting Procedure, a simplified example for one pixel is presented in Figure 5.3. In this example, the pixel resolutions over stages are the same as one used in our study. The description is as follows: The retrieved rain rate R 1 in stage 1 is refined in stage 2. Only crossed pixels among those 4 pixels are refined by the new solution from the lower to the higher. In stage 2, the retrieved values (R 2 , R 3 , R 4 ) are obtained for the crossed part of the area while uncrossed pixel R 1 keeps the same value from the stage 1. Similarly, each pixel at 1 × 1 km 2 can be split into 4 smaller pixels with a resolution of 0.5 × 0.5 km 2 . For example, one can see the indices of the rainfall in the apriori and solution vector. The uncrossed pixels keep the same value obtained in the previous solution.

Another important point is that if we change the resolution from the lower to the higher, the pixel sizes change and the retrieved values at lower resolution is placed into the pixels at the higher resolution. In other words, the uncrossed pixels at 1 × 1 and 0.5 × 0.5 km 2 keep the pixel values retrieved at 2 × 2 and 1 × 1 km 2 , respectively.

Rec .

Tr.

R 1 R 3 R 4 R 2 R 1 R 1 R 1 R 1 R 1 R 0 Tr. Rec . R 1 R 1 R 1 R 1 R 1 R 1 R 1 R 1 R 3 R 3 R 3 R 4 R 2 R 2 R 2 R 3 R 3 R 3 R 4 R 2 R 2 R 2 R 3 R 4 R 4 R 4 R 2 R 14 R 15 R 11 R 12 R 8
Rec .

Tr. rainfall vector.

Sensitivity analysis

The sensitivity analysis gives a clear insight into understanding an interrelation between the algorithm parameters. The importance of the sensitivity to the algorithm parameters is to test the influence of different values on the solution vector (5.3). ' In this section, the sensitivity analysis to the parameters of the algorithm is discussed. First, subsection 5.7.1 describes the sensitivity analysis protocol. Next, subsection 5.7.2 addresses the influence of the apriori parameters (i.e.

k and ∆) on the rainfall retrieval. Then, subsection 5.7.3 presents the influence of the model error magnitudes on the overall retrieval performance.

Sensitivity analysis protocol

Overall structure of the sensitivity analysis test to the algorithm parameters is given in Figure 5.4. Here, level 1 and 2 represent the values of sensitivity parameters to be tested which are apriori rainfall choice, i.e. k, and decorrelation distance, i.e. ∆, respectively. The rainfall retrieval process shown in level 3 is performed according to the following steps:

1. Rainfall map: Extract rainfall map from weather radar. The rainfall map at 0.25 × 0.25 km 2 resolution is obtained from the radar data presented in chapter 4.

Generate data:

The rain attenuation data over the microwave links network are generated using the procedures described in subsection 4.4.3 of chapter 4. Here, the default percentage of α is 5 percent given in the (4.4). The a and b coefficients at 18, 23 and 38 GHz are obtained from (ITU-R, 2005) which will also be applied in the rainfall retrieval, i.e. following step.

Rainfall retrieval:

The Grid Nesting Procedure (section 5.5) is applied to carry out rainfall retrieval.

Comparison:

The retrieved rainfall map is compared with observed radar rainfall map using NSE metric given in section 5.6.

Combinations of different values of k from 5 to 256 and decorrelation distance, ∆ from 1 km to 15 km are tested in each rain retrieval. For each pair of k and ∆, the level 3 is applied to 207 rainfall maps which represent light rain, shower, organised and unorganised storm events.

The characteristics of the rainfall data set can be found in subsection 4.3.2 (chapter 4). values which indicate that the algorithm in these choices performs best. We tested the sensitivity of the algorithm to k in each rainfall event, separately, i.e. light rain, shower, organised and unorganised storm cases. It seems that the algorithm is quite sensitive to the change in k parameter. That is to say, increase in k values decreases the accuracy of the algorithm regardless the presence of rainfall variability. However, one important feature found in these results is that the higher variability of rainfall is leading to increase ranges of optimum values of k parameter. For instance, one can see the red dots count difference between light rain and unorganised storm.

Sensitivity to the decorrelation distance

The covariance of the apriori rainfall C R 0 consists of two parameters: (i) the standard deviation σ R 0 , in mm.hour -1 , and the decorrelation distance ∆, in km, given in (5.9). The former parameter cannot be considered to be the part of the sensitivity analysis, because it is directly obtained from R 0 . However, the latter is used in the sensitivity test. The objective is to test the sensitivity of the retrieval algorithm to the parameter ∆. Selected values of ∆ are between the range of 1 km and 15 km.

The results of the analysis in Figure 5.5 show that the algorithm is also sensitive to ∆ value. It seems that the influence of ∆ increases with the variability of rainfall. In contrast to k parameter, the optimum choice of ∆ value regardless the rain type is found to be quite strict,

i.e. very few.

The optimum parameter choices of the algorithm for each rainfall type are summarized in

Influence of the model error

The objective in this subsection is to understand the influence of the model error on the performance of the algorithm. The term used for the magnitude of the measurement model error is denoted as β parameter which is the part of σ 2 model shown in (5.13). The error magnitude is expressed as a certain fraction of the generated data A 0 added to the overall variance σ 2 A of apriori covariance C A 0 (5.12).

As previously discussed in subsection 3.3.2 of chapter 3, exact values of a and b coefficients used in the model i.e. A(R) in (5.1), are approximate and unknown. The reason is that they are obtained experimentally on the basis of measured drop size distribution for certain climate conditions. Therefore, testing different percentages of β parameter is an indication to the fact that A(R) employed in the algorithm is inaccurate and can affect the stability of the retrieval algorithm. Here, solid(black), dashed(blue), dotted(green) and dotted-dashed(red) lines represent the error magnitude percentage (β) at 5, 10, 20 and 50 %, respectively. It should be noted that this sensitivity test is carried out based on the optimum parameters of the algorithm previously obtained in Table 5.1. When the NSE metric is considered, i.e. Figure 5.6, we can see that there is a slight influence of the model error on the algorithm with increasing rainfall variability, see shower event for example. However, it seems that the algorithm is stable enough showing almost identical performance in all cases, i.e.

highest (β = 50%) and lowest (β = 5%). The results based on the BIAS metric, i.e. Figure 5.7, show a bit different scenario. That is to say, a negative bias is increasing gradually with model error magnitude. Overall, the algorithm applied with optimal apriori parameters is found to be insensitive to the model error magnitudes.

Evaluation

The objective is to evaluate the capability and limitations of the proposed rainfall retrieval algorithm. First, the evaluation principle and conditions are described in subsection 5.8.1. Then, obtained results are discussed in subsection 5.8.2.

Evaluation principle

As a basis, a series of rainfall events recorded by weather radar are used to evaluate the proposed algorithm. These events consist of 207 rainfall maps grouped by light rain, shower, organised and unorganised storm events representing different spatial and temporal variability.

A detailed description about the characteristics of the rainfall data and weather radar is given in section 4.3 of chapter 4.

The evaluation of the proposed method is performed at 0.5 × 0.5 km 2 resolution based on the optimum parameters of the algorithm given in Table 5.1. First, we discuss the overall comparison between global and local apriori knowledge introduced in the algorithm. In addition, we also analyse the overall influence of the inverse algorithm and the apriori knowledge alone on the retrieval performance. This helps to understand the added value of the apriori knowledge depending on the variability of the rainfall. Second, the proposed algorithm is evaluated over the entire network in the presence of the measurement α and model β error using only two extreme cases: (i) α=5%, β=5%, (ii) α=20% and β=20%. Logically, increasing the value of α parameter should decrease the retrieval performance. Third, the retrieval efficiency of the network system depending on the link density is also evaluated. This is based on the pixel density map given in Figure 4.4 that represents the density level of the microwave links. In addition, overall retrieval performance in the densest and whole area of the network system is presented.

Evaluation results

Regarding the first aspect of our analysis, examples of retrieved rainfall for light rain, shower, organised and unorganised storm case are shown in Figure 5.8. These maps are obtained using the Level 3 which is the retrieval protocol depicted in Figure 5.4. It is clear to see that the spatial variability in all rain types is far better captured when the local apriori choice is applied (shown in the right column). The reason is that the apriori initialization procedure heavily influences the algorithm outcome as it has been proven in the sensitivity analysis test.

Therefore, the algorithm performs better in local than in global apriori initialization. However, we can see that the algorithm is smoothing the solution in all cases compared to observed maps.

Especially, this is clearly visible in shower rain type (depicted in the second row). This can be explained by the following reason: The Grid Nesting Procedure given in section 5.5 is the first cause for that because the retrieved rainfall vector in stage 1 remains as a solution vector without being updated by the new solution. Similarly, the retrieved rain vector in stage 2 remains as a solution vector in stage 3. On the one hand, the application of the Grid Nesting Procedure is decreasing the quality of the retrieval i.e. smoothing. On the other hand, it is increasing the retrieval area.

We have seen that the local apriori is significantly improving the performance of the retrieval algorithm. However, two important issues arise: (i) whether this improvement is coming only from the apriori knowledge, (ii) whether the application of the inverse algorithm adds a value for the rainfall retrieval. To address these issues, we compare the apriori knowledge with the inverse algorithm performance. Figure 5.9 illustrates the results of such comparison in NSE metric. It should be noted that the figures in the 1st, 2nd and 3rd columns are the results of 1st, 2nd and 3rd stages of Grid Nesting Procedure, respectively, see Figure 5.3.

Interestingly, the performance of the inverse algorithm is becoming almost identical as the retrieval resolution decreases from 2 × 2 km 2 to 1 × 1 km 2 . This is due to the fact that the apriori solution in the 2nd and 3rd columns do not really represent the local apriori choice. Rather, they are the results of the algorithm run in the previous stages, accepted as the apriori knowledge for the next stages. In contrast to the local apriori solution, the inverse algorithm has a significant influence on the outcome at 2 × 2 km 2 resolution. Especially, this influence becomes very large with increasing variability of rainfall, see the shower event (second row). That is to say, the more variability of the rainfall is, the larger influence of the inverse algorithm will be. with moderate density label (red line) shows unusual performance in the first few retrieved maps. The explanation to such behaviour is that the number of retrieved pixels to be compared with observed ones is much smaller than the other ones i.e. red or blue regions. One can find the number of pixels and area size in Table 4.2.

The scatter plots of the retrieved rainfall given in Figure 5.11 provide another view on the performance of the proposed algorithm. Regression lines in black and red color belong to the observed maps by weather radar and estimated maps by the algorithm, respectively. There is a systematic underestimation along all events for both whole and high density regions in the study area. However, this underestimation is decreasing in the dense regions of the network, see the figures on the right column. We can see that higher rainfall rates are always smoothed resulting in underestimation. This smoothing is explained by two reasons: (i) the solution vector retrieved in the 1st stage of Grid Nesting Procedure procedure mostly belongs to rural areas. Then, it is not updated by the algorithm in the next stages, (ii) the microwave links with longer lengths are also causing underestimation.

Summary

The objective in this chapter was to assess the feasibility of rainfall mapping by means of inverse method applied to generated signal attenuation data based on existing microwave link network. The nonlinear inverse algorithm is proposed to demonstrate the capability of microwave links for rainfall monitoring. We presented that the apriori knowledge used to initialize the algorithm heavily influences the retrieval accuracy although rainfall variability is high. A series of retrieval tests have been carried out to understand the influence of difference error sources such as measurement, model. A detailed sensitivity analysis to the apriori parameters have also been performed to test the stability of the algorithm. The algorithm applied has been found very stable even the problem to be solved is too under-determined. Various evaluation procedures have been carried out to analyse the performance of the monitoring system based on the proposed pixel density map. One weakness of the algorithm found in the evaluation is systematic underestimation especially in higher rain intensity. However, the proposed algorithm is able to handle with larger model errors.

Achieved results are naturally a subject to further works: (i) testing the algorithm in time,

(ii) a real case study instead of realistic that requires attenuation data directly obtained by cellular phone network companies, (iii) validation of the proposed approach using a large set of rainfall data and testing different network topologies.

Retrieval Model 2. Discrete Mojette

Transform Algorithm

Introduction

This chapter presents the second retrieval model that is used to assess the feasibility of rainfall mapping based on attenuation measurements from commercial microwave links. The retrieval principle is based on a tomographic approach which converts 1D attenuation data to 2D rainfall map.

Most applications based on tomographic approach (medical, non-invasive control of objects) are subject to under-determined problem which means that the number of unknown parameters (denoted as R in the previous chapter) to be estimated are much larger than the available measured data set (denoted as A in the previous chapter as well). In the context of microwave link based rainfall mapping, the introduction of tomographic approach is even more challenging due to for various reasons that have been experienced by [START_REF] Giuli | Tomographic reconstruction of rainfall fields through microwave attenuation measurements[END_REF][START_REF] Zinevich | Estimation of rainfall fields using commercial microwave communication networks of variable density[END_REF][START_REF] Cuccoli | Tomography applied to radiobase network for real time estimation of the rainfall rate fields[END_REF] as discussed in the state of the art chapter 3. First, a geometry of commercial microwave links network is unavoidably arbitrary and operating frequencies in the network are inhomogeneous. To be applicable, a tomographic reconstruction would require a homogeneous network. Secondly, the problem to be solved is nonlinear and under-determined (see, subsection 5.2.2) if we see the monitoring area globally. Therefore, we examine the applicability as well as limitations of a specific discrete tomographic approach for rainfall mapping from generated attenuation data. In particular, we address these issues using adjustment procedures applied to a real network of microwave links presented in the case study.

The retrieval model presented in this chapter uses a discrete tomographic algorithm based on the Mojette Transform (MT), developed by [START_REF] Guédon | Psychovisual Image Coding via an Exact Discrete Radon Transform Algorithm[END_REF][START_REF] Guédon | The Mojette transform: the first ten years[END_REF]. Unlike other methods proposed by [START_REF] Giuli | Tomographic reconstruction of rainfall fields through microwave attenuation measurements[END_REF][START_REF] Cuccoli | Tomography applied to radiobase network for real time estimation of the rainfall rate fields[END_REF], the MT uses completely different way of reconstructing 2D rainfall map from 1D attenuation data relying only on the geometry of the area to be reconstructed and of the microwave links network.

This chapter is organised as follows. First, section 6.2 describes a core principle of tomography. Second, section 6.3 presents a mathematical background of the applied tomographic algorithm, i.e. the MT algorithm, in this study. Then, section 6.4 addresses specific procedures used to adjust the parameters of the MT algorithm by taking into account an arbitrary geometry of the network topology, non-uniform distribution of links frequency and lengths. Finally, section 6.5 discusses the results obtained from the retrieval model. In addition, we also present a comparison between the performance of the tomographic algorithm to be explained in this chapter and the one presented in chapter 5.

What is tomography?

Tomography allows to reconstruct an interior part of an object knowing only an attenuated signal that has passed through it following different paths. Therefore, the tomographic approach is a reconstruction of images from slices. From a series of 2D images acquired through an object, using tomographic reconstruction we can get a 3D reconstruction of this object or from a series of 1D projections of an object we can reconstruct a 2D image of this object. Over the past few decades, the tomographic approach has emerged as a powerful platform in various applications, e.g. computer tomographic scanner in medical imaging, seismic and acoustic retrieval systems in geology and tools in non-destructive controls for manufactured goods as well [START_REF] Grangeat | La Tomographie[END_REF].

An object can be reconstructed from the acquisition of a finite series of its projections (Figure 6.1). For doing so, four elements of the tomography should be defined: object, bin, projection and its direction. The object is also called a region of interest that needs to be reconstructed. Rainfall map is the object in our case. The projection is the attenuation distribution at projection direction, i.e. θ angle. This is obtained by penetrating the object with a number of rays, depicted as black arrows. Each discrete value on this projection is called a bin.

Presentation of the Radon Transform

In 1917, Johann Radon defined the transform named after him [START_REF] Radon | Uber die Bestimmung von Functionen durch ihre Integralwerte langs gewisser Mannigfaltigkeiten[END_REF]. The direct transform describes the projection of an object along lines integral. This projection can be done following various acquisition geometry [START_REF] Defrise | La tomographie[END_REF] such as parallel, cone beam, fan beam, helicoïdal but we will describe and use only the parallel one here.

A 2D Radon transform projects an object (x, y) onto a 1D projection proj(t, θ) following θ projection angle. Values onto the projection depend on the integral of the values through the object along the integral line with the θ direction (Figure 6.1).

The Radon transform proj(t, θ) of function f (x, y) is given in its continuous mode by:

proj(t, θ) = +∞ -∞ +∞ -∞ f (x, y)δ(x cos θ + y sin θ -t) dxdy, (6.1)
where 0 ≤ θ ≤ π is the projection angle with the axis x and -∞ ≤ t ≤ +∞ is the curvilinear coordinate onto the projection with respect to the projection of the origin of the reference system.

FBP: Filtered Back Projection

The tomographic reconstruction is an ill-posed problem, i.e. either it may have no solution or its solution is not unique or its solution does not depend continuously on the data. In this last case, a noise onto the projection makes the inverse reconstruction unstable.

The Filtered Back Projection is one of the most used reconstruction method to retrieve an object from its projections. We will present it in the following explanation given by [START_REF] Servières | Reconstruction tomographique mojette[END_REF].

Given the Radon transform proj(t, θ) of function f (x, y) on eq. 6.1 and its Fourier transform F (λ, µ), the reconstructed object g(x, y) could be written from the inverse Fourier transform of F (λ, µ) as :

g(x, y) = +∞ -∞ +∞ -∞
F (λ, µ) exp(2iπ(λx + µy)) dxdy. (6.

2)

The Central Slice theorem gives us the 1D Fourier transform of a parallel projection of one object is equal to a slice of the 2D Fourier transform of this object:

P (ν, θ) = F (λ, µ) (6.3)
A demonstration is given in Appendix A.

For g(x, y) we obtain:

g(x, y) = +∞ -∞ +∞ -∞ P (ν, θ) exp(2iπ(λx + µy)) dxdy. (6.4)
Changing variables to use a polar coordinate system:

   λ = ν cos θ µ = ν sin θ , (6.5)
gives:

g(x, y) = Taking into account:

p(t, θ) = p(-t, θ + π) = p(t, θ + 2π), (6.7)
and splitting (6.6) on [0, π] and [π, 2π], the reconstructed object is written:

g(x, y) = π 0 +∞ -∞ P (ν, θ) exp(2iπν(x cos θ + y sin θ)) |ν| dνdθ. (6.8)
With t = x cos θ + y sin θ we get :

g(x, y) = π 0 +∞ -∞ P (ν, θ) exp(2iπνt) |ν| dνdθ. (6.9)
|ν| in this equation express a ramp filter onto the projections.

With:

Q(t, θ) = +∞ -∞ P (ν, θ) |ν| exp(2iπνt) dν, (6.10)
we get: .11) In the backprojection step, each projection is filtered by ν and then backprojected onto the image.

g(x, y) = π 0 Q(t, θ) dθ. ( 6
The equations 6.10 and 6.11 are the two steps of the filtered backprojection. In (6.10), a filtering operator is associated to each projection p θ (t) and gives a filtered projection:

pθ (t) = p θ (t) * k(t), (6.12)
with

k(t) = +∞ -∞
|ν| exp(2iπνt) dν. (6.13)

In the backprojection step, each projection is filtered and then backprojected onto the image with (6.11) (Radon dual operator): Cartesian grid and pixels it is considered as Dirac, a continuous object f (x, y) is described as:

g(x, y) = R * pθ (t) = π 0 pθ (t) dθ. ( 6
f (x, y) = +∞ k=-∞ +∞ l=-∞ f (k, l)δ(x -∆ p k)δ(y -∆ p l), (6.15)
where f (k, l) is the object value at the integer point on the grid (k, l) with:

δ(x) =    1 if x = 0 0 otherwise , (6.16)
and ∆ p is the pixel size.

Then, the Radon transform for the previously defined object is:

proj δ (t, θ) = +∞ -∞ +∞ -∞ f (x, y)δ(t + x sin θ -y cos θ) dxdy = +∞ k=-∞ +∞ l=-∞ f (k, l) +∞ -∞ +∞ -∞ δ(x -∆ p k)δ(y -∆ p l)δ(t + x sin θ -y cos θ) dxdy = +∞ k=-∞ +∞ l=-∞ f (k, l)δ(t + ∆ p k sin θ -∆ p l cos θ). (6.17)
(p, q) is the discrete projection direction with p (resp. q) the shift in pixels along x (resp. y). If several (p, q) directions give the same angle (for ex. (1, 1) and (2, 2)) the one with PGCD(p, q) = 1 is chosen [START_REF] Guédon | L'Histoire de la transformée Mojette (la forme et la projection)[END_REF].

The directions (p, q) and (-p, -q) give the same projection then only θ in [0, π[ are kept.

The discrete and continuous projections directions are linked by the following equations:

tan θ = q p , cos θ = p p 2 + q 2 and sin θ = q p 2 + q 2 , ( 6.18)
Then, the equation (6.17) is written [START_REF] Guédon | The Mojette transform: theory and applications[END_REF]:

proj δ (t, θ) = +∞ k=-∞ +∞ l=-∞ f (k, l)δ(t p 2 + q 2 + kq∆ p -pl∆ p ). (6.19)
The projection of the Cartesian grid gives a regular sampling onto the projection with respect to the projection direction (p, q):

t = b∆ p p 2 + q 2 , ( 6.20)
where b is the bin index onto the projection.

Bins number onto the projection depends on the (p, q) projection direction and on the size P × Q of the bonding box of the projected 2D object.

#bins p,q,P,Q = (P -1) |p| + (Q -1) |q| + 1. (6.21)

PGCD(p, q) = 1 means that only integer values of b are corresponding to (k, l) pixels.

Then, we get the Mojette Dirac equation: Mojette transform can also be defined for other pixel models and in higher dimension.

[M δ f ](b, p, q) = proj δ (b, p, q) (6.22) = +∞ k=-∞ +∞ l=-∞ f (k, l)∆(b + kq -pl), (6.23
Interested reader will find details in [START_REF] Servières | Reconstruction tomographique mojette[END_REF]; Guédon, 2013).

Projections directions

The discrete projection directions used to process the Mojette projections are based on the use of Farey-Haros series [START_REF] Svalbe | Farey Sequences and Discrete Radon Transform Projection Angles[END_REF].

The Farey-Haros series of order N (F N ) is the set of all fractions in lowest terms between 0 and ∞ whose denominators do not exceed N , arranged in order of magnitude. For instance, F 4

is composed of { 0 1 , 1 4 , 1 3 , 1 2 , 2 3 , 3 4 , 1 1 , 4 3 , 3 2 , 2 1 , 3 1 , 4
1 , 1 0 }. The Farey-Haros series give a set of discrete angles. With such a series, we have a set of discrete angles on the q p form arranged in order of magnitude between [0, π 2 [. Angles over [0, π[ are obtained by symmetry. Then, Farey-Haros series give the discrete angles (p, q) with irreducible fraction q p (Figure 6.3).

FBPM: Filtered Back Projection Mojette

The Filtered Back Projection Mojette is an adaptation in Mojette geometry of the FBP presented in 6.2.2. It has been defined in [START_REF] Servieres | Some links between continuous and discrete Radon transform[END_REF][START_REF] Subirats | Angular assessment of the Mojette filtered back projection[END_REF]. It is divided in two steps, a filtering step then a backprojection step.

k 0 filter

The filter used to do a FBP Mojette reconstruction is based on a Spline 0 pixel model. The continuous exact Spline 0 FBP filter was derived in [START_REF] Guédon | Bandlimited and Haar Filtered Back-Projection Reconstuction[END_REF]:

k 0 (t, θ) = 1 π sin(2θ) ln t 2 - ∆ 2 p 4 (1 + sin(2θ)) t 2 - ∆ 2 p 4 (1 -sin(2θ)) , (6.24)
for t = 0 and θ ∈]0, π 4 ] with ∆ p the pixel size, and

k 0 (t, 0) = -1 π 2∆ p 4t 2 -1 . ( 6
k 0 (b, p, q) = p 2 + q 2 2πpq ln b 2 -( p+q 2 ) 2 b 2 -( p-q 2 ) 2 , for b = 0 and (p, q) = (1, 0) (6.26)
with

k 0 (b, 1, 0) = -1 π 2s 4b 2 s 2 -1 and k 0 (0, p, q) = p 2 + q 2 πpq ln p + q p -q . ( 6.27)
As stated in [START_REF] Servières | The Mojette transform: Discrete angles for tomography[END_REF], this filter can not be implemented in a straightforward manner in the spatial domain. It exhibits discontinuities at the points where the projection of the pixels (described by a trapezoidal shape) is continuous but not differentiable. The Ramp filter acting as a derivative operator, the values at these points have to be computed using a Dirichlet condition.

Mojette Back Projection

The Mojette backprojector is defined as the dual of the Mojette projector. For a Dirac pixel model, the Mojette exact backprojector is defined as [START_REF] Servières | A Discrete Tomography Approach to PET Reconstruction[END_REF]:

[M * proj δ ](k, l) = +∞ i=-∞ +∞ j=-∞ δ(k -i)δ(l -j) +∞ b=-∞
proj δ (b, p, q)∆(b + qi -pj). (6.28)

Reconstruction of rain map using Filtered Back Projection Mojette

In order to reconstruct 2D rainfall map from the 1D vector of attenuation data we use the FBP Mojette algorithm, however the adjustment procedures are applied before implementing the algorithm.

Usually, the FBP Mojette algorithm is applicable only under the following conditions:

1. Ray frequencies are the same and constant along all the projection directions (see, subsection 6.3.2).

2. Geometry of the data acquisition system is regular and predefined as given in Fig- In order to find a compromise between 3 conditions of the FBP Mojette algorithm and those 3 existing challenges in microwave links we do adjustment procedures to the algorithm properties (i.e. the way of taking projection, projection direction, angles and interpolation of the projected bins).

Homogenizing the attenuation data set

As previously shown in chapter 4, the attenuation data generated over the network has been obtained at 18, 23 and 38 GHz using the k-R relation. However, this data set is inhomogeneous and cannot be applied directly in the tomography context. The reason is that the k-R relation coefficients obtained at these frequencies (except 38 GHz) are greater than 1 indicating that the problem type is nonlinear. Hence, the problem is nonlinearity and inhomogeneity of the measured data set which is not compatible with above defined applicability of the tomography conditions. To overcome this issue, the data transformation procedure is applied in order to homogenize the generated attenuation data. The purpose of this procedure is to express the respectively. If we consider that the average rain rate, R, along this link is the same for both frequencies we can express the equation (4.3) as follows:

A 1 = L * a 1 * R b 1 (6.29) A 2 = L * a 2 * R b 2 (6.30)
Defining R from both equations leads a relation: Using (6.31), we transform all the generated attenuation data at 18 and 38 GHz into the equivalent one at 23 GHz. However, to make sure if the transformation procedure is 'consistent' enough we test the scatter plot between transformed and original values in noise free conditions.

A 2 = L * a 2 * ( A 1 L * a 1 ) b2 b1
This helps to understand whether the transformed data set at 23 GHz can be equivalent to one that has been generated at 18, 23 and 38 GHz. Figure 6.4 shows such consistency test performed in four rain events. It turns out that the correlation between those variables is almost 1 which justifies that the transformation approach is reasonable. Based on this evidence, we can assume that the transformed data is equivalent to the generated one in chapter 4. In the second step, we divided the study region into 'reconstructible' sub-grids so that each region would be reconstructed independently with respect to the attenuation measured through it. Only zones totally crossed by the microwave link were chosen for the reconstruction process.

In Figure 6.5, the 'reconstructible' sub-grids are depicted as black rectangles with different sizes over the study area. In our case, we selected 21 valid sub-grids for the reconstruction process. However, the number of sub-grids are reduced to only one sub-grid (i.e. sub-grid 11), in the evaluation of the tomographic algorithm. So, two areas i.e. whole area and sub-grid 11 is evaluated to see the algorithm performance difference between those zones. It is worth mentioning that the number of sub-grids depend on the crossed number of links in the study 114CHAPTER 6. RETRIEVAL MODEL 2. DISCRETE MOJETTE TRANSFORM ALGORITHM area. We consider the sub-grid region to be 'reconstructible' if at least two links cross the chosen sub-grid and this was performed manually.

Taking projections from the microwave link

A set of projections M p,q with discrete Mojette directions were found for the MT acquisition in each grid. We have chosen discrete (p, q) projections directions fitting the best microwave link directions and minimizing the number of bins. These possible directions depend on the 'reconstructible' sub-grid size. However, only very few bins onto each projection were filled. The Figure 6.6 shows an example of taking projections for one link (red solid line in the black rectangle) from the 4 th sub-grid at two projection directions P roj 1 , P roj 2 . This sub-grid is the one of the chosen sub-grids in Figure 6.5 Below, the steps for taking projections from one link are given:

1. Define projection directions P roj(p, q) using Farey-Haros series (subsection 6.3.2);

2. Find all possible projections P roj(p, q) that are nearly orthogonal to the link L: P = {P p,q ⊥ L} > A * , where A * is the orthogonality level;

3. Choose only P with minimum number of bins in the set of the projections (P p,q ): P p,q = argmin{N bin (P p,q )} where, N bin -number of bins in (p, q) projection. The following steps are computed to perform the angular interpolation technique on to the empty bins:

1. Find the central point (C x , C y ) of the projection area, i.e. microwave link area. Then, it is possible to assume that there is a line (denoted as D) parallel to the projection (p3, q3) which passes through the central point (C x , C y ). The equation of this line is {(C x -q 3 × t; C y + p 3 × t)|t ∈ }.

2. Take projection of bin b 3 on A point: b 3 = -q 3 × C x + q 3 2 × t + p 3 × C y + p 3 2 × t (6.32)

3. Project A point on (p 1 , q 1 ) and (p 2 , q 2 ) and compute the crossed points. Here, the ith 116CHAPTER 6. RETRIEVAL MODEL 2. DISCRETE MOJETTE TRANSFORM ALGORITHM projected coordinate is computed as follows: t i = -q i × C x + q i × q 3 × t + p i × C y + p i × p 3 × t (6.33)

In our example, two points t 1 and t 2 are computed using the formula above.

4. Interpolate the b 3 using the crossed points t i on each projection. This is implemented using the values of t 1 and t 2 and the corresponding angles α 1 and α 2 :

b 3 = α 1 × t 1 + α 2 × t 2 α 1 + α 2 (6.34)
A detailed explanation about the interpolation method can be found in [START_REF] Servières | Interpolation method for the Mojette transform[END_REF]. Using these four steps we interpolate all the empty bins based on the neighbouring bin values. This process is repeated for all the selected sub-grids in the study area.

Summary of the reconstruction protocol

The reconstruction of rainfall maps is performed as follows:

1. Generate the attenuation data using the protocol given in subsection 4.4.3;

2. Transform the generated attenuation data at 23 GHz using the Equation 6.31;

3. Apply the adjustment procedures:

-Choose valid sub-grids at 2 km × 2 km resolution (subsubsection 6.4.2.1); -Get projection (subsubsection 6.4.2.2); -Interpolate (subsubsection 6.4.2.3); 4. Reconstruct the attenuation map using the FBP Mojette algorithm (subsection 6.3.3).

The reconstructed images in all sub-grids are assembled then normalized with the maximum attenuation value in each sub-grid; The reconstruction of rainfall map can be performed using these 6 steps. It should be noted that only some parts of the radar rainfall map (shown in black rectangles in Figure 6.8)

are compared with reconstructed rainfall maps since the algorithm is able to retrieve rainfall sults of the test give better understanding of the reconstruction capability as well as limitations of the MT algorithm over different geometry of the network.

Surprisingly, the densest regions are heavily decreasing the accuracy of the MT algorithm showing RMSE error a bit less than 600 mm.h -1 for i.e. sub-grid 11 and a bit more than 1000 mm.h -1 for sub-grid 12. This result is coming from overfeeding too many empty bins onto projections in the interpolation procedure. Very long links in sub-grid 12 can be a driving source for that because they produce large attenuation values, especially in high rain rates. Therefore, the overestimation grows with increasing rain intensity as has been shown in Figure 6.9. The same pattern applies to the sub-grid 11. Although, the link lengths in this region are very short it may lead to overfeeding the empty bins during the interpolation procedure.

Interestingly, the MT algorithm performs better in very sparse regions of the network which have a few links and the density is very low. We can see that the lowest RMSE error over all events is in sub-grid 19 which is the second sparsest among the considered regions, for example, showing less than 130 mm.h -1 in organised storm event. Next sparse region, i.e. sub-grid 15, has almost identical trend in comparison with sub-grid 19. This result also confirms that the interpolation plays a vital role in the reconstruction accuracy of the MT algorithm.

Overall, the capability of tomographic approach based on the MT algorithm clearly represents a poor performance due to the following reasons: (i) The accuracy of the tomographic algorithm used in this chapter strongly depends on a number of projections which are not sufficient along the microwave links, (ii) The 'reconstructible' sub-grid choice which was manually fixed in this study.

Summary

The objective in this chapter was to demonstrate the feasibility of rainfall mapping by means of tomographic reconstruction applied to signal attenuation data generated based on real existing microwave links network. The attempt in this present chapter was to explore the capability as well as limitations of the tomographic approach on the basis of the Mojette Transform algorithm. A series of reconstruction tests have been carried out in 4 rain events with strong spatial and temporal variability. It is worth noting that the transform-based tomographic algorithm was applied to inhomogeneous microwave links network topology, length and frequency for the first time. The proposed algorithm seems to be quite stable for capturing a spatial variability of rainfall over less dense network. 122CHAPTER 6. RETRIEVAL MODEL 2. DISCRETE MOJETTE TRANSFORM ALGORITHM

In this chapter, we have addressed three important issues faced in the context of applicability of the tomography approach to rainfall retrieval by microwave links. These issues are related to (i) adjustment procedures of the tomographic algorithm (ii) Inhomogeneous link frequency and lengths, (iii) Non-linearity of the reconstruction problem. However, a further study with more focus on the following issues is suggested: (i) improvement on the 'reconstructible' sub-grid choice, (ii) improvement on the way of taking projection, (iii) A sensitivity analysis to the application conditions such as 'orthogonality level' and quantification of associated error sources of the model parameters needs to be investigated. Further, a definition of a priori knowledge about the state of the problem needs to be sought since the rainfall retrieval problem to be solved is highly ill-conditioned in the tomography context.

Conclusions and Perspectives

Conclusions

The main objective of this research was to assess the feasibility of retrieving rainfall maps using signal attenuation data from commercial microwave links that are operated by cellular network companies. For that purpose, we adopted a simulation framework applied to a real case study.

We recall two challenges stated at the beginning of the study: i) To generate rain attenuation data along microwave links, ii) To retrieve rainfall map using the generated rain attenuation data.

We addressed the first challenge using three chapters i.e. chapter 2, chapter 3 and chapter 4.

In chapter 2, we presented a theoretical background on rain attenuation along the microwave links. In particular, we showed steps for computing electromagnetic signal scattering as well as absorption by rain drop on the basis of Mie scattering theory and drop size distribution models. Then, these steps were used to demonstrate the influence of drop size on microwave signal at frequencies at 18, 23 and 38 GHz used by microwave links of cellular network in this study. We also demonstrated a dependence of the specific rain attenuation on frequency and rain intensity.

In chapter 3, the state of the art in rainfall measurement using microwave link(s) was reviewed in detail. In this chapter, we analysed all existing problems and scientific progress regarding the applications of microwave links in rainfall measurement which were grouped by three most important issues. The conclusions drawn from these issues are summarised as follows:

Issue 1. Measurement errors sources that might occur along a single or multiple links:

Reviewed papers were classified by two types of errors: instrumental (i.e. quantization, link length, frequency, antenna wetting) and environmental(drop size distribution, temperature, rainfall variability, gas and water vapour effects). The rain attenuation data presented in chapter 4 was generated on the basis of general conclusions drawn on these error types. In chapter 4, the objective was to present a case study that uses data sets i.e. weather radar rainfall maps and cellular network. In addition, we also provided a procedure for generating rain attenuation data based on these given data sets. For this, we applied the k-R relation (given in chapter 2) to generate rain attenuation data along each microwave link of the cellular network. The measurement error sources caused by drop size distribution variability and signal quantization were simulated based on the conclusions drawn in chapter 3. Obtained results in this chapter were used to represent realistic signal attenuation data that could be obtained from microwave antenna stations of cellular phone companies.

To address the second challenge of the simulation framework, we developed two retrieval algorithms in chapter 5 and chapter 6 both of which were used to reconstruct rainfall map from the rain attenuation data generated in chapter 4. The conclusions drawn from these chapters are as follows:

In chapter 5, we proposed the first retrieval algorithm which uses the principle of inverse problems. The retrieved rainfall maps by this algorithm were compared with weather radar rainfall data. The conclusion based on more than 200 retrieval tests shows that the present algorithm is capable of capturing spatial variability of rainfall even in a large magnitude of measurement (i.e. upto 20 %) and model errors (i.e. upto 50%). However, one limitation was found that the algorithm systematically underestimates the rainfall especially in higher rain intensity.

In chapter 6, we proposed the second retrieval algorithm which is based on the principle of tomography. The attempt in this chapter was more focused on exploring capabilities as well as limitations of the tomographic approach. We addressed three important issues which were found challenging in the application of tomographic algorithms for the main objective of this research: i) adjustment procedures on the tomographic algorithm parameters ii) Inhomogeneous link frequency and lengths, iii) Nonlinearity of the retrieval problem. The proposed algorithm was tested applying the same rainfall data as the one used in chapter 5. It was found that the retrieval performance is not as expected resulting very large errors compared to the first algorithm.

Perspectives

This research opens several new directions for the future study to better understand rainfall retrieval capability of commercial microwave links.

Regarding our simulation framework, a further research is required to improve the quality of generated rain attenuation data in order to reflect a real nature of signal data that can be measured at microwave antenna stations. Measurement error sources added in rain attenuation generation (see chapter 4) were simulated based on a simplified assumption i.e random variable with variance equal to 5% and 20% percent of the total attenuation using Gaussian normal distribution. This could be improved by formulating a new structure which can include link characteristics (essentially frequency and length) and climate conditions such as temperature, humidity, gas, water vapour effects (see chapter 3). This error structure can be developed on the basis of a methodology for modelling error sources and uncertainty proposed by [START_REF] Zinevich | Prediction of rainfall intensity measurement errors using commercial microwave communication links[END_REF].

Regarding the inverse algorithm presented in chapter 5, a further study could be to assess the long-term effects of the cellular network for rainfall monitoring based on longer data set, e.g. one week instantaneous rainfall provided by weather radar. In addition, model uncertainty structure introduced in the algorithm could be changed using the same recommendation given in the previous paragraph. This helps to get a deeper understanding of capabilities as well as

Bahtiyor ZOHIDOV

Reconstitution des champs de pluvieux en utilisant des mesures d'atténuation sur les liaisons hyper-fréquences commerciales. Une étude de modélisation de faisabilité.

Retrieval of rainfall fields using signal attenuation measurements from commercial microwave links. A modeling feasibility study.

Mots-clés

La pluie, Liaison Hyper-Fréquence, Atténuation, Méthode d'Inverse, Tomographie.

Abstract

The objective in this study is to assess the feasibility of retrieval of rainfall fields in urban areas using commercial microwave links, with a special focus given on cellular network, in a simulation framework applied to a realistic case study. The

Key Words

Rainfall, Microwave link, Attenuation, Inverse method, Tomography.

Figure 1 . 1 -

 11 Figure 1.1 -Rain attenuation between Transmitter and Receiver station of a microwave link. (This figure was obtained from the source 1 ).

Figure 2 . 1 -

 21 Figure 2.1 -Electromagnetic scattering in a spherical particle. Modified version of figure by[START_REF] Hahn | Light scattering theory[END_REF]

  where, notations a n and b n are the complex Mie coefficients which are obtained from matching the boundary conditions at the surface of the sphere. They are expressed in terms of spherical Bessel functions of the first and second kind (shown in Figure 2.2 and Figure 2.3, respectively)

Figure 2 . 2 -

 22 Figure 2.2 -Spherical Bessel functions of the first kind.

Figure 2 . 3 -

 23 Figure 2.3 -Spherical Bessel functions of the second kind.

Figure 2 .

 2 4 illustrates the electromagnetic signal scattering dependence on different sizes of rain drop that has been calculated for frequencies ranging from 10 GHz to 100 GHz. The figure on the left is the real part whereas the figure on the right is the imaginary part of the forward scattering function. The imaginary part is directly proportional to the microwave signal attenuation.

Figure 2 . 4 -

 24 Figure 2.4 -Real (on the left) and imaginary (on the right) parts of the complex forward scattering functions depending on the frequency range.

2 . 3 . 4 .

 234 Drop radius: 0 to 0.35 cm; Temperature at 0 o C; Drop shape is spherical.

Figure 2 .

 2 Figure 2.5 shows the relation between extinction, scattering and absorption with dependence on the rain drop radius size at 18, 23, 38 GHz. As can be seen that the extinction cross section values increase as the frequency changes over ascending order with dependence on the drop radius range. Both extinction and absorption efficiencies are similarly increasing according to the frequency ranges. The drop radius below about 0.175 cm gradually influences on scattering amplitude which indicates that attenuation is considerably large in average-sized rain drops rather than large rain drops, especially, in between 0.2 cm and 0.35 mm. According to the Mie cross section calculations it can be seen that maximum extinction efficiency is observed when the drop radius reaches the values between 0.2 cm and 0.3 cm at 38 GHz. Comparing values of the extinction cross section in Figure2.5 (figures on the right column), we can see that

Figure 2 . 5 -

 25 Figure 2.5 -Mie Efficiency coefficients (on the left) and Cross Sections (on the right). Extinction, scattering and Absorption at 18, 23, 38 GHz (top-bottom order).

Figure 2 . 6 .

 26 Figure 2.6. (Pruppacher and Pitter, 1971) investigated the shape of raindrop in the framework of modelling. They analysed variation of drop deformation with drop size. A diameter of the drop can only exist not greater than around 7 mm (Laws and Parsons, 1943; Pruppacher and Beard, 1970; Pruppacher and Pitter, 1971). Extremely large drops are often split into new small sized drops during falling towards the ground.

  DSD is typically described using a distribution function N (D) that provides N (D)dD the mean number of drops per unit of air volume with diameters between D and D + dD. So, let us see the formulation and characteristics of these two (Gamma and Marshal Palmer) models in detail in the next paragraph.

Figure 2 . 6 -

 26 Figure 2.6 -A) Raindrops are not tear-shaped; B) Very small raindrops are almost spherical in shape; C) Larger raindrops become flattened at the bottom, like that of a hamburger bun, due to air resistance; D) Large raindrops have a large amount of air resistance, which makes them begin to become unstable; E) Very large raindrops split into smaller raindrops due to air resistance 3 .

Figure 2 . 8 -

 28 Figure 2.8 -Exponential DSD function at 5 mm.hour -1 (Table2.2).

Figure 2 . 9 -

 29 Figure 2.9 -Specific rain attenuation as a function of frequency between 1-100 GHz.

Figure 2 . 10 -

 210 Figure 2.10 -Specific attenuation plot at 18, 23 and 38 GHz using MP and Gamma DSD at temperature 0 o C. Note that the parameters for plotting MP and Gamma DSD were derived from Table2.2 and the one demonstrated by[START_REF] Zhang | Project IEEE 802.16 Broadband Wireless Access Working Group[END_REF], respectively.

3. 2 3 . 2 . 1

 2321 Recent advances in microwave link based rainfall measurement Rainfall measurement along a microwave link The beginning of the current research goes back to late 1960s. Over the past four decades there has been an extensive research in microwave link based systems (telecommunication, cellular network, broadcasting, terrestrial and satellite) that focused on understanding the relationship between signal attenuation and rainfall intensity in different climate regions of the world (Hogg, 1968; Atlas and Ulbrich, 1977; Olsen et al., 1978; Jameson, 1991; Zhang and Moayeri, 1999). In the same line of these investigations, Das et al., 2010 established the rain attenuation model at frequencies from 10 to 100 GHz in the case of tropical India. These achievements improved the understanding of the rain attenuation variability and designing the microwave link frequency ranges depending on the climate. However, these findings were only focused on establishing a quality communication in different regions of the world.

42CHAPTER 3 .

 3 STATE OF THE ART: RAINFALL MEASUREMENT USING MICROWAVE LINKS data set of drop size distribution (about one year) in the case of the Netherlands. Thus, the link estimates were found to be almost identical to those measured by rain gauges. The true path averaged rainfall intensity along the link was obtained by the records of 7 rain gauges which were closely placed along the link. A related point to consider is that the authors give more emphasis on steps of pre-processing the signal data by taking into account various attenuation effects which are different from rain. These are related to antenna wetting after the rain event or signal attenuation before the rain event[START_REF] Minda | High temporal resolution path-average rain gauge with 50-GHz band microwave[END_REF]; Leijnse et al., 2008; Leijnse et al., 2010; Zinevich et al., 2010). (Leijnse et al., 2007a) suggested that using the same link it is possible to measure evaporation as well. However, the studies discussed above account for only research dedicated links. In an operational framework, (Messer et al., 2006) demonstrated the feasibility of rainfall monitoring on the basis of commercial microwave links employed by digital fixed radio systems. The measurement process is the same as previously mentioned studies, but using the signal attenuation data recorded by cellular networks. The rainfall intensity estimates by cellular-link were found consistent with those obtained by rain gauge and weather radar. Several experiments in very different contexts have confirmed those achievements. For example, (Leijnse et al., 2007c) conducted a similar experiment applied to the case study in the Netherlands territory. In their experimental set-up, two cellular links, both operating at 38 GHz, were used. These links consist of the same receiver station at which the signal levels are recorded. The objective in this experiment was to measure the path average rainfall intensity along the link based on the signal attenuation data obtained by those two links. For this purpose, the authors used k-R relation which relates a rainfall intensity (R, mm per hour) to a specific attenuation (k, dB per km), see the details about the k-R relation in chapter 2. Drop size distribution

2 . 2 .(

 22 [START_REF] Doumounia | Rainfall monitoring based on microwave links from cellular telecommunication networks: First results from a West African test bed[END_REF] obtained encouraging results in Sahelian West Africa test bed using 29 km long cellular microwave link at 7 GHz. Rain attenuation data was obtained using the difference between transmitted and received signal levels which have been recorded every second with a power resolution of 1 dB during monsoon period. Authors employed the k-R relation to convert signal attenuation data to the path average rainfall intensity along the link following the same principle used in(Leijnse et al., 2007c). Since only rain induced attenuation data is useful for rainfall measurement, the method proposed by (Schleiss and Berne, 2010) was applied to remove attenuation in dry period that does not belong to rain. The measured rainfall intensity by the cellular link shows a good performance when compared with rain gauge data at 5 min time interval; however, measuring rainfall in this time step along longer links is not representative due to the high variability of rainfall. Interestingly, it is worth mentioning that even though low frequencies are considered to be less sensitive to rainfall, the link tested in this experiment is still capable of measuring rainfall rate in real time. Based on the results, the authors state that the measurement by cellular links are more in accordance to the result obtained by rain gauge data than the ones estimated by satellite.Recent experiments[START_REF] Overeem | Measuring urban rainfall using microwave links from commercial cellular communication networks[END_REF][START_REF] Rayitsfeld | Comparison of two methodologies for long term rainfall monitoring using a commercial microwave communication system[END_REF] which group longer data sets in terms of time length recorded by larger number of microwave links and validation procedures based on rainfall measurement devices (rain gauges, disdrometers) provide robust assessments of these initial findings.

  al., 2000; Rahimi et al., 2003; Upton et al., 2005; Minda and Nakamura, 2005; Leijnse et al., 2008; Fenicia et al., 2012) and commercial (Leijnse et al., 2007c; Schleiss and Berne, 2010; Zinevich et al., 2010; Overeem et al., 2011) microwave links.Error sources due to gas and water vapour effects are often called attenuation baseline. The influence of this error depends on atmospheric conditions such as air temperature and humidity level in a particular region.(Leijnse et al., 2007b) state that the microwave link can estimate path-averaged rainfall at a high resolution if the attenuation baseline of the microwave link 46CHAPTER 3. STATE OF THE ART: RAINFALL MEASUREMENT USING MICROWAVE LINKS signal is perfectly known. The main challenge is to identify and remove this effect from the measured signal attenuation along the link. It is obvious that the attenuation baseline does not represent the rain induced attenuation. A simple procedure for identifying this error source is to measure the signal attenuation along the link before and after the rain event which are referred to as dry and wet periods, respectively[START_REF] Rahimi | Use of dual-frequency microwave links for measuring path-averaged rainfall[END_REF][START_REF] Minda | High temporal resolution path-average rain gauge with 50-GHz band microwave[END_REF][START_REF] Upton | Microwave links: The future for urban rainfall measurement?[END_REF]. The attenuation baseline could be used as a threshold to separate the total rain attenuation measured in the latter period. Another way is to employ a weather radar or a network of rain gauges closely located along the link so that dry/wet periods can be detected(Upton et al., 2005; Minda and Nakamura, 2005; Cummings et al., 2009; Overeem et al., 2011; Fenicia et al., 2012). It should be noted that the procedure to estimate the attenuation baseline proposed in all studies differ from each other in terms of underlying assumptions about the nature of this error source. For example, the methods developed by (Upton et al., 2005; Schleiss and Berne, 2010) are based on the assumption of a constant baseline attenuation that the signal level in dry period is constant in all rainy events. However, (Fenicia et al., 2012) argued that the assumption of a constant baseline attenuation seems to be unrealistic due to the signal variability in realtime conditions. The authors proposed two baseline estimation methods, namely constant and variable. Comparisons between those models showed that the variable baseline attenuation has a better performance to detect the wet/dry periods. According to (Zinevich et al., 2010) the attenuation baseline is the major source of error compared to other effects such as quantization and drop size distribution variability. A remaining challenge in identifying the baseline errors is due to low-intensity wet periods in which the signal level is almost identical to one in dry periods. Another source of error is the attenuation due to antenna wetting which is also called wet antenna. The attenuation due to wet antenna arises from the remaining water films on the antenna radome after (or during) the rain event. Previous studies (Minda and Nakamura, 2005; Leijnse et al., 2007b; Leijnse et al., 2007c; Leijnse et al., 2008) applied correction methods for removing the attenuation caused by antenna wetting. (Minda and Nakamura, 2005) corrected the attenuation due to antenna wetting including two other effects, namely temperature and water vapour. The wet antenna correction procedure was applied during as well as after the rain event. (Leijnse et al., 2008) quantified the errors caused by antenna wetting using an experimental set-up located in the Netherlands. In their study, the signal attenuation data by (Leijnse et al., 2007b) were used to investigate the error source caused by antenna wetting.The procedure for correcting the wet antenna attenuation was computed by subtracting the expected path-integrated attenuation based on the rain gauges measurements from the measured total attenuation. The authors demonstrated the relation between wet antenna attenuation and frequency for different rain intensities for a single antenna. They found that increasing rain rate is responsible for higher wet antenna attenuation as the frequency increases from 5 to 50 GHz.The influence of the wet antenna attenuation is very large for shorter links with low frequency compared to the measured signal attenuation along the path. However, it was shown that the wet antenna attenuation does not depend on frequency in the range of 17 and 23 GHz. The authors state that the effects of wet antenna and quantization errors can be reduced if a correct temporal sampling strategy is selected.[START_REF] Zinevich | Prediction of rainfall intensity measurement errors using commercial microwave communication links[END_REF] carried out a comprehensive investigation in which the errors and uncertainties due to the wet antenna, baseline, quantization, drop size distribution variability can be modelled. They applied 21 vertically polarized links at 18-23 GHz with maximum link length of 7.26 km in order to independently quantify the impact of each error source on the rainfall measurement accuracy. Their conclusion regarding the wet antenna attenuation is that the error can increase with the link length because of the increased spatial variability.

(

  Leijnse et al., 2007c) found that a and b coefficients of the k-R relation are relatively insensitive to the type of rainfall and the drop size distribution for the Netherlands case. However, the conclusion drawn in this experiment is based on the: (i) microwave links at 38 GHz (ii) the k-R relation coefficients was computed using drop size distribution data measured by[START_REF] Wessels | Metingen van regendruppels te De Bilt[END_REF]) (iii) signal attenuation is caused by spherical rain drops.(Leijnse et al., 2007b) investigated the effect of temperature on the measurement model.The k-R relation for different drop size distributions was found to be linear at 27 GHz. It was reported that the the k-R relation coefficients are not influenced by the temperature variations between -10 to 40 o C. In addition, this effect was also found to be minor for single frequency links above 10 GHz. (Leijnse et al., 2008) quantified the effects of the link frequency and length, non-linearity of the R-k relation, different sampling strategies, and two other factors (quantization and wet antenna) which have been discussed in subsection 3.3.1. The microwave signal attenuation along the links lengths from 0.24 to 10.08 km with frequencies from 5 to 100 GHz was simulated using the rain rate maps obtained from X-band weather radar data recorded over more than 1.5 years. These findings indicate that the errors and uncertainties related to a nonlinearity of the R-k relation increase with link lengths at lower frequencies, e.g. below 10 GHz. It seems possible that these results are due to the spatial variability of rain rate along the links with long lengths. (Leijnse et al., 2010) examined the errors and uncertainties caused by the spatial variability of rainfall as well as drop size distribution as a function of link length and frequency. The simulated microwave link signals using the drop size distribution data sets recorded over a long period (nearly 2 years) were applied to estimate the path averaged rain rate along the link. The authors found that, for single-frequency with single polarization links, the errors and uncertainties associated with spatial variability of drop size distribution are limited (mean bias error 10% and uncertainties of 50%) at frequencies between 20 and 40 with lengths from 2 to 10 km. The least amount of errors and uncertainties were found at frequencies around 30 GHz where the linearity of the k-R relation is dominant. 50CHAPTER 3. STATE OF THE ART: RAINFALL MEASUREMENT USING MICROWAVE LINKS (Zinevich et al., 2010) modelled all possible error sources, i.e. the instrumental and environmental (i.e. related to the measurement model) that might occur during the rainfall measurement along microwave links. Although this study has been discussed in subsection 3.3.1, in terms of the measurement model related error source, the authors reported that the errors caused by drop size distribution variability also increases because of lowering frequency band from 22

(

  [START_REF] Zinevich | Estimation of rainfall fields using commercial microwave communication networks of variable density[END_REF]) applied a tomography algorithm based on the SIRTsimultaneous iterative reconstruction technique described by[START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF] to map rainfall fields.The algorithm was tested in a simulation framework over a network of 249 real commercial microwave links of lengths 0.5-27 km, operating at 8-38 GHz during 2 hour-period of intense rain. The authors addressed the important issues raised by[START_REF] Giuli | Tomographic reconstruction of rainfall fields through microwave attenuation measurements[END_REF] such as the nonlinearity of the k-R relation and irregular geometry of the network in terms of link lengths and location. However, the proposed algorithm does not take into account of the modelling error sources that could arise from the k-R relation in the reconstruction process which was a subject to a further investigation.[START_REF] Zinevich | Frontal rainfall observation by a commercial microwave communication network[END_REF] combined a space-time advection model of rainfall and nonlinear extended Kalman filter[START_REF] Welch | An introduction to the kalman filter[END_REF] to perform 52CHAPTER 3. STATE OF THE ART: RAINFALL MEASUREMENT USING MICROWAVE LINKS the rainfall mapping from path integrated attenuation measurements. The applied model takes the advantages of the motion of the rainfall field in space and time.[START_REF] Goldshtein | Rain rate estimation using measurements from commercial telecommunications links[END_REF]

(

  [START_REF] Gosset | Improving Rainfall Measurement in gauge poor regions thanks to mobile telecommunication networks[END_REF] emphasize the potential interest of radio links for mapping rainfall in areas with a poor or non-existent coverage of classical (rain gauge and weather radar) measurement devices including urban areas.

Figure 4 . 1 )

 41 . This city has a temperate and oceanic climate with no significantly variable topography. Both rain attenuation measurement and rainfall retrieval experiments are conducted in the chosen part of the study area with a size of 40 × 40 km 2 depicted in yellow square which corresponds to Nantes city center with a large number of links (Figure 4.1). This city is well-equipped with commercial microwave links that are used for broadcasting, wireless internet service, and cellular phone communication systems. We selected only cellular phone network which groups 52 microwave links at 18 GHz, 102 at 23 GHz and 103 at 38 GHz (256 microwave links in total) with both horizontal and vertical polarization. Their spatial coverage is almost 1368 km 2 area.

Figure 4 . 1 -

 41 Figure 4.1 -Cellular phone network located in Nantes city. The zoomed-in part of the figure (yellow square) represents the location of microwave links employed cellular network operating at 18, 23 and 38 GHz. Two yellow arrows are also depicted as an example to indicate the transmitter stations. The chosen area size is 40 × 40 km 2 .

Figure 4 . 2

 42 which shows three histograms for the links at 18, 23 and 38 GHz separately. It turns out that almost 95 % of the links at 38 GHz are shorter than 5 km while the same percentage shows about 11 and 10 km for the links at 18 and 23 GHz, respectively. This is a clear indication to the fact that microwave link lengths at higher frequencies are most likely expected to be shorter and vice versa.

Figure 4 . 2 -

 42 Figure 4.2 -Histograms of microwave links length at 18, 23 and 38 GHz. In the figure, 'All links' represents a histogram for the lengths of all the microwave links regardless of their frequencies.

Figure 4 .

 4 4 as an example.

  e. D vector. A greater value of D indicates that a pixel is more informative and well-sampled by the link, otherwise, less informative. The closer the value of D to zero, the lower density level is in the network.

Figure 4 . 3 -

 43 Figure 4.3 -Histogram of the pixel density map.

Figure 4 .

 4 Figure 4.3 shows a histogram of the pixel density computed using the microwave links. In this figure, the distribution of D is highly right-skewed indicating that the average value of D is

Figure 4 .

 4 Figure 4.4 shows the pixel density map obtained at 2 × 2 km 2 resolution over the entire network. The classified density levels are depicted in grey, red and blue color represent low, moderate and high density regions of the network, respectively. It should be noted that these

Figure 4 . 4 -

 44 Figure 4.4 -Pixel density map at 2 × 2km 2 resolution. The colors blue, red, grey represent the high, moderate, low density areas, respectively.

4. 3

 3 Presentation of weather radar rain events 4.3.1 Weather radar Rainfall maps used in this study were collected between the year of 2009 and 2012 by Cband weather radar of Treillières, operated by Météo-France. The weather radar is located about 10 km north of the center of Nantes city. In Figure 4.1, its location is marked as C-band. The characteristics of this radar are given in Table4.3. The measured reflectivity profiles in a polar coordinate system are projected every 5 min onto a 128 km × 128 km Cartesian grid with a spatial resolution of 0.25× 0.25 km 2 . Then, the instantaneous reflectivity profiles are converted into rainfall maps using the Z -R relationship by Marshall-Palmer (i.e. Z = 200 R1.6 , with Z and R being the radar reflectivity in mm 6 .m -3 and the rainfall rate in mm.h -1 , respectively).

Figure 4 . 5 -

 45 Figure 4.5 -Example of weather radar rainfall maps with a spatial (0.25×0.25 km 2 ) and temporal (5 min) resolutions.

4. 4 . 1 . 1

 411 Main assumptionsLet our study area be described in two dimensional (2D) discrete space, Figure4.6. Then, it is possible to assume that each pixel has a constant rainfall rate, R 1 , ..., R 6 in mm.hour -1 (see, the part of the figure on the right). In this 2D discrete space, the microwave link propagates a signal only in the form of straight line with length of L km, between transmitter and receiver, and passes through some parts of the discrete area (R 1 , R 2 , R 3 , R 4 ) as shown in Figure4.6. Each pixel represents the value of constant rainfall rate. Besides, the impact of antenna altitudes and signal direction (either from transmitter or receiver) on the signal transmission trajectory are neglected since it is only valid in three dimensional space. Based on such assumptions, the problem of generating rain attenuation data along the link is formulated under the following conditions:1. A transmitted signal is attenuated by rain rate pixels that are only crossed by the microwave link and these pixel values are constant, see the zoomed-in part of Figure4.6.2. Each pixel size in the study area is equal to the resolution of rainfall map by weather radar and this is 0.25 km × 0.25 km.3. Attenuation measurement model and its parameters are perfectly known and it can be expressed as k -R relation, see the formula (2.26) given in chapter 2.

4 .

 4 Rain attenuation measurement along the link is affected by different sources of errors which will be discussed in subsection 4.4.2.

Figure 4 . 6 -

 46 Figure 4.6 -Example of discretizing the study area. The figure on the left represents the discretized area over the cellular networks at 18, 23 and 38 GHz. The figure on the right shows an example of intersected pixels by one link with length L.

- 1 ;

 1 l ji -the length of the crossed part of the jth link in ith pixel, km. As discussed in chapter 2, the a and b coefficients are constant values which depend on signal frequency, polarization, climate temperature and drop size distribution. In this study, these coefficients were obtained by the recommendations of International Telecommunication Union (ITU-R, 2005) for three frequencies, namely 18, 23 and 38 GHz. The term ε is the measurement error source which we will define in subsection 4.4.2. Applying the Equation4.3 the rain attenuation data over the cellular network can be generated.4.4.2 Measurement error sourcesSources of error in rain attenuation measurement process are related to k -R relation and instrumental impairments. Based on the conclusions from the literature review discussed in chapter 3, two major sources of error, namely the influence of drop size distribution and signal quantization are added to compute path integrated rain attenuation data along the link in Equation4.3. Since the quantifying the error sources precisely is beyond the scope of this study, we consider a simple approach in order to simulate the nature of the real signal attenuation caused by rain.

Figure 4 .

 4 7 shows an example of generated rain attenuation data for 6 chosen links at 18, 23 and 38 GHz with different lengths located at different places of the study area. The x and y axes represent time slots of each event, e.g. 113 in light rain, and path integrated attenuation computed along the given link, respectively. The characteristics of these chosen links are given in

Figure 4 .

 4 7 belong to shorter (∼ 2 km) and longer (∼ 10 km) links, respectively. Expectedly, it turns out that an overall trend of the attenuation data measured along the longer links seems to be quite high compared to the shorter links. Beside, this trend goes up with increasing variability and intensity of the rainfall, see the figure difference between light rain and shower event for example.

Figure 4 .

 4 Figure 4.8 illustrates all the generated attenuation data during lightrain, shower, organised and unorganised storm events. The y axis represents the average values of path integrated

Figure 4 . 7 -

 47 Figure 4.7 -Example of generated path integrated rain attenuation along the sample links at 5 % error magnitude with 0.1 dB quantization. Note that the link lengths on the left and right column are ∼ 2 km and ∼ 10 km, respectively. Each row represents the type of rain event (from top to bottom direction): light rain, shower, organised and unorganised storms.

Figure 4 . 8 -

 48 Figure 4.8 -Boxplot of overall generated path integrated rain attenuation data set. Note that y axis represents the mean values of path integrated attenuation computed for each rain event.

Figure 4 .

 4 Figure 4.9 illustrates the average standard deviation of measurement error sources quantified with a magnitude of α = 5% and quantization level at 0.1 dB (see subsection 4.4.2). In general, these results reflect the nature of the error sources discussed in experimental studies (Leijnse et al., 2007b; Leijnse et al., 2010; Zinevich et al., 2010). It is clearly visible that the standard deviation value gradually increases with growing link lengths and frequencies. This

subsection 3 . 3 .

 33 2 in chapter 3). In fact, the error magnitude should decrease as the frequency increases as confirmed both experimentally and theoretically by studies (Atlas and Ulbrich, 1977; Leijnse et al., 2007b; Leijnse et al., 2008; Leijnse et al., 2010; Zinevich et al., 2010).

Figure 4 . 9 -

 49 Figure 4.9 -Box plot: Average standard deviation of measurement error source at 5% magnitude and 0.1 dB quantization. Note that the values of standard deviation in y axis were the mean values of all generated standard deviations in each event: light rain, shower, organised and unroganised storm.

Figure 4 . 10 -

 410 Figure 4.10 -Box plot: Average standard deviation of measurement error source at 20% magnitude and 0.1 dB quantization. Note that the values of standard deviation in y axis were the mean values of all generated standard deviations in each event: light rain, shower, organised and unroganised storm.

  A : The problem is under-determined. The suburb areas or rural regions of the microwave link network where the link density is very low only enabling limited information about rainfall, see the sparse part of the microwave links in Figure 4.1.

Figure 5 . 1 -

 51 Figure 5.1 -Diagram of iterative algorithm used to minimize the likelihood φ given in (5.1).

( b ) 2 .

 b2 Compute the mean value R of these rain rates obtained by the formula (5.5): term m is the number of microwave links. (c) Assign mean rain rate value R to ith pixel: R 0i = R (5.7) Local apriori: The a priori rainfall vector (R 0 ) values are obtained by k number of links which are the closest to the center of the ith pixel denoted as R 0i in Figure 5.2.

Figure 5 . 3 -

 53 Figure 5.3 -An example of Grid Nesting Procedure applied to one pixel crossed by a link. The left and the right columns represent the apriori rainfall and the retrieved rain rate by the algorithm, respectively. The rows (top-down direction) are the stage 1, 2 and 3 in Grid Nesting Procedure. The line (in blue) is a sample microwave link and its labels, 'Tr.' and 'Rec.', stand for Transmitter and Receiver, respectively.

Figure 5 . 4 -

 54 Figure 5.4 -Diagram of overall structure of the sensitivity analysis test.

Figure 5 . 5 -

 55 Figure 5.5 -Influence of the apriori parameters on the efficiency of the algorithm. The red dots in the figure indicate the location of the optimal parameters of both a priori and decorrelation distance. Note that, the NSE values are averaged over all maps for each event.

Figure 5 . 6 -

 56 Figure 5.6 -Sensitivity of the algorithm to the model error β in NSE metric. Note that this test is in the presence of the measurement error with magnitude of α = 5%. The evaluation is performed for 4 types of rainfall: Light rain, shower, unorganised and organised storm (clockwise direction). The retrieval resolution is 0.5 × 0.5km 2 .

Figure 5 . 7 -

 57 Figure 5.7 -Sensitivity of the algorithm to the model error β in BIAS metric. The description of Figure 5.6 also applies to this figure.

Figure 5 . 6 and

 56 Figure 5.6 and Figure 5.7 illustrate the results of the sensitivity of the algorithm to the influence of the model error in NSE and BIAS metrics, respectively. Here, solid(black), dashed(blue),

Figure 5 . 8 -

 58 Figure 5.8 -Retrieved rain map in 4 rain events: (a) lightrain, (b) shower, (c) organised and (d) unorganised storm. The retrieval resolution is 0.5 × 0.5 km 2 . The first column is the observed rain map by weather radar. The 2nd and 3rd columns are the rain maps retrieved using 'global' and 'local' apriori choices, respectively.

Figure 5 . 9 -

 59 Figure 5.9 -Comparison between inverse model and a priori knowledge at 2 × 2, 1 × 1 and 0.5 × 0.5km 2 resolutions which belong to 1 st 2 nd and 3 rd columns, respectively. Red solid lines indicate the NSE values using only apriori knowledge before applying the inverse algorithm. Blue solid lines represent the NSE values obtained after applying the inverse algorithm.

Figure 5 .

 5 Figure 5.10 gives overall retrieval efficiency of the algorithm in NSE metric depending on the density of microwave links. This is performed based on the pixel density map (Figure 4.4) obtained in subsection 4.2.2. It is important to note that the line colors in Figure 5.10 are the same as in Figure 4.4 except yellow color which is labelled as Whole density indicating the evaluation result without density map. It seems that the accuracy of the algorithm in the denser part of the network i.e. blue line is quite high compared to the sparser part. That is to say, the rainfall retrieved in the city center is more accurate than in the suburb area. However, regions

Figure 5 . 10 -

 510 Figure 5.10 -The retrieval efficiency in NSE at at 0.5 × 0.5 km 2 resolution in different parts of the monitoring area classified as high, moderate, low and whole density. The evaluation test is in the presence of magnitude of measurement (α = 5%), model (β = 5%) and quantization ( = 0.1 dB) errors. The figures in the clock-wise direction represent light rain, shower, unorganized and organized storm. Note that map numbers in x axis are ordered to illustrate NSE values in ascending order.

Figure 5 . 11 -

 511 Figure 5.11 -Scatter plot: The measurement, quantization and model errors are at the magnitude of α = 5%, ∆ = 0.1 dB, β = 5%. Right and left hand side of the figure show the scatterplots computed for whole and high density regions of the study area, see Figure 4.3. Note that the linear fit has been calculated over the whole dataset for each event at 0.5 × 0.5 km 2 .

Figure 6 . 1 -

 61 Figure 6.1 -Tomographic projection following θ angle (from Servières, 2005).

P

  (ν, θ) exp(2iπν(x(cos θ) + y(sin θ))) |ν| dνdθ.(6.6)

  (b, p, q) ∈ Z 3 .

Figure 6 .

 6 Figure 6.2 gives an example of the Dirac Mojette projections of a 4 × 4 image for 3 projections directions.

108CHAPTER 6 .Figure 6 . 2 -

 662 Figure 6.2 -Mojette Dirac transform of a 4 × 4 image with (p, q) ∈ [(1, 0), (1, 1), (-2, 1)] (from Servières, 2005).

Figure 6 . 3 -

 63 Figure 6.3 -Farey-Haros series of order 10 and its symmetrics (from Servières, 2005).

. 25 )

 25 Discretizing this filter in the Mojette geometry gives[START_REF] Servières | The Mojette transform: Discrete angles for tomography[END_REF]:

ure 6. 1 . 3 .

 13 Measurement model that relates the ray attenuation and object values to be reconstructed is completely linear. However, we can see exactly the opposite in the case of commercial microwave links that:(i) the geometry of the network topology is irregular, (ii) the operating frequencies and the link lengths are different and inhomogeneous (Figure4.1), (iii) the measurement model that relates signal attenuation, A in dB, to rain intensity, R in mm.h -1 , is nonlinear, i.e A -R model given in equation (4.3).

6. 4 .

 4 RECONSTRUCTION OF RAIN MAP USING FILTERED BACK PROJECTION MOJETTE111 generated attenuation data with nonlinear k-R relation by an "equivalent" attenuation associated with a quasi-linear k-R relation. Hence, the result of the transformation procedure should give 'homogenized' and 'linear' data which is equivalent to the one generated in chapter 4 . The data transformation procedure is explained as follows. Let us assume that one microwave link operates at both 18 and 23 GHz as well. Sometimes, we can see this situation in the given network that two microwave links with the same length overlaying each other (Figure 4.1). We denote the rain attenuation along the link, L at 18 and 23 GHz as A 1 and A 2 ,

(6. 31 )

 31 The formula(6.31) transforms the attenuation measured at 18 GHz into the equivalent attenuation at 23 GHz. The same transformation can be done for 38 GHz by only changing the values of A 1 , a 1 and b 1 . Here, the attenuation data is supposed to be homogeneous at 23 GHz. Because the A-R relation at this frequency is approximately linear showing its exponent b 2 being equal to 1.05 which is nearly 1.

Figure 6 . 4 -

 64 Figure 6.4 -Scatter plot between the original data at 18, 23 and 38 GHz and transformed data at 23 GHz.

Figure 6 . 5 -

 65 Figure 6.5 -Selected 21 sub-grids in black rectangles with different sizes for projection. Note that the area size of each sub-grid has been chosen visually based on the density of the links at 2 × 2 km 2 resolution.

Figure 6 . 6 -

 66 Figure 6.6 -Link projection in 4 th sub-grid. Here, the projections from one link (red segment line) in 2 directions are labelled as P roj 1 and P roj 2 . The term α 1 and α 2 are the angles between projections and link.

6. 4 .

 4 RECONSTRUCTION OF RAIN MAP USING FILTERED BACK PROJECTION MOJETTE115 6.4.2.3 Angular interpolation An angular interpolation technique proposed by (Servières et al., 2005) was applied to fill empty bins in each projection based on nearby values. Figure 6.7 exemplifies an implementation of the angular interpolation technique. The main idea here is to fill the bin value b 3 located in the projection (p3, q3) by using the values of b 1 and b 2 located in neighbouring projections (p1, q1) and (p2, q2), respectively.

Figure 6 . 7 -

 67 Figure 6.7 -Example scheme of angular interpolation. Here, two bins (b 1 and b 2 ) are interpolated onto the bin b 3 . (C x , C y ) represents the central point of the pixel. (Servières et al., 2006)

5 . 6 .

 56 Convert the reconstructed attenuation map into rainfall rate field using the inverse form of the k-R relation at 23 GHz; Compare the reconstructed images with rainfall maps by weather radar at 2 km × 2 km resolution.

Issue 2 .

 2 Integration of a single or multiple links with traditional sensors to refine rainfall measurement information: Recent advances in this category were outlined for radar and rain gauges separately. Issue 3. Development of spatial rainfall retrieval algorithms: Existing algorithms that can be used to reconstruct rainfall map using signal attenuation data from an existing network of commercial microwave links into rainfall maps were discussed. Algorithms include: Interpolation (Kriging, Inverse Distance Weighting), tomography. Very few studies have dealt with this issue. Our contribution has been found in this growing area of the research.

RésuméL

  'objectif de cette étude est d'évaluer la faisabilité de la reconstitution des champs de précipitations dans les zones urbaines à l'aide de liaisons hyper-fréquence commerciales. Dans un cadre de simulations appliquées à une étude de cas réaliste, une attention particulière est portée au réseau de téléphonie mobile et à la ville de Nantes en France. Comme base, nous utilisons un ensemble de données composé de 207 images radar météo enregistrées par Météo-France (bande C) à haute resolution spatiale (250m x 250m) et temporelle (5 min). Nous avons généré des données d'atténuation de pluie le long de ces liaisons qui sont utilisés pour représenter le signal réel. Ces données d'atténuation peuvent être enregistrées dans les stations d'antenne à hyper-fréquence. Deux algorithmes de reconstitution, à savoir l'algorithme inverse et la tomographie, sont proposés pour reconstruire des champs de pluie sur la base des données d'atténuation de pluie générées. Pour valider les performances de la méthode, nous réalisons des séries de tests de reconstitution de pluie pour des pluies légères, des averses, des orages organisés et non organisés avec des variabilités spatiales et temporelles. Les comparaisons entre les algorithmes inverse et tomographique sont réalisées également pour plus de 40 champs de précipitations et les 4 types d'événements afin de définir les limites et les capacités des deux algorithmes.

  study domain is the city of Nantes, France. As a basis, we use a data set consisting of 207 weather radar images recorded by Météo-France C band at high spatial (250m x 250m) and temporal (5 min) resolutions. We generated rain attenuation data along those links that are used to represent real signal attenuation data that can be recorded at microwave antenna stations. Two retrieval algorithms, namely inverse and tomography, are proposed to retrieve rainfall fields based on the generated rain attenuation data. To validate the system performance, we carry out series of rain retrieval tests in light rain, shower, organized and unorganized storm events with high spatial and temporal variability. The comparisons between the inverse and the tomographic algorithms are also made over 40 rainfall fields in 4 events in order to define limitations and capabilities of both algorithms.

Table 2 .

 2

1 -Gamma model parameter values N 0 and Λ, (source: Fiser, 2010). Rain type N 0 , mm -4 m -3 Λ, mm -1 µ Convective 6.29 × 10 5 R -0.416 8.35R -0.185 3 Stratiform 2.57 × 10 4 R 0.012 5.5R -0.129 3 Figure 2.7 -Gamma DSD plot in convective and stratiform rain, (Table 2.1).

Table 2 . 2 -

 22 Exponential DSD parameters: N 0 and λ in different rain types[START_REF] Joss | The variation of rain drop-size distributions at Locarno, paper presented at[END_REF]).

	(source: Fiser, 2010).
	Rain type	N 0 , mm -1 m -3	Λ, mm -1
	Average	8, 000	4.1R -0.21
	Drizzle	30, 000	5.7R -0.21
	Widespread	7, 000	4.1R -0.21
	Thunderstorm	1, 400	3R -0.21
	2.4 Rainfall intensity
	Rainfall intensity is measured directly by raingauges, distrometers or indirectly by weather
	radar with taking into account drop size distribution and raindrop velocity parameters. If the
	effects of wind (notably updrafts and downdrafts), turbulence, and raindrop interaction are ne-
	glected, the (stationary) rain rate, denoted as R in mm.hour -1 , is related to the rain-drop size
	distribution N (D) according to:

  , 12.5, 25, 50, 100, 150 mm.hour -1 . The specific attenuation values as a function of frequency are depicted in Figure 2.9. The figures a and b represent Marshal Palmer (2.19) and Gamma DSD models (2.18), respectively.

	R) -drop size distribution as a function of drop diameter and rain rate, defined in the
	expression (2.18) and (2.19) of section 2.3, mm -4 ; D max , D min -maximum and minimum rain
	drop diameter, mm; D -Drop diameter, mm; R -Rain intensity, mm.hour -1 , defined in the
	expression (2.20) of section 2.4; the term e is constant (Euler's number) equal to approximately
	2.718.
	Due to its simplicity, the following empirical relation is most commonly used to compute
	the specific rain attenuation:
	k = aR b	(2.26)
	Where, a and b -the power law coefficients. The expression (2.26) is often called k-R relation
	(Atlas and Ulbrich, 1977; Olsen et al., 1978). The power law coefficients (i.e. a and b) depend
	on frequency, polarization, temperature and DSD model. These coefficients are obtained exper-
	imentally on the basis of k (the equation 2.25) and R (the equation 2.20) both of which depend
	on the DSD recorded by a distrometer in certain climate conditions. It is important to note that
	the output of the equation (2.26) is almost the same as that of (2.25). However, these coefficients
	are climate dependent variables which mean that the approximated values may differ from one
	region to another. Therefore, the main challenge is to accurately measure the DSD data in order
	to compute those a and b values.
	Many studies (Olsen et al., 1978; Zhang and Moayeri, 1999; Das et al., 2010) established
	the k-R relation for different climate regions of the world. If there is no available measured DSD
	data, International Telecommunication Union (ITU-R, 2005) provides standard coefficients for
	a wide range of frequencies from 1 to 1000 GHz which can be used globally.
	2.6.1 Specific attenuation as a function of frequency
	In order to better understand the specific rain attenuation dependence on signal frequency
	magnitude we compute the formula Equation 2.25 as a function of frequency. We choose dif-
	ferent rainfall intensities of 2.5, 5

  Some of the issues addressed in recent studies include the errors and uncertainties related to a measurement model which is used to convert signal attenuation data into rain intensity.Due to its simplicity, the k-R relation is often used as a measurement model, see more details in section 2.6 of chapter 2. The a and b coefficients of this relation are computed on the basis of drop size distribution measured by distrometer[START_REF] Bringi | Raindrop size distribution in different climatic regimes from disdrometer and dual-polarized radar analysis[END_REF] and the extinction cross section for a given frequency[START_REF] Oguchi | Electromagnetic wave propagation and scattering in rain and other hydrometeors[END_REF]. As discussed in chapter 2, these coefficients essentially depend on 3 factors: (i) the link frequency and polarization, (ii) the temperature of a rain drop, (iii) the variability of drop size distribution. Such factors are the main source of errors that can lead to under/over estimation of the actual rain rate[START_REF] Atlas | Path-and area-integrated rainfall measurement by microwave attenuation in the 1-3 cm band[END_REF]

Leijnse et al., 2007b;[START_REF] Leijnse | Errors and uncertainties in microwave link rainfall estimation explored using drop size measurements and highresolution radar data[END_REF][START_REF] Zinevich | Prediction of rainfall intensity measurement errors using commercial microwave communication links[END_REF]

. So, the issue is to quantify the influence of each factor on the measurement model accuracy.

Table 4

 4 .1. The rows count, mean, min, max represent the number of microwave links, average, minimum and maximum lengths, respectively. The percentiles 25%, 50%, 75% show the value below which the percentage of the link lengths may be found. Except count, all rows have the same unit in km. Maximum length of the microwave links at 18, 23 and 38 GHz is 12.6, 16.8 and 9.1 km, respectively. At all frequencies, the minimum lengths are almost the same and less than 1 km.

Table 4 . 1 -

 41 Statistical description of the microwave links lengths. 18 GHz 23 GHz 38 GHz Total

	count	51	102	103	256
	mean	5.3	5.2	2.6	4.2
	min	0.8	0.5	0.3	0.3
	25%	2.3	2.9	1.5	2.1
	50%	5.2	4.8	2.5	3.6
	75%	7.6	7.1	3.7	5.2
	max	12.6	16.8	9.1	16.8

Table 4 . 2 -

 42 Characteristics of the link density in the study area.

	Number of pixels at resolutions

Table 4 .

 4 2 gives the overall summary of the pixel density map. In addition to the classified density regions, we also denote 'Whole zone' to express an entire region without pixel density map. Number of pixels at different resolution shown in the table correspond to the classified density levels in blue, red and grey colors. However, it is important to note that the density level range in D depends on the pixel size. Therefore, it can be different for another type of microwave link topology.

Table 4 . 3 -

 43 Characteristics of the Treillières radar.

	Frequency band	C
	Latitude	47 o 20 12 N
	Longitude	1 o 39 09 N
	Antenna altitude	80 m
	Beam width	1.25 o at -3 dB
	Plan Position Indicator	0.4 o , 0.8 o , 0.5 o
	Temporal resolution	5 min
	Spatial resolution	0.25 × 0.25 km 2

Table 4 . 4 -

 44 Rainfall event periods selected for the case study.

	Event type	Duration Number of maps	Period
	Light rain	03h05	37	Sep-2009
	Light rain	06h20	76	Jun-2012
	Shower	01h30	18	Jun-2012
	Shower	40 min	8	Jul-2012
	Organised storm	02h50	34	Jun-2012
	Unorganised storm	45 min	19	May-2009
	Unorganised storm	75 min	15	Jun-2009

Table 4 .

 4

			5 -Summary statistics of rain events.
		Light rain Shower Organised storm Unorganised storm Total
	count	113	26	34	34	207
	mean	1.8	4.1	4.0	6.3	3.0
	std	1.6	8.2	6.8	10.8	5.8
	q1	0.8	0.4	0.7	1.1	0.8
	q2	1.4	1.2	2.2	2.9	1.4
	q3	2.2	3.9	4.5	7.0	2.9
	max	64.6	133.2	133.2	133.2	133.2
	4.4 Generating rain attenuation data

4.4.1 Formulation of the rain attenuation generation

Table 4 .

 4 6.

Table 4 . 6 -

 46 Sample links used for visualizing the generated signal

	Link ID Frequency, GHz Length, km
	35	18	1.92
	43	18	9.89
	58	23	1.73
	129	23	9.9
	168	38	1.96
	203	38	9.06
	It is worth mentioning that the left and right columns of

Table 5 .

 5 [START_REF] Atlas | Path-and area-integrated rainfall measurement by microwave attenuation in the 1-3 cm band[END_REF]. When the values shown in this table are used the algorithm performs best. It should be mentioned that optimum choices also include all the combinations of parameters represented in red dots in Figure5.5.

Table 5 . 1 -

 51 Optimum parameters of the model obtained in the sensitivity analysis.

		Light rain Shower	Organised storm	Unorganised storm
	Decorrelation distance, ∆ (km)	3	3	3	2
	Closest links number, k	5	10	7	5

Table 5 . 2 -

 52 Overall statistical performance of the retrieval algorithm at 0.5×0.5 km 2 resolution.

	Rain type	NSE	α = 5%, β =5% ρ BIAS RMSE NSE	α= 20%, β = 20% ρ BIAS RMSE
	Light rain	0.58 0.76 0.04	0.76	0.54 0.74 0.02	0.80
	Shower	0.45 0.68 0.17	3.96	0.44 0.66 0.05	4.02
	Organised storm	0.44 0.66 0.06	3.90	0.42 0.65 0.01	3.96
	Unorganised storm 0.55 0.75 0.52	6.70	0.53 0.74 0.41	6.82

Table 5 .

 5 2 illustrates the overall statistical performance of the retrieval algorithm in the presence of measurement and model error for two extreme cases: 5% and 20%. It is worth mentioning that the values in this table are obtained by averaging all metric values in each rain type. Surprisingly, bias seems to be counter-intuitive showing the values less in 20% case than in 5%. However, the algorithm performance decreases in all statistical metrics as the variability of rainfall increases (from light rain to unorganised storm). Overall, the algorithm is stable enough indicating negligible difference between those extreme cases.

  .14) 106CHAPTER 6. RETRIEVAL MODEL 2. DISCRETE MOJETTE TRANSFORM ALGORITHM 6.3 A Discrete Tomographic Method: the Mojette transform 6.3.1 A discrete version of Radon transform: the Dirac Mojette projections The Mojette projection is a discrete version of the Radon transform defined about twenty years ago by Jeanpierre Guédon and others (Guédon et al., 1995; Guédon, 2013). With a

https://en.wikipedia.org/wiki/Drop_(liquid)

support the hypothesis that microwave links are capable of providing space-time rainfall information based on the signal attenuation measurement. It is worth mentioning that not a single study has failed to prove this hypothesis. However, three most important scientific challenges from the studies discussed so far which are still subject to further investigations: (i) Quantifying/modelling/estimating error sources regarding microwave link based rainfall measurement;ii) Integration of a single or multiple links with traditional sensors to refine rainfall measurement information; (iii) the feasibility of rainfall mapping based on the attenuation measurement from commercial microwave links, in particular, the development of algorithms. In this regard, the latter subject is the main objective of the present study. Overall, the reviewed works help conduct our feasibility study in the framework of simulation that will be presented in the next chapter.

(a) Light rain, Time: 00h30, 17.09.2009 (b) Shower, Time: 20h45, 20.06.2012 (c) Organized storm, Time: 05h25, 07.06.2012 (d) Unorganized storm, Time: 18h10, 11.05.2009

(a) Light rain event (b) Shower event (c) Organised storm event (d) Unorganised storm event

Acknowledgements

this project. He is very generous and kind man. This space is not enough to describe his support in the last three years. I would like to thank Erasmus Mundus of European Union for funding this research and the staff of ARCADE project, especially Mme Francesca Chicco at the Department of International Affairs at Politecnico di Torino. and Prof. Benis Fouad at International Relations Office

List of Figures

Metrics

The evaluation of the proposed algorithm is performed using the following statistical metrics:

1. Nash-Sutcliffe Efficiency (NSE) gives an overall assessment between the observed and retrieved rainfall map:

The range of NSE is (-∞; 1]. The more close NSE value to 1 is, the more accurate the model retrieval will be; 2. Pearson's correlation coefficient (ρ) measures the strength and direction of a linear relationship between observed and retrieved rainfall map: 3. Bias estimates the difference between observed and retrieved rainfall map:

The bias metric ranges between [-∞; ∞]. The estimator is said to be unbiased if bias = 0, positive and negative bias if bias > 0 and bias < 0, respectively.

4.

Root Mean Square Error (RMSE) gives the difference between observed and retrieved rainfall map. In other words, it demonstrates the standard deviation of the residuals between observed and retrieved rainfall map.

The value of RMSE is always positive. The model accuracy is perfect if RMSE = 0.

In all metrics mentioned above, R is observed rainfall vector, R is retrieved rainfall vector, R i is mean value of observed rainfall vector, n is the number of values in the only in 'reconstructible' sub-grids. To compute the accuracy of the reconstruction algorithm the equations 5.14 -5.17 are used.

Results and discussions

In this section, the results obtained by the MT algorithm are presented. In particular, in order to understand its capability as well as limitations of the proposed algorithm we discuss two important aspects: (i) General reconstruction performance over an entire region and comparison with the Inverse algorithm presented in the previous chapter, (ii) The influence of the microwave links network topology on the algorithm accuracy.

To evaluate the performance of the MT algorithm, we carried out 40 reconstruction tests using the protocol explained in subsection 6.4.3. The reconstruction tests are grouped by event type, i.e. 10 maps for each light rain, shower, organised and unorganised storm events which have also been used for validation purpose in chapter 5. All these tests have been performed in the presence of 5% magnitude of measurement and 0.1 dB quantization errors.

Regarding the first aspect of our analysis, the example of the reconstructed rain map for each event is illustrated in Figure 6.9. The reconstructed area consists of 106 pixels covering 448 km 2 surface of the study area with a resolution of 2 × 2 km 2 . Since the resolutions of rain maps by weather radar and inverse model are 0.25 × 0.25 km 2 and 0.5 × 0.5 km 2 , respectively, they have been averaged at 2 × 2 km 2 resolution to be able to compare with the results of the Mojette Transform algorithm. Then, the averaged maps by weather radar have been considered

to be a ground-truth. We can see that the performance of the MT algorithm is poor compared to the Inverse algorithm although it is capturing the spatial structure of the rainfall. Perhaps, the most striking and unexpected finding is the systematic overestimation in all events. However, the MT algorithm is only failing if the rainfall intensity is very high, especially more visible in shower event (see, the figure on the right in the 2nd row in top-bottom direction).

Another important finding is an increase in overestimation with respect to higher rain rates resulting in larger error in all events. A possible explanation for this may be the lack of adequate adjustment procedures, especially, choosing the 'reconstructible' sub-grids since it was selected manually in this study. In addition to that, interpolating the empty bins by their neighbourhood values can be another indication, perhaps a major factor to such outcomes. This indication will be analysed in detail in Figure 6.10. same trend with increasing variability in rain events. Besides, a similar performance can be found in Bias and RMSE as well. One interesting finding is that the MT algorithm shows the minimum Bias in all events, resulting in 270.02 mm.h -1 (see, shower event), while the Inverse algorithm has its maximum Bias for the same case. Overall trend seems to suggest that the MT algorithm can perform well if the rainfall variability is very high while this is the opposite case for the Inverse algorithm. To further investigate this behaviour, we address the second aspect of our analysis which was stated in the beginning of the section. Regarding the second aspect stated above, Figure 6.10 illustrates spatial statistics in RMSE metric that has been computed for testing the influence of network geometry on the reconstruction accuracy of the MT algorithm. We chose the RMSE to carry out the test because the measurement unit in this metric is the same as rain intensity which is mm.h -1 . In this test, only 5 specific locations have been selected which are depicted in the legend of Figure 6.10. These locations represent sub-grid region numbered in Figure 6.5. The sub-grids 11, 12 and 20 are the densest parts of the network area while the sub-grids 15 and 19 are very sparse regions. The re-limitations of the algorithm. Since this algorithm heavily relies on apriori knowledge, optionally, different types of apriori knowledge can be introduced to see a greater degree of accuracy on this matter.

Various issues on the tomographic algorithm proposed in chapter 6 need to be addressed: i) improvement on the 'reconstructible' sub-grid choice, ii) improvement on the way of taking projection, iii) A sensitivity analysis to the application conditions such as 'orthogonality level' and quantification of associated error sources of the model parameters needs to be investigated.

Furthermore, it would be interesting to test these two proposed retrieval algorithms in different cities of the world based on a real case study that requires signal attenuation data directly obtained from cellular phone companies. However, this idea will naturally bring a challenge that signal data recorded at antenna stations must be cleansed before using for rainfall retrieval.

Signal attenuation unrelated to rainfall includes baseline effects (i.e. gas, water vapour, temperature, wind) and antenna wetting after the rainy period. Algorithms developed by ( Global perspective on the main objective of this research could be to combine commercial microwave links with weather radar to improve the accuracy of real-time rainfall monitoring over urban areas. This will reduce a lot of money for hydrology applications.

A Central slice theorem

The Central Slice theorem gives us the 1D Fourier transform of a parallel projection of one object is equal to a slice of the 2D Fourier transform of this object.

Given p(t, θ) the projections of the f function defined for t ∈ R and θ ∈ [0, π[: (A.4)

The centrale slice theorem comes from writing the 2D Fourier transform of the f function in a

Cartesian coordinate system:

f (x, y)e -2iπ(xU +yV ) dxdy.