IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 2, FEBRUARY 2007 463
The Regularized Iteratively Reweighted
MAD Method for Change Detection
in Multi- and Hyperspectral Data
Allan Aasbjerg Nielsen
Abstract—This paper describes new extensions to the previously
published multivariate alteration detection (MAD) method for
change detection in bi-temporal, multi- and hypervariate data
such as remote sensing imagery. Much like boosting methods
often applied in data mining work, the iteratively reweighted
(IR) MAD method in a series of iterations places increasing
focus on “difficult” observations, here observations whose change
status over time is uncertain. The MAD method is based on the
established technique of canonical correlation analysis: for the
multivariate data acquired at two points in time and covering
the same geographical region, we calculate the canonical variates
and subtract them from each other. These orthogonal differences
contain maximum information on joint change in all variables
(spectral bands). The change detected in this fashion is invariant to
separate linear (affine) transformations in the originally measured
variables at the two points in time, such as 1) changes in gain and
offset in the measuring device used to acquire the data, 2) data
normalization or calibration schemes that are linear (affine) in
the gray values of the original variables, or 3) orthogonal or
other affine transformations, such as principal component (PC)
or maximum autocorrelation factor (MAF) transformations. The
IR-MAD method first calculates ordinary canonical and original
MAD variates. In the following iterations we apply different
weights to the observations, large weights being assigned to
observations that show little change, i.e., for which the sum of
squared, standardized MAD variates is small, and small weights
being assigned to observations for which the sum is large. Like the
original MAD method, the iterative extension is invariant to linear
(affine) transformations of the original variables. To stabilize
solutions to the (IR-)MAD problem, some form of regularization
may be needed. This is especially useful for work on hyperspectral
data. This paper describes ordinary two-set canonical correlation
analysis, the MAD transformation, the iterative extension, and
three regularization schemes. A simple case with real Landsat
Thematic Mapper (TM) data at one point in time and (partly)
constructed data at the other point in time that demonstrates the
superiority of the iterative scheme over the original MAD method
is shown. Also, examples with SPOT High Resolution Visible data
from an agricultural region in Kenya, and hyperspectral airborne
HyMap data from a small rural area in southeastern Germany are
given. The latter case demonstrates the need for regularization.
Index Terms—Canonical correlation analysis (CCA), iteratively
reweighted multivariate alteration detection (IR-MAD), MAD
transformation, regularization or penalization, remote sensing.
Manuscript received March 18, 2005; revised July 12, 2006. This work was
supported in part by the EU funded Network of Excellence Global Monitoring
for Security and Stability, GMOSS. The associate editor coordinating the review
of this manuscript and approving it for publication was Dr. Jacques Blanc Talon.
The author is with Informatics and Mathematical Modelling, Technical Uni-
versity of Denmark, DK-2800 Kgs. Lyngby, Denmark (e-mail: aa@imm.dtu.
dk).
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/TIP.2006.888195
I. I
NTRODUCTION
T
HIS paper deals with detection of nontrivial change
in multi- and hypervariate, bi-temporal data. The term
“nontrivial” here means nonaffine change between two points
in time. This means that changes—due to, for instance, 1) an
additive shift in mean level (offset) or a multiplicative shift in
calibration of a measuring device (gain), 2) data normaliza-
tion or calibration schemes that are linear (affine) in the gray
values of the original variables, or 3) orthogonal or other affine
transformations such as principal component (PC) or max-
imum autocorrelation factor (MAF) transformations—are not
detected. This invariance is an enormous advantage over most
other multivariate change detection schemes published, see [1]
for an early survey and [2] for a more recent one. For more
recent work on temporal dynamics in remote sensing image
data including change detection, see, for example, [3]–[5].
The method described here which is called iteratively
reweighted multivariate alteration detection (IR-MAD) is a new
extension to the previously published multivariate alteration
detection (MAD) method [6]–[9] which, in turn, is based on
the established multivariate statistical technique canonical
correlation analysis (CCA), first described by Hotelling in 1936
[10]. Inspired by boosting methods often applied in data mining
work [11] and by [12], iteratively reweighted MAD in a series
of iterations places increasing focus on “difficult” observations;
in a change detection setting, “difficult” observations are the
ones whose change status over time is uncertain. This is done
by calculating a measure of no change based on the sum of
squared, standardized MAD variates in each iteration. This
measure is then used as a weighting function for the calculation
of the statistics used to calculate the MAD transformation in the
next iteration. The idea in using such a scheme is to establish
an increasingly better no-change background against which
to detect change. Other types of robustification of the change
detection method are briefly mentioned.
To prevent a change detection method from detecting unin-
teresting change due to noise or arbitrary spurious differences,
this paper also describes the application of regularization
(also known as penalization). Regularization in the form of
smoothing of the (IR) CCA/MAD solution is described in some
detail. The paper also mentions the exploitation of the affine
transformation invariance of the MAD method as a regularizing
measure, and a combination of the two types of regularization.
Regularization may be important especially when change
detection is applied to hyperspectral data.
1057-7149/$25.00 © 2006 IEEE
464 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 2, FEBRUARY 2007
Geometrical and other corrections required in order to carry
out change studies are not dealt with here. For special prob-
lems with high spatial resolution, oblique viewing data, see, e.g.,
[13] and [14] on change detection in IKONOS data and [15] on
change detection in QuickBird data.
Methods such as the ones described in this paper are well
suited for inclusion in image processing packages and in geo-
graphical information systems (GIS).
Section II introduces multivariate change detection, Sec-
tion III very briefly describes canonical correlation analysis,
and Section IV defines the MAD transformation with sub-
sections on both the suggested iterative re-weighting and
regularization schemes. Section V gives three data examples
with Landsat Thematic Mapper data (partly constructed data),
SPOT High Resolution Visible data, and hyperspectral HyMap
data. Section VI gives conclusion and dicusses future work. An
Appendix gives more detail on selected aspects of canonical
correlation analysis including regularization.
II. M
ULTIVARIATE
CHANGE
DETECTION
When we analyze changes in panchromatic image data with
additive noise taken at different points in time, it is customary
to calculate the difference between two images. The idea is, of
course, that areas which exhibit no or small changes have zero
or low absolute values and areas with large changes have large
absolute values in the difference image. If we have two mul-
tivariate images with variables at a given location written as
vectors (without loss of generality we assume that the expec-
tation values
), and
, where is the number of spectral bands, then
a simple spectral change detection transformation is the vector
of band-wise differences also known as the change vector
(1)
In general, simple differences make sense only if the data are
normalized to a common zero and scale or calibrated over time.
If our image data have (many) more than three spectral bands,
it is difficult to visualize change in all bands simultaneously.
To overcome this problem and to concentrate information on
change, linear transformations of the image data that optimize
some measure of change (also termed a design criterion) can be
considered. A linear transformation that will maximize a mea-
sure of change in the simple multispectral difference image is
one that maximizes deviations from no change, for instance the
variance
(2)
Areas in the image data with high absolute values of
are maximum change areas. A multiplication of vector by a
constant
will multiply the variance by . Therefore, we must
make a choice concerning
. A natural choice is to request that
is a unit vector, . Maximizing the variance in (2)
under the constraint amounts to finding principal components
of the simple difference images. Principal components analysis
was developed by Hotelling in 1933 [16] based on a technique
described by Pearson in 1901.
A more parameter rich measure of change that allows dif-
ferent coefficients for
and and different numbers of spec-
tral bands in the two sets,
and , respectively, , are
linear combinations
(3)
(4)
and the difference between them
. This measure in
principle also accounts for situations where the spectral bands
are not the same but cover different spectral regions, for in-
stance if one set of data comes from the Landsat MultiSpectral
Scanner (MSS) and the other set comes from the Landsat The-
matic Mapper (TM) or from the SPOT High Resolution Visible
(HRV) which may be valuable in historical change studies. In
this case, one must be more cautious when interpreting the mul-
tivariate difference as multivariate change.
To find
and , [17] uses principal components (PC) anal-
ysis on
and considered as one concatenated vector vari-
able; [18] applies PC analysis to simple difference images as de-
scribed above. This approach requires normalized or calibrated
data and results depend on the scale at which the individual vari-
ables are measured (for instance it depends on gain settings of
a measuring device). Also, it forces the two sets of variables to
have the same coefficients (with opposite signs), and it does not
allow for the case where the two sets of images have different
numbers of spectral bands.
Other change detection schemes based on simple difference
images include factor analysis and maximum autocorrelation
factor (MAF) analysis [19]–[21].
[12] deals with (iterated) PC analysis of the same variable
at the two points in time and consider the second PC as a (mar-
ginal) change detector for that variable. [12] also introduces spa-
tial measures such as inverse local variance weighting in statis-
tics calculation and Markov random field modelling of the prob-
ability of change (versus no change).
Another approach is to define a set of
and simultaneously
in the fashion described below. Again, let us maximize the vari-
ance, this time
. A multiplication of and
by a constant will multiply the variance by . Therefore, we
must make choices concerning
and , and natural choices in
this case are requesting unit variance of
and , see
Section III and the Appendix. The criterion then is maximize
with . With this
choice, we have
(5)
We request that
and are positively correlated, see
the next section on canonical correlation analysis. Therefore,
determining the difference between linear combinations with
maximum variance corresponds to determining linear combi-
nations with minimum (non-negative) correlation. Determina-
tion of linear combinations with extreme correlations brings the
theory of canonical correlation analysis to mind.
NIELSEN: REGULARIZED ITERATIVELY REWEIGHTED MAD METHOD 465
III. CANONICAL CORRELATION
ANALYSIS
Canonical correlation analysis investigates the relationship
between two groups of several variables. It finds two sets of
linear combinations of the original variables, one for each group.
The first two linear combinations are the ones with the largest
correlation. This correlation is called the first canonical correla-
tion and the two linear combinations are called the first canon-
ical variates. The second two linear combinations are the ones
with the largest correlation subject to the condition that they
are orthogonal to the first canonical variates. This correlation is
called the second canonical correlation and the two linear com-
binations are called the second canonical variates. Higher order
canonical correlations and canonical variates are defined simi-
larly.
Since we are looking for canonical variates that are as sim-
ilar as possible as measured by correlation, we request positive
canonical correlations.
If we denote the variance-covariance matrix, also known as
the dispersion (matrix), of the one set of variables
,
the dispersion of the other set of variables
, the co-
variance between them
, and the canonical correlation
, we get (see the Appendix)
(6)
(7)
or in terms of Rayleigh quotients
(8)
i.e., we find the desired projections for
by considering the
mutually orthogonal (also known as conjugate) eigenvectors
corresponding to the eigenvalues
of with respect to . Similarly, we find the
desired projections for
by considering the conjugate
eigenvectors
of with respect to
corresponding to the same eigenvalues .
This technique was first described in [10] and a treatment is
given in most textbooks on multivariate statistics (good refer-
ences are [22] and [23]).
Multiset canonical correlation analysis where we investigate
the relationship between more than two groups of several vari-
ables first introduced in [24], [25] is described and applied to
remote sensing data in [6] and [26]. Nonlinear (two- and mul-
tiset) canonical correlation analysis is dealt with in [27]–[30].
IV. MAD T
RANSFORMATION
Inspired by Sections II and III, we define the multivariate al-
teration detection (MAD) transformation as
.
.
.
(9)
where
and are the defining coefficients from a standard
canonical correlation analysis. To maximize variance in (5),
we must minimize
; therefore, we have reversed the order of
the differences between the canonical variates in (9) so MAD
variate 1 is the difference between the highest order canonical
variates, MAD variate 2 is the difference between the second
highest order canonical variates, etc.
The dispersion matrix of the MAD variates is
(10)
where
is the unit matrix and is a matrix containing
the ascendingly sorted canonical correlations on the diagonal
and zeros off the diagonal so the MAD variates are orthogonal
with variance
(11)
The MAD transformation has the very important property
that if we consider linear combinations of two sets
(of vari-
ables) and
(of variables, ) that are positively correlated
then the
th difference shows maximum variance among such
variables. The
th difference shows maximum variance
subject to the constraint that this difference is uncorrelated with
the previous
ones. In this way, we sequentially extract uncorre-
lated difference images where each new image shows maximum
difference (change) under the constraint of being uncorrelated
with the previous ones. If
, then the projection of on
the eigenvectors corresponding to the eigenvalues 0 will be in-
dependent of
. That part may be considered the extreme case
of multivariate change detection.
As opposed to the principal components of simple differ-
ences, the MAD variates are invariant to affine transformations
(including linear scaling), which means that they are sensitive
to neither, for example, changes in gain settings and offset in a
measuring device, nor to linear (affine) radiometric and atmo-
spheric correction schemes.
Because the MAD variates are linear combinations of the
measured variables, they will have approximately a Gaussian
distribution because of the Central Limit Theorem, see, e.g.,
[31]. In addition, if there is no change at pixel
, then the th
MAD value, MAD
, has mean 0. Assuming also independence
of the orthogonal MAD variates we may expect that the sum
of the squared MAD variates for pixel
after standardization
to unit variance approximately follows a
distribution with
degrees of freedom, i.e., approximately
(12)
The standardization should ideally be done by means of the stan-
dard deviation of the no-change observations. This standard de-
viation can be estimated by means of expectation–maximiza-
tion (EM) based methods for determining thresholds for differ-
entiating between change and no change in the difference im-
ages, and for estimating the variance-covariance structure of the
no-change observations [32]–[37]. Below, we use the standard
deviation for all observations for simplicity. Provided that the
proportion of changed pixels is small, this will have minimal
effect.
Equation (12) can be used to assign labels “change” or
“no-change” to each observation by means of percentiles in the
466 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 16, NO. 2, FEBRUARY 2007
distribution. We may choose to assign the label “change”
to observations with
values greater than, say, the 99% per-
centile and similarly the label “no-change” to observations with
values smaller than, say, the 1% percentile. Since the MAD
transformation is invariant to linear (affine) transformations
these no-change observations are suitable for carrying out an
automated normalization between the two points in time. This
is described in detail in [38].
The spatial aspect introduced in change detection in [20] can
be applied here also by postprocessing the MAD variates with
the (change strength weighted) MAF transformation, see [8].
The spatial aspect is dealt with elegantly in a Markov random
field setting in [33] and [34].
The main feature of the MAD method is the transformation
from a space where the originally measured variables are or-
dered by wavelength into a feature space where the transformed,
orthogonal variables are ordered by similarity (as measured
by linear correlation). This latter ordering is considered to
be more relevant for change detection purposes. Differences
between corresponding pairs of variables in this latter space,
i.e., differences between the canonical variates, give us the or-
thogonal MAD variates which can be considered as generalized
difference images well suited for change detection.
A. Iteratively Reweighted MAD, IR-MAD
Inspired by boosting methods often used in data mining [11]
and by [12], the idea in iteratively reweighted (IR) MAD is
simply in an iterative scheme to put high weights on observa-
tions that exhibit little change over time. This can be done in
several ways. We start with the original MAD transformation,
i.e., we assign the same weight (= 1) to all pixels. A natural
choice is to weight pixel
in the next iteration by , which
is a measure of no change, namely the probability of finding a
greater value of the
value in (12)
(13)
This weight enters into the calculation of mean values, variances
and covariances (
is the number of pixels)
(14)
for the mean value of
, and
(15)
for the covariance between
and (if , we get the
variance of
).
Iterations are performed until the largest absolute change in
the canonical correlations becomes smaller than some preset
small value
, e.g., . This weighting scheme maps the
weights applied to the interval
and avoids very high
weights. Unlike boosting methods, weights from the early
iterations are not used in this scheme, only the weights from the
final iteration are used so the “committee” and voting scheme
often involved in boosting are not used here.
Of course, other reweighting schemes for example leading to
robust estimation [39], [40] of the variance–covariance structure
of the data could be used. A limited number of tests indicates
that the iterated scheme suggested in this section performs better
than robust estimation; see also [41].
B. Regularized IR-MAD
If we have many (correlated) variables, the solutions to the
coupled generalized eigenvalue problems in (6) and (7) may be-
come unstable due to (near) singular variance–covariance ma-
trices causing small changes in the data to lead to dramatically
different solutions. A possible solution to such (near) singularity
problems in hyperspectral data change detection may be regu-
larization (also known as penalization) where, inspired by ridge
regression described in [42], we add
to and
to in (6) and (7). are (small) non-negative numbers that
can be chosen subjectively or estimated from the data, see the
Appendix. This was first described in the CCA context with
as the identity matrix in [43]. [44] penalizes high local variation
using a second order derivative-type
. To obtain a continuous
and differentiable second order derivative, [45] in a functional
setting suggests a fourth order derivative-type
.
Since
and here are the same type of data, we choose
and . We choose in
to penalize curvature of the elements in and considered as
functions of wavelength. Choosing the usual discrete approxi-
mation to the second order derivative we penalize
and
where is with (typically )
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
(16)
leading to
.
.
.
.
.
.
.
.
.
(17)
(which is penta-diagonal and rank
); see also [36],
[46]–[50], and the Appendix.
Alternatively, the invariance of the MAD variates to linear
(affine) transformations of the original variables can be ex-
ploited. Possible (near) singularities may also be remedied
by means of principal component analysis (PCA), maximum
autocorrelation factor (MAF), projection pursuit (PP) analysis
or other dimensionality reducing projections applied to the
variables at the two points in time separately before doing
canonical correlation and MAD analysis. This approach is used
in [49] and [51].
NIELSEN: REGULARIZED ITERATIVELY REWEIGHTED MAD METHOD 467
Fig. 1. Landsat TM data from June 6, 1986, covering a forested region in
Northern Sweden, spectral bands 1, 2, 3, 4, 5, and 7 row-wise.
If regularization is needed or desired, one may use either the
former, the latter or a combined scheme. The ordering of the
projection variates in the dimensionality reducing regulariza-
tion scheme is by some projection index (such as variance, au-
tocorrelation, deviation from normality or other) rather than by
wavelength. This ordering makes penalizing for example curva-
ture un-natural. So the two regularization schemes do not readily
combine.
Inspired by [52], we apply the dimensionality or feature re-
duction scheme above to adjacent, nonoverlapping groups of
spectral bands. For example, we may replace bands 1, 2, and
3 with one projection, 4, 5, and 6 with another, etc. In this way,
we reduce the dimensionality of the data (in the example by
a factor of three) while retaining the main spectral features of
the original data and their order. This preservation of order fa-
cilitates the application of further regularization by penalizing
for example curvature as described in the former regularization
scheme above.
If we use this combined regularization scheme, the general
transformation invariance may be lost depending on the choice
of dimensionality reduction scheme. If we choose the MAF
transformation, we retain the invariance to any transformations
that are linear (or affine) in the individual original variables.
V. E
XAMPLES
The examples include a partly constructed case with Landsat
Thematic Mapper (TM) data from a forested region in Northern
Fig. 2. Constructed Landsat TM data covering a forested region in Northern
Sweden: the 512
2
128 leftmost part of the image consists of data from June
27, 1988, padded into the Landsat TM data from June 6, 1986, spectral bands
1, 2, 3, 4, 5, and 7 row-wise.
Sweden, a case with SPOT High Resolution Visible (HRV) data
covering an agricultural region in tropical Kenya, as well as a
case with hyperspectral HyMap data from a small rural area in
southeastern Germany.
All images in this paper are stretched linearly between mean
and
three standard deviations unless otherwise stated.
A. Partly Constructed Landsat TM Data, Northern Sweden
This case compares results from the original MAD method
with those from the new iterated scheme where data at one point
in time are constructed so that we know where change did not
occur. Data at the first point in time are 512
512 25 m 25
m Landsat Thematic Mapper (TM) spectral bands 1, 2, 3, 4, 5,
and 7 from June 6, 1986, covering a forested region in Northern
Sweden. Data at the second point in time are 512
128 Landsat
TM (same bands) from June 27, 1988, covering the same re-
gion padded into the leftmost part of the 1986 image. Hence,
by construction there is no change in the rightmost 512
384
part of the image. In this simple case, band-wise differences
will give the desired zero change but as mentioned in general
simple differences make sense only if the data are normalized
to a common zero and scale or calibrated over time. Fig. 1 shows
the measured 1986 data and Fig. 2 shows the partly constructed
1988 data.
Fig. 3 shows the original MAD variates and Fig. 4 shows the
IR-MAD variates after 30 iterations (and convergence). Visual
