5. Polarimetric SAR Change Detection¶
The /home/imagery directory contains the polametric SAR data and is shared with the host. In the present case there are 3 Radarsat-2 quadpol images in SLC (single-look complex) format. Acquistion times range from April 26, 2010 (20100426) to July 7, 2010 (20100707):
!ls -l /home/imagery | grep "_SLC$"
The images are level one SLC (single look complex format). For example, below are the contents of the image directory corresponding to acquistion date 20100426 (April 26, 2010). The four polarization combinations HH, HV,VH and VV are are stored as complex numbers in GeoTiff format:
ls -l /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC
The subdirectory polsarpro/T3 contains the polarimetric coherency matrix elements generated from the polarization combinations. This can be done with the PolSARpro provided as freeware by the European Space Agency (ESA) or with the Sentinel-1 Toolbox.The PolSARpro-type image files are in ENVI format:
ls -l /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC/polsarpro/T3
We'll return to these images later. But first we'll review the statistical properties if polarimetric SAR data:
Statistical Properties of polSAR Images¶
Vector and matrix representations¶
A fully polarimetric SAR measures a $2\times 2$ scattering matrix $S$ at each resolution cell on the ground. The scattering matrix relates the incident and the backscattered electric fields $E^i$ and $E^b$ according to
$$ \pmatrix{E_h^b \cr E_v^b} =\pmatrix{S_{hh} & S_{hv}\cr S_{vh} & S_{vv}}\pmatrix{E_h^i \cr E_v^i}. $$Here $E_h^{i(b)}$ and $E_v^{i(b)}$ denote the horizontal and vertical components of the incident (backscattered) oscillating electric fields directly at the target. These can be deduced from the transmitted and received radar signals via the so-called far field approximations. If both horizontally and vertically polarized radar pulses are emitted and discriminated then they determine, from the above Equation, the four complex scattering matrix elements. The per-pixel polarimetric information in the scattering matrix $S$, under the assumption of reciprocity ($S_{hv} = S_{vh}$), can then be expressed as a three-component complex vector
$$ s = \pmatrix{S_{hh}\cr \sqrt{2}S_{hv}\cr S_{vv}}, $$where the $\sqrt{2}$ ensures that the total intensity (received signal power) is consistent. It is essentially these vectors which are provided in the SLC level one data discussed above. The total intensity is referred to as the span and is the complex inner product of the vector $s$,
$$ {\rm span} = s^\top s = |S_{hh}|^2 + 2|S_{hv}|^2 + |S_{vv}|^2. $$This is a real number and the corresponding gray-scale image is called the span image. The observation vector $s$ can be shown to be a realization of a complex multivariate normal random variable. An equivalent and often preferred representation is in terms of the coherency vector
$$ t = {1\over\sqrt{2}}\pmatrix{S_{hh} + S_{vv}\cr S_{hh} - S_{vv} \cr 2S_{hv}}. $$This so-called Pauli decomposition allows for a more physical interpretation of the backscattered signal:
$t_1 =$ single bounce (flat surface)
$t_2 =$ double bounce (buildings)
$t_3 =$ volume scattering (forest canopy)
The span images are the same in both cases, i.e.,
$$ s^\top s = t^\top t $$The polarimetric signals can also be represented by taking the complex outer product of $s$ with itself:
$$ C = s s^\top = \pmatrix{ |S_{hh}|^2 & \sqrt{2}S_{hh}S_{hv}^* & S_{hh}S_{vv}^* \cr \sqrt{2}S_{hv}S_{hh}^* & 2|S_{hv}|^2 & \sqrt{2}S_{hv}S_{vv}^* \cr S_{vv}S_{hh}^* & \sqrt{2}S_{vv}S_{hv}^* & |S_{vv}|^2 }. $$The diagonal elements of $C$ are real numbers, with span $= {\rm tr}(C)$, and the off-diagonal elements are complex. This matrix representation contains all of the information in the polarized signal.
Alternatively one can work with
$$ T = t t^\top $$Multi-looking¶
The matrix $C$ (or $T$) can be averaged over the number of looks (number of adjacent cells used to average out the effect of speckle) to give an estimate of the covariance matrix of each multi-look pixel:
$$ \bar{C} ={1\over m}\sum_{\nu=1}^m s(\nu) s(\nu)^\top = \langle s s^\top \rangle = \pmatrix{ \langle |S_{hh}|^2\rangle & \langle\sqrt{2}S_{hh}S_{hv}^*\rangle & \langle S_{hh}S_{vv}^*\rangle \cr \langle\sqrt{2} S_{hv}S_{hh}^*\rangle & \langle 2|S_{hv}|^2\rangle & \langle\sqrt{2}S_{hv}S_{vv}^*\rangle \cr \langle S_{vv}S_{hh}^*\rangle & \langle\sqrt{2}S_{vv}S_{hv}^*\rangle & \langle |S_{vv}|^2\rangle }, $$where $m$ is the number of looks. This matrix (or rather the equivalent coherency matrix $ \bar T = \langle t t^\top \rangle$) is generated by the Sentinel toolbox and stored in the subdirectory polsarpto/T3 which we listed in a previous cell above. Rewriting the first of the above equalities,
$$ m\bar{C} = \sum_{\nu=1}^m s(\nu) s(\nu)^\top =: x. $$This quantity $x$ is a realization of a complex sample matrix which is known to have a complex Wishart distribution with $N\times N$ covariance matrix $\Sigma$ (here $N=3$) and $m$ degrees of freedom:
$$ p_{W_c}( x) ={|x|^{(m-N)}\exp(-{\rm tr}(\Sigma^{-1} x)) \over \pi^{N(N-1)/2}|\Sigma|^{m}\prod_{i=1}^N\Gamma(m+1-i)},\quad m \ge N, $$provided that the vectors $s(\nu)$ are independent and drawn from a complex multivariate normal distribution.
Equivalent number of looks¶
However when multi-look averaging takes place, the observation vectors $s(\nu)$ are not completely independent and will generally be correlated somewhat. In order to account for this, the complex Wishart distribution is often parameterized with ENL (rather than $m$) degrees of freedom, where ENL is the so-called equivalent number of looks. This quantity can be estimated from the image itself. For exampe, if
$$ z = \langle |S_{hh}|^2\rangle $$represents a pixel intensity in the first band of the covariance image, ideally we have
$$ {\rm var}(z) = {\langle z\rangle^2 \over m} $$where $m$ is the number of looks. So we use
$$ {\rm ENL} \approx {\langle z\rangle^2\over {\rm var}(z)}. $$Dual and single polarimetric imagery¶
The scattering vector given above corresponds to so-called full, or quad polarimetric SAR. Satellite-based SAR sensors often operate in reduced, power-saving polarization modes, emitting only one polarization and receiving two (dual polarization) or one (single polarization). The look-averaged covariance matrices are reduced in dimension correspondingly. For instance for dual polarization with horizontal transmission and horizontal and vertical reception,
$$ \bar{C} = \pmatrix{ \langle |S_{hh}|^2\rangle & \langle S_{hh}S_{hv}^*\rangle \cr \langle S_{hv}S_{hh}^*\rangle & \langle |S_{hv}|^2\rangle }, $$and, for single polarization and horizontal transmission/reception, we get simply the intensity image
$$ \bar{I} = \langle |S_{hh}|^2\rangle \quad {\rm or} \quad \bar{I} = \langle |S_{vv}|^2\rangle. $$In the the dual pol case, the observations are complex Wishart distributed with $N=2$, in the single pol case they are gamma distributed.
Hypothesis testing¶
Some terminology:
A statistical hypothesis is a conjecture about the distributions of one or more random variables. It might, for instance, be an assertion about the mean of a distribution, or about the equivalence of the variances of two different distributions. One distinguishes between simple hypotheses, for which the distributions are completely specified, for example:
The mean of a normal distribution with variance $\sigma^2$ is $\mu=0$,
and composite hypotheses, for which this is not the case, e.g.:
The mean is $\mu\ge 0$.
In order to test such assertions on the basis of samples of the distributions involved, it is also necessary to formulate alternative hypotheses. To distinguish these from the original assertions, the latter are traditionally called null hypotheses.
Thus we might be interested in testing the simple null hypothesis $\mu = 0$ against the composite alternative hypothesis $\mu\ne 0$.
An appropriate sample function for deciding whether or not to reject the null hypothesis in favor of its alternative is referred to as a test statistic, often denoted by the symbol $Q$.
An appropriate test procedure will partition the possible realizations of the test statistic into two subsets: an acceptance region for the null hypothesis and a rejection region. The latter is customarily referred to as the critical region.
Referring to the null hypothesis as $H_0$, there are two kinds of errors which can arise from any test procedure:
$H_0$ may be rejected when in fact it is true. This is called an error of the first kind and the probability that it will occur is denoted $\alpha$.
$H_0$ may be accepted when in fact it is false, which is called an error of the second kind with probability of occurrence $\beta$.
The probability of obtaining a value of the test statistic within the critical region when $H_0$ is true is thus $\alpha$.
The probability $\alpha$ is also referred to as the level of significance of the test.
It is generally the case that the lower the value of $\alpha$, the higher is the probability $\beta$ of making a second kind error. Traditionally, significance levels of 0.01 or 0.05 are used.
The Likelihood Ratio Test¶
Consider a vector random sample $z(\nu)$, $\nu=1\dots m$, from a multivariate population whose density function is $p(z\mid \theta)$. The likelihood function for the sample is
$$ L(\theta) = \prod_{\nu=1}^m p(z(\nu)\mid \theta). $$The parameter set $\hat\theta$ which maximizes the likelihood function, i.e., which gives the largest value for all of the realizations, is called the {\it maximum likelihood estimate} of $\theta$.
For example, for normally distributed random vectors $Z$, maximum likelihood parameter estimators for $\theta = \{\mu,\Sigma\}$ turn out to be
$$ \bar{Z} = {1\over m}\sum_{\nu=1}^m Z(\nu) $$$$ S = {1\over m}\sum_{\nu=1}^m( Z(\nu) - \bar Z ) (Z(\nu) - \bar Z)^\top. $$Let $\omega\in R^q$ be space of all possible values of the parameter set $\theta$ and $\omega_0\in R^r$, $r\le q$, be a subset of that space. The likelihood ratio test (LRT) for the null hypothesis $\theta\in\omega_0$ against the alternative $\theta\in\omega-\omega_0$ has the critical region
$$ Q = {\max_{\theta\in \omega_0} L(\theta) \over \max_{\theta\in\omega} L(\theta)} \le k. $$This definition simply reflects the fact that, if $H_0$ is true, the maximum likelihood for $\theta$ when restricted to $\omega_0$ should be close to the maximum likelihood for $\theta$ without that restriction. Therefore, if the likelihood ratio is small, (less than or equal to some small value $k$), then $H_0$ should be rejected.
Often the distribution of $Q$ is unknown, so the critical region may be hard to define. However it can be shown that:
For each $\theta\in\omega_0$ (under the null hypothesis), $-2\log Q$ has an asymptotic chi-square distribution with $q-r$ degrees of freedom as $m\to\infty$.
We will see an example of this below.
Change Detection¶
Bitemporal data¶
The following is discussion is based on Conradsen et al (2004).
As we have seen, we can represent a pixel vector in an $m$ look-averaged polarimetric SAR image in covariance matrix format by $\bar C$, where
$$ m\bar C = x = \sum_{\nu=1}^m s(\nu) s(\nu)^\top $$is a realization of a random matrix $X$ with a complex Wishart distribution. This distribution is described by the $3\times 3$ complex covariance matrix $\Sigma$.
In order to derive a change detection procedure for two polarimetric SAR images, we formulate a statistical test. For each pixel in the two images, define the null (or no-change) simple hypothesis as
$$ H_0:\quad \Sigma_1 = \Sigma_2 = \Sigma, $$and the alternative composite hypothesis as
$$ H_1:\quad \Sigma_1 \ne \Sigma_2. $$Under $H_0$ the maximum likelihood for $\Sigma$ can be shown to be given by
$$ L(\hat\Sigma) = { |x_1|^{m-3}|x_2|^{m-3}\exp(-2m\cdot{\rm tr}(I)) \over \left({1\over 2m}\right)^{3\cdot 2m}| x_1+ x_2|^{2m}\Gamma_3(m)^2 }, $$where $I$ is the $3\times 3$ identity matrix and ${\rm tr}(I)=3$. Under $H_1$ the maximum likelihood for $\Sigma_1$ and $\Sigma_2$ is
$$ L(\hat\Sigma_1,\hat\Sigma_2) = { |x_1|^{m-3}|x_2|^{m-3}\exp(-2m\cdot{\rm tr}(I)) \over \left({1\over m}\right)^{3m}\left({1\over m}\right)^{3m} |x_1|^m |x_2|^m\Gamma_3(m)^2 } $$Then the likelihood ratio test has the critical region for rejection of the no-change hypothesis
$$ Q = {L(\hat\Sigma) \over L(\hat\Sigma_1,\hat\Sigma_2) } = 2^{6m}{ |x_1|^m |x_2|^m \over |x_1 + x_2|^{2m} } \le k. $$For a large number of looks $m$, we get the following asymptotic distribution for the test statistic: As $m\to\infty$, the quantity $$ -2\log Q = -2m( 6\log 2 + \log|x_1| + \log|x_2| -2\log|x_1+x_2| ) $$ is a realization of a chi-square random variable with 9 degrees of freedom:
$$ {\rm Pr}(-2\log Q\le z) \simeq P_{\chi^2;9}(z) = P_{\chi^2;N^2}(z). $$This follows from the fact that the number of parameters required to specify a $3\times 3$ complex covariance matrix is 9: three for the real diagonal elements and $2\times 3$ for the three complex elements above the diagonal. Thus the parameter space $\omega$ for $(\Sigma_1,\Sigma_2)$ has dimension $q=2\times 9$, and the subspace $\omega_0$ for the simple null hypothesis has dimension $r=9$, so $q-r=9$.
For smaller $m$, one can derive (Conradsen et al (2004)) the following approximation for the statistical distribution of the test statistic $Q$:
$$ {\rm prob}(-2\rho\log Q\le z) \simeq P_{\chi^2;N^2}(z) + \omega_2\left[ P_{\chi^2;N^2+4}(z) - P_{\chi^2;N^2}(z) \right], $$where $P_{\chi^2;m}(z)$ is the chi square distribution wth m degrees of freedom, $$ \rho = 1 - {2N^2-1\over 6N}\cdot{3\over 2 m}, $$
and
$$ \omega_2 = - {N^2\over 4}\cdot\left(1-{1\over \rho}\right)^2 + {N^2(N^2-1)\over 24 \rho^2}\cdot{7\over 4m^2}. $$Note that, as $m\to\infty$, $\rho\to1$ and $\omega_2\to 0$.
In practice we choose a significance level $\alpha$, e.g., $\alpha = 0.01$, and decision threshold $z$ such that
$$ {\rm prob}(-2\rho\log Q\le z) = 1-\alpha $$and interpret all pixels with larger values of $-2\rho\log Q$ as change.
Multitemporal data¶
The preceding discussion generalizes in a straightforward way to a time series of $k$ images (Conradsen et al (2015)). In the equation for the test statistic $Q$ above, the numerator consists of a product $k$ determinants $|x_1|\cdot|x_2|\cdot\dots|x_k|$ and the denominator similarly the determinant of the sum of the $k$ observations $|x_1+x_2+\dots x_k|$. The multitemporal test is referred to as the omnibus test and will in general be more powerful (have a higher detection probability) for the same significance level than a bitemporal test just involving the first and last images.
A Processing Chain¶
The change detection method implies the following processing sequence in order to generate a change map from a time series of polarimetric SAR images provided at the single look complex (SLC) processing level:
First of all the multi-look polarimetric SAR images in covariance or coherency matrix format are generated from the SLC data with PolSARpro (or, if prefered, from the Sentinel 1 Toolbox also available from the European Space agency). Presently this must be done outside of the Docker container (and IPython) since the ESA software is only available in the form of a graphics interface. The coherenecy matrix has the same eigenvalues and hence the same determinant as the covariance matrix, so that the hypothesis test described above can be applied unchanged to either format. The rest of the processing takes place in the IPython notebook.
The matrix images are imported by MapReady for georeferencing and, if desired, terrain correction. The bash script /home/mapready.sh in the home directory automates the procedure. MapReady will output the geocoded covariance/coherency matrix image in the form of co-registered GeoTiff files, one for each diagonal matrix element and two (real and imaginary parts) for each off-diagonal component. A python script /home/ingest.py is called automatically to combine these files to a single multi-band image in floating point format.
The ENL (equivalent number of looks) can (optionally) be estimated with the script enlml.py. A multivariate estimator is used as described by Anfinsen et al. (2009).
The bitemporal change detection algorithm is invoked with the bash script /home/wishart.sh. This script calls the Python programs /home/register.py to co-register the two images and then /home/wishart.py to perform the pixel-wise hypothesis tests. The test statistics $-2\rho\log Q$ and change probabilities ${\rm prob}(-2\rho\log Q\le z)$ are written to a two-band GeoTiff file. Additionally a change map showing changes at significance level 0.01 in red overlayed onto the span image is writen to a three-band (RGB) GeoTiff file. The multitemporal algorithm is invoked similarly with /home/omnibus.sh, which calls /home/register.py to co-register all of the images with the first in the series, and then /home/omnibus.py to perform the multitemporal omnibus test.
An Example: Radarsat-2 Quadpol Images¶
Returning now to the Radarsat-2 image acquired April 26, 2010, we will geocode it with MapReady (step 2 in the above processing chain). The bash script /home/mapready takes two arguments, the acquisition date in the format yyymmdd and the sensor (presently rs2quad or tsxdual):
!/home/mapready.sh 20100426 rs2quad
We see that the multi-look images were created from PolSARpro data with $4\times 3 = 12$ looks. This corresponds to a square pixel size of close to $20.5\times 20.5$ meters. The combined coherency matrix image at this resolution is stored in geotiff format.
Before we can display the geocoded image, we have to enable Matplotlib functionality within the notebook with the so-called magic command
%matplotlib inline
The command below will generate an RGB color composite of the three diagonal elements:
run /home/dispms -p [1,6,9] -d [300,300,1000,1000] -f /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3//polSAR.tif
The scene above was acquired over the city of Bonn, Germany (upper right hand corner with the Rhine river). The gray, featurless areas are mixed forest.
To check the number of looks, We will run the /home/enlml.py script (step 3 in the processing chain) on a spatial subset which includes a lot forested land cover. The spatial subset is entered with the -d flag as in the /home/dispms script.
Here we choose -d [800,400,500,500]:
run /home/enlml -d [800,400,500,500] /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3//polSAR.tif
There are two modes (maxima) at about 6 and 12 looks. Here is the ENL image
run /home/dispms -e 2 -f /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR_enl.tif
The bright areas correspond to mixed forest (around 12 equivalent looks) and allow for the best estimate of the ENL, see Anfinsen et al. (2009). So we can conclude that the nominal and eqivalent number of looks are the same and equal to 12.
Next we geocode two more Radarsat-2 images, the ones acquired on May 20, 2010 (20100520) and June 7, 2010 (20100707):
!/home/mapready.sh 20100520 rs2quad
!/home/mapready.sh 20100707 rs2quad
Now let's perform the last processing step, polSAR change detection. The bitemporal bash script /home/wishart.sh needs five input parameters, the two acquistion times in yyymmdd format, a spatial subset, and the two ENL values (which in principle can be different). We choose first of all the April and July images:
!/home/wishart.sh 20100426 20100707 [300,300,1000,1000] 12.0 12.0
Here is the 1% significance level change map image generated by the above script:
run /home/dispms -e 4 -p [1,2,3] -f /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/wishart(20100426-20100707)_cmap.tif
Over the interval late April - early July separating the two acquisitions the significant changes are mostly in the agricultural areas near center right and upper left. Barge movements on the Rhine river are evident. Zooming in on the upper left hand corner we can see a flooded sand quarry pit with two dredging arms that are in continual motion, giving rise to significant change signals.
run /home/dispms -e 4 -p [1,2,3] -d [100,0,400,400] -f /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/wishart(20100426-20100707)_cmap.tif
Finally, we'll run the omnibus test statistic on all three images. The script requires the acquisition times, the spatial subset, the number of looks and the desired significance level:
!/home/omnibus.sh 20100426 20100520 20100707 [300,300,1000,1000] 12 0.01
We will display the bitemporal and multitemporal change maps side-by-side:
run /home/dispms -e 4 -p [1,2,3] -E 4 -P [1,2,3] -d [100,0,400,400] -D [100,0,400,400] \
-f /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/wishart(20100426-20100707)_cmap.tif \
-F /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/omnibus(20100426-1-20100707)_cmap.tif
The multitemporal test (right hand image) detects more changes than the bitemporal test, at the same 1% significance level.