This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING 1
Haze Detection and Removal in Remotely
Sensed Multispectral Imagery
Aliaksei Makarau, Rudolf Richter, Rupert Müller, and Peter Reinartz
Abstract—Haze degrades optical data and reduces the accuracy
of data interpretation. Haze detection and removal is a challenging
and important task for optical multispectral data correction. This
paper presents an empirical and automatic method for inhomoge-
neous haze detection and removal in medium- and high-resolution
satellite optical multispectral images. The dark-object subtraction
method is further developed to calculate a haze thickness map,
allowing a spectrally consistent haze removal on calibrated and
uncalibrated satellite multispectral data. Rare scenes with a uni-
form and highly reflecting landcover result in limitations of the
method. Evaluation on hazy multispectral data (Landsat 8 OLI
and WorldView-2) and a comparison to haze-free reference data
illustrate the spectral consistency after haze removal.
Index Terms—Haze removal, Landsat 8 OLI, spectral consis-
tency, WorldView-2.
I. INTRODUCTION
H
AZINESS is a common artifact in optical remotely
sensed data originating from fractions of water vapor,
ice, fog, sand, dust, smoke, or other small particles in the at-
mosphere. Haze transparency leaves an opportunity for optical
imagery restoration; nevertheless, for an inhomogeneous and
structured haziness (see the example shown in Fig. 1), a precise
detection is necessary, and a haze thickness map (HTM) has to
be calculated. An efficient method for multi- and hyperspectral
remotely sensed optical data dehazing handling all types of
haziness and haze thickness is still a challenge.
Richter [1] developed a dehazing method accurately dealing
with haze transition (haze–nonhaze) regions. Clear, hazy, and
cloud regions are separated and stored. A haze boundary region
is introduced to generate a smoother transition from haze to
clear areas. The algorithm matches the histogram of the haze
areas to the histogram of the clear part leading to dehazing.
Moro et al. [2] presented a framework to reduce haze and
calculate some vegetation indices. Haze removal is performed
based on the haze optimized transform (HOT; see original
work of [3]) technique. The improved HOT consists of the
following: determination of a man-made feature mask, masking
of estimated haze for water bodies and man-made features, in-
terpolation of masked-out areas, and final haze subtraction from
the image. Dark-object subtraction (DOS) is then performed
Manuscript received August 17, 2013; revised October 2, 2013 and
October 31, 2013; accepted November 21, 2013.
The authors are with the German Aerospace Center (DLR), 82234 Wessling,
Germany (e-mail: aliaksei.makarau@dlr.de).
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.2013.2293662
on the dehazed scene in order to estimate approximate soil
reflectances necessary to calculate the vegetation indices.
Liu et al. [4] developed a technique to remove spatially vary-
ing haze contamination, comprising three steps: haze detection,
haze perfection, and haze removal. The background suppressed
haze thickness index (BSHTI) is used to indicate relative haze
thickness. The virtual cloud point method based on BSHTI
is used in haze removal. Human intervention is necessary to
outline thick haze and clear regions. Here, 76 paired (dehazed
and haze-free) regions were compared, resulting in a correlation
coefficient higher than 0.7.
He et al. [5] proposed an efficient haze removal and depth
map estimation method for outdoor colored RGB images. The
method is based on dark channel prior calculation (a search
of local patches in the image containing pixels with very low
intensity). The hazy image dark channel values are mainly
contributed by the airlight (path radiance), and the dark pixels
can provide an estimation of the haze transmission. A haze
imaging model allows us to recover a high-quality haze-free
image and a corresponding depth map.
Chavez [6], [7] proposed an improvement for DOS technique
to correct optical data for atmospheric scattering. The method
assumes to find a pixel or several pixels with a very low
reflectance. Due to the haziness in the scene, these pixels are not
completely dark, and the digital numbers (DNs) in the pixels
are nonzero. The values of the DNs are assumed to be the
haze thickness in the image. Assuming a constant haze over
the image, a subtraction of the DN value from the whole image
allows us to perform dehazing.
In this paper, we present a further development of the DOS
technique by Chavez [6]. Instead of searching only several dark
objects in the whole scene, it is proposed to search dark objects
locally in the whole image. A local search of the dark objects
allows us to construct an HTM. Assuming an additive model
of the haze influence, a subtraction of the HTM from a hazy
image allows us to restore the haze-free signal at the sensor. The
dehazed data can be later employed for atmospheric correction
[8] and for further processing like landcover classification [9]–
[11] or change detection [12].
This paper is organized as follows. In Section II, the new
haze removal algorithm is presented. The method of the HTM
calculation as well as a correction for large highly reflective ob-
jects is given. Haze thickness estimation per band is developed,
allowing multispectral image dehazing. Section III contains a
description of hazy and haze-free data, data properties, and ex-
amples of dehazed data evaluation. The conclusion and possible
further developments of the method are given at the end of this
paper.
0196-2892 © 2013
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.
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
2 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING
Fig. 1. Typical types of haziness in optical imagery. (a) Uniform haze (WorldView-2). (b) Moderately thick haze (WorldView-2). (c) Moderately thin haze
(WorldView-2). (d) Thick structured haze (AVNIR).
II. PROPOSED METHOD FOR HAZE REMOVAL
A. Haze Model
Small field-of-view optical systems are employed to acquire
medium, high, or very high spatial resolution multispectral data.
Then, the distance from any object in the scene to the sensor
can be assumed as constant. Our simplified model (assumption
is verified in the experiments) can be described as follows (see
the additive model by Chavez [6]):
L
sensor
= L
0
+ HR (1)
where L
sensor
is the acquired radiance, L
0
is the sum of the path
radiance and the surface reflected radiance, and HR is the haze
contribution.
Assuming a linear nature of the DN to radiance conversion,
the equation is adapted to band DN
DN
sensor
i
(x, y)=DN
i
(x, y)+HR
i
(x, y) (2)
where (x, y) are the coordinates of a pixel in an optical image,
i is the band number, DN
sensor
i
(x, y) is the acquired DN, and
DN
i
(x, y) is the DN without the influence of haze. The term
HR
i
(x, y) is dependent on the haze thickness and can be
defined as HTM HTM
i
(x, y).
B. HTM
Since the HTM is employed as an additive term (2), a
subtraction of the HTM from an optical image allows us to
restore the haze-free image. The task is to calculate HTM
i
(x, y)
for each band i.
Haze thickness is usually variable over a scene [e.g.,
Fig. 1(b)–(d)], and the HTM has to be calculated by searching
dark pixels (shaded objects [6]) over all of the image. A medium
to small ground sampling distance (GSD; 30–2 m) of a sensor
allows us to record shadows from objects. The shading can
be caused by a varying relief in the scene, natural or man-
made objects. These pixels have a very low DN value (small
ground reflectance+path radiance), and with the presence of
haze, the recorded DN value is also nonzero and higher than
the corresponding haze-free contribution. These nonzero DN
values can be employed to estimate the HTM of the image. It
should be noted that the HTM includes the thickness of the haze
and the aerosol thickness. The compensation for the subtracted
clear scene aerosol thickness is given in the dehazing section
(Section II-E).
The dark pixels are searched using a local nonoverlapping
window w (i.e., area w × w pixels). Pixelwise search for dark
pixels can result in an occasional selection of pixels over
nonshaded or bright objects, and these pixels cannot be used
to estimate haze thickness (the HTM will contain errors). The
window is employed to increase the chance to locate dark
pixels over shaded regions. The other pixels in the window
are assumed to have the same haze thickness as the thickness
estimated in the pixel with the minimal value. A larger window
allows easier detection of dark pixels, but the border as well
as the structure (density in each pixel or local area) of haze
can be smoothed, leading to an underestimation of the HTM
(Fig. 2). A smaller window allows better estimation of haze
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
MAKARAU et al.: HAZE DETECTION AND REMOVAL IN REMOTELY SENSED MULTISPECTRAL IMAGERY 3
Fig. 2. HTM preciseness dependence on the nonoverlapping window size. (a) Example of an image (Landsat 8 OLI, RGB true color combination) with
inhomogeneous haze (thick in the middle but moderately thin over all of the subimage). (b) HTM calculated using the window size w =5. (c) HTM calculated
using the window size w =21.
Fig. 3. Example of the HTM calculation for an area containing large objects with moderately high reflectance (soils). The HTM is corrected for the regions with
the bright objects. (a) Landsat 8 OLI RGB true color band combination (haze noticeable at the upper part of the image). (b) Mask for large bright objects.(c)HTM
(large bright objects can be confused with haze). (d) Corrected HTM (the bright object pixels are interpolated by triangulation).
structure together with a possibility that outside the haze pure
dark pixels are not found (and the pixels with a moderately high
intensity are selected). The size of the window is set according
to the GSD of the sensor: the smaller the GSD, the smaller
size of the window can be set, and a more precise HTM can
be estimated. Usually, the size of the window is constant for
a selected sensor. By the experimental analysis, it is found that
the minimum window size should be three times higher than the
sensor GSD, e.g., Landsat 8 OLI (30-m GSD for multispectral
data) should have a minimum window size of 90 m × 90 m, and
WorldView-2 (2-m GSD for multispectral data) should have a
minimum window size of 6 m × 6 m, respectively. The dark
pixels selected using the window (the size is w × w) are stored
in a matrix with the size w times smaller than the original hazy
image. Then, the matrix is smoothed by a median filter (usually
3 × 3) and interpolated (cubic interpolation) up to the size of
the original band. Therefore, the possible haze overestimation
is suppressed.
The presence of large bright objects can make the selection
of dark pixels difficult using a small window (see Fig. 3).
Such objects can be agricultural fields, regions with bright
sand, uniform snow cover, etc. These objects are labeled by a
segmentation method (in this work, a mean-shift segmentation
is employed [13], but any other segmentation methods can be
also used). The HTM values in the regions of bright objects are
interpolated using a triangulation. In an extreme case having
a uniform landcover in the whole scene with a high reflectance
(e.g., snow), the HTM might be impossible to estimate. Usually,
such cases are rare in remotely sensed data and not of interest
for the use of this methodology.
The HTM calculation depends on the employed spectral band
for the dark pixel search and the size of the nonoverlapping
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
4 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING
Fig. 4. AVNIR-2 data dehazing example (subset). The employment of band 1 (0.42–0.50 μm) for the HTM calculation can lead to overdehazing in the same
band. (a) Original band 1. (b) Dehazed band 1 (overdehazing is notable). (c) Dehazed band 1 (HTM was calculated using a new band produced by extrapolation
of AVNIR-2 band 1 and band 2).
window
HTM(x, y)=DARK_PIX_SEARCH (Band
b
(x, y),w)
(3)
where Band
b
(x, y) is the employed image (usually a blue
band), w is the size of the nonoverlapping window, and
HTM(x, y) is the HTM (unitless). In our experiments, w
usually takes values ranging from 3 to 9. With w =3, rea-
sonable results on medium- and high-resolution data (Land-
sat 7, Landsat 8 OLI, Advanced Visible and Near Infrared
Radiometer (AVNIR), AVIRIS, WorldView-2, and RapidEye)
were achieved.
It is necessary to label hazy and haze-free regions. The
information in haze-free regions is further employed in the
dehazing and compensating of clear scene aerosol. In order
to find haze-free regions, the haze mask is created. The haze
mask is calculated by the generation of an additional HTM
with moderately large window size (GSD ∗ 20, w =21) and a
thresholding of the HTM. The HTM generated with the large
window size is used only to label hazy and haze-free regions.
The coefficient of 20 is estimated experimentally. The thresh-
old value is equal to the mean value (water area under haze
is segmented with the employment of the near infrared band
and excluded), and reasonable results were achieved for all
processed data. The haze mask as well as the water mask does
not necessarily have to be precise, and a mislabeling of objects
in the scene does not influence dehazing (the haze mask is used
only to calculate a relative haze thickness coefficient, and this
is robust to mislabeled outliers). The employment of such large
windows allows us to find dark pixels and estimate the HTM
irrespectively of large bright objects.
In multispectral data, the search of dark pixels should be
performed in a band with a minimal ground reflectance and
a maximal attenuation by haze. A spectral band in the blue
spectral region (0.37–0.49 μm) is most suitable for this purpose.
Spectral bands in the red or near infrared spectral regions ex-
hibit higher ground reflectance for land and are less influenced
by haziness; therefore, they are not suitable in detecting and
properly estimating haze thickness.
C. Band Extrapolation
The employment of a blue band from a multispectral image
for the HTM generation can lead to overdehazing of this band.
Fig. 4 illustrates such a situation on AVNIR-2 data (see [14]). A
per-band analysis should be carried out [compare Fig. 4(a) and
(b)]. In order to overcome this problem, a new synthetic band
is created by a linear extrapolation of two bands of the image
cube:
Band
extrapol
(x, y)
= EXTRAPOL(Band
m
(x, y),Band
n
(x, y)) (4)
where Band
extrapol
(x, y) is the extrapolated band and
Band
m
(x, y) and Band
n
(x, y) are two bands from the data
cube used for an extrapolated band creation. An employment
This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.
MAKARAU et al.: HAZE DETECTION AND REMOVAL IN REMOTELY SENSED MULTISPECTRAL IMAGERY 5
Fig. 5. Estimation of the HTM in Landsat 8 OLI data. (a) RGB true color composite (bands 4, 3, and 2). (b) Band 1 (0.4430 μm). (c) Extrapolated band
Band
extrapol
(λ =0.4233; extrapolated using Landsat 8 OLI bands 1 and 2). (d) HTM calculated on the extrapolated band in (c).
Fig. 6. Example of varying haze thickness in WorldView-2 spectral channels.
(a) Band 1 (0.4283 μm). (b) Band 4 (0.6080 μm). (c) Band 7 (0.8277 μm).
of the extrapolated band yields better dehazing results [see
Fig. 4(a) and (c)]. Usually, the first two bands are selected
(in the case of Landsat 8 OLI, Band
m
is 0.4430 μm, and
Band
n
is 0.4825 μm band). Fig. 5 illustrates an example of
an extrapolated band: band 1 at 0.4430 μm [Fig. 5(b)] and an
extrapolated band (Fig. 5(c); λ =0.4232 μm). The reflectance
of the surface in the extrapolated (a lower wavelength) band is
less, haze thickness is higher in comparison to the VNIR region
bands, and the HTM can be estimated more precisely. In the
case of noise in the interpolated band, a median filtering can be
applied.
D. Haze Thickness Per Band
In order to dehaze a multispectral image, a subtraction of the
HTM has to be done for each spectral band. The derived haze
thickness is different in bands and usually decreases from the
shortest to the longest wavelength band (Fig. 6). Subtraction of
a constant HTM from all of the bands can lead to overdehazing
in the red and near infrared bands and a loss of spectral
properties in the dehazed data. To preserve spectral consistency,
a band-specific HTM (HTM
i
) should be calculated. Taking into
account that spectral channels in the far short-wave infrared
range are only marginally influenced by haze, the dehazing is
performed for the visible, near infrared, and short-wave infrared
bands, terminating the dehazing wavelength at 2.2 μm.
The thickness of the haze in each band (HTM
i
) is estimated
relative to the HTM calculated using the image Band
b
(x, y)
[see (3)]. A temporary HTM HTM
i
(x, y) is calculated for each
band i, and the regression coefficient (slope of the fitted line
in linear regression) (k
i
) of HTM
i
(x, y) versus HTM(x, y) is
stored in an array ( K) of regression coefficients. The regression
coefficient is calculated using the pixel values in the segmented
hazy regions (haze mask is employed). According to haze
thickness, the regression coefficients are expected to have a
decreasing trend for the bands from the shortest to the longest
wavelength. A check if K has a meaningful decreasing trend
and a correction of the trend are performed. Each k
i+1
should
be less than k
i
and (k
i+1
− k
i
) ≤ 0.1. The initial start value
of 0.1 is an empirical threshold, and it was calculated exper-
imentally on multispectral data. It might be changed during
the later iteration (7). A linear scaling of K into the range [1,
0] is performed [the maximal value corresponds to maximal
HTM density in the first band with shortest wavelength, and
the minimal value corresponds to haze thickness in the longest
wavelength band at 2.2 μm (haze is not visible)]. Therefore,
a band-specific haze thickness is subtracted from the original
data. Note on the linear scaling of K: in principle, a nonlinear
scaling with wavelength could also be used, requiring scene-
derived parameters that influence the radiative transfer, e.g., the
haze particle size distribution and scattering phase functions.
This information cannot usually be retrieved. However, we
checked the potential influence of the phase function: in the
relevant spectral region for haze removal (400–900 nm), the
change of haze scattering functions with wavelength is negli-
gible. This behavior was checked with the MODTRAN [15]
haze/fog models. Summarizing, HTM
i
is calculated as follows:
HTM
i
(x, y)=DARK_PIX_SEARCH (Band
i
(x, y),w)
(5)
where the nonoverlapping window size is w =21
k
i
= SLOPE (HTM(x, y), HTM
i
(x, y)) (6)
where k
i
is the regression coefficient (slope of the fitted line)
for band i and K is the array of the coefficients (k
i
∈ K).
HTM(x, y) is the independent variable, and HTM
i
(x, y) is the
dependent variable. Note that the HTM(x, y) and HTM
i
(x, y)
values are taken only in the segmented hazy regions (haze mask
is employed)
K = LINEAR_SCALE (K, [1, 0]) for λ =[λ
blue
, 2.2 μm].
(7)
![Fig. 7. Landsat 8 OLI band 3 (0.525–0.600 μm) dehazing (subset). (a) Original band 3. (b) Dehazed band 3. (c) Histogram of the original band 3 (dotted), dehazed (dashed), and the compensated [aerosol thickness compensation (10)] band (solid line). Note that negative values can appear (regions of shadows, vegetation, deep water, and dark objects) since the HTM is calculated using the extrapolated band and the HTM has higher thickness values than the original bands.](/figures/fig-7-landsat-8-oli-band-3-0-525-0-600-mm-dehazing-subset-a-2jgzsqwn.webp)
