Downscaling and Forecasting of Evapotranspiration ... - IEEE Xplore

12 downloads 520 Views 3MB Size Report
9, SEPTEMBER 2008. Downscaling and Forecasting of Evapotranspiration. Using a Synthetic Model of Wavelets and. Support Vector Machines. Yasir H. Kaheil ...
2692

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 9, SEPTEMBER 2008

Downscaling and Forecasting of Evapotranspiration Using a Synthetic Model of Wavelets and Support Vector Machines Yasir H. Kaheil, Enrique Rosero, M. Kashif Gill, Mac McKee, and Luis A. Bastidas

Abstract—Providing reliable forecasts of evapotranspiration (ET) at farm level is a key element toward efficient water management in irrigated basins. This paper presents an algorithm that provides a means to downscale and forecast dependent variables such as ET images. Using the discrete wavelet transform (DWT) and support vector machines (SVMs), the algorithm finds multiple relationships between inputs and outputs at all different spatial scales and uses these relationships to predict the output at the finest resolution. Decomposing and reconstructing processes are done by using 2-D DWT with basis functions that suit the physics of the property in question. Two-dimensional DWT for one level will result in one datum image (low–low-pass filter image) and three detail images (low–high, high–low, and high–high). The underlying relationship between the input variables and the output are learned by training an SVM on the datum images at the resolution of the output. The SVM is then applied on the detailed images to produce the detailed images of the output, which are needed to help downscale the output image to a higher resolution. In addition to being downscaled, the output image can be shifted ahead in time, providing a means for the algorithm to be used for forecasting. The algorithm has been applied on two case studies, one in Bondville, IL, where the results have been validated against AmeriFlux observations, and another in the Sevier River Basin, UT. Index Terms—Data management, learning systems, multiresolution techniques, water resources, wavelet transforms.

I. I NTRODUCTION

E

VAPOTRANSPIRATION (ET) is the total water loss from the land surface to the atmosphere, both evaporated from the surface and transpired by the vegetation. ET links the water Manuscript received March 8, 2007; revised November 19, 2007. This work was supported in part by the Utah Water Research Laboratory, by the Utah Center for Water Resources Research, and by the Utah State University Research Foundation. The work of E. Rosero was supported by the National Oceanic and Atmospheric Administration under Grant NA07OAR4310216. Y. H. Kaheil is with the Utah Water Research Laboratory and Department of Civil and Environmental Engineering, Utah State University, Logan, UT 84322 USA, and also with the International Research Institute for Climate and Society at Columbia University, New York, NY 10027 USA (e-mail: [email protected] ). E. Rosero is with the Department of Geological Sciences, John A. and Katherine G. Jackson School of Geosciences, University of Texas at Austin, Austin, TX 78712 USA (e-mail: [email protected]). M. K. Gill is with the Pacific Northwest National Laboratory, Richland, WA 99352 USA (e-mail: [email protected]). M. McKee and L. A. Bastidas are with the Utah Water Research Laboratory and Department of Civil and Environmental Engineering, Utah State University, Logan, UT 84322 USA. Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2008.919819

cycle and the energy balance together and influences carbon and nutrient cycles. For monitoring mass and energy exchanges between the land surface and the atmosphere, providing boundary conditions for weather and climate models, and for agricultural applications such as irrigation scheduling, detailed and accurate knowledge of ET is critical. Methods for quantifying ET based on micrometeorological point measurements (Bowen ratio and eddy correlation system) [1] or area-averaged estimates based on turbulence theory (scanning Raman lidar and large aperture scintillometer) [2] are still unable to provide spatially resolved ET at large scales, mostly because of the heterogeneity of the land surface and the dynamic nature of the heat transfer processes [1], [3], [4]. Driven by meteorological data, ET can be calculated by using conventional empirical or semiempirical methods such as the Penman–Monteith equation [5], the Priestly–Taylor equation [6], and the crop coefficient method (FAO 56) [7]. On a watershed scale, these process-based models provide good estimates of potential ET. However, in order to calculate the real rate of water consumption (actual ET), detailed knowledge of the complex spatial distribution of soil thermal and physical properties and their relationship with soil moisture, weather and environmental conditions, and crop variety and developing stage is needed [1], [3]. Such information is available from remote sensors over various temporal and spatial scales. Current satellite-based methods for estimating ET are indirect and rely on the assumption that a combination of vegetation and bare soil exists over a given area and that the evapotranspiring vegetation has a different thermal signature than the bare soil [8]. The following two categories of remote sensing methods for the estimation of ET exist: 1) residual methods of the energy budget that combine physical modules and empirical relationships to calculate ET by subtracting sensible heat flux from net radiation, such as two-source energy balance [9], surface energy balance algorithm for land [10], or surface energy balance system [11], and 2) methods that use visible and near infrared sensors to extract a vegetation index (VI) and the surface radiative temperature to make an inference about its corresponding skin temperature (Ts) [8]. By using the VI-Ts scatter plot, Nemani and Running [12] quantified changes in land-surface evaporation resistance. Moran et al. [13] explored the geometry of the VI-Ts and presented an algorithm to estimate crop water deficit index. Gillies et al. [14] used an inverted land-surface model (LSM) to extract moisture content information from the VI-Ts

0196-2892/$25.00 © 2008 IEEE

KAHEIL et al.: DOWNSCALING AND FORECASTING OF EVAPOTRANSPIRATION

2693

distributions. Building upon that line of research, Nishida et al. [15] presented the Moderate Resolution Imaging Spectroradiometer (MODIS) MOD16 algorithm. Via isolation of the end members of a VI over a scene with a large enough number of pixels to allow the estimation of the endpoint Ts signature, it globally derives an evaporative fraction at 1-km resolution every eight days [15]. In theory, two central advantages of models based on or ingesting remote sensing imagery over purely process-based models are the following: 1) broad spatial coverage and regular temporal sampling of satellite remote sensing and 2) waterconstraining variables (such as stomatal conductance), and their spatial and temporal distributions are no longer needed. Hence, at regional to continental scales, actual ET could be accurately predicted by remote sensing models. In a recent application, Yang et al. [16] investigated the use of an inductive machine learning technique or support vector machine (SVM) in modeling ET over the continental U.S. The SVM was trained on ET measurements of 25 AmeriFlux [17] sites to reproduce distributed ET, given remotely sensed land surface temperature, enhanced VI, land cover, and short wave radiation at an 8-km resolution and eight-day averages. They showed that the SVM outperformed neural networks and multiple regressions and that it could reproduce seasonal and regional patterns of ET. However, they concluded that VI at finer temporal resolution and short wave data at a finer scale could substantially enhance their estimates [16]. Courault et al. [3] note that obtaining satellite imagery at high spatial and temporal resolutions is the main problem for routine monitoring (i.e., daily to weekly) of surface energy fluxes. Furthermore, particularly during the growing season, frequent data acquisitions are needed for proper crop monitoring. Meteorological satellites (e.g., the Advanced Very High Resolution Radiometer and Geostationary Operational Environmental Satellite) offer the necessary frequency of observations (daily and hourly), but the spatial resolution (∼1–5 km) is too coarse to distinguish ET from each field (∼101–102 m) because such estimates would be representative of averages over ∼1–25-km2 areas, as pointed out by Kustas et al. [18]. On the other hand, pattern rich information in the visible and near infrared wavelengths used to compute VI is available at resolutions that are an order of magnitude higher, i.e., from the Land remote-sensing satellite Thematic Mapper (Landsat TM) (120 m in thermal infrared (TIR) for Landsat-5 and 60 m for Landsat-7) and the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) (15 m in three visible near-infrared bands and 90 m in the TIR band), but it is temporally sparse (every ∼16 days). Downscaling is required to transfer large-scale observations to a finer scale using a disaggregation procedure, which requires physical relationships that govern the processes to be used to decompose the aggregated value into a pattern at a smaller scale. McCabe et al. [19] point out that unlike attempts that made use of geostatistical techniques and interpolation procedures to determine spatial distributions of ET, the use of spatially dense remote sensing information has received relatively little attention and rather disappointing results (e.g., [20]) despite the great potential it offers. Moreover, they recognize that to expand our knowledge of changes in surface fluxes throughout time, an approach that

makes use of intermittent snapshots of the surface condition would prove very useful. Due to their coarse spatial resolution, current ET estimates are of little use for analyzing the spatial distribution of water use in irrigation networks. To better address the problem of reservoir and canal operations in irrigated river basins, the amount of water that should be released from an upstream reservoir and/or diverted into an irrigation canal to serve a certain field (proportional to the ET that the crop in this field is going to need in a few days) needs to be known in advance. The main goal behind this paper is to provide short-term ET forecasts at the farm scale. This paper presents an approach that uses available coarse resolution ET products, other satellite remote sensing products at different scales, and spatially distributed meteorological forecasts to provide downscaled ET forecasts at a finer spatial resolution. Forecasting ET can be a straightforward process if a single crop at a given crop stage is considered. When dealing with multiple crops in the same coarse ET footprint, this process becomes more complex. However, because the crop growth stage is partially reflected by multispectral images that are available at fine spatial resolutions, these images, along with other inputs, make this process plausible. Additionally, because the challenge is to provide those estimates shifted in time (which corresponds to the travel time of the water since it was released from a reservoir into a river and then into a canal; this could take a few days depending on the size of the river basin and the irrigation delivery system), this paper discusses the skill that the model has in forecasting ET. II. M ATERIALS AND M ETHODS A. SVMs Regression techniques are used to relate remotely sensed information with biophysical variables according to empirical models (e.g., linear, exponential, and polynomial) fitted to training samples. Because those models are too simple to adequately capture the relationship and the spatiotemporal variability of the quantity, neural-network-based techniques for nonlinear regression have been used. Several drawbacks of the neural networks, specifically the complexity of their implementation, risk of overfitting, and degraded performance with sparse data, have favored the use of SVMs [21], [22] in a variety of applications [23]–[26]. For instance, Bruzzone and Melgani [26] showed that SVMs achieve the same accuracy as neural networks but with reduced computational cost and that SVMs outperform neural networks with regard to their ability to generalize for sparse noisy inputs. Developed in the early 1990s for applications in classification, SVM for regression was developed by Vapnik [21] to extend his work. SVMs consist of a quadratic programming problem that can be efficiently solved and for which there is a guaranteed global extremum [22]. While working with statistical learning tools, the ultimate goal is to estimate a functional dependence f (x) between inputs {x1 , x2 , . . . , xl }, which are taken from x ∈ RK , and corresponding outputs {y1 , y2 , . . . , yl } with y ∈ R taken from a set L of independent and identically distributed

2694

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 9, SEPTEMBER 2008

observations. Hence, f (x) is found by minimizing the following regularized functional:  1 w2 + C (ξi + ξi∗ ) 2 i=1 L

Minimize

Subject to

 K  L    wj xji − b ≤ ε + ξi  yi −   j=1 i=1 K  L   wj xji + b − yi ≤ ε + ξi∗    j=1 i=1  ξi , ξi∗ ≥0

(1)

where f (x) = w, x + b with w ∈ X, b ∈ R f (x) =

K 

wj xj + b.

(2)

j=1

Also, w, x denotes the dot product of w and x, K is the input dimension, and b is the bias. The aforementioned formulation makes the solution sparse because it ignores errors that are less than ε. The first term in the objective function is the regularization term. It reduces the complexity of the function and makes it as flat as possible. The second term is the ε-insensitive loss function. Also, shown here are the slack variables ξi and ξi∗ that determine the degree to which samples will be penalized when errors are larger than ε. The value of the capacity constant C determines the tradeoff between the complexity of function f and the amount to which deviations that are larger than ε are tolerated. Equation (1) is usually solved in dual form by employing Lagrange multipliers. Writing (1) in dual form and solving by differentiating with respect to primal variables results in the following: Maximize : W (α∗ , α) = −ε

L 

(αi + αi∗ ) +

i=1



L 

yi (αi − αi∗ )

i=1

L 

  1 (αi − αi∗ ) αj − αj∗ k(xi , xj ) 2 i,j=1

subject to constraints L 

(αi∗ − αi ) = 0,

0 ≤ αi , αi∗ ≤ C

(3)

i=1

and the approximating function f (x) =

N 

(αi∗ − αi ) k(x, xi ) + b

(4)

i=1

where α∗ , α are Lagrange multipliers and k(x, xi ) is the kernel that replaces the dot products of the input vectors. The points/samples with nonvanishing coefficients are called support vectors, (xi in the aforementioned formulation); hence, the name is support vector machines. Basically, the support vectors are the sample points that support the “decision surface”

or “hyperplane” that best fits the data according to the criteria specified previously. Therefore, the SVM regression estimation steps are the following: 1) selection of a suitable kernel and kernel parameter γ; 2) specifying the ε parameter; and 3) specifying the capacity C. Kernels are one of the building blocks of SVMs. Unlike neural networks, which try to define complex combinations of functions in the input feature space, kernel methods map the complex nonlinear data into high-dimensional feature spaces (this is called the “kernel trick”) and then use linear functions to create decision hyperplanes or a regression function [22]. Choosing a kernel for the data projection in SVM is analogous to the problem of choosing a suitable architecture in neural networks [24], [27]. SVMs have many attractive features that make them an accurate, robust, and efficient method for nonlinear regression [26], including: 1) their strong connection with statistical learning theory [21] that guarantees better solutions in terms of generalization capabilities; 2) a simple architectural design (only three parameters); 3) guaranteed global optimality as opposed to gradient-based training methods for neural networks that might converge to a local optima; 4) a training process that is significantly faster than neural networks and requires, on average, less computational cost; 5) an ability to handle highdimensional input spaces that makes them particularly well suited to work with hyperspectral data; 6) effective measures to avoid overfitting by controlling the hyperplane margin; and 7) a reduced sensitivity to noisy inputs and sparse training samples [24]–[26]. B. Wavelet Decomposition Wavelet transformation decomposes a signal into a scale frequency space, allowing the determination of the relative contributions of the different spatial scales present within an image [28]. If a process exhibits self-similar scaling behavior, the wavelet coefficients obtained through a wavelet transform preserve that self-similarity. If a process does not show self-similar behavior, this methodology still permits analysis of the actual scaling behavior through a multiresolution framework. Wavelet decompositions are powerful tools in analyzing the variation in signal properties across different resolutions of geophysically distributed variables such as precipitation and soil moisture [29]–[31]. The wavelet transform overcomes the inability of the Fourier transform to represent a signal in the time and frequency domain at the same time by using a fully scalable modulated window that is shifted along the signal. For every position, the spectrum is calculated. After repeating the process, each time with a different window size, the result is a collection of time–frequency representations of the signal, all with different resolutions. Data are separated into multiresolution components; each is studied with a resolution that matches its scale [32]. The high-resolution components capture the fine scale features in the signal, whereas the low-resolution components capture the coarse scale features in the signal. Wavelet analysis represents any arbitrary (nonlinear) function by a linear combination of a set of wavelets or alternative basis functions, making them very suitable to be used both as an integration kernel for

KAHEIL et al.: DOWNSCALING AND FORECASTING OF EVAPOTRANSPIRATION

2695

analysis to extract information about the process and as a basis for representation or characterization of the processes [31]. A good introduction on multiresolution analysis is presented in [33]. 1) 1-D Wavelet Decomposition and Reconstruction: Following the notation in [31], the general form of the wavelet transform of a function f (t) with finite energy is defined as the√integral transform with a family of functions ψλ,t (u) = 1/ λ · ψ((u − t)/λ) and is given as follows: ∞ f (u)ψλ,t (u)du,

W f (λ, t) = −∞ ∞

= −∞

1 f (u) √ ψ λ



u−t λ

λ>0 du

(5) Fig. 1. DWT of a 1-D signal using Haar basis function.

where λ is a scale parameter, t is a location parameter, and the functions ψλ,t (u) are called wavelets. The inverse wavelet transfer is similarly given as follows: (u)dλdt (6) f (u) = W f (λ, t)ψλ,t (u) ψλ,t

is now the inverse of the basis function. The where wavelet transform can also be interpreted by using a time-scale transform given as follows:

∞ √ t λf (λu)ψ u − W f (λ, t) = du. λ

(7)

−∞

The mapping f (t) → f (λt) is achieved by contracting f (t) when λ > 1 and magnifying when λ < 1. This means that when the scale grows, a contracted version of the function is seen through a fixed size filter and vice versa. The scale factor λ therefore has the interpretation of scale in maps. The discrete wavelet transform (DWT) employs the convolution of the input data set with the low- and high-pass filters of the wavelet and the downsampling of the output [32]. DWT is the transformation of information from a fine scale to a coarser scale by extracting information that describes the fine scale variability and the coarse scale smoothness. It is given as follows: {sm } = [L]{sm+1 } {dm } = [H]{sm+1 }

(8)

where s represents smooth coefficients, d represents detail coefficients, m is the level, and L and H are the convolution matrices based on the wavelet basis function. Higher values of m signify finer scales of information. The complete wavelet transform is the process that recursively applies the aforementioned equation from the finest to coarsest scale. This describes a scale-to-scale extraction of the variability information at each scale. The smooth coefficients generated at each scale are used for the extraction in the next coarser scale. Because DWT signals are analyzed over a discrete set of scales, commonly powers of two, the most common implementation of the DWT is the dyadic-filter tree [32].

The inverse DWT (IDWT) is similarly implemented via a recursive recombination of smooth and detail information from the coarsest to finest scale; it is given as follows: {sm+1 } = [L] {sm } + [H] {dm }.

(9)

The matrices L and H are now the transposes of the L and H matrices, respectively. This idea is illustrated by using the flow diagram shown in Fig. 1. Let us say that a signal at resolution m = 3 is represented as s = {s1 , s2 , s3 , s4 }. The wavelet basis function in this case is a two-coefficient Haar function, which is one of the simplest wavelets, which resembles a step function. 2) 2-D Wavelet Decomposition and Reconstruction: The properties that make wavelets attractive for analyzing 1-D data sets also hold for images, matrices, and other 2-D data sets. In order to apply wavelets to images, an extension is made to the DWT to obtain the 2-D DWT. The 2-D DWT computes the coefficients of the 2-D wavelet series expansion for an m × n image Im,n . The 2-D DWT maps the image Im,n to an m × n matrix of wavelet coefficients wm,n . The 2-D DWT is implemented by an extension of Mallat’s forward pyramid algorithm [28]. This algorithm consists of passing the 1-D lowand high-pass filters through the rows of the image data set, while retaining every other column. The same filters are then passed through the columns of the resulting data set, while retaining every other row. Two sets of wavelet coefficients are obtained, approximation [low–low (LL)] and detail [high–low (HL), low–high (LH), and high–high (HH)] signals of the original data. In theory, the decomposition can be extended to p = log2 (N ) levels, where N is the length of the input signal. In practice, the maximum level also depends on the choice of the mother wavelet [34] and there is not a unique optimal level of decomposition for a particular resolution ratio [35]. Fig. 2 shows an output of a three-level 2-D DWT indicating the various wavelet output matrices. The matrix (LL3) located in the lower left-hand corner shows the smooth wavelet coefficients that approximate the original image. The matrices (LH3, LH2, and LH1) located along the xaxis capture the vertical edges of the image. The matrices (HL3, HL2, and HL1) located along the y-axis capture the horizontal

2696

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 9, SEPTEMBER 2008

Fig. 2. Wavelet output matrices for a three-level 2-D DWT.

edges. The matrices (HH3, HH2, and HH1) located along the diagonal capture the diagonal edges. Analogous to reversing the DWT to obtain an original signal from its wavelet coefficients, the 2-D DWT algorithm can be reversed to obtain an original image from its wavelet coefficients. This process is called the 2-D IDWT. DWTs have proven very effective for the analysis of remotely sensed imagery [36], mainly in data compression [37], edge detection, denoising, and multiscale data fusion [33]. C. Algorithm This downscaling/forecasting algorithm builds multiple relationships between inputs and outputs at different spatial scales. These relationships are then used to downscale and forecast the output at the finest scale. As a conservative assumption, this model assumes that the output image is available at the coarsest resolution. All inputs are upscaled to the coarsest output resolution. Upscaling is carried out by using 2-D DWT with the basis functions suiting the property in physical terms. Twodimensional wavelet decomposition for one level will result in one datum image (low–low-pass filter image or LL) and three detailed images (i.e., LH, HL, and HH). Once the inputs are available at the spatial resolution of the output, an SVM can be employed to learn the underlying physics between the inputs and the output using a random subset of pixels. The outcome of this SVM will be the LL image for this particular resolution. Because some inputs are available at higher spatial resolutions, these “leftover” inputs can still be upscaled to the output resolution and another SVM can be implemented to learn the relationship between these leftover inputs and the output. This SVM will be applied on the three high-pass components of these leftover inputs. The result of this SVM will be inherently biased due to the convolution processes performed at the decomposition step. This bias is linear, and it could be corrected. The linear bias corrector could be obtained at a coarser resolution where the three output detailed images are available. Once corrected for linear bias, the result of the SVM will be the three high-pass components. These can be used along with the datum image established in the first SVM to reconstruct the output at the next finer spatial resolution. The algorithm continues in this manner until all the inputs at higher spatial resolutions are consumed. Similar to the (downscaling only) explanation, the SVM can be trained against an output image ahead of time, assuming all inputs are at time t and the output is at time t + n, where n is the number of time steps ahead. This provides the framework

for the algorithm, serving not only as a means for downscaling but also for forecasting. The algorithm will have three parameters per SVM machine. Fig. 3 shows the scheme of this model. The triangles pointing up represent the 2-D wavelet decomposition operations, whereas those pointing down represent 2-D wavelet reconstruction operations. The dashed lines connecting the abbreviation “SVM” refer to “operational mode,” whereas the solid lines refer to “training mode.” Encircled images represent observed images. In this schema, there are two inputs and three different resolutions. One input is observed at fine resolution, whereas the other is observed at a medium resolution. The output is observed at a coarse resolution. D. Algorithmic Steps Notation: 1) O: denotes output. 2) I: denotes input. 3) {}: denotes set, e.g., {I} is a set of input images. 4) x: length of each side of the image, assuming the image to be square. 5) r: resolution index. Subscripts: 1) 0, i, i + 1, and n: denote the resolution of an image. Zero is the resolution of the observed output image, whereas n is the resolution of the highest resolution observed input. 2) |T =t : denotes the time step at which the image is produced. Superscripts: 1) LL: low–low-pass filter image. 2) HL, LH, and HH: high–low, low–high, and high–high filter images. 3) H-: all high-pass filter images, i.e., HL, LH, and HH. Format: 1) Bold: denotes observed. 2) Italic bold: denotes calculated. Objective: Downscale and forecast output image O0 |T =0 to On |T =τ using inputs {I}0,1,...,n |T =0 , where τ is the shift in time step. Given: 1) O0 |T =0,τ : a coarse image observed output at resolution 0 and time 0, τ . 2) {I}0,1,...,n |T =0 : a set of input images observed at different spatial resolutions and time 0. Each input in this set is only observed at one resolution. Also, these inputs are known to be significant predictors of the predictant O. Algorithm: 1) Initialization: set r = 0 2) Decomposition: Decompose {I}r+1,...,n |T =0 (1 to n − r levels) using 2-D DWT to produce {ILL }r |T =0 and {IH− }r |T =0 . 3) SVM training, calibration, and testing: a) Using a uniform sampler, sample three sets: {tr}, {ca}, and {ts}, where {tr} ⊂ {1, 2, . . . , x2r }, {ca} ⊂ {1, 2, . . . , x2r }, and {ts} ⊂ {1, 2, . . . , x2r }

KAHEIL et al.: DOWNSCALING AND FORECASTING OF EVAPOTRANSPIRATION

Fig. 3.

2697

Flowchart of the downscaling algorithm.

and {tr} ∩ {ca} = Φ, {tr} ∩ {ts} = Φ, and {ts} ∩ {ca} = Φ. b) Train a radial kernel SVM on the {tr}-indexed pixels of {ILL }r |T =0 versus the corresponding pixels of Or |T =0 or Or |T =τ (in case of forecasting), calibrate the SVM parameters on set {ca}, and test the performance on set {ts}. c) At resolution r = 0 only, concatenate {I}0 |T =0 with {ILL }r |T =0 , and using the result as input set and O0 |T =τ as output, train an SVM machine on {tr}indexed pixels, calibrate on {ca}, and test on set {ts}. This machine is used in operational mode to produce O0 |T =τ . 4) SVM application: a) Using the SVM machine from 3)-b), transform {IH− }r |T =0 onto the output space. Because the results of that process will be linearly biased, a linear bias remover is needed. This linear relationship can be obtained at the next coarser resolution where OH− r−1 are available. b) Using the linear bias remover, correct the output of the SVM machine from 3)-a) to obtain OH− r |T =τ . 5) Reconstruction: Use the inverse 2-D discrete wavelet analysis with Or |T =τ from step 3)-c) (as LL) along with OH− r |T =τ (as LH, HL, and HH) from step 4)-b) to produce Or+1 |T =τ . If r > 0, image Or |T =τ will be available from the previous iterations, i.e., step 3)-c) does not apply. 6) Repetition: Set r = r + 1 and repeat steps 2)–6). If r > n, halt the algorithm.

E. Operational Strategy To enable the forecasting functionality of the algorithm, multiple pretrained models from one or many previous crop season(s) should be employed. Given a current set of input images that reflect the crop growth stage, a distance metric could be implemented to a previous set of input images to know which pretrained model should be used. Given that weather forecasts are incorporated as inputs, the model should be able to account for the weather variability as well as for the new time step. Fig. 4 shows a flowchart of the model’s operational strategy. Forecasting ET/latent heat (LE) would not be possible if such a strategy is not implemented. III. C ASE S TUDIES The aforementioned algorithm is applied to two case studies. Because ET MODIS (MOD16) is not released yet, the first case study is applied toward downscaling and forecasting the photosynthesis (PNS) MODIS product, MOD17A2, in the Sevier River Basin, UT. The second case study deals with downscaling and forecasting the LE model output in Bondville, IL. The results in this case study have been validated with the AmeriFlux [16] data over different snapshots in the irrigation season. Because the MOD16 product [15] was not released at the time of this research, LE simulations from the Noah LSM [38] run offline at a 1-km resolution through the Land Information System (LIS) [39], which is the sister project of the Global Land Data Assimilation System (GLDAS) [40], were used.

2698

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 9, SEPTEMBER 2008

TABLE I INPUT IMAGES USED IN THE FIRST CASE STUDY: THE S EVIER R IVER B ASIN

Fig. 4. Flowchart of the model’s operational strategy. A pretrained model’s flowchart is shown in Fig. 3.

the basin and much of the real-time database utilized in this research is available at http://www.sevierriver.org [43]. In this application, the algorithm was used to downscale the PNS MODIS product (PNS-MOD17). The inputs used in this application are shown in Table I. The Landsat TM multispectral image, which is available at 15-m resolution, does not satisfy the relationship s × 2n = S, where s is the original spatial resolution (15 m), S is the upscaled spatial resolution (250 or 1000 m), and n is an integer. Therefore, the image had to be resampled at resolution s = 15.625 m, in which case, n will be four (for S = 250) or six (for S = 1000). PNS was forecasted eight days ahead using all inputs at, or upscaled to, a resolution of 1000 m. Then, it was downscaled from 1000 to 250 m using the inputs available at, or upscaled to, 250 m. Finally, it was downscaled to 15.625 m using the inputs that are available at a 15.625-m resolution. Fig. 6 shows the results of the model at the three benchmark spatial resolutions. Pending the release of MOD16, the algorithm could be applied to forecast ET in the Sevier River Basin. B. Bondville, IL

Fig. 5. Sevier River Basin, UT map.

A. Sevier River Basin, UT The Sevier River Basin, which is a closed system in rural south central Utah, is one of the state’s major drainages, encompassing 12.5% of the state’s total area. The Sevier River Basin, shown in Fig. 5, has five subwatersheds and is divided into two major divisions, which are the upper and lower basins, for the administration of water rights. Average annual precipitation ranges from 6.4 to 13.0 in in the valleys; the growing season ranges from 60 to 178 days [41], [42]. Most of the surface water runoff comes from snowmelt during the spring and early summer months. The primary use of water in the basin is for irrigation. The average annual amount of water diverted for cropland irrigation is 903 500 acre-feet. Of this amount, approximately 135 000 acre-feet is pumped from groundwater. The irrigation season in the basin generally extends from April to the end of October. About 40% of the diversions are return flows from upstream use [42]. More detailed information about

The Bondville site (40◦ 0 21.96 N, 88◦ 17 30.72 W, 300 masl) is located near Champaign, IL, in an agricultural farm that rotates annually between corn (C4) and soybeans (C3) [44] (Fig. 7). The site is sponsored by the Global Energy and Water Cycle Experiment (GEWEX) Americas Prediction Project, which is maintained by the National Oceanic and Atmospheric Administration (NOAA) Atmospheric Turbulence and Diffusion Division (T. Meyers, PI), and is one of the reference sites of the Coordinated Enhanced Observing Period initiative of GEWEX. The mean annual precipitation is 582 mm, and the mean temperature is 11.5 ◦ C. The climate regime is temperate continental. Soils are moderately well-drained silt loams. The height of the canopy is about 3.00 and 0.9 m for corn and soybeans, respectively. The Noah LSM [39] was used to provide LE simulations. The Noah LSM is a stand-alone 1-D column model, which can be executed in either coupled or uncoupled mode. Noah is the land module of the operational prediction system of the U.S. NOAA. It simulates soil moisture (both liquid and frozen), soil temperature, skin temperature, snow depth, snow water equivalent, canopy water content, and the land surface energy and water fluxes. The model applies finite-difference spatial discretization methods and a Crank–Nicholson time-integration

KAHEIL et al.: DOWNSCALING AND FORECASTING OF EVAPOTRANSPIRATION

2699

Fig. 6. Net PNS forecasting/downscaling units are in kilograms of carbon per meter square. (a) PNS at 1000-m resolution forecasted to t + 8 days, (b) PNS at 250-m resolution forecasted to t + 8 days, and (c) PNS at 15.625-m resolution forecasted to t + 8 days. TABLE II INPUT IMAGES USED IN THE SECOND CASE STUDY: BONDVILLE, IL

Fig. 7.

Bondville, IL map.

2700

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 9, SEPTEMBER 2008

Fig. 8. Robust variogram plots and index of dispersion statistics corresponding to input layers used for downscaling in Bondville case study. r and s2 /m are resolution and index of dispersion, respectively. Refer to Table II for more details.

scheme to numerically integrate the governing equations of the physical processes of the soil-vegetation-snowpack medium. The model was run offline at 0.010 resolution through the LIS [39], which is the sister project of GLDAS [40]. Because GLDAS ingests a variety of satellite and ground-based observational products to parameterize, force, and constrain LSMs, it presents an optimal way to showcase not only the capabilities of the methodology exposed in this paper but also its suitability of being incorporated into operational mode. The Noah LSM was run at 30-min time steps for the January 1 to December 1, 2006 period. Output was stored every 6 h at 00:00, 06:00, 12:00, and 18:00 coordinated universal time (UTC). The 13 × 13 domain of the simulation has its northeast corner in 40.07 N, 88.23 W and the southwest in 39.95 N, 88.34 W. The model output was then projected on the universal traverse mercator to match the input data and was later resampled at 960-m square cell resolution. The reprojected resampled domain’s northeast corner is 40.035 N, 88.24 W, and its southwest corner is 39.947 N, 88.35 W. The AmeriFlux tower is located at 40.01 N, 88.29 W. Among the data that the AmeriFlux tower provides is LE flux, which is used in this case study for validation purposes. LE flux is continuously monitored by using an open path infrared gas analyzer in a tower at a 10-m height. In our algorithm, the ASTER L1B data were used as inputs at two different time steps. The intention was to run the algorithm in different time steps throughout the irrigation season to show the algorithm performance over the different stages of the crop. However, we only ran the algorithm over two time steps where the ASTER images had 0% cloud cover. For more on the ASTER L1B data, the reader is referred to

http://edcdaac.usgs.gov/aster/ast_l1b.asp. Table II provides a summary of the inputs used in this application. It is worth noting that when forecasting LE across the irrigation season, distributed weather forecasts could be used as inputs. However, because only two time steps were used in this application, both of which were in the summer time, there was no significant need to incorporate such inputs. The algorithm uses two sets of input images corresponding to a defined current point in time and a preceding point with an offset number of days corresponding to the satellite return period. Assuming that the current time is t, the input sets of images used to run the algorithm correspond to t − r and t to forecast the LE image corresponding to t + f, where r and f are the satellite return period and number of forecasting days, respectively. The algorithm was used to downscale LE images from 960to 15-m resolution and forecast five days in advance for July 16 as the current input images and six days in advance for August 1 as the current input images. The model is set to forecast LE images at UTC 18:00. Fig. 8 shows robust variogram plots as well as the index of dispersion statistics for each layer of input used for downscaling. The algorithm bracketed the AmeriFlux LE readings at the two time steps near the tower location. The AmeriFlux LE reading on July 21 at UTC 18:00 could be interpolated to 586 W/m2 , whereas the model output records a mean value of 463 with a standard deviation of 146 within a 24-pixel neighborhood. On August 7, the AmeriFlux value was 86–153 and the model output was 96.3 with a standard deviation of

KAHEIL et al.: DOWNSCALING AND FORECASTING OF EVAPOTRANSPIRATION

2701

Fig. 9. LE forecasting (July 16–21)/downscaling (960–15 m) units are in watts per meter square. (a) LE at 15-m resolution forecasted to t + 5 days, (b) LE at 60-m resolution, and (c) LE at 240-m resolution. High-pass filter images at 240 m: (d) LH, (e) HL, and (f) HH.

43 within a 24-pixel neighborhood. The prediction capability of the SVM at the coarsest resolution varied from an rmse of 6.3 for July 21 to 3.3 for August 7. Figs. 9 and 10 show the results from the two time periods. The high-pass filter images are shown at a resolution of 240 m to illustrate an example of the ability of the model to perform the downscaling operation at two completely different scales of values.

In addition to forecasting at t + 5 in the month of July, we have evaluated the model forecasting capability by running the model to produce t + 1, t + 2, t + 3, and t + 10. Fig. 11 shows the model output corresponding to these dates. Fig. 12 shows the performance of the SVM machine forecasting at t + 10. Although the rmse reported for these two time series was ∼14 W/m2 , the pattern is still captured.

2702

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 9, SEPTEMBER 2008

Fig. 10. LE forecasting (August 1–7)/downscaling (960–15 m) units are in watts per meter square. (a) LE at 15-m resolution forecasted to t + 6 days, (b) LE at 60-m resolution, and (c) LE at 240-m resolution. High-pass filter images at 240 m: (d) LH, (e) HL, and (f) HH.

Fig. 13 shows the model performance corroborated by the AmeriFlux data, assuming its average error to be 25% [41]. It is seen in this figure that the model loses its capability to produce a good LE forecast as we move ahead in time. We think that this happens due to the fact that the images used for the downscaling no longer reflect the stage of the crops on the ground. Table III summarizes the model performance,

given the number of forecasting days ahead, the rmse of the SVM at the coarse resolution (datum image), and the model output versus the ground reading at the AmeriFlux tower location. For an overall view of the model output across different spatial resolutions, Fig. 14 shows the histograms of the LE output images. The mean and standard deviations are also shown.

KAHEIL et al.: DOWNSCALING AND FORECASTING OF EVAPOTRANSPIRATION

2703

Fig. 11. Model output at resolution of 15 m corresponding to (a) t+1 (July 17, 2006), (b) t+2 (July 18, 2006), (c) t+3 (July 19, 2006), and t+10 (July 26, 2006).

Notice that the mean does not change because the model uses the arithmetic averaging wavelet basis function. The standard deviation decreases as the resolution goes coarser. We realize that validating at one point is by no means sufficient to make a statement about the model performance in terms of capturing the spatial variability. Fig. 15, however, shows that a robust variogram plot of the output at a given spatial resolution looks similar to those shown in Fig. 8 at the same resolution. We also believe that due to the great difference in magnitude between the pixel values of LE and the input layers, a comparison between the index of dispersion statistics is not as meaningful. In terms of the model forecasting ability, we believe that the two sets of fine input images corresponding to two points in time reflect the crop stage; thus, the farther we move away from

the most recent point in time, the less informative these images become. IV. D ISCUSSION AND C ONCLUSION Both at the farm and the irrigation project level, an accurate estimation of ET and detailed knowledge of its spatial variability is required for adequate water management and crop monitoring. SVMs have been proposed as accurate, efficient, and robust methods for nonlinear regression, particularly suited for large input spaces [25], [26]. Among the attractive features of SVMs are better generalization capabilities (performance on previously unseen data), low sensitivity to small training samples and noisy data [26], and effective avoidance of overfitting [25].

2704

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 9, SEPTEMBER 2008

Fig. 12. SVM output versus testing set pixel values at coarse resolution (960 m) forecasting at t + 10. (a) Time series plot. (b) Scatter plot. TABLE III INPUT IMAGES USED IN THE MODEL FORECASTING CAPABILITY FOR JULY 2006

Fig. 13. Model forecasting ability versus the number of forecasting days ahead.

Decomposing and reconstructing processes are done by using 2-D DWT with basis functions that suit the physics of the property in question. Two-dimensional DWT for one level results in one datum image (LL) and three detailed images (LH, HL, and HH). The underlying physics between the input variables and the output are learned by the SVM. In this application, a synthetic algorithm of wavelets and SVMs was devised to simultaneously perform downscaling and forecasting. The algorithm ingests the available satellite data at different resolutions to produce an output at the finest spatial resolution. Due to its flexible design, this algorithm could be implemented on a variety of applications. The algorithm is applied to two different case studies and validated for the second one at the AmeriFlux tower location. As pointed out

by Anderson et al. [46], the evaluation of the results of remote sensing modeling is hampered by “large spatial and temporal gaps due to cloud cover and infrequent image availability as governed by the satellite overpass schedule.” For this reason, a comprehensive corroboration of the model results in both the temporal and spatial scales has not been made possible. Further testing of the forecasting capability of the algorithm shows that it is only suited for short-term forecasting of ET. Nonetheless, the results are encouraging and more applications are being sought as more satellite products are available. ACKNOWLEDGMENT The authors would like to thank the Utah Water Research Laboratory, the Utah Center for Water Resources Research, and the Utah State University Research Foundation for their support of this research. The authors would also like to thank Dr. F. Yang at the University of Wisconsin and Dr. M. White at Utah State University for their help in providing the data, as well as Dr. L. G. de Goncalves at the NASA Goddard Space Flight Center for his constructive advice and help in running the LIS. The authors would also like to thank Dr. A. F. Khalil at Columbia University, L. E. Gulden at the University of Texas

KAHEIL et al.: DOWNSCALING AND FORECASTING OF EVAPOTRANSPIRATION

2705

Fig. 14. Output histograms at different spatial resolutions.

Fig. 15. Robust variogram plots and index of dispersion statistics corresponding to output at different resolutions. r and s2 /m are resolution and index of dispersion, respectively. Compare with Fig. 8.

2706

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 46, NO. 9, SEPTEMBER 2008

at Austin, and Dr. B. T. Neilson and Dr. C. M. U. Neale, both at Utah State University, for their insightful comments that helped to improve this paper. The authors would also like to thank the anonymous reviewers of their original manuscript for their extremely valuable comments. R EFERENCES [1] W. J. Shuttleworth, “Putting the ‘vap’ into evaporation,” Hydrol. Earth Syst. Sci., vol. 11, no. 1, pp. 210–244, Jan. 2007. [2] R. J. Hill, R. A. Bohlander, S. F. Clifford, R. W. McMillan, J. T. Priestley, and W. P. Schoenfeld, “Turbulence-induced millimeter-wave scintillation compared with micrometeorological measurements,” IEEE Trans. Geosci. Remote Sens., vol. 26, no. 3, pp. 330–342, May 1988. [3] D. Courault, B. Seguin, and A. Olioso, “Review on estimation of evapotranspiration from remote sensing data: From empirical to numerical modeling approaches,” Irrig. Drain. Syst., vol. 19, no. 3/4, pp. 223–249, Nov. 2005. [4] J. S. Famiglietti and E. F. Wood, “Effects of spatial variability and scale on areally averaged evapotranspiration,” Water Resour. Res., vol. 31, no. 3, pp. 699–712, 1995. [5] J. L. Monteith, “Evaporation and environment,” in Proc. 19th Symp. Soc. Exp. Biol., 1965, vol. 19, pp. 205–234. [6] C. H. B. Priestley and R. J. Taylor, “On the assessment of surface heat flux and evaporation using large-scale parameters,” Mon. Weather Rev., vol. 100, no. 2, pp. 81–92, 1972. [7] R. G. Allen, L. S. Pereira, D. Raes, and M. Smith, Crop Evapotranspiration—Guidelines for Computing Crop Water Requirements, vol. 56. Rome, Italy: FAO, 1998. [8] L. Jiang and S. Islam, “An intercomparison of regional latent heat flux estimation using remote sensing data,” Int. J. Remote Sens., vol. 24, no. 11, pp. 2221–2236, Jun. 2003. [9] J. M. Norman, W. P. Kustas, and K. S. Humes, “Source approach for estimating soil and vegetation energy fluxes in observations of directional radiometric surface temperature,” Agric. For. Meteorol., vol. 77, no. 3/4, pp. 263–293, Dec. 1995. [10] W. G. M. Bastiaanssen, M. Menenti, R. A. Feddes, and A. A. M. Holtslag, “A remote sensing surface energy balance algorithm for land (SEBAL). 1. Formulation,” J. Hydrol., vol. 212/213, pp. 198–212, Dec. 1998. [11] Z. Su, “The surface energy balance system (SEBS) for estimation of turbulent heat fluxes,” Hydrol. Earth Syst. Sci., vol. 6, no. 1, pp. 85–99, 2002. [12] R. R. Nemani and S. W. Running, “Estimation of regional surface resistance to evapotranspiration from NDVI and thermal-IR AVHRR data,” J. Appl. Meteorol., vol. 28, no. 4, pp. 276–284, Apr. 1989. [13] M. Moran, T. Clarke, U. Inoue, and A. Vidal, “Estimating crop water deficit using the relation between surface-air temperature and spectral vegetation index,” Remote Sens. Environ., vol. 49, no. 3, pp. 246–263, Sep. 1994. [14] R. R. Gillies, T. N. Carlson, J. Cui, W. P. Kustas, and K. S. Humes, “A verification of the ‘triangle’ method for obtaining surface soil water content and energy fluxes from remote measurements of the normalized diference vegetation index (NDVI) and surface radiant temperature,” Int. J. Remote Sens., vol. 18, no. 15, pp. 3145–3166, Oct. 1997. [15] K. Nishida, R. R. Nemani, J. M. Glassy, and S. W. Running, “Development of an evapotranspiration index from Aqua/MODIS for monitoring surface moisture status,” IEEE Trans. Geosci. Remote Sens., vol. 41, no. 2, pp. 493–501, Feb. 2003. [16] F. Yang, M. White, A. Michaelis, K. Ichii, H. Hashimoto, P. Votava, A. X. Zhu, and R. R. Nemani, “Prediction of continental-scale evapotranspiration by combining MODIS and AmeriFlux data through support vector machine,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 11, pp. 3452–3461, Nov. 2006. [17] D. Baldocchi, E. Falge, L. Gu, R. Olson, D. Hollinger, S. Running, P. Anthoni, C. Bernhofer, K. Davis, R. Evans, J. Fuentes, A. Goldstein, G. Katul, B. E. Law, X. Lee, Y. Malhi, T. Meyers, W. Munger, W. Oechel, K. T. Paw, U. K. Pilegaard, H. P. Schmid, R. Valentini, S. Verma, T. Vesala, K. Wilson, and S. Wofsy, “FLUXNET: A new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities,” Bull. Amer. Meteorol. Soc., vol. 82, no. 11, pp. 2415–2434, Nov. 2001. [18] W. P. Kustas, J. N. Norman, M. C. Anderson, and A. N. French, “Estimating subpixel surface temperatures and energy fluxes from the vegetation index-radiometric temperature relationship,” Remote Sens. Environ., vol. 85, no. 4, pp. 429–440, Jun. 2003.

[19] M. F. McCabe, J. D. Kalma, and S. W. Franks, “Spatial and temporal patterns of land surface fluxes from remotely sensed surface temperatures within an uncertainty modeling framework,” Hydrol. Earth Syst. Sci., vol. 9, no. 5, pp. 467–480, 2005. [20] Y. Chemin and K. Honda, “Spatiotemporal fusion of rice actual evapotranspiration with genetic algorithms and an agrohydrological model,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 11, pp. 3462–3469, Nov. 2006. [21] V. N. Vapnik, The Nature of Statistical Learning Theory. New York: Springer-Verlag, 1998. [22] B. Scholkopf and A. Smola, Learning With Kernels. Support Vector Machines, Reguralization, Optimization and Beyond. Cambridge, MA: MIT Press, 2002. [23] M. Brown, H. G. Lewis, and S. R. Gunn, “Linear spectral mixture models and support vector machines for remote sensing,” IEEE Trans. Geosci. Remote Sens., vol. 38, no. 5, pp. 2346–2360, Sep. 2000. [24] E. J. Kwiatkowska and G. Fargion, “Application of machine-learning techniques toward the creation of a consistent and calibrated global chlorophyll concentration baseline dataset using remotely sensed ocean color data,” IEEE Trans. Geosci. Remote Sens., vol. 41, no. 12, pp. 2844–2860, Dec. 2003. [25] G. Camps-Valls, L. Gomez-Chova, J. Calpe, E. Soria, J. D. Martin, L. Alonso, and J. Moreno, “Robust support vector method for hyperspectral data classification and knowledge discovery,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 7, pp. 1530–1542, Jul. 2004. [26] L. Bruzzone and F. Melgani, “Robust multiple estimator systems for the analysis of biophysical parameters from remotely sensed data,” IEEE Trans. Geosci. Remote Sens., vol. 43, no. 1, pp. 159–174, Jan. 2005. [27] T. Asefa, M. W. Kemblowski, G. Urroz, M. McKee, and A. Khalil, “Support vectors-based groundwater head observation networks design,” Water Resour. Res., vol. 40, no. 11, W11 509, Nov. 2004. DOI: 0.1029/ 2004WR003304. [28] S. Mallat, “A theory for multiresolution signal decomposition: The wavelet representation,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 11, no. 7, pp. 674–693, Jul. 1989. [29] P. Kumar and E. Foufoula-Georgiou, “A multicomponent decomposition of spatial rainfall fields. 1. Segregation of large- and small-scale features using wavelet transforms,” Water Resour. Res., vol. 29, no. 8, pp. 2515– 2532, 1993. [30] P. Kumar and E. Foufoula-Georgiou, “Wavelet analysis for geophysical applications,” Rev. Geophys., vol. 35, no. 4, pp. 385–412, 1997. [31] Z. Hu, Y. Chen, and S. Islam, “Multiscaling properties of soil moisture images and decomposition of large- and small-scale features using wavelet transforms,” Int. J. Remote Sens., vol. 19, no. 13, pp. 2451–2467, Sep. 1998. [32] I. Daubechies, Ten Lectures on Wavelets, ser. CBMS-NSF Regional Conference Series Applied Mathematics, vol. 61. Philadelphia, PA: SIAM, 1992. [33] B. Aiazzi, L. Alparone, S. Baronti, and A. Garzelli, “Context-driven fusion of high spatial and spectral resolution images based on oversampled multiresolution analysis,” IEEE Trans. Geosci. Remote Sens., vol. 40, no. 10, pp. 2300–2312, Oct. 2002. [34] L. Mann Bruce, C. Morgan, and S. Larsen, “Automated detection of subpixel hyperspectral targets with continuous and discrete wavelet transforms,” IEEE Trans. Geosci. Remote Sens., vol. 39, no. 10, pp. 2217– 2226, Oct. 2001. [35] P. S. Pradhan, R. L. King, N. H. Younan, and D. W. Holcomb, “Estimation of the number of decomposition levels for a wavelet-based multiresolution multisensor image fusion,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 12, pp. 3674–3686, Dec. 2006. [36] N. K. Howard, “Multiscale analysis of landscape data sets from northern Ghana: Wavelets and pattern metrics,” Ph.D. dissertation, Univ. Bonn, Bonn, Germany, 2005. Ecology and Development Series No. 31. [37] P. Luigi-Dragotti, G. Poggi, and A. R. P. Ragozini, “Compression of multispectral images by three-dimensional SPIHT algorithm,” IEEE Trans. Geosci. Remote Sens., vol. 38, no. 1, pp. 416–428, Jan. 2000. [38] M. B. Ek, K. E. Mitchell, Y. Lin, E. Rogers, P. Grunmann, V. Koren, G. Gayno, and J. D. Tarpley, “Implementation of Noah land surface model advances in the national centers for environmental prediction operational mesoscale Eta model,” J. Geophys. Res., vol. 108, no. D22, 8851, Nov. 2003. DOI: 10.1029/2002JD003296. [39] S. V. Kumar, C. D. Peters-Lidard, Y. Tian, P. R. Houser, J. Geiger, S. Olden, L. Lighty, J. L. Eastman, B. Doty, P. Dirmeyer, J. Adams, K. Mitchell, E. F. Wood, and J. Sheffield, “Land information system: An interoperable framework for high resolution land surface modeling,” Environ. Model. Softw., vol. 21, no. 10, pp. 1402–1415, Oct. 2006. [40] M. Rodell, P. R. Houser, U. Jambor, J. Gottschalk, K. Mitchell, C.-J. Meng, K. Arsenault, B. Consgrove, J. Radakovich, M. Bosilovich, J. K. Entin, J. P. Walker, D. Lohmann, and D. Toll, “The global land data

KAHEIL et al.: DOWNSCALING AND FORECASTING OF EVAPOTRANSPIRATION

2707

assimilation system,” Bull. Amer. Meteorol. Soc., vol. 85, no. 3, pp. 381– 394, Mar. 2004. K. Zeller, G. Zimmerman, T. Henn, E. Donev, D. Denny, and J. Welker, “Analysis of inadvertent microprocessor lag time on Eddy covariance results,” J. Appl. Meteorol., vol. 40, no. 9, pp. 1640–1646, Sep. 2001. DOI: 10.1175/1520-0450. Utah Board of Water Resources, Utah’s Water Resources: Planning for the Future, Div. Water Resources Publication, Salt Lake City, UT, 2001. B. Berger, R. Hansen, and A. Hilton, “Using the world-wide-web as a support system to enhance water management,” in Proc. 18th ICID Congr., 53rd IEC Meeting, Montreál, QC, Canada, 2002. A. Khalil, M. McKee, M. Kemblowski, and T. Asefa, “Sparse Bayesian learning machine for real-time management of reservoir releases,” Water Resour. Res., vol. 41, no. 11, W11 401, Nov. 2005. DOI: 10.1029/ 2004WR003891. T. P. Meyers and S. E. Hollinger, “An assessment of storage terms in the surface energy balance of maize and soybean,” Agric. For. Meteorol., vol. 125, no. 1/2, pp. 105–115, Sep. 2004. M. C. Anderson, J. M. Norman, J. R. Mecikalski, J. A. Otkin, and W. P. Kustas, “A climatological study of evapotranspiration and moisture stress across the continental United States based on thermal remote sensing: 1. Model formulation,” J. Geophys. Res., vol. 112, D10 117, May 2007. DOI: 10.1029/2006JD007506.

M. Kashif Gill received the B.S. degree in civil engineering from the University of Engineering and Technology, Lahore, Pakistan, in 2002 and the M.S. and Ph.D. degrees in civil engineering from Utah State University, Logan, in 2004 and 2006, respectively. He is currently a Postdoctoral Research Associate with the Pacific Northwest National Laboratory, Richland, WA. His research interests involve stateof-the-art data-driven tools, lumped and distributed hydrologic modeling, climate impacts, evolutionary computation, data assimilation/fusion, etc.

[41]

[42] [43] [44]

[45] [46]

Yasir H. Kaheil received the Ph.D. degree in civil and environmental engineering from the Utah State University, Logan, in 2007, where he focused on developing Bayesian and evolutionary model calibration algorithms as well as scale reconciliation, data mining, and data fusion techniques for applications in hydrology. He is currently a Postdoctoral Research Scientist at the International Research Institute for Climate and Society at Columbia University, New York, NY, where he is focusing on hydrologic modeling and downscaling of climate model output for water resources planning and management. He is also with the Utah Water Research Laboratory and Department of Civil and Environmental Engineering, Utah State University.

Enrique Rosero received the B.S. degree in civil engineering from Escuela Politecnica Nacional, Quito, Ecuador, in 2001 and the M.S. degree in environmental fluid mechanics from the University of Karlsruhe, Karlsruhe, Germany, in 2003. He is currently working toward the Ph.D. degree in geosciences in the Department of Geological Sciences, John A. and Katherine G. Jackson School of Geosciences, University of Texas, Austin, under the supervision of Dr. Z.-L. Yang. He combines data assimilation and parameter estimation techniques within a multimodel ensemble framework to improve probabilistic flood forecasting. His research interests include high-resolution land-surface modeling, uncertainty estimation, statistical pattern recognition, and machine learning. Mr. Rosero is the recipient of the National Weather Service/Office of Hydrologic Development Graduate Fellowship.

Mac McKee received the B.S. degree in philosophy, the M.S. degree in systems ecology, and the Ph.D. degree in civil and environmental engineering (water resources and hydrology) from Utah State University, Logan, in 1972, 1979, and 1985, respectively. He is currently the Director of the Utah Water Research Laboratory, Utah State University (USU). He teaches graduate and undergraduate courses in water resources systems analysis, water resources engineering, and water resources planning and management in the Department of Civil and Environmental Engineering, USU. He has been active in domestic and international water resources planning and management projects for the past 33 years. His diverse experience has included assignments in flood control planning and design in the Philippines, the development of an environmental baseline document for Uzbekistan, water quality management in Samoa, integrated river basin planning and management in India, and the development of a comprehensive water resource master plan and integrated water management plans for the West Bank and Gaza Strip. On many assignments, he had the chief responsibility for providing decision-relevant information on contentious water development issues at local, interstate, and international levels. Examples include the design of plans for managing interstate rivers in India and the formulation and implementation of a planning framework for the development of Palestinian water resources in the West Bank and Gaza Strip. His current research is focused on the improvement of irrigation canal and reservoir operation efficiency through the use of statistical learning theory and remote sensing.

Luis A. Bastidas received the B.S. degree in civil engineering from Escuela Politécnica Nacional, Quito, Ecuador, in 1985, the M.Sc. degree in hydraulic engineering from the University of Newcastle upon Tyne, Newcastle, U.K., in 1990, and the Ph.D. degree in hydrology from the University of Arizona, Tucson, in 1998. He was a Research Faculty with the University of Arizona from 1998 to 2002. Since 2002, he has been with the Department of Civil and Environmental Engineering and the Utah Water Research Laboratory, Utah State University, Logan. His research includes hydrological and hydrometeorological modeling, system and statistical techniques applied to hydrologic modeling, remote sensing in hydrological modeling, inverse modeling, model uncertainty analysis, and model performance evaluation. Dr. Bastidas is in various national and international science committees.