Radar ranging¶
%matplotlib inline
from numpy import *
import matplotlib.pyplot as plt
def chirp(t,t0):
result = 0.0*t
idx = array(range(2000))+t0
tt = t[idx] - t0
result[idx] = sin(2*math.pi*2e-3*(tt+1e-3*tt**2))
return result
t = array(range(5000))
plt.plot(t,chirp(t,400)+9)
plt.plot(t,chirp(t,800)+6)
plt.plot(t,chirp(t,1400)+3)
signal = chirp(t,400)+chirp(t,800)+chirp(t,1400)
kernel = chirp(t,0)[:2000]
kernel = kernel[::-1]
plt.plot(t,signal)
plt.plot(0.003*convolve(signal,kernel,\
mode='same')-5)
plt.xlabel('Time')
plt.ylim((-8,12))
plt.show()
The power received by a radar transmitting/receiving antenna reflected from a distributed (as opposed to point) target a distance $D$ from the antenna is given by
$$ P_R = {P_TG_TG_R\lambda^2\sigma^o\Delta_a\Delta_r\over (4\pi)^3D^4}\ [W],\quad (1) $$where $P_T$ is the transmitted power $[W\cdot m^{-2}]$, $\lambda$ is the operating wavelength~$[m]$, $G_T(G_R)$ is the transmitting (receiving) antenna gain, $\Delta_a(\Delta_r)$ is the azimuth (ground range) resolution $[m]$ and $\sigma^o$ is the unitless scattering coefficient (referred to as the {\it radar cross section}) of the target surface. The scattering coefficient is related to the (bio)physical properties of the surface being irradiated, notably its water content.
Polarimetric SAR¶
We will be concerned primarily with fully and partially polarimetric SAR data. A full, or quad, 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
$$ E^b = S E^i \quad\hbox{or}\quad \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}.\quad (2) $$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 Equation (2), the four complex scattering matrix elements.
Complex numbers¶
A complex number is an expression of the form $z=a+ib$, where $i^2=-1$. The real part of $z$ is $a$ and the imaginary part is $b$. The number $z$ can be represented as a point or vector in the complex plane with the real part along the $x$-axis and the imaginary part along the $y$-axis. Complex number addition then corresponds to addition of two-dimensional vectors:
$$ (a_1 + i b_1) + (a_2+i b_2) = (a_1+a_2) + i(b_1+b_2). $$Multiplication is also straightforward, e.g.,
$$ (a_1 + i b_1)(a_2+i b_2) = a_1a_2 - b_1b_2 +i (a_1b_2-a_2b_1). $$The complex conjugate of $z$ is $z^*=a-i b$. Thus
$$ z^*z = (a-i b)(a+i b) = a^2+b^2 = |z|^2, $$where $|z|=\sqrt{a^2+b^2}$ is the magnitude of the complex number $z$. If $\theta$ is the angle that $z$ makes with the real axis, then
$$ z = a + i b = |z|\cos(\theta) + i |z|\sin(\theta). $$From the well-known Euler's Theorem we can write this as
$$ z = |z|e^{i\theta}. $$Complex numbers provide a convenient representation for the amplitude $E$ of an electric field:
$$ E = |E|cos(\omega t + \phi) = {\rm Re}\left(|E|e^{i (\omega t +\phi)}\right), $$where $\omega = 2\pi f$ and $\phi$ are the angular frequency and phase of the radiation and Re denotes real part. It is usually convenient to work exclusively with complex amplitudes $|E|e^{i (\omega t +\phi)}$, bearing in mind that only the real part is physically significant. When the oscillating electric fields are described by complex numbers in this way, the scattering matrix elements are also complex.
A full polarimetric, or quad polarimetric SAR image then consists of four complex bands
$$ s_{hh}, s_{hv}, s_{vh}, s_{vv}, $$one for each pixel-wise determination of an element of the scattering matrix. So-called reciprocity, which normally applies to natural targets, implies that
$$ s_{hv} = s_{vh} $$so we might represent a measured quad pol pixel by the complex vector
$$ s = \pmatrix{s_{hh} \cr \sqrt{2}s_{hv} \cr s_{vv} }. $$The squared amplitudes of the scattering coefficients $|s_{hh}|^2, 2|s_{hv}|^2, |s_{vv}|^2$ constitute the radar cross sections for each polarization combination. These in turn replace $\sigma^o$ in Equation (1) and determine the received power in each polarization channel.
Speckle statistics¶
(snap)
The intensity of the backscattered signal in a polarimetric band is exponentially distributed. The argument runs as follows:
Even when the properties of the reflecting surface are uniform, random irregularities on the wavelength scale within the imaged pixel will act like a very large number of incremental scatterers, leading to a random superposition of small scattering amplitudes. These in turn combine coherently to form the backscattered signal and are responsible for speckle. A simple model for this effect is to write the backscattered amplitude in the form
$$ E^b = {E\over\sqrt{n}}\sum_{k=1}^n e^{i\phi_k}. $$Here we assume that the only effect of the incremental scatterers is to introduce local random phase shifts $\phi_k$, $k=1\dots n$, in an overall scattered signal amplitude $E$. The phase shifts may be considered to be uniformly distributed on the interval $[-\pi,\pi]$, so that we have
$$ |E^b|^2 = {|E|^2\over n}\sum_k\sum_\ell e^{i(\phi_k-\phi_\ell)} ={|E|^2\over n}\sum_k e^{i(\phi_k-\phi_k)} = |E|^2. $$We can express the real and imaginary parts in the form
$$ {\rm Re}(E^b) = {E\over\sqrt{n}}\sum_k \cos\phi_k = {E\over\sqrt{n}} X\ $$$$ {\rm Im}(E^b) = {E\over\sqrt{n}}\sum_k \sin\phi_k = {E\over\sqrt{n}} Y, $$where
$$ X=\sum_k\cos\phi_k, \quad Y=\sum_k\sin\phi_k. $$Now the intensity of the backscattered signal can be written
$$ |E^b|^2 = [{\rm Re}(E)]^2 + [{\rm Im}(E)]^2 = E^2{1\over n}(X^2+Y^2). $$We can think of $X$ and $Y$ as random variables. Since the sines and cosines contributing to $X$ and $Y$ are distributed symmetrically about zero, the means vanish:
$$ \langle X\rangle = \langle Y\rangle = 0. $$Furthermore,
$$ {\rm var}(X) = \langle X^2\rangle - \langle X\rangle^2 = \sum_k\sum_{\ell}\langle\cos\phi_k\cos\phi_\ell\rangle $$$$ = \sum_k \langle\cos^2\phi_k\rangle = \sum_k\langle{1\over 2}(\cos 2\phi_k+1)\rangle = {n\over 2}, $$with the same result for ${\rm var}(Y)$. Moreover, $X$ and $Y$ are uncorrelated:
$$ {\rm cov}(X,Y) = \langle XY\rangle -\langle X\rangle\langle Y\rangle = \sum_k\sum_{\ell}\langle\cos\phi_k\sin\phi_\ell\rangle = 0. $$Both $X$ and $Y$ are superpositions of a large number of identically distributed random variables ($\cos\phi_k$ and $\sin\phi_k$). According to the Central Limit Theorem they independently and normally distributed with zero mean and variance $\sigma^2=n/2$. The joint probability density of $X$ and $Y$ is therefore the bivariate distribution with
$$ \mu = \pmatrix{0\cr0},\quad \Sigma = {n\over 2}\pmatrix{1 & 0\cr 0 & 1} $$that is,
$$ p(x,y) = {1\over 2\pi (n/2)}\exp\left(-{1\over 2(n/2)}(x^2+y^2)\right). $$With $z = x+i y$, this can be written in the form
$$ p(z) = {1\over\pi n}\exp\left(-z^*\left({1\over n}\right)z\right). $$This shows that the complex random variable $Z= X + i Y$, and hence the backscattered amplitude
$$ E^b = {E\over\sqrt{n}} (X + i Y) = {E\over\sqrt{n}}Z, $$has a complex normal distribution.
The sum $U$ of the squares of the standardized random variables $X$ and $Y$,
$$ U = {X^2\over n/2} + {Y^2\over n/2} = {2\over n}(X^2+Y^2), $$is chi-square distributed with 2 degrees of freedom, that is,
$$ p_u(u) = p_{\chi^2;2}(u) = {1\over 2} e^{-u/2}. $$This is the exponential probability density with parameter $\beta=1/2$. The backscattered signal intensity is then
$$ |E^b|^2 = E^2 U/2. $$A measurement of a radar cross section, say $|s_{hh}|^2$, derived from the emitted and received polarized signal intensities will have the same exponential distribution. In other words, if $z$ is the measured intensity, $x = |s_{hh}|^2$ is the underlying signal and $u$ is the speckle, then
$$ z = x u/2. $$Applying the rule
$$ p_z(z) = p_u(u)\left |{du\over dz}\right | = {1\over 2} e^{-u/2}{2\over x} = {1\over x} e^{-z/x} $$Thus, speckle behaves as an exponentially distributed noise with mean and variance equal to the underlying signal strength.
Multi-look data¶
The measured amplitude signal for a quad pol image pixel is
$$ s = \pmatrix{s_{hh} \cr \sqrt{2}s_{hv} \cr s_{vv} }. $$The total instensity is referred to as the span,
$$ {\rm span} = s^+ s = |s_{hh}|^2 + 2|s_{hv}|^2 + |s_{vv}|^2, $$and the corresponding gray-scale image as the span image.
Multi-look processing essentially corresponds to the averaging of neighbourhood pixels, with the objective of reducing speckle and compressing the data. If m measurements such as
$$ z(\nu) = |s_{hh}|_\nu^2 $$are averaged to give
$$ \bar z = {1\over m}\sum_{\nu=1}^m z(\nu) $$then the sum
$$ m\bar z = \sum_{\nu=1}^m z(\nu) $$is gamma dstributed with $\alpha=m$ and $\beta = x$. So its mean and variance are
$$ \alpha\beta = mx,\quad \alpha\beta^2= mx^2. $$Then the mean of the multi-look pixel $\bar z$ is
$$ \langle\bar z\rangle = {1\over m}mx = x $$and its variance is
$$ {\rm var} (\bar z) = {\rm var}\left({1\over m}\sum_{\nu=1}^m z(\nu)\right) = {1\over m^2}mx^2 = {x^2\over m}.\quad (1) $$We see that the variance of the look-averaged image decreases inversely as the number of looks.
In real SAR imagery, the neighborhood pixel intensities contributing to the look average will be correlated to some extent. This is accounted for by defining an equivalent number of looks (ENL) whose definition is motivated by Equation (1), that is, by solving it for $m$:
$$ {\rm ENL} = m = {x^2\over {\rm var}(\bar z)}. $$The ENL can be estimated in homogeneous regions of an image, where the contribution of speckle to image statistics may be expected to dominate.
Polarimetric matrix representations and pre-processing¶
Another representation of the polarimetric signal is the complex covariance matrix $c$, obtained by taking the complex outer product of $s$ with itself:
$$ c = s s^+ = \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, and their sum is the span. The off-diagonal elements are complex. The covariance matrix representation contains all of the information in the polarized signal.
It can also be averaged over the number of looks to give
$$ \bar{c} ={1\over m}\sum_{\nu=1}^m s(\nu)s(\nu)^+ = \langle s s^+ \rangle $$or
$$ \bar{c}= \pmatrix{ \langle |s_{hh}|^2\rangle & \sqrt{2}\langle s_{hh}s_{hv}^*\rangle & \langle s_{hh}s_{vv}^*\rangle \cr \sqrt{2}\langle s_{hv}s_{hh}^*\rangle & \langle 2|s_{hv}|^2\rangle & \sqrt{2}\langle s_{hv}s_{vv}^*\rangle \cr \langle s_{vv}s_{hh}^*\rangle & \sqrt{2}\langle s_{vv}s_{hv}^*\rangle & \langle |s_{vv}|^2\rangle }. $$Rewriting the first equation above,
$$ m\bar{c} = \sum_{\nu=1}^m s(\nu)s(\nu)^+. $$This is seen to be a realization of a complex sample matrix with a complex Wishart distribution with $3\times 3$ covariance matrix $\Sigma$ and $m$ degrees of freedom, provided that the vectors $s(\nu)$ are independent and drawn from a complex multivariate normal distribution.As already mentioned, the $s(\nu)$ will generally be correlated. In order to account for this, the complex Wishart distribution is often parameterized with ENL (rather than $m$) degrees of freedom. The complex Wishart distribution of $m\bar c$ will be exploited later for change detection.
Another representation is the coherency vector
$$ t = \pmatrix{{1\over \sqrt{2}}(s_{hh}+s_{vv}) \cr {1\over\sqrt{2}}(s_{hh}-s_{vv}) \cr \sqrt{2}s_{hv} }, $$and the corresponding look-averaged coherency matrix
$$ \bar{k} = {1\over m}\sum_{\nu=1}^m t(\nu)t(\nu)^+ = \langle t t^+\rangle. $$$m\bar k$ is also complex Wishart distributed.
(snap)
ls -l /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC
cd /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC
ls -l polsarpro/T3
!/home/mapready.sh 20100426 rs2quad
run /home/dispms -f /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif -p [1,6,9]
The gamma MAP filter¶
We have seen that an $m$ look-averaged SAR intensity image is described by the random variable $Z$ with mean $\langle Z\rangle=x$ and variance ${\rm var}(Z) = x^2/m$, where $x$ is the underlying signal. Since $Z$ is gamma-distributed, its parameters $\alpha$ and $\beta$ must satisfy
$$ \alpha\beta = x,\quad \alpha\beta^2 = x^2/m, $$from which follows $\alpha= m$ and $\beta = x/m$. The probability density function of $Z$ for given value of $x$ is therefore
$$ p(z\mid x) = {1\over (x/m)^m\Gamma(m)}z^{m-1}e^{-zm/x}.\quad (1) $$Let
$$ z = xv. $$Then $v$ has the probability distribution
$$ p_v(v) = p_z(z)\cdot\left|{dz\over dv}\right| = {1\over (x/m)^m\Gamma(m)}x^{m-1}v^{m-1}e^{-mv}\cdot x = {m^m\over \Gamma(m)}v^{m-1}e^{-vm}. $$which is the gamma distribution with mean 1 and variance $1/m$. Because of this special multiplicative noise nature of speckle, conventional smoothing filters of the kind used for multispectral imagery are not particularly suitable as an aid to SAR image interpretation.
According to Bayes' Theorem, the posterior probability for signal $x$, given observation $z$, is
$$ {\rm Pr}(x\mid z) = { p(z\mid x){\rm Pr}(x)\over p(g) },\quad (2) $$where $p(z\mid x)$ is given by Equation (1). This formulation allows us to include prior knowledge of the signal statistics (or texture) if available. An empirical statistical model for $x$ is suggested by measurements of backscatter from ocean waves:
$$ {\rm Pr}(x) \sim \left({\alpha\over\mu}\right)^\alpha {x^{\alpha-1}\over\Gamma(\alpha)}e^{-\alpha x/\mu}. \quad (3) $$This is just the gamma probability density with $\beta=\mu/\alpha$, and hence with mean $\alpha\beta= \mu$ and variance
$$ {\rm var}(x) = \alpha\beta^2 = \mu^2/\alpha. $$The parameters $\mu$ and $\alpha$ must be estimated from the image intensity data.
Combining Equations (1), (2) and (3) and neglecting the total probability $p(g)$, which doesn't depend on $x$, the posterior probability for $x$ given measurement $z$ is
$$ {\rm Pr}(x\mid z) \sim {1\over (x/m)^m\Gamma(m)}z^{m-1}e^{-zm/x}\left({\alpha\over\mu}\right)^\alpha {x^{\alpha-1}\over\Gamma(\alpha)}e^{-\alpha x/\mu} =: L $$We get the maximum a posteriori (MAP) value for $x$ given the observed pixel intensity $z$ by maximizing $\log L$ with respect to $x$:
$$ {d\over dx}\log L = -m/x + mz/ x^2 + (\alpha-1)/x - \alpha/\mu = 0. $$This leads to a quadratic equation for the most probable signal intensity $x$,
$$ {\alpha\over\mu}x^2 + (m+1-\alpha)x - mz = 0, $$where the parameters $\mu$ and $\alpha$ are estimated locally in an $n\times n$ window:
$$ \mu \approx \langle z\rangle_n $$and
$$ \alpha \approx {\mu^2\over {\rm var}(x)}. $$We can get ${\rm var}(x)$ as follows. Recall that $z = xv$ where $\langle v\rangle = 1$ and ${\rm var}(v) = 1/m$. Since signal $x$ and speckle noise $v$ are independent,
$$ \langle z^2\rangle_n = \langle x^2\rangle_n\langle v^2\rangle_n $$or
$$ {\rm var}(z)+\langle z\rangle_n^2 = ({\rm var}(x) + \langle x\rangle_n^2)({\rm var}(v) + \langle v\rangle_n^2). $$Solving for ${\rm var}(x)$ we obtain, with $\langle x\rangle_n = \langle z\rangle_n =\mu$,
$$ {\rm var}(x) = {{\rm var}(z) - {1\over m}\mu^2 \over 1 + {1\over m}}. $$run /home/gamma_filter -h
Directional filter
run /home/gamma_filter -d [400,400,400,400] /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif 12
run /home/gamma_filter -p -d [400,400,400,400] /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif 12
run /home/dispms -p [1,6,9] -P [1,2,3] -d [400,400,400,400] \
-f /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif \
-F /home/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR_gamma.tif
