IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. 14, NO.

12, DECEMBER 2017 2255

Haze Correction for Contrast-Based

Multispectral Pansharpening
Simone Lolli , Luciano Alparone , Andrea Garzelli , Senior Member, IEEE,
and Gemine Vivone , Member, IEEE

Abstract— In this letter, we show that pansharpening of analysis (MRA) [1]. The Pan image is preliminarily histogram-
visible/near-infrared (VNIR) bands takes advantage from a cor- matched, that is, radiometrically transformed by constant gain
rection of the path-radiance term introduced by the atmosphere and offset in such a way that its low-pass version exhibits mean
during the fusion process. This holds whenever the fusion and variance the same as the component that shall be replaced,
mechanism emulates the radiative transfer model ruling the for CS methods, or the MS band that shall be sharpened, for
acquisition of the Earth’s surface from space, that is, for methods
MRA methods. The injection model rules the combination of
exploiting a contrast-based injection model of spatial details
extracted from the panchromatic (Pan) image into the inter- the low-pass MS image with the spatial detail extracted from
polated multispectral (MS) bands. Such methods are high-pass Pan. Such a model is stated between each of the resampled
modulation (HPM), Brovey transform, synthetic variable ratio MS bands and the low-pass version of the Pan image.
(SVR), University of New Brunswick pansharp, smoothing filter- The most popular injection models are: 1) the projection
based intensity modulation, and spectral distortion minimization. model, which may be derived from the Gram-Schmidt (GS)
The path radiance should be estimated and subtracted from orthogonalization procedure, representing the basis of the
each band before the product by Pan is accomplished and added GS spectral sharpening [2] and of the context-based deci-
back after. Both empirical and model-based estimation techniques sion (CBD) [3] and 2) the multiplicative or contrast-based
of MS path radiances are compared within the framework of or modulation-based model, which is the basis of such
optimized SVR and HPM algorithms. Simulations carried out
on QuickBird and IKONOS data highlight that haze correction
techniques as high-pass modulation (HPM) [4], Brovey
of MS before fusion is always beneficial, especially on vegetated transform (BT) [5], synthetic variable ratio (SVR) [6],
areas and in terms of spectral quality. University of New Brunswick (UNB) pansharp [7], smoothing
filter-based intensity modulation (SFIM) [8], and spectral
Index Terms— Haze, image fusion, multispectral (MS) distortion minimizing (SDM) injection model [9]. Unlike the
pansharpening, path radiance, radiative transfer model, remote projection model, which may be either global, as for GS, or
sensing.
local, as for CBD, the contrast-based model is inherently local,
I. I NTRODUCTION because the injection gain changes at each pixel [10].
Although considerations on atmospheric effects were
P ANSHARPENING techniques take advantage of the com-
plementary spatial and spectral resolutions of multi-
spectral (MS) and panchromatic (Pan) data to synthesize a
already present in SVR [6] and unspecified empirical adjust-
ments in the baseline of UNB pansharp [7], the paper that
unique product that exhibits as many spectral bands as the introduced SFIM [8] first gave an interpretation of the mul-
MS image, each with the same spatial resolution as the tiplicative injection model in terms of the radiative transfer
Pan image. model ruling the acquisition of an MS image from a real-
After the MS bands have been interpolated and coregistered world scene [11]: a low spatial resolution spectral reflectance,
to the Pan image, spatial details are extracted from Pan preliminarily estimated from the MS bands, is sharpened via
and added to the MS bands according to a certain injection multiplication by the high spatial resolution Pan image.
model. The detail extraction step can follow the spectral Currently, very few authors have ever considered, and
approach, originally known as component substitution (CS), never explicitly, the path radiance of each MS band, which
or the spatial approach, which may rely on multiresolution is an undesired energy scattered by different atmospheric
constituents that reaches the aperture of the instrument without
Manuscript received July 11, 2017; revised September 12, 2017; accepted being reflected by the Earth’s surface. The path radiance,
September 27, 2017. Date of publication October 31, 2017; date of current which appears as a haze in a true-color display, should be
version December 4, 2017. (Corresponding author: Luciano Alparone.) estimated and subtracted from each band before modulation
S. Lolli is with the Institute of Methodologies for Environmental Analysis,
National Research Council, 85050 Tito Scalo, Italy, and also with NASA
and reinserted later, to restore the unbiased sharpened image.
Goddard Space Flight Center Joint Center for Earth Systems Technology, In this letter, several methods for estimating the path
Greenbelt, MD 20771 USA (e-mail: simone.lolli@cnr.it; slolli@umbc.edu). radiance are presented; noteworthy, the Fu–Liou–Gu (FLG)
L. Alparone is with the Department of Information Engineering, University radiative transfer model [12] has been considered. Experiments
of Florence, 50139 Florence, Italy (e-mail: luciano.alparone@unifi.it). on QuickBird and IKONOS images, in which the bandwidth
A. Garzelli is with the Department of Information Engineering
and Mathematics, University of Siena, 53100 Siena, Italy (e-mail: of Pan encompasses the four spectral bands, demonstrate that
andrea.garzelli@unisi.it). the advantages of removing the estimated path radiances for
G. Vivone is with the Department of Information Engineering, Electrical calculating the injection model are more consistent in terms
Engineering and Applied Mathematics, University of Salerno, 84084 Fisciano, of spectral quality (color hues) with respect to spatial quality
Italy (e-mail: gvivone@unisa.it).
Color versions of one or more of the figures in this letter are available
(geometric sharpness). Significant improvements are notable
online at http://ieeexplore.ieee.org. on vegetated areas and less relevant on other landscapes.
Digital Object Identifier 10.1109/LGRS.2017.2761021 A combination of empirical and statistical path-radiance
1545-598X © 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
2256 IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. 14, NO. 12, DECEMBER 2017

estimators attains the performance of the model and of an Under the hypothesis of substitution of a single component
exhaustive search performed at a degraded spatial scale. that is a linear combination of the input bands, the fusion
process can be obtained without the explicit calculation of the
II. R ADIATIVE T RANSFER M ODEL forward and backward transformations, but through a proper
injection scheme [1], thereby leading to fast implementations
The radiative transfer model relates the at-sensor spectral
radiance to the surface reflectance, top-of-atmosphere (TOA) of CS methods, whose general formulation is given by
solar irradiance, a.k.a. solar constant, atmospheric transmit- M̂k = M̃k + Gk · (P − I L ) , k = 1, . . . , N (3)
tances and upward scattered radiance, a.k.a. path radiance [11]
in which k is the band index, G = [G1 , . . . , Gk , . . . , G N ] the
ρ(λ)τu (λ)(E s (λ) cos(θ )τd (λ) + E d (λ))
L(λ) = 2 ·π
+ L P (λ) array of injection gains, while the intensity, I L , is defined as
dES
(1) 
N
IL = wi · M̃i (4)
in which λ is the wavelength of the electromagnetic i=1
radiation [μm]; L(λ) is the at-sensor spectral radiance in which the weight vector w = [w1 , . . . , wi , . . . , w N ]
[W·m−2 ·sr−1 ·μm−1 ]; ρ(λ) is the surface reflectance [unitless]; is the first row of the forward transformation matrix.
τu (λ) is the upward transmittance of atmosphere [unitless]; Alternatively [2], a multivariate linear regression is exploited
E s (λ) is the mean TOA solar irradiance [W·m−2 ·μm−1 ]; θ is to model the relationship between the low-pass-filtered
the solar zenith angle [degrees]; τd (λ) is the downward Pan, P L , and the interpolated MS bands
transmittance of atmosphere [unitless]; E d (λ) is the diffuse
irradiance at the surface [W · m−2 · μm−1 ]; dES is the Earth– 
N
Sun distance [astronomical units] (dES ≈ 1); and L P (λ) is the PL = ŵi · M̃i +  = Î L +  (5)
upward scattered radiance [W · m−2 · sr−1 · μm−1 ]. i=1
Estimation of surface reflectance requires preliminary cor-
rection of the offset (path radiance) of the kth spectral band, in which Î L is the optimal intensity component and  is the
L P (k), corresponding to a certain wavelength interval, and minimum-variance space-varying residue. The set of space-
then rescaling of all pixels by the product of the atmospheric constant optimal weights {ŵk }k=1,...,N is calculated as the
upward transmittance, τu (k), and by the total solar irradiance minimum mean square error (MMSE) solution of (5).
measured in the kth spectral interval of the instrument, namely, In (3), defining space varying injection gains, G, such that
E t (k)  (E s (k) · cos(θ ) · τd (k) + E d (k))/dES
2
M̃k
Gk = , k = 1, . . . , N (6)
(L(k) − L P (k)) · π IL
ρ(k) = , k = 1, . . . , N (2)
τu (k) · E t (k) yields
in which ρ(k)/π is the average of a Lambertian bidirectional M̃k P
reflectance distribution function, with maximum ρ(k). M̂k = M̃k + · (P − I L ) = M̃k · , k = 1, . . . , N (7)
In applications concerning different acquisition dates, e.g., IL IL
change detection [13] and multitemporal pansharpening [14], which, in the case of spectral weights all equal to 1/N, is the
sunlight and atmospheric corrections should be performed widely known BT pansharpening method [5]. An evolution of
according to (2). MS pansharpening, which produces a BT was SVR [6], which used a modified atmospheric model
sharp MS image having the same format as the original that accepts target reflectance and relative spectral response
MS image [15], generally does not require corrections, unless and measured the representative reflection spectra of the five
a multiplicative detail-injection model is exploited [6], [8]. land-cover classes of urban, soil, water, trees, and grass. The
parameters {wk } were then obtained through regression analy-
III. S PECTRAL AND S PATIAL F USION M ETHODS sis of the values simulated using the atmospheric model and
The math notation used hereafter is detailed in the follow- the five classes. After construction of I L , a linear histogram
ing. Vectors are indicated in bold lowercase (e.g., x) with the adjustment had to be performed to force the Pan image to
i th element indicated as x i . 2-D and 3-D arrays are expressed match the I L image in order to eliminate atmospheric and
in bold uppercase (e.g., X). An MS image M = {Mk }k=1,...,N illumination differences. UNB pansharp [7] exploits a multi-
is a 3-D array composed by N bands indexed by the subscript variate regression of original Pan to interpolated MS bands
k = 1, . . . , N; accordingly, Mk indicates the kth band of M. to yield the set of {wk },k=1,...,N . Thus, BT, SVR, and UNB
The Pan image is a 2-D matrix and will be indicated as P pansharp fit (3) with the choice of injection gains (6).
and its histogram-matched version with P̌. Also, M̃k and M̂k
indicate interpolated and sharpened MS bands, respectively. B. Spatial or Multiresolution Analysis Methods
Unlike conventional matrix product and ratio, such operations The spatial approach relies on the injection of high-pass
are intended as product and ratio of terms of the same positions spatial details of Pan into the resampled MS bands [16].
within the array. The most general MRA-based fusion may be stated as

A. Spectral or Component Substitution Methods M̂k = M̃k + Gk · (P − P L ) , k = 1, . . . , N. (8)

The class of CS, or spectral, techniques is based on the According to (8), the different approaches and methods
projection of the MS image into another vector space assuming belonging to this class are uniquely characterized by the low-
that the forward transformation splits the spatial structure and pass filter employed for obtaining the image P L and by the
the spectral diversity into separate components. set of space-varying injection gains, {Gk }k=1,...,N .
LOLLI et al.: HAZE CORRECTION FOR CONTRAST-BASED MS PANSHARPENING 2257

The contrast-based version of MRA pansharpening is V. E XPERIMENTAL R ESULTS

A. Methods
M̃k P
M̂k = M̃k + · (P − P L ) = M̃k · , k = 1, . . . , N. (9) Path-radiance correction has been considered for two opti-
PL PL mized contrast-based methods: one relying on CS and another
Equation (9) accommodates HPM [4], SFIM [8], and SDM [9], on MRA [15]. The two methods with path-radiance correc-
which differ from one another by the filter to achieve P L . tion are labeled CSw/PRC (10) and MRAw/PRC (13). The
two versions without path-radiance correction, CSw/oPRC
IV. C ONTRAST-BASED F USION W ITH H AZE C ORRECTION and MRAw/oPRC, are given by (10) and (13), respectively,
with L P (k) = 0, ∀k. The method labeled Exp denotes plain
In this section, path radiance correction is introduced interpolation without injection of details.
in (7) and (9) in order to produce accurate estimates of low
spatial resolution spectral reflectance, which is the key to
contrast-based pansharpening, both CS and MRA. B. Data Sets
Starting from the general model of contrast-based 1) QuickBird-Trento: A QuickBird image has been
CS fusion (7), the haze-corrected version is given by acquired over an industrial area near Trento, in Italy,
in October 2005. The spatial sampling interval (SSI) of
P̌Î∗ the geo-coded products is 2.8 m for MS (blue (B),
M̂k = (M̃k − L P (k)) · ∗
L
+ L P (k), k = 1, . . . , N (10) green (G), red (R), and near infrared (NIR)) and 0.7 m
Î L for Pan, respectively. The data format is spectral radiance
in which [W · m−2 · sr−1 · μm−1 ] with 11 b.
∗   2) IKONOS–Toulouse: An IKONOS image has been
Î L = ŵi · (M̃i − L P (i )) = Î L − ŵi · L P (i ) (11) acquired on the urban area of Toulouse, in France, on
i i May 15, 2000. The SSI is 4 m for MS (blue, green, red,
and NIR) and 1 m for Pan, respectively. The format of the
and P̌Î∗ is P histogram-matched to Î∗L
L data is spectral radiance with 11-b word-length.
σÎ∗
P̌Î∗  (P − μP ) · L
+ μÎ∗ . (12) C. Assessments
L σP L L
Since the path radiance arguably does not depend on the
μ and σ denote the mean and square root of variance, spatial resolution, quality evaluations have been carried out at
respectively. a spatial scale four times greater than that of the original Pan,
Starting from the general model of contrast-based checking the synthesis property of Wald’s protocol [17]. The
MRA fusion (9), the haze-corrected version is given by availability of reference originals entails the use of widespread
(k) vector (dis)similarity scores: spectral angle mapper (SAM),
P̌ erreur relative globale adimensionnelle de synthèse (ERGAS),
M̂k = (M̃k − L P (k)) · (k)
+ L P (k), k = 1, . . . , N (13)
P̌ L and Q4 [3], the four-band extension of the universal image
quality index, widely referred to as Q.
in which P̌(k) is P histogram-matched to M̃k − L P (k) and P̌(k)
L
its low-pass-filtered version D. Estimation of Path Radiances
σM̃k The haze-corrected versions of CS (10) and MRA (13)
P̌(k)  (P − μP ) · + μM̃k − L P (k), k = 1, . . . , N.
σP L require path radiances estimated for each band. Several
(14) approaches are feasible. The path radiance is arguably a
fraction of the minimum of spectral radiance attained over
Equations (10) and (13) may be easily explained by watch- the scene. If the scene is large enough and hence statistically
ing (2), in which the radiative transfer model is inverted consistent, setting the actual minimum equal to the 1 percentile
to yield the spectral reflectance product. Accordingly, the of the histogram ensures robustness to noise and outliers.
reflectance is given by the spectral radiance minus the path The path radiance of the B-band may be approximated by the
radiance (offset) divided by a gain that is the product of 1-p-tile. Once the L P (B) is known, the intercept of the
upward atmospheric transmittance by total solar irradiance. G-to-[B–L P (B)] scatterplot yields an estimate of L P (G);
Each pixel of M̃k is a sample of spectral radiance. The analogously for the R channel, L P (R) may be estimated
offset is L P (k) and is assumed to be constant over the from the R-to-[G–L P (G)] scatterplot. This empirical/statistical
scene. Each pixel either of the MMSE intensity calculated approach holds only for the visible bands. For calculating the
from the path radiance-free MS bands, Î∗L , or of the low- path radiance of NIR, which is uncorrelated with the visible
resolution Pan image histogram-matched to the kth MS band, bands, the scatterplot method fails. A reasonable physical
(k) approximation is that the L P of NIR is set equal to zero.
P̌ L , measures the solar irradiance reported at the aperture
of the instrument, that is, the denominator of (2). In this Also a modeling of the atmosphere was considered. The
way, upward atmospheric transmittance and solar irradiance FLG radiative transfer model [12] requires acquisition year,
are eliminated [6], [8]. Once a map of low spatial resolution month, day, local time, longitude, latitude, and possibly type
spectral reflectance is obtained, it is sharpened by multiply- of landscape (urban or rural) for setting aerosols [18]. Such
ing its pixels by a high spatial resolution irradiance, that a model directly yields values of path radiance in predefined
is, either Pan histogram matched to the path radiance-free bands, roughly corresponding to those of MS scanners, like
MMSE intensity or Pan histogram matched to the kth path Landsat 8 OLI. With small adjustments to fit the R and NIR
radiance-free interpolated MS band. bands of QuickBird and IKONOS, it was found that for the
2258 IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. 14, NO. 12, DECEMBER 2017

Fig. 1. True-color details of QuickBird (256 × 256, at 2.8 m). (a) 11.2-m degraded MS expanded to 2.8 m. (b) 0.7-m Pan degraded to 2.8 m.
(c) Reference 2.8-m MS. (d) CSw/PRC. (e) CSw/oPRC. (f) MRAw/PRC. (g) MRAw/oPRC.

TABLE I The benefits of the path-radiance correction are significant

S CORES OF D EGRADED Q UICK B IRD T RENTO D ATA for all three indexes (SAM, ERGAS, and Q4), but partic-
ularly remarkable for the spectral angle error. Interestingly,
while there may be differences in performance between
CS and MRA without correction, the corrected version per-
forms approximately the same, for all indexes.
This is another clue that CS and MRA fusion methods,
once they have been optimized to the physical instrument and
imaging model, should be identical in performance.
The QuickBird images at 2.8-m scale are portrayed
in Fig. 1. What immediately stands out is that the synthesis
of vegetated areas is much more accurate for the haze-
two data sets, the modeled path radiance is well approximated corrected methods. Without correction [Fig. 1(e) and (g)], the
by 95% of the 1-p-tile of B, 45% of the 1-p-tile of G, texture of the canopies, which appears in Pan [Fig. 1(b)],
40% of the 1-p-tile of R, and 5% of the 1-p-tile of NIR. being originated from the NIR wavelengths, but is much less
The results reported in Section V-E are the best attainable notable in the bands of the reference ground truth covering
varying with the estimated path radiances. An exhaustive the visible spectrum [Fig. 1(c)], is transplanted in the fusion
search at steps of one [W · m−2 · sr−1 · μm−1 ] was per- products and originates an unlikely over-enhancement, which
formed for each data set: the optimal path radiances are is responsible for the remarked loss of performance. The
those that optimize the three fusion scores, in average, at visual quality of nonvegetated areas is generally good for all
the degraded spatial scale, i.e., when the ground truth is methods. While the difference between corrected and non-
available. With modeled path radiances, the performance is corrected methods is substantial, the difference between the
about 0.1 % lower, thereby indicating that: 1) the radiative CS and MRA approaches, both with or without correction, is
transfer model actually rules the performance of contrast- imperceivable.
based MS pansharpening and 2) the accuracy of path radi- 2) IKONOS–Toulouse: The IKONOS image portraying the
ance estimation is not crucial for the fusion performance, at city center of Toulouse exhibits a very low amount of vegetated
least if the format of the pansharpened product is a spectral areas. This explains the fusion scores in Table II. Notwith-
radiance. standing the increment in performance due to the correction of
haze is moderate, in average, all corrected methods are better
E. Fusion Simulations than their basic versions. SAM benefits more than ERGAS
1) QuickBird–Trento: Table I reports the scores achieved by and Q4 from the correction, especially for CS, in which
the two CS and MRA contrast-based methods, with and with- SAM passes from 4.84° to 3.01°.
out path-radiance correction. The interpolated low-resolution Also fusion tests at full scale (1 m) have been performed.
MS and the reference high-resolution MS (ground truth) are Small icons of a vegetated square are shown in Fig. 2.
included in the comparison as Exp and Ref, respectively. Notwithstanding a reference is missing, the difference between
LOLLI et al.: HAZE CORRECTION FOR CONTRAST-BASED MS PANSHARPENING 2259

TABLE II of its histogram is consistent to a model based path-radiance

S CORES OF D EGRADED IKONOS T OULOUSE D ATA estimation and to an exhaustive search of the set of L P (k) that
maximizes fusion performance at a degraded scale. Both visual
and objective assessments highlight improvements, mostly in
terms of spectral angle and in vegetated areas. The procedure
may need adjustments for WorldView-2 data, in which the two
outermost bands are not encompassed by Pan.

R EFERENCES
[1] L. Alparone, B. Aiazzi, S. Baronti, and A. Garzelli, Remote Sensing
Image Fusion. Boca Raton, FL, USA: CRC Press, 2015.
[2] B. Aiazzi, S. Baronti, M. Selva, and L. Alparone, “Enhanced
Gram–Schmidt spectral sharpening based on multivariate regression of
MS and Pan data,” in Proc. IGARSS, Jul./Aug. 2006, pp. 3806–3809.
[3] L. Alparone, L. Wald, J. Chanussot, C. Thomas, P. Gamba, and
L. M. Bruce, “Comparison of pansharpening algorithms: Outcome of the
2006 GRS-S data-fusion contest,” IEEE Trans. Geosci. Remote Sens.,
vol. 45, no. 10, pp. 3012–3021, Oct. 2007.
[4] R. A. Schowengerdt, Remote Sensing: Models and Methods for Image
Processing, 2nd ed. Orlando, FL, USA: Academic, 1997.
[5] A. R. Gillespie, A. B. Kahle, and R. E. Walker, “Color enhancement of
highly correlated images. II—Channel ratio and ‘chromaticity’ transfor-
mation techniques,” Remote Sens. Environ., vol. 22, no. 3, pp. 343–365,
Aug. 1987.
[6] C. K. Munechika, J. S. Warnick, C. Salvaggio, and J. R. Schott, “Reso-
lution enhancement of multispectral image data to improve classification
accuracy,” Photogramm. Eng. Remote Sens., vol. 59, no. 1, pp. 67–72,
Jan. 1993.
[7] Y. Zhang, “A new merging method and its spectral and spatial effects,”
Int. J. Remote Sens., vol. 20, no. 10, pp. 2003–2014, Jul. 1999.
[8] J. G. Liu, “Smoothing Filter-based Intensity Modulation: A spectral
preserve image fusion technique for improving spatial details,” Int. J.
Remote Sens., vol. 21, no. 18, pp. 3461–3472, Dec. 2000.
[9] L. Alparone, B. Aiazzi, S. Baronti, and A. Garzelli, “Sharpening of very
high resolution images with spectral distortion minimization,” in Proc.
IEEE IGARSS, Jul. 2003, pp. 458–460.
[10] G. Vivone, R. Restaino, M. Dalla Mura, G. Licciardi, and J. Chanussot,
“Contrast and error-based fusion schemes for multispectral image
pansharpening,” IEEE Geosci. Remote Sens. Lett., vol. 11, no. 5,
pp. 930–934, May 2014.
Fig. 2. True-color details of IKONOS at 1 m and related spatial distortion [11] F. Pacifici, N. Longbotham, and W. J. Emery, “The importance of
values, D S . (a) CSw/PRC. (b) CSw/oPRC. (c) MRAw/PRC. (d) MRAw/oPRC. physical quantities for the analysis of multitemporal and multiangular
optical very high spatial resolution images,” IEEE Trans. Geosci. Remote
Sens., vol. 52, no. 10, pp. 6241–6256, Oct. 2014.
[12] Q. Fu and K. N. Liou, “On the correlated k-distribution method for
corrected and uncorrected methods stands out. The haze radiative transfer in nonhomogeneous atmospheres,” J. Atmos. Sci,
correction, most relevant in the blue band, prevents textures vol. 49, no. 22, pp. 2139–2156, Nov. 1992.
of tree canopies, collected by the broadband Pan in the NIR [13] F. Bovolo, L. Bruzzone, L. Capobianco, A. Garzelli, S. Marchesi, and
wavelengths, from being injected in the visible bands, thereby F. Nencini, “Analysis of the effects of pansharpening in change detection
originating a blueish texture. Values of the full-scale spatial on VHR images,” IEEE Geosci. Remote Sens. Lett., vol. 7, no. 1,
pp. 53–57, Jan. 2010.
distortion index, D S , of the QNR protocol (see [17]) highlight [14] B. Aiazzi, L. Alparone, S. Baronti, R. Carlà, A. Garzelli, and L. Santurri,
that the sharpness of the uncorrected versions is unlikely. “Sensitivity of pansharpening methods to temporal and instrumental
changes between multispectral and panchromatic data sets,” IEEE Trans.
VI. C ONCLUSION Geosci. Remote Sens., vol. 55, no. 1, pp. 308–319, Jan. 2017.
[15] L. Alparone, A. Garzelli, and G. Vivone, “Intersensor statistical match-
In this letter, we pointed out that the step of MS path- ing for pansharpening: Theoretical issues and practical solutions,” IEEE
radiance estimation and correction is the key to attain Trans. Geosci. Remote Sens., vol. 55, no. 8, pp. 4682–4695, Aug. 2017.
improved performance from pansharpening methods based on [16] L. Alparone, S. Baronti, B. Aiazzi, and A. Garzelli, “Spatial methods
spatial modulation of pixel spectra. Whenever the bandwidth for multispectral pansharpening: Multiresolution analysis demystified,”
IEEE Trans. Geosci. Remote Sens., vol. 54, no. 5, pp. 2563–2576,
of Pan encompasses the MS bands, an optimized intensity May 2016.
component can be achieved through multivariate regression of [17] F. Palsson, J. R. Sveinsson, M. O. Ulfarsson, and J. A. Benediktsson,
low-pass-filtered Pan to MS and optimal injection of spatial “Quantitative quality evaluation of pansharpened imagery: Consistency
details can be achieved through haze estimation and removal versus synthesis,” IEEE Trans. Geosci. Remote Sens., vol. 54, no. 3,
pp. 1247–1259, Mar. 2016.
from MS bands. The path radiance of each band, L P (k), [18] S. K. Pani et al., “Assessment of aerosol optical property and radiative
is estimated, subtracted before the spatial modulation, and effect for the layer decoupling cases over the northern South China
reinserted after fusion. The empirical trick of setting the path Sea during the 7-SEAS/Dongsha experiment,” J. Geophys. Res.-Atmos.,
radiance of the kth band equal to a fraction of the 1 percentile vol. 121, no. 9, pp. 4894–4906, 2016.