Surface-water quality and flow Modeling Interest Group
|
1Hydrologist USGS, Water Resources Division 6520 Mercantile Way, Suite 5 Lansing, MI 48910 Internet: djholtsc@usgs.gov Phone: (517) 887-8910 FAX: (517) 887-8937 |
2Professor of Electrical Engineering School of Engineering and Computer Science California State University - Fullerton Fullerton, CA 92807 Internet: mgrewal@fullerton.edu Phone: (714) 278-3013 |
Citation:
Holtschlag, D.J. and Grewal, M.S., 1998, Estimating ice-affected streamflow
by extended Kalman filtering: ASCE Journal of Hydrologic Engineering,
3(3), p. 174-181.
The document displayed below is based on the final draft provided to the journal. Minor discrepancies between this document and the published version, therefore, may exist.
This work was done in cooperation with the Cold Regions Research and Engineering Laboratory, U.S. Army Corps of Engineers, 72 Lyme Road, Hanover, NH 03755-1290 (Contact: Kathleen D. White, Research Hydraulic Engineer).
Abstract
Introduction
Formulation of the Problem
Application of Extended Kalman Filtering
Implementation
Temporal Updates
Observational Updates
Available Data
Results
St. John River at Dickey, Maine
Platte River at North Bend, Nebraska
Conclusions
ReferencesUSGS operates a network of about 7,000 continuous-record streamflow-gaging stations nationwide. Although stage (water-surface elevation of the stream) is continually recorded, direct measurements of streamflow and stage are obtained only at about 6-week intervals. The average (monotonically increasing) relation between stage and streamflow during open-water conditions (no ice cover) is referred to as the "stage-discharge rating." Streamflow indicated by the rating for a particular stage is referred to as the "apparent streamflow." Deviations in the average relation between stage and streamflow, which are interpreted from individual direct measurements of streamflow and stage, are described by shifts in the rating. Together with hourly or more frequent measurements of stage, the rating and shifts are used to compute streamflow records for annual publication, following a comprehensive analysis of streams in a network of gaging stations.
During open-water conditions, uncertainty associated with shifts in the rating are usually small with respect to the need for information affecting dam operations. Thus, real-time estimates of streamflow, traditionally determined directly from the rating and possibly the shift defined from the most recent direct measurement, provide sufficient accuracy. However, during ice-affected periods, the accuracy of real-time streamflow estimates are degraded by rapidly varying ice-backwater conditions and large time-varying shifts. In these conditions, the ratio of the true streamflow to the apparent streamflow, referred to as the "streamflow ratio," can vary from 1 during open-water conditions to a value near zero. These shifts create sufficient uncertainty in traditional real-time streamflow estimates that dam operations can be degraded.
Melcher and Walker (1992) evaluated and compared 17 methods for estimating ice-affected streamflow that have traditionally been or could potentially be used in the nationwide streamflow-gaging station network maintained by USGS. The methods were divided into two general categories, subjective and analytical, depending on whether judgement was necessary for method application. They identified two subjective methods that were more accurate than other subjective methods analyzed and approximately as accurate as the best analytical method. Holtschlag (1996) developed a dynamical-systems approach for computing ice-affected streamflow. This approach ranked higher than the 11 analytical methods investigated by Melcher and Walker on the basis of accuracy and feasibility criteria. The difference equation formulated by Holtschlag provided the basis for the filter developed in this report.
Formulation of the Problem
Based upon analysis of historical streamflow characterizations developed by
traditional methods for estimating ice-affected streamflow, the dynamics of
ice effects were classified into three modes. Mode 1 dynamics are associated
with the sudden formation/ablation of ice effects as indicated by abrupt
changes in apparent streamflow. Mode 2 dynamics are associated with stable
ice conditions and are approximated by a first-order difference equation
relating the streamflow to air temperature. Finally, mode 3 dynamics are
associated with the eventual elimination of ice effects at warm air
temperatures. Each mode corresponds to a unique process (state) model.
The process model for mode 1 dynamics is an algebraic equation of the form
where x1(k) (or x1(k-1)) is the estimated ratio of the actual to apparent streamflow at time indexed by variable k (or k-1). The variable h corresponds to the apparent streamflow. Mode 1 dynamics were applied for either of two conditions: (1) if the one-day change in apparent streamflow increased by more than q_dl% and the average air temperature was less then t_lo, or (2) if the one-day change in apparent streamflow decreased by more than q_dl% when the streamflow ratio was less than 1. Because of the lack of derivative information in equation 1, parameters q_dl and t_lo were included among five threshold parameters estimated outside the extended Kalman filter.
The process model for mode 2 dynamics is a first-order difference equation driven by air temperature, u. The form of the process model for mode 2 dynamics is
which indicates that at times of ice effect and constant air temperature
x5, the streamflow ratio (x1) is in equilibrium about a
nominal value x2. Changes from the equilibrium value are described
by a difference equation that includes an autoregressive component with
parameter x3 and a forcing function term associated with daily air
temperature. Air temperatures that vary from a nominal value of x5
change the streamflow ratio from its nominal value of x2 at a rate
of x4. This form of a difference equation is nonlinear in
parameters because rate parameters (x3 and x4) and
offset parameters (x2 and x5) are estimated
simultaneously. Prior information on the distribution of streamflow ratios
during periods of ice effect was used in hope of facilitating the solution of
this inherently difficult estimation problem. Air temperature,
, indicates an exponentially weighted average of
temperature from the three previous days. An exponential weighting factor,
t_wt, used in computing , was included among
threshold parameters because it did not occur explicitly in equation 2.
The process model for mode 3 dynamics is an equation of the form
The mode 3 dynamic model is applied when the exponentially weighted air temperature value exceeds t_hi. Then, the streamflow ratio increases from its value at time k-1 to a value of one when the exponentially weighted air temperatures equals t_ou. Because t_hi and t_ou occur explicitly in the process model, they could be included as parameters in the state vector. However, because of dangerous conditions for direct measurement of streamflow during mode 3 dynamics and correspondingly limited direct measurement data, t_hi and t_ou were included among threshold parameters estimated outside the filter.
Estimates of streamflow ratio for all dynamic modes were constrained to an interval between a maximum of 1 and a minimum greater than zero. In this analysis, the minimum was the minimum ratio of published to apparent streamflow determined from historical record.
The process models for the three dynamic modes were developed so that the first (only) element in the state vector was the signal element (the estimated streamflow ratio). For this application, change in mode affected only the signal element. In contrast, the state parameter vector was treated as if mode 2 dynamics were always applied. Although this convention can create uncertainty in the state parameter vector, the uncertainty is minimized because updates affecting the parameter vector are computed only for days of direct streamflow measurement. Direct measurements are not generally possible during conditions of mode 1 or 3 dynamics because of unsafe measuring conditions associated with these ice conditions.
Application of Extended Kalman Filtering
Grewal and Andrews (1997), Bar-Shalom and Li (1993), Bozic (1994), Mendel
(1995), and Brown and Hwang (1997) provide detailed information on the
mathematical development and general application of Kalman filtering. A
Kalman filter is an estimator of the state of a dynamic system given
measurements that are related to the state. The extended form of the Kalman
filter was selected because the state vector was formulated to include both a
signal element (the streamflow ratio) and unknown parameters. In addition to
the unknown parameters in the state vector, five threshold parameters are
estimated externally to the extended Kalman filter.
A discrete-time formulation was used for consistency with the availability of hydrologic and climatic data. In this report, the length of the discrete time step was one day, a length that facilitated analysis of extensive historical periods of record.
A discrete-time extended Kalman filter was developed to account for effects of ice on streamflow. The filter consists of two models, a nonlinear process model and a linear measurement model. The general form of the nonlinear process model is
where
x(k) is the state vector. In this paper, the state vector is partitioned into two components. The first component is the streamflow ratio and is referred to as the "state signal element." The second component is referred to as the "state parameter vector." The total number of elements in the state vector is referred to as the "dimension of the state space."
f(x(k-1),k-1) is a nonlinear function of the state at the previous time step plus other information on auxiliary variables available at time k-1.
w(k-1) is a value from a random sequence representing process noise at time k-1. (In the analysis of dynamic systems, noise refers to random inputs that cannot be directly measured or controlled). The sequence w is assumed to be independent and normally distributed with a mean of zero and a covariance structure, Q(k-1), commonly written w~N(0,Q(k-1)). In this application, only the variance of x1(k) was assumed be nonzero; no process noise was associated with the state parameters.
The time-varying linear measurement model is of the form
where
z(k) is the streamflow at time k,
H(k) is the time-varying measurement sensitivity matrix that is represented
by the vector
[h(k) 0 0], where h(k) is the apparent streamflow, and
v(k) is a value from a random sequence representing measurement noise at time k. The sequence v is assumed to be independent and normally distributed, with a mean of zero and variance R(k). The subset of days indexed by k on which direct measurements of streamflow were obtained is denoted as k'. For the purpose of developing an estimate, the variance R(k') was assumed to be proportional to the published streamflow on days of direct streamflow measurement. In open-water conditions, the standard error was assumed to be 2.5% of the published streamflow; during ice-affected conditions, the standard error was assumed to be 8.0% of the published streamflow. On days for which direct measurements were not made, published streamflow data were not used to update the estimate of the streamflow ratio.
Implementation
The extended Kalman filter was implemented by recursively computing daily
updates to the state vector x and the state error covariance matrix P,
given an initial estimate of the state vector x0 and
P0. The update is done in two steps: a temporal update and an
observational update. The temporal update is computed at each time step; the
observational update is computed only on days of direct streamflow
measurement. The magnitude of P increases with temporal updates and decreases
with observational updates. Bootstrap estimates of x0 and
P0 were computed by iteratively replacing
x0(j+1) with xf(j) until
xf(j+1) ~ xf(j), where
x0 and xf indicate the initial and final values for the
state vector for the j or j+1 iteration of the extended Kalman filter,
respectively.
Temporal Updates
The temporal update represents the best linear estimate of the state at time
k on the basis of measurements of z available through k-1. The form of the
general state equation is
In case of mode 1 dynamics, the signal element is computed as follows:
For mode 2 dynamics, the temporal update of the state vector is computed as follows:
Finally, in the case of mode 3 dynamics, the signal element is computed as
for specified threshold parameters t_hi and t_ou.
A temporal update of streamflow is computed by multiplying the apparent streamflow by the estimated streamflow ratio; or, for consistency with the extended Kalman filter notation, by multiplying the time-varying measurement sensitivity matrix, H(k), by the temporal update of the state vector as
A temporal update of the error covariance matrix is computed as
A first-order approximation of the state transition matrix is
The estimate of the transition matrix used in this analysis is
The value of Q was determined such that
where 
is the standardized normal deviate equal
to 1.64.
Observational Updates
Observational updates were computed for days of direct streamflow
measurement. First, the Kalman gain matrix was computed as
Then, the covariance matrix was updated as
The state vector update was computed as
Finally, the observational update was computed as
On days for which direct streamflow measurements were not made, observational updates were set equal to the temporal updates computed for that time step. Thus, estimated values computed by use of the extended Kalman filter were equal to the observational updates on days of direct streamflow measurement and were equal to the temporal updates otherwise. No adjustment was included for uncertainty in the apparent streamflow values.
Available Data
Sites in Maine and Nebraska were selected for initial application of the
filter. The St. John River at Dickey site is in northern Maine, about 160 mi
(miles) north of Bangor. Hydrologic data were obtained from USGS
streamflow-gaging station 01010500. Drainage area at the gage is 2,680
mi2 (square miles); average streamflow from 1946 through 1995 was
4,771 ft3/s (cubic feet per second). Climatic data for the site
was obtained from Fort Kent, Maine, which is about 25 mi east-northeast of
the gaging station. Climatic data for the Fort Kent station are published by
the National Oceanic and Atmospheric Administration (NOAA) as station 2878.
Seventeen years of daily hydrologic and climatic data for 1971-74, 1978-79,
and 1983-93 were selected for analysis; intervening periods were omitted
because of intervals of missing hydrologic or climatic data. The
stage-discharge relation is usually affected by ice from early December
through mid-April. Within the selected periods, 37.1% of the daily values
were ice affected.
The Platte River at North Bend site is in east-central Nebraska, about 45 mi northwest of Omaha. Hydrologic data were obtained from USGS streamflow-gaging station 06796000. Drainage area at the gage is 70,400 mi2; average streamflow from 1949 through 1995 was 4,569 ft3/s. Climatic data for the site were obtained from Fremont, Nebr., which is about 15 mi east of the streamflow-gaging station (written commun., Mat Werner, Climate Resources Specialist, High Plains Climate Center, University of Nebraska-Lincoln, 1996). For this analysis, 27 years of daily streamflow and climatic data were analyzed, including the periods 1965-67, 1969-70, 1972-90, and 1992-94; intervening periods were omitted because of intervals of missing hydrologic or climatic data. Ice effects are common on Platte River in November through March; some severe ice jams occur at the gage in most years. Within the period analyzed, 28.9% of daily streamflow values were affected by ice.
Results
Filter estimates generally refer to estimates at time k based on measurements
up to and including time k. In contrast, forecast estimates refer to
estimates at time k based on data up to, but not including, time k. In this
paper, extended Kalman filter estimates of streamflow are usually forecasts
determined on the basis of the temporal updates. However, on days of direct
measurement, more accurate filter estimates are computed by use of
observational updates.
St. John River at Dickey, Maine
The extended Kalman filter was initialized to St. John River data by manually
adjusting preliminary estimates of threshold parameter values to minimize the
sum of squared errors in the extrapolated streamflow ratio,
x1(-)(k') - x1(k'), where k'
indicates days of direct streamflow measurements. Once satisfactory estimates
of the threshold parameters were obtained, they were fixed (Table 1), and the
filter was run repetitively using the state vector and error covariance
matrix computed on the last iteration of the previous run to initialize the
subsequent run. The filter was run repeatedly until elements in the state
parameter vector were virtually constant. In this process, initial estimates
for the state error covariance matrix converged from an initially specified
diagonal matrix with nonzero components of [0.5 0.5 0.01 10] to the standard
errors listed in Table 1.
Final estimates for the state parameters indicate that mode 2 dynamics are highly autoregressive, as indicated by the parameter x3=0.981, about a streamflow ratio offset of x2=0.544. Streamflow ratios increase at a rate x4=0.000855 oC-1 (degrees Celsius-1) from the temperature offset x5=-3.19 oC. (Note, for consistency with the natural correspondence between the freezing point of water and the zero of the Celsius scale, temperature values are referenced in the International System of Units rather than the English System of Units. To convert from Celsius to Fahrenheit, multiply the value reported in Celsius degrees by 9/5 and add 32.) Although the value for x2 is higher than 0.15 (the mode of the distribution of empirical streamflow ratios) it is physically realizable. Similarly, x5 is consistent with the distribution of air temperatures during periods of ice effects, which generally ranged from -20.5 to 2.5 oC. However, analysis of the state error covariance matrix indicates a large positive correlation (>0.99) between x2 and x5 and a large negative correlation (< -0.93) between x3 and x4. Thus, although the form of the difference equation used to describe mode 2 dynamics resulted in parameters with physically realizable values, the large correlations indicate ambiguity concerning their true values. The high correlations in the state error covariance matrix indicate a potential for reducing the dimension of the state parameter vector without loss of filter accuracy.
Sensitivities for threshold parameters (Table 1) were estimated as the change in the sum of squared errors in the streamflow ratio divided by the change in the corresponding parameter near the selected values. Results of simulations indicate that filter computations were most sensitive to changes in the q_dl parameter and least sensitive to the t_lo parameter. Formal optimization of the threshold parameters could lead to further improvement in filter performance.
One measure of the estimation accuracy of the filter is the relation between temporal updates of streamflow, z1(-)(k'), and published daily streamflow values, z(k'), where k' indexes days of direct streamflow measurement (Fig. 1). Although this value is modified by observational updates to more precisely estimate streamflow on days of direct measurement, z1(-)(k') provides a conservative indication of the filter accuracy. The temporal update is considered a conservative indicator of accuracy because the variance of the estimation increases monotonically with time from the previous measurement and the length of the forecast lead is maximum just before the observational update. Results for St. John River indicate that the correlation between log-transformed values of z1(-)(k') and z(k') is 0.777 based on 40 days of ice-affected measurements and 0.998 based on 63 days of open-water measurements. Time between measurements at St. John River used in this analysis averaged 8.6 weeks.
Figure 1. Relation between published streamflow and temporal updates of streamflow on selected days of direct measurement at St. John River at Dickey, Maine, from 1971 through 1993.
Another measure of filter accuracy is the relation between published and estimated streamflows during periods of ice effects. The relation between published and estimated values is linear in the logarithm of transformed values (Fig. 2). Uncertainty occurs in both the published and estimated values. Published streamflow values during periods of ice effects are typically rated "fair" or "poor" on a subjective basis. A rating of "fair" implies that about 95% of the daily values are within 15% of the true value; a "poor" rating indicates that daily streamflow values have less than "fair" accuracy (Novak, 1985). Discrepancies between published and estimated values during periods of ice effects (Fig. 3) were computed as
Analysis of the distribution of discrepancies shows that 87.2% of the time, the absolute value of elements in the e sequence is less than 8%; 96.6% of the time, it is less than 15%.
Figure 2. Relation between published and estimated streamflow during selected ice-affected periods at St. John River at Dickey, Maine, from 1971 through 1993.
Figure 3. Distribution of discrepancies between logarithms of published and estimated streamflow during selected ice-affected periods at St. John River at Dickey, Maine, from 1971 through 1993. Logarithms are in base 10.
Streamflow and climatic data for the St. John River at Dickey, Maine site from November 1986 through April 1987 are shown on figure 4. Upper and lower estimations were computed by adjusting the variance of Q to 0.0035 so that interval formed between the upper and lower estimations about the temporal updates at k' contained the published values 90% of the time.
Figure 4. Streamflow and climatic data for the St. John River at Dickey, Maine, November 1986 through April 1987.
Results for Platte River data also indicate that mode 2 dynamics are highly autoregressive, as indicated by the parameter x3=0.990. Streamflow ratios increase at a rate x4=0.000939 oC-1 about a temperature offset x5=-9.37 oC, a lower temperature than that estimated for St. John River in Maine. Unfortunately, the estimated streamflow ratio offset of x2=-0.068 is not physically realizable. Analysis of the state error covariance matrix indicates a maximum positive correlation of 0.81 between x3 and x5 and a maximum negative correlation of -0.74 between x3 and x4. The magnitudes of these correlations are not thought to be sufficient to cause significant degradation of parameter estimates. However, given the small magnitude of the estimated x2 value, in future applications it may be possible to eliminate (set to zero) the streamflow ratio offset from the difference equation for mode 2 dynamics. Such an elimination would reduce the dimension of the state space, which would also likely reduce parameter ambiguity caused by high correlations in the state error covariance matrix. Suboptimal values for the threshold parameters also present a possible explanation for the discrepancy between the estimated value of x2 and the conceptualized value at the mode (0.30) of the empirical distribution of discharge ratios.
Sensitivities for threshold parameters (Table 2) were estimated as the change in the sum of squared errors in the streamflow ratio estimate divided by the change in the corresponding parameter near the selected values. Results of simulations indicate that estimates are most sensitive to changes in the q_dl parameter and least sensitive to the t_lo parameter. Again, formal optimization of the threshold parameters could lead to further improvement in filter performance.
The temporal updates of streamflow on days of streamflow measurements compare closely with published daily mean streamflows (Fig. 5). Results indicate that the correlation between log-transformed values of z1(-)(k') and z(k') is 0.864 based on 87 days of ice-affected measurements and 0.997 based on 345 days of open-water measurements. Time between measurements at Platte River used in this analysis averaged 3.2 weeks.
Figure 5. Relation between published streamflow and temporal updates of streamflow on selected days of direct measurement at Platte River at North Bend, Nebr., from 1965 through 1994.
The relation between published and estimated streamflows at Platte River during periods of ice effects is linear and unbiased in the logarithms of streamflow (Fig. 6). Analysis of the distribution of discrepancies between published and estimated values (Fig. 7) during periods of ice effects by use of equation 19 indicates that 90.7% of the time, the absolute value of elements in the e sequence is less than 8%, and that 97.7% of the time, it is less than 15%.
Figure 6. Relation between published and estimated streamflow during selected ice-affected periods at Platte River at North Bend, Nebr., from 1965 through 1994.
Figure 7. Distribution of discrepancies between logarithms of published and estimated streamflow during selected ice-affect periods at Platte River at North Bend, Nebr., from 1965 through 1994. Logarithms are in base 10.
Streamflow and climatic data for the Platte River at North Bend, Nebr. site from November 1977 through April 1978 are shown in Fig. 8. Upper and lower limits for estimates were computed by adjusting the variance of Q to 0.0039 so that the interval formed by the upper and lower estimate about the temporal update at k' included the published values 90% of the time.
Figure 8. Streamflow and climatic data for the Plate River at North Bend, Nebr., November 1977 through April 1978.
Three dynamic modes of ice effects were identified on the basis of historical interpretations of ice-affected streamflow records. The first mode is associated with rapid transitions in ice effects related to ice formation and ablation processes, as indicated by abrupt changes in apparent streamflow. The second mode is associated with gradual changes in ice effects with changes in air temperature. The third mode is associated with elimination of ice effects at warm air temperatures. Equations for these dynamic modes were developed within a discrete-time extended Kalman filter to estimate daily values and uncertainties of ice-affected streamflow.
The filter consists of two models, a nonlinear process model and a time-varying linear measurement model. For the dominant, mode 2 dynamics, the process model computes the ratio of actual to apparent streamflow (streamflow ratio) by use of a first-order difference equation driven by daily air-temperature values and four filter parameters. During periods of mode 1 or mode 3 dynamics, the difference equation is replaced by an algebraic expression to estimate the streamflow ratio by use of five threshold parameters.
The measurement model estimates a daily mean streamflow on the basis of the computed streamflow ratio and the apparent streamflow. On days of direct streamflow measurement, the estimate is computed and the covariance matrix is updated to account for the measurement information. The utility of the filter was evaluated by application to historical data from two long-term streamflow-gaging stations.
The filters developed in this paper were stable, and parameters converged for both stations, which allowed estimation of ice-affected streamflows. Results for the gaging station at St. John River at Dickey, Maine, indicate that, during periods of ice effects, logarithms of estimated streamflow values were within 8% of the logarithms of published values 87.2% of the time and within 15% 96.6% of the time. Results for the gaging station at Platte River at North Bend, Nebr. indicate that logarithms of estimated streamflow values were within 8% of the logarithms of published daily values 90.7% of the time and within 15% 97.7% of the time. The correlation between temporal updates and published values of streamflow on days of direct streamflow measurement are 0.777 and 0.998 for data from St. John River at Dickey, Maine, and 0.864 and 0.997 for data from Platte River at North Bend, Nebr., respectively. Analysis of the state error covariance matrix for both the St. John River and the Platte River both indicates that the number of parameters associated with the difference equation formulated for mode 2 dynamics could possibly be reduced.
The extended Kalman filter developed in this paper provides a basis for
estimating ice-affected streamflow at other gaging stations by adjusting
filter parameters to site-specific conditions. The filter can be used to
estimate daily mean streamflow during periods of ice effects by use of
real-time climatic and hydrologic data.
References
Bar-Shalom, Y., and Li, X.-R. (1993). Estimation and tracking:
principles, techniques, and software. Artech House, Boston, Mass.
Bozic, S.M. (1994). Digital and Kalman filtering, 2nd Ed. Halsted Press, New York.
Brown, R.G., and Hwang, P.Y.C. (1997). Introduction to random signals and applied Kalman filtering, 3rd Ed. John Wiley, New York.
Grewal, M. S., and Andrews, A.P. (1997). Kalman filtering theory and practice, 4th Ed. Prentice Hall Information and System Sciences Series, Englewood Cliffs, N.J.
Holtschlag, D.J. (1996). "A dynamical-systems approach for computing ice-affected streamflow" Water-Supply Paper 2473 U.S. Geological Survey.
Melcher, N.B., and Walker, J.F. (1992). "Evaluation of selected methods for determining streamflow during periods of ice effect" Water-Supply Paper 2378 U.S. Geological Survey.
Mendel, J.M. (1995). Lessons in estimation theory for signal processing, communications, and control. Prentice Hall PTR, Englewood Cliffs, N.J.
Novak, C.E. (1985). "WRD data reports preparation guide" Open-File Rep. 85-480 U.S. Geological Survey.
Back to the SMIG Features Page
Home | Mailing List | Features | Conferences | Classes | Reading | Model Archives | Feedback