A multi-spectral spatial convolution approach of rainfall
forecasting using weather satellite imagery
Chiang Wei
a, Wei-Chun Hung
b, Ke-Sheng Cheng
b,*aExperimental Forest Administration, National Taiwan University, Nantou, Taiwan, ROC
bDepartment of Bioenvironmental Systems Engineering, National Taiwan University, Taipei 106, Taiwan, ROC
Received 4 August 2004; received in revised form 4 August 2005; accepted 4 August 2005
Abstract
Flood forecasting has long been a major topic of hydrologic research. Recent events and studies indicate that the success of flood forecasting in Taiwan depends heavily on the accuracy of real-time rainfall forecasting. In this study, we demonstrate a multi-spectral spatial convolution approach for real-time rainfall forecasting using geostationary weather satellite images. The approach incorporates cloud-top temperatures of three infrared channels in a spatial convolution context. It not only characterizes the input–output relation-ship between cloud-top temperature and rainfall at the ground level, but also is more consistent with physical and remote sensing prin-ciples than single-pixel matches. Point rainfall measurements at raingauge sites are up-scaled to pixel-average-rainfall by block kriging, then related to multi-spectral cloud-top temperatures derived from Geostationary Meteorological Satellite images by spatial convolution. The kernel function of the multispectral spatial convolution equation is solved by the least squares method. Through a cross-validation procedure, we demonstrate that the proposed approach is capable of achieving high accuracy for 1- to 3-h-lead pixel-average-rainfall forecasting.
2005 COSPAR. Published by Elsevier Ltd. All rights reserved.
Keywords: Rainfall forecasting; Remote sensing; Spatial convolution
1. Introduction
Taiwan is located at the center of the western Pacific Rim and is particularly vulnerable to threat by typhoons. On average, there are 3.5 typhoons passing through Tai-wan annually. Approximately 79% of these typhoons occur in the period from July to September. Typhoons often draw huge amounts and high intensity rainfall and may sometimes result in high casualties and severe property damage. For example, a ferocious typhoon Nari passed through and ravaged Northern Taiwan in September, 2002. During its passage, more than 700 mm of rainfall was recorded near the capital city Taipei. Overbank flood flow caused extensive inundation and tremendous property damage in the city.
The water resources agency (WRA) has a flood forecast-ing and warnforecast-ing system in operation. The system comprises three major components: a rainfall forecasting system using weather radar, a hydrological rainfall-runoff model, and a three-dimensional inundation model. The three compo-nents work in series to yield a projected inundation map of the study area. Among the three components, rainfall forecasting is considered the most difficult, owing to high spatial and temporal variations of rainfall. Furthermore, the extremely rugged terrain in upstream mountainous areas, which causes beam blockage and ground clutter, of-ten sabotage successful utilization of radar images for rain-fall forecasting. In contrast, the weather satellite provides images of cloud-top temperature and water vapor, which are not affected by land surface features.
The success of flood forecasting in Taiwan depends heavily on accurate real-time rainfall forecasting since time-of-concentrations of major watersheds in Taiwan
0273-1177/$30 2005 COSPAR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.asr.2005.08.017
*
Corresponding author.
E-mail address:[email protected](K.-S. Cheng).
are relatively short (usually a few hours only) and flow at the watershed outlet quickly responds to rainfall in the up-stream area. Unless accurate rainfall forecasts can be made, there will not be enough time to issue flood warnings and evacuate people who live in flood-prone areas. There-fore, the objective of this study is to develop a real-time rainfall forecasting approach which can be integrated into a flood forecasting system, using weather satellite imagery. Estimation of rainfall using weather satellite images has a long history of nearly four decades. Thermal infrared images are most widely used for such applications (Scofield, 1987; Griffith et al., 1978; Woodley et al., 1979; Negri and Adler, 1987; DÕsouza et al., 1990; Martin et al., 1990; Ba and Nicholson, 1998; Ba and Gruber, 2001; Mishra and Sharma, 2001). Techniques coupling multi-sensors images (for example, thermal infrared and space-borne microwave images) have also been developed (Kummerow and Giglio, 1995; Turk et al., 1998; Berg et al., 1999; Xu et al., 1999; Kazumasa and Liu, 2000; Pra-bhakara et al., 2000; Miller et al., 2001; Todd et al., 2001; Ferreira et al., 2001; Toshiaki et al., 2001). Since rainfall forecasting must be conducted in real-time, many multi-sensor techniques that utilize space-borne microwave images are not applicable and we develop a rainfall fore-casting algorithm using only geostationary weather satellite images. It is also worthy to mention that if the space-borne microwave images are used only to calibrate the thermal infrared data, then the delay of the microwave data is not an issue in applying the combined IR/microwave algo-rithms in real time, unless other data such as rain gauge data are also considered. Examples of such algorithms in-cludeTurk et al. (1998) and Sorooshian et al. (2000). 2. Study area and data
The Danshui River watershed, which covers 2700 km2 drainage area in northern Taiwan, was selected for this study. It is composed of three major tributaries – the Hsindien River (900 km2), Dahan River (1200 km2) and Keelung River (600 km2). Terrain elevation in the wa-tershed ranges from near mean sea level at the outlet to over 3500 m in headwater mountainous area.Fig. 1shows the location and drainage systems of the Danshui River watershed. GMS images and hourly rainfall data from a network of 37 raingauge stations (see Fig. 1) for 12 ty-phoon events with durations ranging from 10 to 57 h (see Table 1) were collected. GMS has three thermal infrared channels (IR1–IR3) and one panchromatic (PAN) channel. The spectral ranges and spatial resolu-tions of GMS images are shown inTable 2. Radiances re-ceived at the satellite sensors in the spectral range of IR1 and IR2 channels (10.5–12.5 lm) are mainly emitted by clouds and are dependent on the cloud-top temperature. Radiances received at the sensors in the IR3 and PAN spectral ranges, respectively, characterize water vapor in the upper atmosphere and visible albedo of the cloud (Conway, 1997).
3. Methodology
Our methodology of real-time rainfall forecasting is composed of three major components: (1) estimation of pixel-average-rainfall (PAR) by block kriging, (2) forecast-ing of pixel-average-rainfall usforecast-ing weather satellite imagery by a multi-spectral spatial convolution (MSSC) approach, and (3) updating MSSC kernel function by Kalman filter-ing. Because of the page limit, this paper presents only
Fig. 1. The study area and raingauge locations.
Table 1
Duration of typhoon events used in this study
Name Nominal date Duration (h)
Winnie 19/08/1997 18 Amber 29/08/1997 24 Zeb 15/10/1998 35 Babs 25/10/1998 48 Yanni 26/09/1998 57 Maggie 08/06/1999 10 Dan 03/10/1999 16 Kaitak 09/07/2000 20 Bilis 22/08/2000 25 Prapiroon 28/08/2000 28 Bopha 09/09/2000 21 Xangsane 31/10/2000 41 Table 2
Spectral range and spatial resolution of GMS images Channel Spectral range
(lm)
Spatial resolution (km)
Related information
IR1 10.5–11.5 5 Cloud-top temperature
IR2 11.5–12.5 5 Cloud-top temperature
IR3 6.5–7.0 5 Water vapor
results of the first two components. Implementation of the Kalman filtering technique for MSSC kernel function up-date and real-time rainfall forecasting will be presented in a future paper.
3.1. Estimation of pixel-average-rainfall by block kriging An important difference between rainfall measured by raingauges and estimated from GMS satellite images is their spatial scales, i.e. the spatial area over which rain-fall amounts are collected. Rainrain-fall measurements by raingauges are considered point measurements, whereas GMS estimates of rainfall rate are areal average (5 km· 5 km) rainfall. Also, satellite cloud images are ac-quired instantaneously while rain gauges give accumu-lated amounts of rainfall that are discrete in time. This further complicates the efforts to relate the two. In order to establish a relationship between ground-based rainfall measurements and satellite-based rainfall estimates, we first discretized the study area into a grid mesh with 5-km intervals so that each grid cell corresponds to a pixel on GMS satellite images. Then, 1- and 3-h cumulative pixel-average-rainfall (hereafter referred to as 1-h-PAR and 3-h-PAR) corresponding to raingauges were esti-mated by block kriging using point rainfall measured at 37 raingauges.
Consider point rainfall at location x as a random vari-able denoted by Z(x) and spatial variation of point rainfall as a random field {ZX(x), x2 X} where X represents the
spatial extent of rainfall field. Rainfall measurements {z(xi), i = 1, 2, . . . , n} are observed by raingauges at
loca-tions xi, i = 1, 2, . . . , n, and average rainfall within an area
of V centered at x0, denoted by zV(x0), can be estimated by
a linear estimator zVðx0Þ ¼
Xn i¼1
kizðxiÞ. ð1Þ
Block kriging estimator z
Vðx0Þ is a linear, unbiased
esti-mator of zV(x0) and has a minimum variance of estimation
error. The following block kriging system of equations is used to solve the linear weights kiÕs:
Pn b¼1 kbcðxa xbÞ þ l ¼ cðV ; xaÞ ða ¼ 1; 2; . . . ; nÞ; Pn b¼1 kb¼ 1:0; 8 > > > < > > > : ð2Þ
where l is a Lagrange multiplier, cðV ; xaÞ ¼V1RVcðx; xaÞ dx,
and c(xi,xj) is the variogram of the random field ZX(x) and
is defined as cðxi; xjÞ ¼ cðjxi xjjÞ ¼12Ef½ZðxiÞ ZðxjÞ2g.
Variogram is a function that characterizes the spatial variation structure of a random field and can be estimated using observed data {z(xi), i = 1, 2, . . . , n},
c _ ðhÞ ¼ 1 2jN ðhÞj X NðhÞ ½zðxiÞ zðxjÞ 2 ; h >0; ð3Þ where N(h)”{(xi xj):jxi xjj = h; i, j = 1, 2, . . . , n}. The
function c_ðhÞ is often termed the experimental variogram. For detailed descriptions of spatial estimation by block kri-ging readers are referred toJournel and Huijbregts (1978) and Chile`s and Delfiner (1999).
3.2. Multi-spectral spatial convolution for pixel-average-rainfall forecasting
Both cloud-top temperature (CTT) observed by weather satellites and rainfall on the ground surface exhibit signifi-cant spatial and temporal variations and can be character-ized using the concept of random fields. Let T(x, y) and R(x, y), respectively, represent the random fields of cloud-top temperature and ground-surface rainfall at image loca-tion (x, y). Transformaloca-tion of cloud-top temperature to rainfall can be considered as a linear system with T(x, y) and R(x, y) being the input and output functions, respec-tively. The following equation of spatial convolution inte-gral is introduced to relate the two random fields:
Rðx; yÞ ¼ Z xþ‘ x‘ Z yþ‘ y‘ Tðx0y0Þf ðx x0; y y0Þ dx0dy0; ð4Þ
where the extent of convolution (also seen as the extent of influence by unit input) is (2‘ + 1) by (2‘ + 1), and f is the kernel function. The above integration equation only works for continuous random fields; however, both CTT and PAR are discrete random fields with 5 km· 5 km sup-port and the following discrete form of spatial convolution should be used: Rðx; yÞ ¼ X xþ‘ x0¼x‘ Xyþ‘ y0¼y‘ Tðx0; y0Þf ðx x0; y y0Þ ¼X N i¼1
Tði; x; yÞf ði; x; yÞ; N ¼ ð2‘ þ 1Þ2. ð5Þ R(x, y) represents pixel-average rainfall corresponding to a pixel at image location (x, y). It is calculated by spatial con-volution of cloud-top temperature of surrounding pixels T(x0, y0) (nine pixels in our study using ‘ = 1 and
N = (2‘ + 1)2). In order to simplify the expression we adopt T(i; x, y) to represent the cloud-top temperatures of nine pixels surrounding the central pixel at image location (x, y). The index i is used merely to identify the ith pixel (i = 1, 2, . . . , N) surrounding the central pixel at image location (x, y). Similarly, f(i; x, y) is the value of kernel func-tion assigned to the cloud-top temperature of the ith pixel surrounding the central pixel at (x, y). The kernel function acts as a moving window of weights corresponding to cloud-top temperatures T(x, y). The rationale of using spa-tial convolution equation is twofold (seeFig. 2): (1) sensor onboard the satellite receives not only radiances from the target cloud pixel (path of solid line in Fig. 2) but also atmospheric-scattered radiances from adjacent pixels (path of dashed line in Fig. 2), and (2) precipitation within a cloud pixel may result in rainfall over an area at the ground
level much larger than the corresponding ground cell. An-other reason to use spatial convolution is that there may be errors in the navigation of the satellite data, or displace-ments between the cloud tops and the associated precipita-tion at ground level due to wind shear. Both processes can be characterized by point-spread-functions.
Suppose that there are n pixel-average rainfall in the study area, we then establish n spatial convolution equa-tions and form the following matrix equation:
Rð1Þ Rð2Þ .. . RðnÞ 2 6 6 6 6 4 3 7 7 7 7 5¼ Tð1; 1Þ Tð1; 2Þ T ð1; N Þ Tð2; 1Þ Tð2; 2Þ T ð2; N Þ .. . .. . . . . .. . Tðn; 1Þ Tðn; 2Þ . . . Tðn; N Þ 2 6 6 6 6 4 3 7 7 7 7 5 fð1Þ fð2Þ .. . fðN Þ 2 6 6 6 6 4 3 7 7 7 7 5. ð6Þ Or, R = T . . . F. T(j, i) (j = 1, . . . ,n; i = 1, . . . ,N) in Eq.(6)
represents the cloud-top temperature of the ith surround-ing pixel associated with the jth pixel-average rainfall (R(j)). The kernel function F can be solved by the least-squared method
F ¼ ðT0TÞ1T0R ð7Þ Eqs.(4)–(6) only consider CTT of single spectral chan-nel, if CTT from k infrared channels are considered, Eq.
(6) can be extended as Rm¼ Xk i¼1 wiRiðmÞ ¼ Xk i¼1 wi XN j¼1 TijðmÞfij ! ¼X k i¼1 XN j¼1 TijðmÞðwifijÞ ¼ Xk i¼1 XN j¼1 TijðmÞfij0; m¼ 1; 2; ; n; ð8Þ
where Rmrepresents the multi-spectral spatial convolution
estimate of the mth pixel-average-rainfall. Previous re-searches suggested using CTT threshold of 253 K for dis-crimination of rain and no-rain pixels (Griffith et al., 1978; Woodley et al., 1979; Negri and Adler, 1987; Adler and Negri, 1988). Therefore, effective temperatures (CTT minus 253 K) are used in the spatial convolution equation (Eqs.(6) and (8)). All cloud pixels with CTT higher than
Fig. 2. Processes of radiances received at sensor and rainfall reaching ground surface.
Fig. 3. Comparison of block kriging estimates and MSSC predictions of 1-h-PAR using various numbers of spectral channels (model building).
253 K correspond to zero effective temperatures, and there-fore, produce no rainfall. Also, readers are reminded that in order to achieve 3-h-lead rainfall forecasting 3-h-PAR (3-h cumulative pixel-average-rainfall starting at time t) to-gether with GMS CTT of time t must be used in Eqs.(6) and (8). The kernel function (F in Eq. (7) and f0
ij in Eq.
(8)) is updated at every real-time step using the Kalman fil-tering algorithm by taking into account the rainfall predic-tion errors of the previous time step.
4. Results and discussion
4.1. 1-h-PAR accuracy assessment in model building Block kriging estimates of 1-h-PAR corresponding to 37 raingauges were compared against predicted 1-h-PAR using GMS images by the MSSC approach.Fig. 3 demon-strates that block kriging estimates and MSSC predictions of 1-h-PAR are highly correlated and inclusion of more spectral channels yield better results. Fig. 3 accounts for a total of 11,174 data points which correspond to PAR of 12 typhoon events with durations ranging from 10 to 57 h.
4.2. Cross-validation of 1- and 3-h-PAR forecast
Eqs. (6) and (8)which involve block kriging PAR esti-mates corresponding to available raingauges and the kernel function f, is solved by least-squared method. Accuracy assessment in Section4.1is therefore, analogue to compar-ison of measurements and their regression estimates. Such comparison, generally, only indicates good data-fitting and will not serve to demonstrate the accuracy of forecasts. In this study, we adopted a cross-validation approach of accu-racy assessment for 1- and 3-h- PAR forecasts.
From a total of n block kriging PAR estimates, we first remove one PAR estimate, say corresponding to the first raingauge, and solve Eqs.(6) or (8)for the kernel function f using the remaining (n 1) PAR estimates. We then cal-culate the PAR forecast corresponding to the first rainga-uge. The same procedure is repeated again by replacing the removed PAR estimate with another one until all of the n PAR estimates have been chosen.
We conducted cross-validation for 1- and 3-h-PAR fore-casts (1- and 3-h-lead forecasting) at pixels corresponding to individual raingauge sites for individual events. Figs. 4 and 5, respectively, show examples of cross validation for
1- and 3-h-PAR forecasts corresponding to four raingauges (Guandu, Chung-Cheng Bridge, Dabao and Gauyi) during Typhoon Winnie occurred in 1997. It can be seen clearly that MSSC forecasts of 1- and 3-h-PAR forecasts are com-patible with block kriging estimates.
In addition to cross-validation, the 3-h-PAR forecasts are also compared with forecasts by simple persistence (i.e. assuming that rainfall amount from the previous 3 h will fall again during the next 3 h). Even though data from many typhoons and raingauges were used in our study, we only show comparisons between 3-h-PAR fore-casts (by spatial convolution and simple persistence meth-ods) and block-kriging PAR (which are considered as observed rainfall) for Chung-Cheng Bridge and Gauyi stations during Typhoon Winnie. As can be seen in
Fig. 6, forecasts by spatial convolution tend to follow the variation trend of block-kriging PAR. Root mean
square errors of the spatial convolution method are 24.6 and 20.4 mm for Chung-Cheng Bridge and Gauyi sta-tions, respectively. As can be expected, the simple persis-tence forecasts follow the variation trend very well, however, it cannot yield forecasts at locations, where no raingauge stations exist. In contrast, the proposed spatial convolution approach solves for the kernel function and then extends forecasts to large cloud-covered area with no riangauge stations.
5. Conclusion
We have developed a multi-spectral spatial convolution approach for real-time rainfall forecasting using geosta-tionary weather satellite images. The approach incorpo-rates cloud-top temperatures of three infrared channels in a spatial convolution context. It not only characterizes
the input–output relationship of cloud-top temperature and rainfall on the ground level but also is more consistent with physical and remote sensing principles than single-pix-el matches. Through cross-validation, the MSSC demon-strates its capability of achieving high accuracy for 1- to 3-h-lead pixel-average-rainfall forecasting.
Acknowledgements
The authors are thankful to Dr. T.K. Chiu of Central Weather Bureau for offering comments and suggestions. We are also grateful for support by the Council of Agricul-ture and staff of the Soil and Water Conservation Bureau, especially, the section chief Mr. J.L. Wang, and Mr. C.S. Yen.
References
Adler, R.F., Negri, A.J. A satellite infrared technique to estimate tropical convective and stratiform rainfall. J. Appl. Meteorol. 27, 30–51, 1988. Ba, M.B., Nicholson, S.E. Analysis of convective activity and its relationship to the rainfall over the Rift Valley Lakes of east Africa
during 1983–90 using METEOSAT infrared channel. J. Appl. Mete-orol. 37, 1250–1264, 1998.
Ba, M.B., Gruber, A. GOES multispectral rainfall algorithm (GMSRA). J. Appl. Meteorol. 40, 1500–1514, 2001.
Berg, W., Bates, J.J., Jackson, D.L. Analysis of upper-tropospheric water vapor brightness temperatures from SSM/T2, HIRS, and GMS-5 VISSR. J. Appl. Meteorol. 38, 580–595, 1999.
Chile`s, J.P., Delfiner, P. Geostatistics – Modeling Spatial Uncertainty. John Wiley and Sons Inc., New York, 1999, pp. 449–451.
Conway, E.D. An Introduction to Satellite Image Interpretation. The Johns Hopkins University Press, Baltimore, MD, 1997.
DÕsouza, G., Barrett, E.C., Power, C.H. Satellite rainfall estimation techniques using visible and infrared imagery. Remote Sens. Revs 4 (2), 379–414, 1990.
Ferreira, E., Amayenc, P., Oury, S., et al. Study and tests of improved rain from the TRMM precipitation radar. J. Appl. Meteorol. 40, 1878– 1899, 2001.
Griffith, C.G., Woodley, W.L., Grube, P.G., et al. Rain estimation from geosynchronous satellite imagery – visible and near Infrared studies. Mon. Weather Revs 106, 1153–1171, 1978.
Journel, A.G., Huijbregts, C.J. Mining Geostatistics. Academic Press, London, 1978.
Kazumasa, A., Liu, G. Passive microwave precipitation retrievals using TMI during the Baiu period of 1998. Part I: algorithm description and validation. J. Appl. Meteorol. 39, 2024–2037, 2000.
Kummerow, C., Giglio, L. A method for combing passive microwave and infrared rainfall observations. J. Atmos. Ocean. Technol. 12, 33–45, 1995.
Martin, D.W., Goodman, B., Schmit, T.J., et al. Estimates of daily rainfall over the Amazon basin. J. Geophys. Res. 95 (D10), 17043– 17050, 1990.
Miller, S.W., Arkin, P.A., Joyce, R.A. Combined microwave/infrared rain rate algorithm. Int. J. Remote Sens. 22 (17), 3285–3307, 2001. Mishra, J.K., Sharma, O.P. Cloud top temperature based precipitation
intensity estimation using INSAT-1D data. Int. J. Remote Sens. 22 (6), 969–985, 2001.
Negri, A.J., Adler, R.F. Infrared and visible satellite rain estimation Part I: a rrid cell approach. J. Clim. Appl. Meteorol. 26, 1553–1564, 1987.
Prabhakara, C., Iacovazzi Jr., R., Weinman, J.A., et al. A TRMM microwave radiometer rain rate estimation method with convective and stratiform discrimination. J. Meteorol. Soc. Jpn 78 (3), 241–258, 2000.
Scofield, R.A. The NESDIS operational convective precipitation estima-tion technique. Mon. Weather Rev. 115, 1773–1792, 1987.
Sorooshian, S., Hsu, K.L., Gao, X., et al. Evaluation of PERSIANN system satellite-based estimates of tropical rainfall. Bull. Am. Mete-orol. Soc. 81 (9), 2035–2046, 2000.
Todd, M.C., Kidd, C., Kniveton, D., et al. A combined infrared and passive microwave technique for estimation of small-scale rainfall. J. Atmos. Ocean. Technol. 18 (5), 742–755, 2001.
Toshiaki, K., Kawanishi, T., Kuroiwa, H., et al. Development of precipitation radar onboard the tropical rainfall measuring mission (TRMM) satellite. IEEE Trans. Geosci. Remote Sens. 39 (1), 102–116, 2001.
Turk, F.J., Marzano, F.S., Smith, E.A., et al. Using coincident SSM/I and infrared geostationary satellite data for rapid updates of rainfall, in: Geoscience and Remote Sensing Symposium Proceedings, IGARSSÕ98, pp. 150–152, 1998.
Woodley, W.L., Griffith, C.G., Griffin, J.S., et al. The inference of GATE convective rainfall from SMS-1 imagery. J. Appl. Meteorol. 19, 388– 408, 1979.
Xu, L., Gao, X., Sorooshian, S., et al. A microwave infrared threshold technique to improve the GOES precipitation index. J. Appl. Mete-orol. 38, 569–579, 1999.
Fig. 6. Comparison of 3-h-PAR forecasts by spatial convolution, simple persistence and block-kriging PAR during Typhoon Winnie at Chung-Cheng Bridge and Gauyi raingauge stations.
