Adjustment of state space models in view of area rainfall estimation

TL;DR: In this paper, state space models and the Kalman filter were used to merge weather radar and rain gauge measurements in order to improve area rainfall estimates, and the importance of the parameters estimation method to the accuracy of mean area precipitation estimates was investigated.

Abstract: This paper uses state space models and the Kalman filter to merge weather radar and rain gauge measurements in order to improve area rainfall estimates. Particular attention is given to the estimation of state space model parameters because precipitation data clearly deviates from the normal distribution, and the commonly used maximum likelihood method is difficult to apply and does not perform well. This work is based on 17 storms occurring between September 1998 and November 2000 in an area including part of the Alenquer river hydrographical basin. Based on these data, the work aims to investigate the importance of the parameters estimation method to the accuracy of mean area precipitation estimates. It was possible to conclude that the distribution-free estimation methods produce, in general, better mean area rainfall estimates than the maximum likelihood. Copyright © 2010 John Wiley & Sons, Ltd.

Summary

1. INTRODUCTION

• Many problems in meteorology and hydrology need an accurate measurement of the total rainfall amount in a certain area.
• Alternatively, weather radar may outline accurate rainfall isopleths , but their point estimates are not so good due to errors of either a meteorological or instrumental nature.
• Krajewski (1987) and Severino and Alpuim (2005) apply an optimal interpolation method based on Kriging and CoKriging.
• There are several ways of designing the state space model that relates the gauges observations, the radar observations and their bias.

M. COSTA AND T. ALPUIM 2

• This is called the measurement equation whereas the equation governing the time variation of the calibration factor is called the transition equation.
• Equation (1) describes the relationship between a rain gauge and a radar cell in the same site, but it can be used to produce a mean-field bias, replacing Gt and Rt by the vectors comprising all the gauges in a certain area and the corresponding radar cells.
• The estimation problem has been addressed since the first works on radar calibration appeared.
• Ahnert et al. (1986) propose intuitively based estimation methods, but nothing is proved about their statistical properties.
• In many applications, the state space models’ parameters are estimated by maximum gaussian likelihood via the Newton-Raphson method (Harvey, 1996) or, more often, the EM algorithm (Shumway and Stoffer, 1982).

M. COSTA AND T. ALPUIM 4

• There are five rain gauges located in the area being studied, namely, Abrigada (A), Olhalvo (O), Penedos (P), Meca (M) and Merceana (Mr), which correspond to a reasonably high density of gauge/20km2 .
• Brown et al. (2001) consider a circular region with an approximate density of one gauge per 500 km2.
• Furthermore, because the original radar data corresponded to 10-minute periods, the authors used the average of the six available measurements in each hour as the hourly data in each cell of 2kmx2km.
• The rain gauges in Abrigada and Olhalvo are used to calibrate the radar whereas the remaining three gauges are used to compute the mean area precipitation, usually known as the “ground truth”, in order to assess the quality of the adjusted models.
• The gauges used for calibration are not used for assessment, since the two processes should be independent.

3. THE KALMAN FILTER ALGORITHM

• The Kalman filter technique, Kalman (1960), has been used in many different scientific areas to describe the evolution of dynamic systems.
• The first one considers that in each site where a rain gauge is available the relationship between the gauge, the corresponding radar cell measurement and the bias is described by the set of two equations Later, a field of calibration factors is calculated for each radar cell, with the help of interpolation methods.
• It has fewer parameters to estimate and may also produce good results when applied to small areas or to storms where the precipitation pattern is spatially homogeneous.
• After estimation of the parameters, the Kalman filter may be applied to produce an estimator of the state variable bt, at each time t.

4. ESTIMATION OF THE PARAMETERS

• There are several methods that may be used to estimate the state space models parameters.
• Meteorological and hydrological data are, generally, not normally distributed.
• This fact, together with the complexity of the ML method for the models in use, shows the need to look for alternative estimation procedures.
• The authors present two different distribution-free methods which are based on the generalized method of moments.

M. COSTA AND T. ALPUIM 6

• Thus, the log-likelihood is given by the sum of the loglikelihoods for each storm.
• The use of numerical methods to maximize this likelihood function may be a difficult and complex task, very often without satisfactory results.
• This problem may occur either because the numerical iterative techniques do not converge or because of the existence of multiple critical points.

4.2 Distribution-free estimation for the parameters – method M1

• Costa and Alpuim (2010) propose consistent distribution-free estimators for the parameters of the univariate model (2)-(3) based on the generalized method of moments.
• These estimators have an explicit analytical expression and do not assume any specific distribution for the errors.
• Considering the multi-factor model as described in equations (2)-(3), the method estimates first the mean of the calibration factor, ][ tbE=µ , and then the autoregressive coefficient φ.
• Finally, note that the estimators for the bias mean and autocovariance are averages calculated over all pairs of observations for which the radar is not null.
• The maximum likelihood method may not converge at all and produces estimates falling outside the parameters space more frequently than this method, especially for small values of the autoregressive coefficient.

M. COSTA AND T. ALPUIM 8

• Alpuim (1999) suggested noise variance estimators for univariate and multivariate state space models, designated in the text by D(k).
• The authors need, however, to estimate the autoregressive parameter of the state equation.
• Thus, in order to apply this method, the authors combine the D(k) estimators of noise variances with the mean and autoregressive parameters estimators as in method M1.
• Usually the indices are taken as k =1 and 2=l , but if these values produce estimates outside the parameters space, other values should be used.
• They tend, however, to produce estimates outside the parameters space more often than method M1.

5. ASSESSMENT OF THE RADAR-RAIN GAUGE ADJUSTMENT

• This section describes how to assess the several models’ performance according to two criteria: point and mean area precipitation.
• First, radar and gauge measurements have different natures: while radar represents an area mean value (radar cell), gauges measure precipitation in a point (gauge location).
• This is considered to give the best quality pattern because, due to its nature, the true value of mean area precipitation cannot be directly measured.
• The “ground-truth” is evaluated using interpolation methods in a reasonable thin grid of rain gauges.
• In this work the authors used also the Thiessen polygons and the inverse square distance methods.

6. RESULTS AND CONCLUSIONS

• Rain gauges located at Abrigada and Olhalvo were used for the estimation and calibration processes while the remaining gauges were kept for the validation process.
• Also, from inspection of figure 3 the authors may conclude that the ratios G/R at Olhalvo have a large variability, larger than at Abrigada, and this characteristic is reflected in the distribution-free estimates: method M1 produces estimates, both for the mean and the error variance, larger at Olhalvo than at Abrigada.
• The analysis of the isopleths of accumulated adjusted radar estimates provides a good appreciation of the final effect of the calibration process.
• --- FIGURE 6 --- Figure 6: Accumulated precipitation, in mm, for the storm of November 2 of 2000 for the radar non calibrated estimates and for the calibrated estimates considering the multi-factor model and the three methods of estimation, M1, M2 and ML.
• For the multi-factor model, the two different interpolation methods were compared.

Full text

For Peer Review
Adjustment of state space models in view of area rainfall
estimation
Journal:
Environmetrics
Manuscript ID:
env-09-0008.R1
Wiley - Manuscript type:
Research Article
Date Submitted by the
Author:
24-Feb-2010
Complete List of Authors:
Costa, Marco; Universidade de Aveiro, Escola Superior de
Tecnologia e Gestão de Águeda
Alpuim, Teresa; Universidade de Lisboa, Departamento de
Estatística e Investigação Operacional da Faculdade de Ciências
Keywords:
Kalman filter, state space model, parameters estimation, area
rainfall estimates
John Wiley & Sons
Environmetrics

For Peer Review
Adjustment of state space models in view of area rainfall estimation
Marco Costa
1,
*
†
and Teresa Alpuim
2
1
Escola Superior de Tecnologia e Gestão de Águeda, Universidade de Aveiro, Apartado 473, 3750-127 Águeda, Portugal
2
Departamento de Estatística e Investigação Operacional, Faculdade de Ciências, Universidade de Lisboa, Edifício C6,
Campo Grande 1749-016 Lisboa, Portugal
SUMMARY
This paper uses state space models and the Kalman filter to merge weather radar and rain gauge measurements
in order to improve area rainfall estimates. Particular attention is given to the estimation of state space model
parameters because precipitation data clearly deviates from the normal distribution, and the commonly used
maximum likelihood method is difficult to apply and does not perform well. This work is based on 17 storms
occurring between September 1998 and November 2000 in an area including part of the Alenquer river
hydrographical basin. Based on these data, the work aims to investigate the importance of the parameters
estimation method to the accuracy of mean area precipitation estimates. It was possible to conclude that the
distribution-free estimation methods produce, in general, better mean area rainfall estimates than the maximum
likelihood.
KEY WORDS: Kalman filter; state space model; parameters estimation; area rainfall estimates.
1. INTRODUCTION
Many problems in meteorology and hydrology need an accurate measurement of the total
rainfall amount in a certain area. For this purpose, the best results are achieved using both rain
gauges and weather radar measurements. However, although rain gauges may provide good
point rainfall estimates, they fail to depict the rainfall spatial distribution. Alternatively,
weather radar may outline accurate rainfall isopleths (Figure 1), but their point estimates are
not so good due to errors of either a meteorological or instrumental nature. This problem can
be reduced to some extent if the radar is carefully calibrated. There are several ways to
combine the rain gauges and radar estimates taking into consideration the different but
complementary nature of the two sensors. Krajewski (1987) and Severino and Alpuim (2005)
apply an optimal interpolation method based on Kriging and CoKriging. Calheiros and
Zawadzki (1987) and Rosenfeld et al. (1993) use a probability matching method. However,
the more commonly used method to merge the radar and rain gauge fields is to relate the two
types of measurements through a state space model and consequent application of the Kalman
filter.
There are several ways of designing the state space model that relates the gauges
observations, the radar observations and their bias. A pioneering work by Ahnert et al. (1986)
suggests the relationship
*
Correspondence to: Marco Costa, Escola Superior de Tecnologia e Gestão de Águeda, Universidade de Aveiro, Apartado 473, 3750-127
Águeda, Portugal.
†
E-mail: marco@ua.pt
Page 1 of 22
John Wiley & Sons
Environmetrics
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

For Peer Review
M. COSTA AND T. ALPUIM 2
G
t
=
b
t
R
t
+
e
t
, (1)
where G
t
and R
t
represent the gauges and radar measurements respectively, and b
t
stands for
the bias or dynamic calibration factor. The measurement error series e
t
is a white noise
sequence with variance
σ
e
2
. This is called the measurement equation whereas the equation
governing the time variation of the calibration factor is called the transition equation. Ahnert
et al. (1986) propose that the bias should follow a random walk as it is equally likely to
increase or decrease. Alpuim and Barbosa (1999), however, achieve better results in the
calibration process using a more general first-order autoregressive, AR(1), process for the
transition equation.
--- FIGURE 1 ---
Figure 1: Image of the radar surface rainfall intensity field.
Equation (1) describes the relationship between a rain gauge and a radar cell in the same site,
but it can be used to produce a mean-field bias, replacing G
t
and R
t
by the vectors comprising
all the gauges in a certain area and the corresponding radar cells. In this approach Lin and
Krajewski (1991) consider a different type of state space model defining the measurement
equation as
z
t
=
b
t
+
e
t
where z
t
represents the observed mean bias, that is, the ratio of the sum of all rain gauge data
available at time t, over the sum of the corresponding radar data. The study of this model was
continued in other works with minor variations in the way the mean-field bias is defined, as
may be seen in Anagnostou et al. (1998), Anagnostou and Krajewski (1999) or Chumchean et
al. (2004).
Brown et al. (2001) propose another type of state space model for the combination of the two
types of data. These authors make the fundamental assumption that the relationship between
gauge and radar reflectance measurements can be described by a power law. Thus, they
conclude that there is a linear relationship between Y
t
=
logG
t
and u
t
=
log R
t
that adjusts well
to the data if the intercept varies in time as an AR(1) process, that is,
Y
t
=
a
t
+
bu
+
Z
t
,
where b is constant and Z
t
is a gaussian white noise process with null mean and variance
σ
Z
2
.
Without being exhaustive, the models described in the previous paragraphs show that state
space models have been extensively used in radar precipitation measurement. The Kalman
filter iterative algorithm applies to these models to estimate the bias as its orthogonal
projection on the gauge measurements up to time t. This estimator is the minimum mean
square error linear predictor which, under the assumption of normality, corresponds to the
minimum variance estimator.
In this work we assume that relationship (1) is true and that the calibration factor b
t
follows an
AR(1) stationary process. Among other possible state space models, we consider this one
because it is simple to deal with and reflects the known type of dynamics relating gauge and
radar measurements. We focus on the statistical adjustment of models, giving particular
attention to the parameter estimation process. The values for the noise variances used in the
Kalman filter have a strong influence on the quality of the bias and area rainfall estimates.
Page 2 of 22
John Wiley & Sons
Environmetrics
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

For Peer Review
Calibration of Radar Precipitation Estimates 3
The estimation problem has been addressed since the first works on radar calibration
appeared. Ahnert et al. (1986) propose intuitively based estimation methods, but nothing is
proved about their statistical properties. Lin and Krajewski (1991) consider an adaptive error
estimation parameter algorithm applied to their alternative state space formulation, which
Anagnostou et al. (1998) compare with the maximum likelihood applied to the same
equations. In many applications, the state space models’ parameters are estimated by
maximum gaussian likelihood via the Newton-Raphson method (Harvey, 1996) or, more
often, the EM algorithm (Shumway and Stoffer, 1982). However, precipitation data may
deviate considerably from the normal distribution and these methods may lead to poor
estimates, fact which is also recognized by Anagnostou and Krajewski (1999). Furthermore,
the loglikelihood function of state space models with time varying coefficients - as is the case
- may have a complex shape which makes it very difficult, or even impossible, to reach the
global maximum.
Hence, there is a serious need to find estimation methods with good statistical properties, easy
to apply so that they are suitable for real time radar estimation rainfall and flexible enough to
adjust to the specific characteristics of rainfall measurements. The main objective of this
paper is to use parameter estimators based on the generalized method of moments (GMM) in
the radar calibration process. We consider two different estimators of this type, as proposed
by Costa and Alpuim (2010) and by Alpuim (1999) and compare them with the maximum
likelihood method. The GMM estimators have good statistical properties, and we generalize
this method to the case of a multivariate state space model that produces a mean field bias
based on several rain gauges in a certain area.
The data available for this study correspond to 17 storms which occurred between September
of 1998 and November of 2000 in an area of high hydrometeorological interest located
around 40 km north of Lisbon. We use this data to investigate the impact of parameters
estimation methods on the accuracy of area precipitation estimates. The authors conclude that
GMM estimators have a better performance in mean area rainfall estimation. It is interesting
to compare this conclusion with the results in Anagnostou et al. (1998), where adaptive
recursive methods involving error correlations compare favorably with maximum likelihood,
although applied to different measurement and transition equations. Besides the good results
in practice, the proposed GMM estimators are a good alternative to maximum likelihood
because they have an analytical expression and are much easier to calculate.
2. DESCRIPTION OF THE DATA
Table 1 describes the data available for this study. They consist of 17 storms between
September of 1998 and November of 2000 in a 10x10 km
2
area, including the Alenquer River
basin, located around 40 km north of Lisbon and between 31 to 44 km distance from the
weather radar in Cruz do Leão.
--- TABLE 1 ---
--- FIGURE 2 ---
Figure 2: Location of the five rain gauges in the portuguese system of coordinates and in the grid of 25 radar cells used in this
work. P – Penedos de Alenquer, A – Abrigada, Mr – Merceana, O – Olhalvo and M – Meca.
Page 3 of 22
John Wiley & Sons
Environmetrics
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

For Peer Review
M. COSTA AND T. ALPUIM 4
There are five rain gauges located in the area being studied, namely, Abrigada (A), Olhalvo
(O), Penedos (P), Meca (M) and Merceana (Mr), which correspond to a reasonably high
density of gauge/20km
2
(see Figure 2). The study area was chosen as the smallest squared
grid of radar cells including the five rain gauges available. This choice maximizes the rain
gauges density and decreases the errors associated with the interpolation methods. Compared
to other works, this area may be considered to have a high density of gauges. For example,
Lebel (1999) refers to the network used in the experimental pioneering work by Huff (1970)
of 49 gauges spread over a 1000 km
2
area, as a dense network. Brown et al. (2001) consider a
circular region with an approximate density of one gauge per 500 km
2
. The area being studied
has the highest gauge density under the radar umbrella. This fact, associated with a very low
concentration time (about 3h), which makes the region particularly subjected to flash floods,
indicates that it is very suitable to conduct radar calibration studies.
The nominal resolution of the radar in Cruz do Leão is 1kmx1km. Nevertheless, as the
resolution commonly used in these type of studies is 2kmx2km, we computed the averages of
each group of four neighbouring cells. Furthermore, because the original radar data
corresponded to 10-minute periods, we used the average of the six available measurements in
each hour as the hourly data in each cell of 2kmx2km.
For each cell of 2km
×
2km we have radar rainfall measurements produced by the recently
installed weather radar in Cruz do Leão, except the storms marked with * in Table 1, where
radar measurements are available only in the five cells with rain gauges. The rain gauges in
Abrigada and Olhalvo are used to calibrate the radar whereas the remaining three gauges are
used to compute the mean area precipitation, usually known as the “ground truth”, in order to
assess the quality of the adjusted models. The gauges used for calibration are not used for
assessment, since the two processes should be independent. The selection of the two locations
for calibration took into consideration the spatial dispersion of the remaining gauges, so that
the latter could represent the whole area and the corresponding “ground-truth”.
3. THE KALMAN FILTER ALGORITHM
The Kalman filter technique, Kalman (1960), has been used in many different scientific areas
to describe the evolution of dynamic systems. The main goal of the algorithm is to find
estimates of unobservable variables based on related observable variables through a set of
equations called a state space model. We now describe the two different formulations that will
be used to calibrate the radar, in section 6.
The first one considers that in each site where a rain gauge is available the relationship
between the gauge, the corresponding radar cell measurement and the bias is described by the
set of two equations
G
t
=
b
t
R
t
+
e
t
(2)
b
t
=
µ
+
φ
b
t −1
−
µ
(
)
+
ε
t
(3)
where the sequence of biases or calibration factors {b
t
} is a stationary AR(1) process, that is,
|
φ
| < 1. Also,
t
G represents the rain gauge measurement,
t
R the radar measurement and
t
e
and
t
ε
are uncorrelated white noise sequences with variances
2
e
σ
and
2
ε
σ
, respectively. We
will refer to equations (2)-(3) as the multi-factor model because it produces one calibration
Page 4 of 22
John Wiley & Sons
Environmetrics
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60