IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING 1
Fusion of MODIS Images UsingKriging With External Drift
Marcio H. Ribeiro Sales, Carlos M. Souza, Jr., and Phaedon C. Kyriakidis
Abstract—The Moderate Resolution Imaging Spectroradiome-ter (MODIS) has been used in several remote sensing studies,including land, ocean, and atmospheric applications. The ad-vantages of this sensor are its high spectral resolution, with 36spectral bands; its high revisiting frequency; and its public domainavailability. The first seven bands of MODIS are in the visible,near-infrared, and mid-infrared spectral regions of the electro-magnetic spectrum which are sensitive to spectral changes dueto deforestation, burned areas, and vegetation regrowth, amongother land-use changes, making near-real-time forest monitoringa suitable application. However, the different spatial resolutionof the spectral bands placed in these spectral regions imposeschallenges to combine them in forest monitoring applications.In this paper, we present an algorithm based on geostatistics todownscale five 500-m MODIS pixel bands to match two 250-mpixel bands. We also discuss the advantages and limitations ofthis method in relation to existing downscaling algorithms. Ourproposed method merges the data to the best spatial resolution andbetter retains the spectral information of the original data.
Index Terms—Forest monitoring, image fusion, kriging with anexternal drift, Moderate Resolution Imaging Spectroradiometer(MODIS), spatial resolution.
I. INTRODUCTION
THE Moderate Resolution Imaging Spectroradiometer(MODIS) measures reflected and emitted energy at 36
bands placed at a wavelength range of 0.405–14.385 μm, withspatial resolutions of 250–1000 m [1], acquired at a near dailybasis. The spectral, spatial, and temporal resolutions of MODISrender it a powerful tool for global and large-scale applicationsof remote sensing. The first seven bands of MODIS cover thewavelength range most related to features of land, cloud, andaerosol properties [2], for example, chlorophyll absorption at
Manuscript received December 21, 2011; revised May 18, 2012; acceptedMay 26, 2012. This work was supported in part by The David and LucilePackard Foundation and in part by the Climate Land Use Alliance. The workof M. H. Ribeiro Sales and P. C. Kyriakidis was supported in part by the U.S.National Geospatial Intelligence Agency.
M. H. Ribeiro Sales is with the Instituto do Homem e Meio Ambienteda Amazônia—Imazon, Belém 66613-397, Brazil, and also with the Pro-duction Ecology and Resource Conservation Graduate School, WageningenUniversity, 6709 PB Wageningen, The Netherlands (e-mail: [email protected]).
C. M. Souza, Jr., is with the Instituto do Homem e Meio Ambienteda Amazônia—Imazon, Belém 66613-397, Brazil (e-mail: [email protected]).
P. C. Kyriakidis is with the Department of Geography, University of Califor-nia at Santa Barbara, Santa Barbara, CA 93106-4060 USA, and also with theDepartment of Geography, University of the Aegean, 81100 Mytilene, Greece(e-mail: [email protected]).
Color versions of one or more of the figures in this paper are available onlineat http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TGRS.2012.2208467
red band (i.e., band 1, 0.620–0.670 μm) and canopy backscat-tering at near infrared (band 2, 0.841–0.876 μm). These bandscan also be combined to retrieve biophysical properties offorests using products such as vegetation indices available at theLand Process Distributed Active Archive Center (LP DAAC)Web site [3]. These features make the first seven MODIS bandsthe most important ones for land-use and land-cover change(LULCC) applications [1]. In Brazil, MODIS data have beenused to detect deforestation at near real time for forest lawenforcement and environmental policy verification [4], [5].
However, a limitation is imposed by the difference in thespatial resolutions of these bands: Bands 1 and 2 are theonly ones available at the 250-m spatial resolution [6] whilebands 3–7 are in the 500-m spatial resolution. This differencein spatial resolutions limits the generation of products basedon MODIS data at the 250-m resolution using the full rangeof spectral information important for forest monitoring. The250-m resolution bands 1 and 2 can, however, be spatiallydegraded to 500 m to create a product at this resolution withseven spectral bands, but the 250-m resolution is preferredbecause many changes in the land cover occur at this spatialscale [7]. Splitting a 500-m pixel into four 250-m pixels wouldresult in a product with seven spectral bands with no necessarilydirect spectral correspondence from the two spatial resolutions.Spectral mixing, for example, is a scale-dependent [8] phe-nomenon that could drastically affect this type of data-fusedproduct. Alternatively, bands 1 and 2 alone have been usedfor land-cover mapping applications [9], [10]. However, theaddition of information on the infrared bands can significantlyimprove the detection and analysis of LULCC [9].
Combining data from different spatial resolutions suchas 250- and 500-m MODIS bands is a typical image fusionproblem. The general objective of image fusion is to combinetwo (or more) images in order to obtain a better image, retainingthe desired characteristic of each of the original images [11].Image downscaling is a special case of image fusion, in whichfine spatial coarser spectral resolution images are combinedwith coarser spatial higher spectral resolution images to yield aspectrally enhanced product at the finer spatial resolution [12].This type of approach is more appropriate than pixel splittingthe 500-m MODIS images because it retains the spatial infor-mation at the finer resolution and the spectral correspondenceof the two spatial-resolution bands is guaranteed. Imagedownscaling has the potential to greatly improve the utility ofMODIS data for several remote sensing applications. For ex-ample, improvements in snow cover mapping have been foundafter the application of a downscaling algorithm to MODISimages [13].
2 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING
Downscaling methods can be classified into two categories[14]: color transformation methods and statistical/mathematicalmethods. The most frequently used color transformationmethod is the intensity–hue–saturation (IHS) method, in whichthree selected bands of a coarse-spatial-resolution multispectralimagery are chosen and converted to red–green–blue valuesand then transformed into their intensity, hue, and saturationcomponents. Then, one of these components is replaced bythe fine-resolution image, and the transformation is reversedto obtain the final downscaled image. The IHS method isrelatively fast and easy to use, and for that, it is a standardfusion method and is implemented in many software packages.However, this method can only deal with a limited number ofbands as it requires the choice of the best three bands for theparticular application and sensor in question to be processed[15]. There are also mathematical methods for selecting bandswith the largest variance [16], but in either case, as not allbands are explicitly used, information is lost, producing spectraldegradation in the downscaled images [17].
Statistical or mathematical methods like principal compo-nent analysis (PCA) and wavelet transforms are preferred bysome because the loss of information and spectral degradationincurred are less than in the IHS-based method. The aforemen-tioned methods are similar to the IHS method in many aspects.For example, in the PCA method, the first principal componentof the coarse spatial multispectral resolution imagery is substi-tuted by the fine-spatial-resolution image, and then, the PCAprocedure is reversed to obtain the downscaled images [18].In the wavelet transform method, wavelet planes are obtainedfrom the multispectral imagery and are replaced by waveletplanes obtained from the fine-resolution image. Subsequently,the inverse of the wavelet transform is applied to the waveletplanes to obtain the downscaled image [17]. The objective is totransfer the high-frequency component from the fine-resolutionimage to the coarse-resolution images. All of these methodsgenerate certain degree of spectral degradation, and because ofthat, methods based on combinations of the above have beendeveloped to overcome their individual problems [19], [20].
Geostatistics is a branch of statistics that explicitly incor-porates the concept of spatial correlation in the analysis ofspatial data [21], [22]. The geostatistical family of estimationmethods called kriging has been mainly used in interpolationapplications [23]–[25] and recently appeared as a potential toolfor image downscaling [12], [26], [27]. Kriging’s main advan-tages over other downscaling techniques are the account of thefollowing: 1) differences between pixel sizes of the two spatialresolutions, by including the sensor’s point spread function(PSF); 2) the spatial correlation structure of the attribute valuesat image pixels [e.g., digital number (DN)]; 3) the preservationof the spectral characteristics of the imagery; and 4) the utiliza-tion of all spectral bands in the downscaling procedure.
Recently, a cokriging approach has been successfully de-veloped and applied for downscaling Landsat-ETM multispec-tral bands [27]. However, the cokriging method requires jointmodeling of autovariograms and cross-variograms among allimage bands, which makes it difficult to automate. Here, wepresent an alternative solution to this problem based on krigingwith external drift (KED). KED is easy to implement and
render operational, as it requires only the estimation of directvariograms (or autovariograms). We then apply our proposedmethodology for increasing the spatial resolution of the 500-mresolution (bands 3–7) of MODIS images to 250-m resolution,using imagery from the Brazilian Amazon region. Next, weevaluate the results visually and quantitatively by investigatingthe spectral distortion of the downscaled images via commonlyused metrics, such as the spectral angle mapper (SAM), theuniversal image quality index (UIQI), the relative dimension-less global error in synthesis (ERGAS), and the spectral infor-mation divergence (SID). We computed the same metrics fromother three existing pansharpening techniques for comparisonin several parts of the image. Finally, we discuss the strengthsand limitations of the proposed KED downscaling algorithm.
II. METHODS
A. Geostatistical Downscaling
In KED downscaling, the fine- and coarse-resolution im-ages are treated as joint realizations of two autocorrelated andcross-correlated random fields, which are referred to as Zv
and ZV , respectively. We will use superscripts k to indicate“source” bands 1 and 2 and l to indicate “target” bands 3to 7. Therefore, Zk
v , k = 1, 2, and ZlV , l = 3, . . . , 7, will be
referred to as the source and target random fields, respectively.Lower case z represents the (pixel) values or realizations ofthese random fields, and s represents the 2-D coordinates ofpixel centers. Thus, zlV (s) and zkv (s
′) represent the outcomesof the random fields Zl
V and Zkv , at pixels centered at s and
s′. Finally, both zlV (s) =∫V zl•(s)PSFV (s)ds and zkv (s
′) =∫v z
k• (s
′)PSFv(s′)ds are assumed to be weighted averages of
their underlying attribute surfaces zl•(s) and zk• (s′) over the
pixel support centered at s and s′, weighted by the PSF of therespective sensors for a given spatial resolution. The term DNwill be used for z for simplicity hereafter. We will estimatetarget DN values of bands 3–7 at the 250-m spatial resolutionusing the available DN values of bands 3–7 at 500-m resolutionand DN values of source bands 1 and 2 at 250-m resolution.
KED downscaling formulates the estimated DN value zlv(s)as a linear combination of N(s) coarse-resolution DN values ofpixels of the same band l
zlv(s) =
N(s)∑n=1
wln(s)z
lV (sn). (1)
In (1), n = 1, . . . , N(s) are local indices of target coarse-resolution pixels V (sn) surrounding a target fine-resolutionpixel v(s). The wl
n(s), n = 1, . . . , N(s), denotes the weightsapplied to coarse-resolution DN values zlV (sn), located in thevicinity of s for estimation of its fine-resolution DN value. Inaddition, the expected value of the random variable Zl
v(s) ismodeled as a function of collocated fine-resolution DN valuesof the K source bands as follows:
E{Zlv(s)
}=
K∑k=0
blk(s)zkv (s), z0v(s) = 1 ∀s. (2)
In (2), the K bands are also known in geostatistical liter-ature as the “drift variables.” Equation (2) is called the KED
RIBEIRO SALES et al.: FUSION OF MODIS IMAGES USING KRIGING WITH EXTERNAL DRIFT 3
regression of DN values of band l on DN values of K sourcebands, and the blk(s) denotes the KED regression coefficients,where s implies that these coefficients are locally estimated foreach fine-resolution-pixel location s.
Based on the aforementioned formulation, KED downscal-ing aims at finding the weights wl
n(s) that yield an unbiasedestimate of zlv(s) and minimize the variance of the estimationerror V ar[Zl
v(s)− Zlv(s)]. Unbiasedness can be achieved by
imposing the following constraints on the weights:
N(s)∑n=1
wln(s)z
kV (sn)=zkv (s), l=3, . . . , 7; k=1, 2. (3)
Equation (3) is verified in Appendix. The KED weightswl
n(s) that minimize the variance of the estimation error, sub-ject to the constraints in (3), are given by system (4), shown atthe bottom of the page, where
ClR(V,V )(sn, sn′) =Cov
(ZlV (sn)−
K∑k=0
blk(sn)ZkV (sn),
ZlV (sn)−
K∑k=0
blk(sn′)ZkV (sn′)
)
=Cov (RV (sn), RV (sn′)) ,n, n′ = 1, . . . , N(s)
are the residual covariances between pixels of the coarse-resolution image at band l surrounding the target finepixel v(s), while Cl
R(V,v)(sn, s) = Cov(RV (sn), Rv(s)), n =
1, . . . , N(s), denotes the residual covariances between thecoarse pixels V (sn) and the fine-resolution pixel v(s). TermsμlK(s), k = 0, . . . ,K, are Lagrange parameters included in the
system to account for the constraints on the weights in (3). Theminimum error variance for each fine-resolution pixel v(s) isobtained as
σle(s)
2 = σlR
2−
⎡⎢⎢⎢⎢⎢⎢⎢⎣
ClR(V,v)(s1, s)
...Cl
R(V,v)
(sN(s), s
)1...
zKv (s)
⎤⎥⎥⎥⎥⎥⎥⎥⎦
T
×
⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣
wl1(s)...
wlN(s)(s)
−μl0(s)...
−μlk(s)
⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦
(5)
where σlR2
is the prior fine-resolution residual variance for bandl. The minimum estimation error variance of (5) provides ameasure of uncertainty resulting from using only the nearby
coarse-pixel DN values in estimating the unknown DN valueat fine pixel s. With no auxiliary bands, system (4) reduces to⎡⎢⎣
wl1(s)...
wlN(s)(s)
⎤⎥⎦
=
⎡⎢⎣
Cl(V,V )(s1, s1) . . . Cl
(V,V )(s1, sN(s))
.... . .
...Cl
(V,V )
(sN(s), s1
). . . Cl
(V,V )
(sN(s), sN(s)
)⎤⎥⎦−1
.
⎡⎢⎣
Cl(V,v)(s1, s)
...Cl
(V,v)
(sN(s), s
)⎤⎥⎦ (6)
and the procedure is called downscaling by simple kriging andis described in detail in [28].
Equation (4) requires knowledge of the residual covariancebetween any two coarse pixels of each target band l. It alsorequires the residual covariance between any coarse and finepixels in the same target band l. These covariances are obtainedby convolution of a punctual residual covariance model Cl
R(•)with the PSF of the sensor over the support(s) involved. In thispaper, the punctual covariance was estimated by deconvolutionof the detrended DN values, as described in the followingapplication section.
B. KED Downscaling Applied to MODIS Imagery
We used the MODIS products MOD09GQK at the 250-mspatial resolution (bands 1 and 2) and MOD09GAV5 at the500-m resolution (bands 3–7). This imagery covers tropicalforests in the Brazilian Amazon and was intentionally chosenfor the analysis due to the high degree of anthropogenic activityin the area, as well as to illustrate the application of themethod for monitoring of deforestation (Fig. 1). The imagehas been preprocessed to correct for atmospheric effects beforeapplication of our algorithm.
Table I shows Pearson’s correlation coefficients betweenthe upscaled fine-spatial-resolution bands (1 and 2) and thecoarse-spatial-resolution bands (3–7). High correlations be-tween bands 1 and 2 and the other bands exist, which indicatesthat these bands carry a significant amount of information onbands 3–7. Under the assumption that all images share the samePSF (see Appendix), we show that the information on fine-resolution bands 1 and 2 can be combined with the information
4 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING
Fig. 1. Study area.
TABLE IPEARSON’S CORRELATION COEFFICIENTS OF BANDS 3–7
WITH UPSCALED BANDS 1 AND 2
on coarse-resolution bands 3–7 in KED to obtain unbiasedestimates of fine-resolution bands 3–7.
The residual covariances ClR(V,V )(si, sj) of the KED system
(4) were obtained by the following procedure:
1) Regression of each target band on bands 1 and 2 of thefine-resolution image, after the latter was upscaled to500-m resolution for resolution compatibility.
2) From this regression model, we obtained five coarse-resolution residual images rlV , l = 1, . . . , 5.
3) Estimation of the residual autocovariance model for eachcoarse-resolution residual image rlV based on the em-pirical variogram calculated from rlV values. Variogramanalysis methodology is described in detail in [18].
4) Punctual residual autocovariance models ClR(•)(h) were
obtained by deconvolution, using the iterative proceduredescribed in [27] (Fig. 2).
5) Using the punctual covariance models ClR(•)(h), the fi-
nal coarse-to-coarse and coarse-to-fine variogram modelswere obtained by convolving the punctual covariancemodel with the PSF adopted for MODIS
ClR(V,V ) =Cl
R(•)(s ∈ V, s′ ∈ V ′) ∗ PSFV (x, y ∈ V )
∗ PSFV (x′, y′ ∈ V ′) (7)
ClR(V,v)(V, v) =Cl
R(•)(s ∈ V, s′ ∈ v) ∗ PSFv(x, y ∈ v)
∗ PSFv(x′, y′ ∈ v′) (8)
where ∗ denotes the convolution operator and PSF =exp((x2 + y2)/2σ2) is the PSF adopted for MODIS im-agery in this work. The PSF of a sensor dictates how
the radiation over an area (specified by the locations xand y) is integrated for generation of a pixel’s DN by thesensor [28]. The variables x and y are coordinate offsetsin relation to the center of a pixel, and σ here is the spreadparameter, which is usually set as the half of the pixelsize.
6) Once the residual autocovariance (coarse resolution) andcross-covariance (coarse-to-fine resolution) models havebeen defined, the kriging weights are calculated using (4)for each fine-resolution pixel v(s). After obtaining theweights, the corresponding KED estimate is computedusing (1).
C. Validation
As only bands 1 and 2 are available at the two spatial res-olutions, we used these bands as reference data for validation.We performed band 1 KED validation by applying KED down-scaling to band 1 using only band 2 as a drift and likewise forband 2 using band 1 as a drift. We compared the estimated andreal images via scatterplots, correlation coefficients, and rootmean square errors. We also evaluated the spectral distortions ofthe resulting images by applying commonly used image qualityevaluation metrics (SAM, ERGAS, SID, and UIQI).
III. RESULTS AND DISCUSSION
The residual variogram model results for γlR(V,V )(h) and
γlR(V,v)(h) are shown in Fig. 3, and the final set of deconvolved
punctual model parameters is shown in Table II. The punctualresidual variogram model has a higher sill than the convolvedvariogram models γl
R(V,V )(h) and γlR(V,v)(h). This is because
convolved variograms are computed via regularization (averag-ing) of point variogram values over an area. Punctual variogrammodels did not vary significantly among bands, except for thesill parameter, which reflects the difference in variance bandDN values in different bands. The range was around 14 km,and the functional form was closely approximated by an ex-ponential model for all bands, except for band 5, which hada range of 8 km and was better approximated by a sphericalmodel. In general, the nugget effect contributions to the residualvariograms were small (again, except for band 5), indicatingsmall subpixel variability (under 500 m) of DN values in thisregion. Band 5 exhibited stripe artifacts, caused by noise in thesensors, and this could explain the difference in behavior fromall other bands. This suggests that variogram results are verysensitive to image artifacts, and such artifacts must be correctedbefore the algorithm is applied.
For estimation, we decided to use only pixels up to thedistance at which 40% of the total sill was reached in eachband, to reduce computation time. The original 500-m MODISbands and their downscaled versions are shown in Fig. 4.Significant visual improvements can be observed in every band.For example, boundaries of deforested areas are better definedin the downscaled images than in the original 500-m resolutionimages, which could simplify interpretation of new smallerdifficult-to-detect areas of deforestation.
RIBEIRO SALES et al.: FUSION OF MODIS IMAGES USING KRIGING WITH EXTERNAL DRIFT 5
Fig. 2. KED downscaling methodology flowchart.
Validation analysis results are shown in Figs. 5 and 6 andTable III. Fig. 5 shows a false-color composite of bands 7,6, and 4 for the entire image, while Fig. 6 shows the samecomposite for subsets randomly selected at different regionsof the scene for more detailed evaluation. Results of otherthree pansharpening methods (PCA, wavelet pansharpening,and Brovey) are also shown for comparison. The visual inspec-tion of the images shows better results for KED downscalingin conserving the original spectral information, as well as indelineating boundaries. The wavelet-derived image seems tohave conserved the spectral information reasonably well butdid not show comparably good boundary delineation. PCA
and Brovey-derived images showed good boundary delineationbut exhibit clear spectral discrepancies when compared to theoriginal image, with the PCA image showing less contrast andthe Brovey image showing excessive contrast.
Table III shows the results of the quantitative comparisonof the downscaled images, via standard metrics such as SAM,ERGAS, SID, and UIQI. Lower scores for SAM and SID, aswell as higher scores for ERGAS and UIQI, indicate betterspectral correspondence between reference and downscaledimages. The quantitative results in general confirm the vi-sual inspection conclusions above: Overall, the wavelet andKED downscaling methods had the best scores, confirming
6 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING
Fig. 3. Empirical 500-m resolution and deconvolved residual point variogrammodels. (Triangles) A 500-m-resolution residual variogram model was com-puted by convolution for comparison.
TABLE IIRESIDUAL POINT VARIOGRAM MODEL PARAMETERS BY MODIS BAND
the efficiency of the two methods in conserving spectral infor-mation. The Brovey transformation method aims at preservingthe relative contributions of the different spectral bands, whichexplains the corresponding extremely low value for SID.
The downscaling procedure greatly increases the value ofMODIS data for remote sensing applications. For example,using the original images, land-cover classification algorithmsfor 250-m spatial resolution are restricted to two bands, thuspotentially limiting the accuracy of classification algorithms.Classification can be performed using all bands for the 500-mresolution, but the resulting images have limited ability toresolve small targets (i.e., < 10 ha). With downscaled images,classification algorithms at 250-m spatial resolution can nowbe developed using all information provided by MODIS data,having the potential to increase classification accuracy andprovide a greater level of detail entailing a more resolvedpicture of targets on the ground [10].
These features are expected to be important for land-covermonitoring applications such as deforestation detection in
Fig. 4. Sample of original and downscaled images.
Amazonian forests. The finer spatial resolution can improveexisting methods of real-time deforestation detection such asthe Real-Time Detection of Deforestation in the Legal Amazon[4] (Portuguese: “Detecção do Desflorestamento da AmazôniaLegal em Tempo Real”) and the Deforestation Alert System(Portuguese: “Sistema de Alerta de Desmatamento”), which areboth based on MODIS data. Using more spectral informationhas the potential to improve the true detection of deforesta-tion and forest degradation at subpixel levels through spectral
RIBEIRO SALES et al.: FUSION OF MODIS IMAGES USING KRIGING WITH EXTERNAL DRIFT 7
Fig. 5. Visual comparison of downscaling results. Red squares are the regions zoomed in, shown in Fig. 6. Axis values represent the Universal TransverseMercator (UTM) coordinates. (a) Original image. (b) KED. (c) PCA. (d) Wavelet. (e) Brovey.
mixture analysis as performed using Landsat imagery [5], whilethe fine spatial resolution can improve the estimates of the forestarea affected by deforestation and degradation. For these sys-tems, the increase in spatial resolution is particularly importantbecause deforestation caused by anthropogenic activity mostlyoccurs within the 250-m resolution [7].
One important note is that, as the method is based ona neighboring averaging, the results are subject to edge oradjacency effects. This can be easily taken into account byworking with overlapping tiles, where the number of overlap-ping pixels needs to be at least half the width of the searchwindow.
8 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING
Fig. 6. Comparison of downscaling results at three randomly selected regions of the study image. Axis values represent the UTM coordinates. (a) Original image.(b) KED. (c) PCA. (d) Wavelet. (e) Brovey.
A major advantage of the proposed method over cokriging,which is another geostatistical-based approach recently appliedtoward the same downscaling objective, is the fact that KEDdownscaling requires parameter estimation for a smaller num-ber of variogram models, thus reducing inference effort andpossibly the amount of human intervention. Using cokriging to
downscale MODIS images would require parameter estimationfor ten additional deregularized cross-semivariograms (five forthe cross-variogram between bands 1 and 3–7 plus five for thecross-variograms between bands 2 and 3–7), which must bejointly modeled via the linear model of coregionalization; thiswould greatly increase manual intervention due to variogram
RIBEIRO SALES et al.: FUSION OF MODIS IMAGES USING KRIGING WITH EXTERNAL DRIFT 9
TABLE IIIQUANTITATIVE COMPARISON OF DOWNSCALING RESULTS
deconvolution (see [29]). The KED procedure, however, worksonly under the assumption that the averaging scheme (i.e., thePSF) used for upscaling the punctual radiation into a pixel valueis invariant for all bands involved in the downscaling process.Otherwise, the method would generate biased coefficients and,consequently, biased estimates [28]. This limitation does notapply to cokriging, which should be the preferred method incases where the objective is to fuse images which bands havedifferent PSFs. In simpler words, the KED method describedhere is appropriate and much easier to implement than cokrig-ing when the images used in the downscaling procedure comefrom sensors with identical PSFs for each band.
The major disadvantage of the proposed method in relationto cokriging is the extended computational cost. While otherpansharpening algorithms can be applied to a MODIS tilein just a few seconds, KED downscaling took c.a. 20 minin a Core i-7 Q840 notebook. This happens because (4) isevaluated locally for prediction at each fine pixel, increasingprocessing time with the size of the neighborhood. In cokriging,the interpolation weights are only a function of the neigh-borhood configuration and the variogram models involved.As that neighborhood configuration typically remains constantin remote sensing applications, cokriging weights need to beevaluated only once and applied verbatim to any predictionpixel. In spatially adaptive downscaling via cokriging (cite thefollowing), the weights change across the image; hence, thecomputational advantage of global cokriging is lost to achievea more realistic (local) variogram model fitting. KED down-scaling in this spatially adaptive case, however, is much moreefficient, since only residual autovariograms are individuallyfitted, instead of jointly permissible autovariograms and cross-variograms for all pairs of bands. Therefore, KED can be apreferred approach, depending on the scale of the applicationand the computational power available, provided that the PSFsof the sensors involved are the same.
IV. CONCLUSION
Our KED approach has improved the spatial informationprovided by the 500-m spatial-resolution bands, generating asharper five-band 250-m-spatial-resolution image. This gain ofdetail could improve precision and detection of classificationmethods and change detection algorithms for near-real-timedeforestation monitoring of the Brazilian Amazon. Also, KEDis easier to automate, since less variogram modeling effort isrequired with regard to cokriging. We recommend the use of thealgorithm when fusing images which bands share the same PSF.Future work should compare the impact of using downscaled
images on the accuracy of land-cover classification and otherremote sensing applications.
APPENDIX
DERIVATION OF UNBIASEDNESS CONSTRAINTS
In this section, we demonstrate how the unbiasedness con-straints on the weights lead to (3). Unbiasedness is ensured byforcing the expectation of the kriging prediction (1) to equal theexpectation of the random variable Zl
v , i.e.,
E
⎧⎨⎩
N(s)∑n=1
wn(s)ZlV (s)
⎫⎬⎭ = E
{Zlv(s)
}.
Therefore, unbiasedness leads to the following conditions:
E
⎧⎨⎩
N(s)∑n=1
wln(s)Z
lV (s)
⎫⎬⎭
= E{Zlv(s)
}2 ⇒
N(s)∑n=1
wln(s)E
{ZlV (s)
}
=
K∑k=0
blkzkv (s) ⇒
N(s)∑n=1
wln(s)
K∑k=0
b∗lk zkV (s)
=
K∑k=0
blkzkv (s) ⇒
K∑k=0
b∗lk
N(s)∑n=1
wln(s)z
kV (s)
=
K∑k=0
blkzkv (s), l = 1, . . . , L, z0V (s), z
0v(s) = 1 ∀s
which lead to (3) if the KED coefficients b∗lk and blk obtained atthe coarse and fine resolutions are the same or, in other words,if the model for (2) is also valid at the coarse resolution. Thiscondition holds if the coefficients blk defined at the punctuallevel do not change between different resolutions. Since a DNvalue for a given band is obtained by the convolution of thepunctual intensity with the sensor’s PSF (here assumed thesame for all bands), we have
zlV (s) =
∫V
PSFV (s)zl•(s)ds
=
∫V
PSFV (s)K∑
k=0
bkzk• (s)ds
=
K∑k=0
bk
∫V
PSFV (s)zk• (s)ds =
K∑k=0
bkzkV (s)
which proves that the regression coefficients of (2) are inde-pendent of support size when the PSF is constant (or at leastsimilar) along all bands.
marcio.sales
Note
Please remove the 2
10 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING
REFERENCES
[1] C. O. Justice, J. R. G. Townshend, E. F. Vermote, E. Masuoka,R. E. Wolfe, N. Saleous, D. P. Roy, and J. T. Morisette, “An overview ofMODIS Land data processing and product status,” Remote Sens. Environ.,vol. 83, no. 1/2, pp. 3–15, Nov. 2002.
[2] X. Zhan, R. A. Sohlberg, J. R. G. Townshend, C. Dimiceli, M. L. Carroll,J. C. Eastman, M. C. Hansen, and R. S. DeFries, “Detection of landcover changes using MODIS 250 m data,” Remote Sens. Environ., vol. 83,no. 1/2, pp. 336–350, Nov. 2002.
[3] U. S. Geological Survey (USGS), The Land Processes Distributed ActiveArchive Center. [Online]. Available: http://lpdaac.usgs.gov/
[4] D. M. Valeriano, Y. S. Shimabukuro, V. Duarte, L. Anderson, F. Espírito-Santo, E. Arai, L. Maurano, R. C. Souza, R. M. Freitas, and L. Aulicino,“Detecção do Desflorestamento da Amazônia Legal em Tempo Real—Projeto DETER,” in Proc. 12th Simpósio Brasileiro de SensoriamentoRemoto., Goiânia, Brazil, 2005, pp. 3403–3409.
[5] C. M. Souza, Jr., S. Hayashi, and A. Veríssimo, “Near Real-Time Defor-estation Detection for Enforcement of Forest Reserves in Mato Grosso,”Instituto do Homem e Meio Ambiente da Amazônia - Imazon, Belém,Brazil, 2009.
[6] J. R. G. Townshend and C. O. Justice, “Towards operational monitoring ofterrestrial systems by moderate-resolution remote sensing,” Remote Sens.Environ., vol. 83, no. 1/2, pp. 351–359, Nov. 2002.
[7] J. R. G. Townshend and C. O. Justice, “Selecting the spatial resolution ofsatellite sensors required for global monitoring of land transformations,”Int. J. Remote Sens., vol. 9, no. 2, pp. 187–236, 1988.
[8] J. B. Adams, D. E. Sabol, V. Kapos, R. A. Filho, D. A. Roberts,M. O. Smith, and A. R. Gillespie, “Classification of multispectral imagesbased on fractions of endmembers: Application to land-cover change inthe Brazilian Amazon,” Remote Sens. Environ., vol. 52, no. 2, pp. 137–154, May 1995.
[9] K. J. Wessels, R. D. Fries, J. Dempewolf, L. O. Anderson, A. J. Hansen,S. L. Powell, and E. F. Moran, “Mapping regional land cover with MODISdata for biological conservation: Examples from the Greater YellowstoneEcosystem, USA and Para State, Brazil,” Remote Sens. Environ., vol. 92,no. 1, pp. 67–83, Jul. 2004.
[10] A. Costa and C. M. S. , Jr, “Comparação entre imagens LANDSAT ETMe MODIS/TERRA para detecção de incrementos de desmatamento naregião do Baixo Acre,” Revista brasileira de cartografia, vol. 57, pp. 93–102, Aug. 2005.
[11] R. H. Rogers and L. Wood, “The history and status of merging multiplesensor data: An overview,” in Proc. ACSM- ASPRS Annu. Conv., ImageProcess. Remote Sens., Denver, CO, 1990, pp. 352–360.
[12] N. Memarsadeghi, J. Le, M. David, M. Mount, and J. Morisette, “A newapproach to image fusion based on cokriging,” in Proc. 8th Int. Conf. Inf.Fusion, 2005, pp. 622–629.
[13] P. Sirguey, R. Mathieu, Y. Arnaud, M. M. Khan, and J. Chanussot, “Im-proving MODIS spatial resolution for snow mapping using wavelet fusionand ARSIS concept,” IEEE Geosci. Remote Sens. Lett., vol. 5, no. 1,pp. 78–82, Jan. 2008.
[14] C. Pohl and J. L. V. Genderen, “Multisensor image fusion in remotesensing: Concepts, methods and application,” Int. J. Remote Sens., vol. 19,no. 5, pp. 823–854, 1998.
[15] C. Sheffield, “Selecting band combinations from multispectral data,” Pho-togram. Eng. Remote Sens., vol. 51, no. 6, pp. 681–687, Jun. 1985.
[16] P. S. Chavez, S. C. Guptill, and J. H. Bowell, “Image processing tech-niques for TM data,” in Proc. 50th Annu. Meeting ASPRS, Washington,DC, 1984.
[17] J. Nunez, X. Otazu, O. Fors, A. Prades, V. Pala, and R. Arbiol,“Multiresolution-based image fusion with additive wavelet decomposi-tion,” IEEE Trans. Geosci. Remote Sens., vol. 37, no. 3, pp. 1204–1211,May 1999.
[18] E. Gringarten and C. Deutsch, “Teacher’s Aide variogram interpretationand modeling,” Math. Geol., vol. 33, no. 4, pp. 507–534, May 2001.
[19] B. Li, J. Wei, S. G. Ungar, S. Mao, and Y. Yasuoka, “Remote sensingimage fusion based on IHS transform, wavelet transform, and HPF,”in Proc. SPIE Image Process. Pattern Recog. Remote Sens., Hangzhou,China, 2003, pp. 25–30.
[20] M. González-Audícana, J. L. Saleta, R. G. Catalán, and R. Garcia, “Fusionof multispectral and panchromatic images using improved IHS and PCAmergers based on wavelet decomposition,” IEEE Trans. Geosci. RemoteSens., vol. 42, no. 6, pp. 1291–1299, Jun. 2004.
[21] E. Isaaks and R. Srivastava, Applied Geostatistics. London, U.K.:Oxford Univ. Press, 1989.
[22] P. Goovaerts, Geostatistics for Natural Resources Evaluation. London,U.K.: Oxford Univ. Press, 1997.
[23] F. Biondi, D. E. Myers, and C. C. Avery, “Geostatistically modeling stemsize and increment in an old-growth forest,” Can. J. Forest Res., vol. 24,no. 7, pp. 1354–1368, Jul. 1994.
[24] N. Nanos and G. Montero, “Spatial prediction of diameter distributionmodels,” Forest Ecol. Manag., vol. 161, no. 1–3, pp. 147–158, May 2002.
[25] M. Sales, C. Souza, P. Kyriakidis, D. Roberts, and E. Vidal, “Improvingspatial distribution estimation of forest biomass with geostatistics: A casestudy for Rondônia, Brazil,” Ecol. Model., vol. 205, no. 1/2, pp. 221–230,2007.
[26] R. Nishii, S. Kusanobu, and S. Tanaka, “Enhancement of low spatialresolution image based on high resolution-bands,” IEEE Trans. Geosci.Remote Sens., vol. 34, no. 5, pp. 1151–1158, Sep. 1996.
[27] E. Pardo-Igúzquiza, M. Chica-Olmo, and P. M. Atkinson, “Downscalingcokriging for image sharpening,” Remote Sens. Environ., vol. 102, no. 2,pp. 86–98, May 2006.
[28] P. Kyriakidis and E.-H. Yoo, “Geostatistical prediction and simulation ofpoint values from areal data,” Geogr. Anal., vol. 37, no. 2, pp. 124–151,Apr. 2005.
[29] E. Pardo-Igúzquiza and P. M. Atkinson, “Modelling the semivariogramsand cross-semivariograms required in downscaling cokriging by nu-merical convolution-deconvolution,” Comput. Geosci., vol. 33, no. 10,pp. 1273–1284, Oct. 2007.
[30] E. Pardo-Iguzquiza, V. F. Rodríguez-Galiano, M. Chica-Olmo, andP. M. Atkinson, “Image fusion by spatially adaptive filtering using down-scaling cokriging,” ISPRS J Photogramm., vol. 66, no. 3, pp. 337–346,May 2011.
Marcio H. Ribeiro Sales received the B.Sc. degree in statistics from theFederal University of Pará (UFPA), Belém, Brazil, in 2002 and the M.A. degreein geography from the University of California, Santa Barbara, in 2010. Heis currently working toward the Ph.D. degree in the Production Ecology andResource Conservation Graduate School, Wageningen University, Wageningen,The Netherlands.
Since 2002, he has been a Researcher with the Instituto do Homem e doMeio Ambiente da Amazonia—Imazon, Belém, where he conducts researchon geostatistical modeling applied to mapping biomass and deforestation riskmodeling. His most recent works include the modeling of uncertainty in CO2
emissions and estimates of the annual risk of deforestation in the Amazon.
Carlos M. Souza, Jr. received the B.S. degree in geology from the FederalUniversity of Pará (UFPA), Belém, Brazil, the M.Sc. degree in soil sciencefrom Pennsylvania State University, University Park, and the Ph.D. degree ingeography from the University of California, Santa Barbara.
Since 1992, he has been with the Instituto do Homem e do Meio Ambienteda Amazonia—Imazon, Belém, a nongovernment organization dedicated topromote sustainable development in the Amazon through research, policyanalysis, dissemination, and capacity building, where he became a SeniorResearcher in 2001 and was the Executive Director from 2005 through 2008.He is currently engaged in the development of open-access geoweb platformsto promote collaborative networking to monitor and protect forests, such asDeforestation Alert System Earth Engine and ImazonGeo. His research focuseson spatial analysis to identify priority areas for conservation in the Amazon,mapping of carbon stocks of forests, and the development of remote sensingtechniques to map and monitor land-cover and land-use changes, particularlyforest degradation caused by selective logging, fires, and fragmentation.
Phaedon C. Kyriakidis received the B.Sc. degree in geology from the AristotleUniversity of Thessaloniki, Thessaloniki, Greece, in 1994 and the Ph.D. degreein geological and environmental sciences (specializing in geostatistics in theEarth Sciences) from Stanford University, Stanford, CA, in 1999.
From 1999 to 2000, he was a Postdoctoral Fellow with the Earth SciencesDivision, Lawrence Berkeley National Laboratory. In 2001, he joined the Fac-ulty of the Department of Geography, University of California, Santa Barbara,where he is currently a Professor of geographic information science. He is alsoa Professor of spatial analysis with the Department of Geography, University ofthe Aegean, Mytilene, Greece. His research interests include geostatistics andspatial analysis, geographic information systems and science, remote sensing,and environmental modeling.
