4. Multispectral change detection¶
To quote a well-known, if now rather dated, review article on change detection (Singh 1989,
The basic premise in using remote sensing data for change detection is that changes in land cover must result in changes in radiance values ... [which] must be large with respect to radiance changes from other factors.
When comparing multispectral images of a given scene taken at different times, it is therefore desirable to correct the pixel intensities as much as possible for uninteresting differences such as those due to solar illumination, atmospheric conditions or sensor calibration. In the case of SAR imagery, solar illumination of course plays no role, but the other considerations are similarly important. If comparison is on a pixel-by-pixel basis, then the images must also be co-registered to high accuracy in order to avoid spurious signals resulting from misaligned pixels.
Algebraic methods¶
A simple way to detect changes in two corrected and coregistered multispectral images, represented by $N$-dimensional random vectors $X$ and $Y$, is to subtract them from each other component by component and then examine the $N$ difference images
$$ Z_i = X_i-Y_i,\quad i=1\dots N. $$Small intensity differences indicate no change, large positive or negative values indicate change, and decision thresholds can be set to define significance. The thresholds are usually expressed in terms of standard deviations from the mean difference value, which is taken to correspond to no change. If the images are uncorrelated, the variances of the difference images are simply
$$ {\rm var}(Z_i) = {\rm var}(X_i) + {\rm var}(Y_i),\quad i=1\dots N, $$or about twice as noisy as the individual image bands.
For SAR imagery, the situation is fundamentally different. The variance of the difference of two uncorrelated $m$-look intensity images is
$$ {\rm var}(X-Y) = {\langle X\rangle^2 + \langle Y\rangle^2 \over m}. $$Simple thresholding of the difference image will yield larger errors for a given change in a bright area (large mean intensity) than in a darker area (small mean intensity). Image ratios are a much better choice for detection of changes in multi-look SAR intensity images as we will see later.
Multivariate Alteration Detection (MAD)¶
We make a linear combination of the intensities for all $N$ bands in the first image thus creating a scalar image characterized by the random variable
$$ U =a^\top X. $$The vector of coefficients $a$ is as yet unspecified.
We do the same for the second image,forming the linear combination $$ V=b^\top Y, $$
and then look at the scalar difference image $U-V$.
We will choose $a$ and $b$ so as to make the images $U$ and $V$ as similar as possible. To do this , we maximize the correlation $\rho$ between the random variables $U$ and $V$:
$$ \rho = {{\rm cov}(U,V)\over \sqrt{{\rm var}(U)}\sqrt{{\rm var}(V)}}.\quad (1) $$Arbitrary multiples of $U$ and $V$ would clearly have the same correlation, so a constraint must be chosen. A convenient one is
$$ {\rm var}(U)={\rm var}V)=1. \quad (2) $$Note that, under this constraint, the variance of the difference image is
$$ {\rm var}(U-V)= {\rm var}(U)+{\rm var}(V)-2{\rm cov}(U,V) = 2(1-\rho).\quad (3) $$The statistical procedure which makes $U$ and $V$ as similar as possible is called canonical correlation analysis.
The bitemporal, multispectral image may be represented by the combined random vector $ \pmatrix{X\cr Y}. $ This random vector has a $2N\times 2N$ covariance matrix which can be written in block form:
$$ \Sigma=\pmatrix{\Sigma_{XX}&\Sigma_{XY}\cr\Sigma_{XY}^\top&\Sigma_{YY}}. $$The variances and covariances in Equations (1) and (2) are then
$$ {\rm var}(U) = a^\top \Sigma_{XX}a,\quad {\rm var}(V)=b^\top \Sigma_{YY}b,\quad {\rm cov}(U,V) = a^\top\Sigma_{XY} b. $$We introduce the Lagrange multipliers ${\nu/2}$ and $\mu/2$ for each of the two constraints, then the problem becomes one of maximizing the unconstrained Lagrange function
$$ L = a^\top \Sigma_{XY} b-{\nu\over 2}( a^\top \Sigma_{XX} a-1)-{\mu\over 2}(b^\top \Sigma_{YY}\ b-1). $$Setting derivatives with respect to $a$ and $b$ equal to zero gives
$$ {\partial L\over \partial a}=\Sigma_{XY} b-\nu\Sigma_{XX} a = 0 $$$$ {\partial L\over \partial b}=\Sigma_{XY}^\top a-\mu\Sigma_{YY} b = 0. $$Multiplying the first of the above equations from the left with the vector $ a^\top$, the second with $ b^\top$, and using the constraints leads immediately to
$$ \nu = \mu = a^\top\Sigma_{XY} b = \rho. $$Therefore we can write
$$ \Sigma_{XY} b - \rho\Sigma_{XX} a = 0 $$$$ \Sigma_{XY}^\top a - \rho\Sigma_{YY} b = 0. $$Next, multiply the first Equations above by $\rho$ and the second from the left by $\Sigma_{YY}^{-1}$. This gives
$$ \rho\Sigma_{XY} b = \rho^2\Sigma_{XX} a $$and
$$ \Sigma_{YY}^{-1}\Sigma_{XY}^\top a = \rho b. $$Combining these two Equations, we obtain the following equation for the transformation coefficient $a$,
$$ \Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{XY}^\top a = \rho^2\Sigma_{XX} a $$Similarly we get for $b$,
$$ \Sigma_{XY}^\top\Sigma_{XX}^{-1}\Sigma_{XY} b = \rho^2\Sigma_{YY} b. $$which can be solved for $a$ and $b$.
Solution of the eigenvalue problems generates new multispectral images
$$ U=(U_1\dots U_N)^\top $$and
$$ V=(V_1\dots V_N)^\top, $$the components of which are called the canonical variates (CVs). The CVs are ordered by similarity (correlation) rather than, as in the original images, by wavelength.
The canonical correlations
$$ \rho_i={\rm corr}(U_i,V_i),\quad i=1\dots N, $$are the square roots of the eigenvalues of the coupled eigenvalue problem:
$$ \rho_i = \sqrt{\lambda_i}. $$The pair $(U_1,V_1)$ is maximally correlated, the pair $(U_2,V_2)$ is maximally correlated subject to being orthogonal to (uncorrelated with) both $U_1$ and $V_1$, and so on.
Taking paired differences (in reverse order) generates a sequence of transformed difference images
$$ M_i = U_{N-i+1}-V_{N-i+1},\quad i=1\dots N, $$referred to as the multivariate alteration detection (MAD) variates.
From Equation (3),
$$ {\rm var}(M_i) = \sigma_{M_i}^2 = 2(1-\rho_{N-i+1}) = 2(1-\sqrt{\lambda_{N-i+1}}).\quad (4) $$
run /home/imad -h
%matplotlib inline
cd /home/imagery
ls -l
The "ordinary MAD transformation", i.e., just one iteration:
run /home/imad -i 1 /home/imagery/20010626 /home/imagery/20010829
run /home/dispms -p [4,5,6] -e 3 -f /home/imagery/MAD(20010626-20010829)
Iterated MAD¶
Imagine two images of the same scene, acquired at different times under similar conditions, but for which no ground reflectance changes have occurred whatsoever.
Then the only differences between them will be due to random effects like instrument noise and atmospheric fluctuation. In such a case we would expect that the histogram of any difference component that we generate will be very nearly Gaussian.
In particular, the MAD variates, being uncorrelated, should follow a multivariate, zero mean normal distribution with diagonal covariance matrix.
Change observations would deviate more or less strongly from such a distribution. We might therefore expect an improvement in the sensitivity of the MAD transformation if we can establish an increasingly better background of no change against which to detect change.
This can again be done in an iteration scheme in which, when calculating the means and covariance matrices for the next iteration of the MAD transformation, observations are weighted by the probability of no change determined on the preceding iteration.
Let the random variable $Z$ represent the sum of the squares of the standardized MAD variates,
$$ Z = \sum_{i=1}^N\left({M_i\over \sigma_{M_i}}\right)^2, $$where $\sigma_{M_i}$ is given by Equation (4). Then, since the no-change observations are expected to be normally distributed and uncorrelated, the random variable $Z$ corresponding to such observations should be chi-square distributed with $N$ degrees of freedom. For each iteration, the observations can thus be given weights determined by the chi-square distribution, namely
$$ {\rm Pr}(\hbox{no change}) = 1-P_{\chi^2;N}(z). $$Here ${\rm Pr}(\hbox{no change})$ is the probability that a sample $z$ drawn from the chi-square distribution could be that large or larger.
Iteration of the MAD transformation continues until there is no significant change in the canonical correlations $\rho_i$, $i=1\dots N$.
run /home/imad /home/imagery/20010626 /home/imagery/20010829
run /home/dispms -e 3 -p [6,5,4] -E 3 -P [7,7,7] \
-f /home/imagery/MAD(20010626-20010829) \
-F /home/imagery/MAD(20010626-20010829)
Clustering changes¶
The image shown on the left is an RGB composite of the three MAD components corresponding to the three most similar canonical variates. The light blue-green (featureless) regions correspond to built-up areas and to forest canopy, neither of which have changed significantly between acquisitions. The significant changes, indicated by the various brightly colored areas, are quite heterogeneous and not easy to interpret, but from their form they are clearly associated with cultivated fields. The main crops in this area are sugar beets, corn and cereal grain. The former two are still maturing in August, whereas the grain fields have been harvested. So we might postulate two main classes of change: maturing crops and harvested crops. Allowing for a no-change class and a "catch-all" class for other kinds of change, we can try to classify the MAD image with the Gaussian mixture clustering algorithm by assuming four classes in all.
run /home/em -p [4,5,6] -K 4 /home/imagery/MAD(20010626-20010829)
run /home/dispms -e 1 -c [1,2,3,4] -f /home/imagery/MAD(20010626-20010829)_em
Harvested crops
run /home/dispms -e 3 -p [4,4,4] -d [0,0,400,400] -f /home/imagery/20010829 \
-C [2] -E 1 -D [0,0,400,400] -F /home/imagery/MAD(20010626-20010829)_em
Automatic radiometric normalizaton¶
Ground reflectance determination from satellite imagery requires, among other things, an atmospheric correction algorithm and the associated atmospheric properties at the time of image acquisition. For many historical satellite scenes such data are not available and even for planned acquisitions they may be difficult to obtain. A relative normalization based on the radiometric information intrinsic to the images themselves is an alternative whenever absolute surface reflectances are not required.
In performing relative radiometric normalization, one usually makes the assumption that the relationship between the at-sensor radiances recorded at two different times from regions of constant reflectance can be approximated by linear functions. The critical aspect is the determination of suitable time-invariant features upon which to base the normalization.
Scatterplot matching¶
With reference to the above Figure, one can derive the following linear transformation, which relates a point $(x,y)$ in the target scatterplot to its transform $(\tilde x,\tilde y)$ in the reference scatterplot:
$$ \tilde x = X_R + (x - X_T){L_R\over L_T}{a_T\over a_R} $$$$ \tilde y = Y_R - (Y_T - y){L_R \over L_T}. $$
IR-MAD Normalization¶
We can use the chi-square image in the IR-MAD transformation to identify invariant pixels. These can be used for radiometric normalization of the two images because canonical correlation anaysis is invariant under linear transformations.
run /home/radcal -h
run /home/radcal /home/imagery/MAD(20010626-20010829)
We compare the August image before and after normalization on absolute scale 0-255 with the July image:
run /home/dispms -e 1 -p [4,5,6] -f /home/imagery/20010829 \
-E 1 -P [4,5,6] -F /home/imagery/20010829_norm
run /home/dispms -e 1 -p [4,5,6] -f /home/imagery/20010626 \
-E 1 -P [4,5,6] -F /home/imagery/20010829_norm
