Theoretical Modeling for Polarimetric Scattering and Information Retrieval of SAR Remote Sensing

1. Introduction

Synthetic aperture radar (SAR) imagery technology is one of most important advances in space-borne microwave remote sensing during recent decades. Completely polarimetric scattering from complex terrain surfaces can be measured. Fully understanding and retrieving information from polarimetric scattering signatures of natural media and SAR images have become a key issue for the SAR remote sensing and its broad applications.

This paper presents some research progress in Fudan University on theoretical modeling of the terrain surface for polarimetric scattering simulation and Mueller matrix solution, mono-static and bistatic SAR image simulation, new parameters for unsupervised surface classification, DEM inversion, change detection from multi-temporal SAR images, and reconstructions of buildings from multi-aspect SAR images, etc. Some applications are discussed.

2. Model of Vegetation Canopy and Mueller Matrix solution

As a polarized wave $$E_{inc}(\chi ,\psi )$$ is incident upon the natural media, the scattering field is written as

$$\left[\begin{array}{l}{E}_{vs}\\ E_{hs}\end{array}\right]=\frac{e^{ikr}}{r}\left[\begin{array}{l}{S}_{vv}{S}_{vh}\\ S_{hv}{S}_{hh}\end{array}\right]\cdot \left[\begin{array}{l}{E}_{vi}\\ E_{hi}\end{array}\right]\equiv \frac{e^{ikr}}{r}\overline{S}\cdot E_{inc}(\chi ,\psi )$$ E1

where 22-D (dimensional) complex scattering amplitude function $$\overline{S}$$ can be measured by the polarimetry. The incident polarization is indicated by the elliptic and orientation angles $$(\chi ,\psi )$$ . Using the Mueller matrix solution of vector radiative transfer equation ( Jin, 1994 ) and Eq. (1), the scattered Stokes vector (four Stokes parameters) can be obtained or measured as (Jin, 2005)

$$I_s(\theta ,\varphi )=\overline{M}(\theta ,\varphi ;\pi -{\theta }_0,{\varphi }_0)\cdot I_i(\chi ,\psi )$$ E2

The Mueller matrix solution of a layer of scatterers model, as shown in Figure 1, is written as

$$\overline{M}(\theta ,\varphi ;\pi -{\theta }_0,{\varphi }_0)={\overline{M}}_{\text{0}}+{\overline{M}}_{\text{1}}+{\overline{M}}_{\text{2}}+{\overline{M}}_{\text{3}}+{\overline{M}}_{\text{4}}$$ E3

$$\begin{array}{c}{\overline{M}}_0=\exp \left[-d\sec \theta {\overline{\kappa }}_e(\theta ,\varphi )\right]\cdot \overline{R}(\theta ,\varphi ;\pi -{\theta }_0,{\varphi }_0)\cdot \exp \left[-d\sec {\theta }_0{\kappa }_e(\pi -{\theta }_0,{\varphi }_0)\right]\\ {\overline{M}}_1=\sec \theta {\displaystyle {\int }_{-d}^{0}d{z}^{\prime }\exp \left[z^{\prime }\sec \theta {\overline{\kappa }}_e(\theta ,\varphi )\right]}\cdot \overline{P}(\theta ,\varphi ;\pi -{\theta }_0,{\varphi }_0)\cdot \exp \left[z^{\prime }\sec {\theta }_0{\overline{\kappa }}_e(\pi -{\theta }_0,{\varphi }_0)\right]\\ {\overline{M}}_2=\sec \theta {\displaystyle {\int }_{-d}^{0}d{z}^{\prime }\exp \left[z^{\prime }\sec \theta {\overline{\kappa }}_e(\theta ,\varphi )\right]\cdot {\displaystyle {\int }_0^{2\pi }d{\varphi }^{\prime }{\displaystyle {\int }_0^{\pi /2}d{\theta }^{\prime }\sin {\theta }^{\prime }}}}\overline{P}(\theta ,\varphi ;{\theta }^{\prime },{\varphi }^{\prime })\\ \cdot \exp \left[(-d-z^{\prime })\sec {\theta }^{\prime }{\overline{\kappa }}_e({\theta }^{\prime },{\varphi }^{\prime })\right]\cdot \overline{R}({\theta }^{\prime },{\varphi }^{\prime };\pi -{\theta }_0,{\varphi }_0)\cdot \exp \left[-d\sec {\theta }_0{\overline{\kappa }}_e(\pi -{\theta }_0,{\varphi }_0)\right]\\ {\overline{M}}_3=\sec \theta \exp \left[-d\sec \theta {\overline{\kappa }}_e(\theta ,\varphi )\right]\cdot {\displaystyle {\int }_0^{2\pi }d{\varphi }^{″}{\displaystyle {\int }_0^{\pi /2}\sin {\theta }^{″}d{\theta }^{″}}}\overline{R}(\theta ,\varphi ;\pi -{\theta }^{″},{\varphi }^{″})\\ \cdot {\displaystyle {\int }_{-d}^{0}d{z}^{\prime }\exp \left[(-d-z^{\prime })\sec {\theta }^{″}{\overline{\kappa }}_e(\pi -{\theta }^{″},{\varphi }^{″})\right]}\overline{P}(\pi -{\theta }^{″},{\varphi }^{″};\pi -{\theta }_0,{\varphi }_0)\cdot \exp \left[z^{\prime }\sec {\theta }_0{\overline{\kappa }}_e(\pi -{\theta }_0,{\varphi }_0)\right]\\ {\overline{M}}_4=\sec \theta \exp \left[-d\sec \theta {\overline{\kappa }}_e(\theta ,\varphi )\right]\cdot {\displaystyle {\int }_0^{2\pi }d{\varphi }^{″}{\displaystyle {\int }_0^{\pi /2}d{\theta }^{″}\sin {\theta }^{″}}}\overline{R}(\theta ,\varphi ;\pi -{\theta }^{″},{\varphi }^{″})\\ \cdot \exp \left[(-d-z^{\prime })\sec {\theta }^{″}{\overline{\kappa }}_e(\pi -{\theta }^{″},{\varphi }^{″})\right]{\displaystyle {\int }_0^{2\pi }d{\varphi }^{\prime }{\displaystyle {\int }_0^{\pi /2}d{\theta }^{\prime }\sin {\theta }^{\prime }\overline{P}(\pi -{\theta }^{″},{\varphi }^{″};{\theta }^{\prime },{\varphi }^{\prime })}}\\ \cdot \exp \left[(-d-z^{\prime })\sec {\theta }^{\prime }{\overline{\kappa }}_e({\theta }^{\prime },{\varphi }^{\prime })\right]\overline{R}({\theta }^{\prime },{\varphi }^{\prime };\pi -{\theta }_0,{\varphi }_0)\exp \left[-d\sec {\theta }_0{\overline{\kappa }}_e(\pi -{\theta }_0,{\varphi }_0)\right]\end{array}$$ E4

where $${\overline{\kappa }}_e,\overline{P},\overline{R}$$ are the extinction matrix, phase matrix of non-spherical scatterers, and reflectivity matrix of underlying rough surface, respectively. Five contributions from scatterers and underlying surface are illustrated in Figure 1. $${\overline{\kappa }}_e,\overline{P}$$ can be expressed by $$S_{pq}{S}_{st}^{*}$$ of $$\overline{S}$$ in Eq. (1). The co-polarized and cross-polarized backscattering coefficients, $${\sigma }_c$$ and $${\sigma }_x$$ , polarization degree $$m_s$$ for scattered Stokes echo with partial polarization and other functions can be numerically calculated ( Jin and Chen, 2002 ).

Eqs. (4a-c) will be applied to numerical simulation of polarimetric scattering from the terrain surface with physical parameters of scatterers (Jin, 2005), as described in next sections.

Figure 1.

A model of non-spherical scatterers for vegetation canopy.

The Mueller matrix is a 4×4-D real matrix with complex eigenvalues and eigenvectors. To be physically realizable, this matrix must satisfy the Stokes criterion together with several other restrictive conditions. Unfortunately, however, these restrictions do not have any direct physical interpretation in terms of the eigenstructure of the Mueller matrix.

The coherency matrix $$\overline{C}$$ is applied to the study of polarimetric scattering of SAR image ( Jin and Cloude, 1994 ). Define the scattering vector as

$$k_P=\frac{1}{2}{[S_{vv}+S_{hh},S_{vv}-S_{hh},S_{vh}+S_{hv},i(S_{vh}-S_{hv})]}^T=\frac{1}{2}{[A,B,C,iD]}^T$$ E5

where A, B, C, D are sequentially defined in Eq. (5), and the subscript P denotes Pauli vectiorization. The coheremcy matrix is defined as

$$\overline{C}=k_P{k}_P^{+}=\frac{1}{4}\left[\begin{array}{cccc}AA*& AB*& AC*& AD*\\ BA*& BB*& BC*& BD*\\ CA*& CB*& CC*& CD*\\ DA*& DB*& DC*& DD*\end{array}\right]$$ E6

The eigenvalue and eigenvector of the coherency matrix are defined as

$$\overline{C}=k_P{k}_P^{+}={\displaystyle \sum _{i=1}^4{\lambda }_i}{k}_{Pi}{k}_{Pi}^{+}$$ E7

All eigenvalues are real and non-negative, and $${\lambda }_1\ge {\lambda }_2\ge {\lambda }_3\ge {\lambda }_4$$ .

The coherency matrix is also related with the Mueller matrix as

$$M_{i+\mathrm{1,}s+1}={\displaystyle \sum _t{\displaystyle \sum _u\frac{1}{2}\text{Trace}({\overline{\sigma }}_i{\overline{\sigma }}_t{\overline{\sigma }}_s{\overline{\sigma }}_u){[\overline{C}]}_{tu}}}$$ E8

where $${\overline{\sigma }}_i,i=\mathrm{0,1,2,3}$$ are the Pauli matrices.

Usually, the case of $$S_{vh}=S_{hv}$$ is simply assumed. The entropy H is defined by the eigenvalues of the coherency matrix as

$$H=-{\displaystyle \sum _{i=1}^3{P}_i{\log }_3{P}_i},P_i={\lambda }_i/({\lambda }_1+{\lambda }_2+{\lambda }_3)$$ E9

The entropy is an important feature since it relates the randomness of scatter media with other physical parameters, such as canopy configuration, biomass, etc.

The coherency matrix is 4×4-D positive semidefinite Hermitian and as such has real non-negative eigenvalues and complex orthogonal eigenvectors. Amplitude and difference of four eigen-values are functionally related with polarimetric scattering process of the terrain canopy. Eigen-analysis of the coherency matrix yields better physical insight into the polarimetric scattering mechanisms than does the Mueller matrix and further, it can be employed to physically identify a Mueller matrix by virtue of the semi-definite nature of the corresponding coherency matrix.

For the first-order solution in small albedo, cross polarization is small and correlations such as the terms $$CA*,CB*,DA*,DB*$$ etc. in Eq. (6) are always very small, and might be neglected in our following derivations. Actually, for weak assumption of azimuthal symmetric orientation, these correlations have been proved to be zero. Applying this approximation to Eq. (6), the eigenvalues of $$\overline{C}$$ have been derived as ( Jin and Cloude, 1994 )

$${\lambda }_{\mathrm{1,2}}=\frac{1}{8}[|A{|}^2+|B{|}^2\pm \sqrt{{(|A{|}^2-|B{|}^2)}^2+4A^{*}BA{B}^{*}}],{\lambda }_3=\frac{1}{4}|C{|}^2$$ E10

It is interesting to see that if $$AB*A*B=AA*BB*$$ at an ordered case, it yields $${\lambda }_2=0$$ .

Substituting A,B,C and D of Eq. (6) into Eq. (10), the eigen-values are derived using the power terms of $$|S_{vv}{|}^2,|S_{hh}{|}^2,|S_{hv}{|}^2$$ as follows,

$${\lambda }_1=\frac{1}{2}(|S_{vv}{|}^2+|S_{hh}{|}^2)-{\lambda }_2,{\lambda }_2=\frac{1}{2}\frac{|S_{vv}{|}^2|S_{hh}{|}^2}{(|S_{vv}{|}^2+|S_{hh}{|}^2)}(1-\Delta ),{\lambda }_3=|S_{hv}{|}^2$$ E11

where the configuration parameter

$$\Delta =\frac{S_{vv}{S}_{hh}{}^{*}{S}_{vv}{}^{*}{S}_{hh}}{|S_{vv}{|}^2|S_{hh}{|}^2}$$ E12

varies from 0 of totally random media to 1 of ordered media, i.e. no-random case.

It can be seen that the first eigen-value $${\lambda }_1$$ now indicates the total vv and hh co-polarized power subtracting their coherence defined by the parameter of Eq. (12), which varies from 1 at ordered case to 0 in totally random medium. The second eigen-value $${\lambda }_2$$ now takes into account for the vv and hh coherent power. As approaches to one, $${\lambda }_2\to \text{ 0}$$ for ordered media, and as approaches to zero, $${\lambda }_2$$ would increase to indicate the loss of the vv and hh power coherency for totally random media. The third eigen-value $${\lambda }_3$$ is now due to depolarization caused by the media randomness.

The (p,q=v,h)-polarized backscattering coefficient is

$${\sigma }_{pq}=4\pi \cos \theta |S_{pq}{|}^2$$ E13

Substituting Eq. (13) into Eqs. (11,12), the functions $$P_i,i=\mathrm{1,2,3}$$ , of Eq. (9) are then derived as

$$P_1=1-P_2-P_3,P_2=X(1-\Delta )(1-\delta ),P_3=\delta (1-\delta )1$$ E14

where

$$X\equiv \frac{{\sigma }_{hh}/{\sigma }_{vv}}{{(1+{\sigma }_{hh}/{\sigma }_{vv})}^2},\delta =\frac{2{\sigma }_{hv}}{{\sigma }_{vv}+{\sigma }_{hh}}$$ E15

Note that all functions $$P_i,i=\mathrm{1,2,3}$$ , can be now calculated by polarized backscattering coefficient $${\sigma }_{pq},p,q=v,h$$ , not by Eq. (7) from full Mueller matrix solution. Then, the entropy H of Eq. (9) is calculated using $$P_i$$ of Eq. (14). It can be seen that the function $$P_1$$ is related with the total vv and hh co-polarized power $${\sigma }_{hh}$$ + $${\sigma }_{vv}$$ ; $$P_2$$ is due to the difference between $${\sigma }_{hh}$$ and $${\sigma }_{vv}$$ ; and $$P_3$$ is due to depolarization $${\sigma }_{hv}$$ . They are modulated by the media configuration and randomness via the parameters $$\Delta$$ and $$\delta$$ . As the medium becomes more random or disordered, $$\delta$$ increases and $$\Delta$$ approaches zero. Vice versa, as the medium becomes more ordered, $$\delta$$ decreases and $$\Delta$$ increases. Thus, the relationship between the entropy H and backscattering measurements $${\sigma }_{pq}(pq=vv,hh,hv)$$ is established.

Define the co-polarized and cross polarized indexes, respectively, as

$$CPI\equiv 10{\log }_{10}({\sigma }_{hh}/{\sigma }_{vv})(dB),XPI\equiv \delta$$ E16

The entropy H is directly related with backscattering indexes CPI and XPI via Eqs. (14~16).

Figure 2(a) presents an AirSAR image of total power of $${\sigma }_{vv},{\sigma }_{hh}$$ and $${\sigma }_{vh}$$ at the L band. The image data from some typical areas such as the lake, an island surface, sparse trees, and thick forest (see the frames in the figure) are chosen, as shown in the discrete points of Figure 2(b), and are compared with theoretical results of H, CPI and XPI. The H is first rigorously calculated by polarimetric data of the Mueller matrix, i.e. the eigenvalues of coherency matrix $$\overline{C}$$ . Using the AirSAR data $${\sigma }_{vv},{\sigma }_{hh},{\sigma }_{vh}$$ and the averaged $$\delta$$ in the respective region, the lines are calculated with appropriate $$\Delta$$ to match the center of discrete data.

Figure 2.

An AirSAR image of total power $${\sigma }_{vv},{\sigma }_{hh},{\sigma }_{vh}$$ at L band and Entropy H vs CPI. Line 1: $$\Delta =0.257$$ , $$\delta =0.176$$ , Line 2: $$\Delta =0.318$$ , $$\delta =0.156$$ , Line 3: $$\Delta =0.365$$ , $$\delta =0.118$$ , Line 4: $$\Delta =0.$$ , $$\delta =0.$$ , Line 5: $$\Delta =0.731$$ , $$\delta =0.0265$$ .

3. SAR Imaging Simulation (MPA)

Given the sophisticated mechanisms of electromagnetic (EM) wave scattering and SAR signal collecting, it is hard to tackle SAR imagery without the aid of forward simulation of scattering based on physical modeling. Simulation of SAR imaging presents a different mode to study inhomogeneous terrain scene under SAR observation. It may evaluate or predict the SAR observation, interpret the physical scattering behavior, and take account of the parameterized characteristics of terrain surface objects. Some approaches to SAR image simulation in previous literature might be classified into two catalogues: one focuses on the statistical characteristics of SAR images, and the other on physical scattering process. Natural scene is more complicated including randomly, inhomogeneously distributed, penetrable or impenetrable objects, such as vegetation canopies, manmade structures and perturbed surface topography. Particularly, volumetric scattering through penetrable scatterers, e.g. timberland forests, crops, green plants etc., plays an essential or dominant role in SAR imagery. It is meaningful to develop computerized simulation of SAR imaging over comprehensive terrain scene with heterogeneous terrain objects.

we present an novel approach of the Mapping and Projection Algorithm (MPA) to polarimetric SAR imaging simulation for comprehensive scenarios, which takes account of scattering, attenuation, shadowing and multiple scattering of spatially distributed volumetric and surface scatterers ( Xu and Jin, 2006 ). In this approach, scattering contributions from scatterers are cumulatively summed on the slant range plane (the mapping plane) following geometric principles, while their extinction and shadowing effects are cumulatively multiplied on the ground range plane (the projection plane). It finally constructs a general imaging picture of the whole terrain scene via mapping and projection operations. A MPA is then devised to speed up the simulation of the whole process of scattering, extinction, mapping and projection in association with a grid partition of the comprehensive terrain scene. Our SAR simulation scheme incorporates polarimetric scattering, attenuation or shadowing of several typical terrain surfaces, as well as the coherent speckle generation.

3.1. The mapping and projection algorithm (MPA)

Since the slant range is much larger than the synthetic aperture, the incidence angle to the same target is nearly invariant during the interval of radar flying over a length of synthetic aperture. Thus in the imaging simulation, the whole scene is divided into cross-track lines along the azimuth dimension, and then the scattering contribution from the terrain objects lying in each line (included in the incidence plane) are calculated sequentially. It finally constructs a scattering picture or a scattering map of the whole scene.

Following from VRT, we derived a general expression of the scattering power at range $r$ of $x$ -position in azimuth dimension as

$$S(x,r)={\displaystyle \underset{{\theta }_0}{\overset{{\theta }_1}{\int }}\exp \left[-{\displaystyle \underset{r_0}{\overset{r}{\int }}\text{d}{r}^{\prime }{\kappa }_e^{+}(x,r^{\prime },\theta )}\right]P(x,r,\theta )\exp \left[-{\displaystyle \underset{r}{\overset{r_0}{\int }}\text{d}{r}^{\prime }{\kappa }_e^{-}(x,r^{\prime },\theta )}\right]r\text{d}x\text{d}r\text{d}\theta }$$ E17

where the phase function $P$ stands for scattering, ${\kappa }_e^{-}(x,r,\theta )$ , ${\kappa }_e^{+}(x,r,\theta )$ for the backward and forward extinction of an differential element $\text{d}v$ at the position $\left(x,r,\theta \right)$ of the imaging space. $\theta$ is incident angle. Integrating $S(x,r)$ over the $n,i$ -th pixel i.e. $x\in \left[x_n,x_{n+1}\right)$ , $s\in \left[s_i,s_{i+1}\right)$ , the corresponding discrete scattering map $S_{n,i}$ can be obtained. Considering polarimetric scattering, the computational discrete form of Eq. (17) for both volume and surface scatterer is given as

$${\overline{S}}_{n,i}={\displaystyle \sum _{m=m_n}^{m_{n+1}-1}{\displaystyle \sum _{p=p_i}^{p_{i+1}-1}{\displaystyle \sum _{q=0}^{Q-1}{\displaystyle \prod _{p^{\prime }=0}^p{\overline{E}}_o^{+}(m,p^{\prime },q)}{\overline{S}}_o(m,p,q){\displaystyle \prod _{p^{\prime }=p}^0{\overline{E}}_o^{-}(m,p^{\prime },q)}}}}$$ E18

$${\overline{S}}_o(m,p,q)=\{\begin{array}{cc}{v}_{ef}{\overline{S}}_{\text{vol}}(m,p,q)& \text{for volume scatterer}\\ a_{ef}{\overline{S}}_{\text{surf}}(m,p,q)& \text{for surface scatterer}\end{array}$$ E19

$${\overline{E}}_o^{\pm }(m,p^{\prime },q)=\{\begin{array}{cc}\exp \left[-d_{ef}{\overline{\kappa }}_e^{\pm }(m,p^{\prime },q)\right]& \text{for volume scatterer}\\ 0& \text{for surface scatterer}\end{array}$$ E20

where $v_{ef}=\gamma v_0$ is the effective volume of volume scatterer, calculated from particle number $\gamma$ and volume $v_0$ ; $a_{ef}=\Delta x\Delta s$ is the effective area of surface scatterer; $d_{ef}$ is the effective penetration depth. All of them take into account the proportion coefficients of each scatterer in the same discrete unit. The subscript $o$ denotes scattering or extinction of a single scatterer. Eq. (18) presents a general description of SAR imaging which deals with single scattering, attenuation or shadowing of all possible objects in the imaging space.

We further devise a mapping and projection algorithm (MPA) to fast compute Eq. (17) The first step is to partition the scene and terrain objects lying in the horizontal plane $xoy$ into multi-grids. Each grid unit contains its corresponding parts of the terrain objects overhead, e.g. ground, vegetation, building, etc.

As shown in Figure 3, two arrays $${\underset{\_}{E}}^{\pm }$$ and $\underset{\_}{S}$ are deployed for temporary storage of scattering and extinction (the under-bar denotes an array). The mapping plane and projection plane are equally divided into discrete cells corresponding to elements of $\underset{\_}{S}$ , $${\underset{\_}{E}}^{\pm }$$ . Attenuation or shadowing effect of each grid unit is counted in the increasing sequence of $y$ axis and cumulatively multiplied into the array $${\underset{\_}{E}}^{\pm }$$ .

Figure 3.

The MPA algorithm for natural media scene.

The MPA calculating sequence of the terrain objects in one grid unit is: first to calculate the scattering power of the current terrain object and accumulate onto the array $\underset{\_}{S}$ , where the attenuation or shadowing from other terrain objects are obtained from array $${\underset{\_}{E}}^{\pm }$$ ; then, to calculate the attenuation or shadowing caused by the current terrain object itself, and cumulatively multiply onto the array $${\underset{\_}{E}}^{\pm }$$ ; ultimately when all grid units of the current line are counted, the scattering map is exactly produced in the array $\underset{\_}{S}$ .

As shown in Figure 3, a vertical segment (length $L$ ) of volume scatterers above the current grid unit is firstly divided into $K$ sub-segments $L_k$ in order to make each sub-segment mapping into one cell ${\text{SI}}_k$ in mapping plane, while its mid-point is projected to cell ${\text{EI}}_k$ in projection plane. After the mapping and projecting, the elements of array $\underset{\_}{S}$ are refreshed as

$$\begin{array}{cc}\underset{\_}{S}[{\text{SI}}_k]:=\underset{\_}{S}[{\text{SI}}_k]+{\underset{\_}{E}}^{+}[{\text{EI}}_k]\cdot {\overline{S}}_o^k\cdot {\underset{\_}{E}}^{-}[{\text{EI}}_k],& k=\mathrm{1,..,}K\end{array}$$ E21

${\overline{S}}_o^k$ is the scattering of the $k$ -th sub-segment. Assuming that volume scatterers always uniformly fill the whole sub-segment, then the effective volume in Eq. (19) is $$v_{ef}=\Delta x\Delta y{L}_k{f}_s$$ . Therefore, we have

$${\overline{S}}_o^k=\Delta x\Delta y{L}_k{f}_s{\overline{S}}_{\text{vol}}=v_{ef}{\overline{S}}_{\text{vol}}$$ E22

where $f_s$ is the fractional volume of scatterers. Otherwise, it is necessary to estimate the total number of scatterers $\gamma$ within the sub-segment, and calculate scattering from Eq. (19). The next step is to refresh the elements of the array $${\underset{\_}{E}}^{\pm }$$ as

$$\{\begin{array}{cc}{\underset{\_}{E}}^{-}[{\text{EI}}_k]:={\overline{E}}_{\text{o}\text{,}k}^{-}\cdot {\underset{\_}{E}}^{-}[{\text{EI}}_k]& k=\text{1}\mathrm{,..,}\text{K}\\ {\underset{\_}{E}}^{+}[{\text{EI}}_k]:={\underset{\_}{E}}^{+}[{\text{EI}}_k]\cdot {\overline{E}}_{\text{o}\text{,}k}^{+}& k=\text{1}\mathrm{,..,}\text{K}\end{array}$$ E23

where ${\overline{E}}_{o,k}^{\pm }=\exp (-d_{ef}{\overline{\kappa }}_e)$ is the attenuation of the $k$ -th sub-segment. Due to discrete calculation for attenuation, the array $${\underset{\_}{E}}^{\pm }$$ should store the averaged attenuation. Thus, the ratio of shadowed area of the sub-segment to the whole projection cell should be taken into account, when counting the contribution of the sub-segment to the mean attenuation. The effective depth in ${\overline{E}}_{o,k}^{\pm }$ is calculated as

$$d_{ef}=\frac{v_{ef}}{R_x{R}_r/\tan \theta }=\frac{\Delta x\Delta y{L}_k}{R_x{R}_r/\tan \theta }=\tan \theta \frac{\Delta x\Delta y{L}_k}{R_x{R}_r}$$ E24

It can be seen from the effective depth $d_{ef}$ that if we redistribute all volume scatterers of the sub-segment and make them cover the whole projection cell under the invariance of density, then the layer depth of redistributed scatterers along the projection line is seen as the effective depth.

In a similar way, larger surface scatterers, like vertical wall surfaces, cliff slopes etc., are segmented according to mapping cells. But, for nearly horizontal surfaces, e.g. ground or roof surface, the facet cut by one grid unit is small enough and is projected into one projection cell, i.e. $K=1$ . Hence, given the normal vector of facet $\widehat{n}$ , the effective area and its scattering in Eq. (21) are calculated in two ways:

$${\overline{S}}_o^k=a_{ef}{\overline{S}}_{\text{surf}}=\{\begin{array}{cc}\Delta x\Delta y{\overline{S}}_{\text{surf}}\cos {\varphi }_2/\cos {\varphi }_1& K=1\\ \Delta x{L}_k{\overline{S}}_{\text{surf}}\cos {\varphi }_2/\sin {\varphi }_1& K1\end{array}$$ E25

where $${\varphi }_1$$ is the angle between $\widehat{n}$ and $z$ axis; ${\varphi }_2$ is the angle between $\widehat{n}$ and $\widehat{i}$ incidence. $L_k$ is the length of sub-segment. Moreover, the array $${\underset{\_}{E}}^{\pm }$$ should be set to zero over all projection interval from the first projection cell to the last one, i.e.

$$\begin{array}{cc}{\underset{\_}{E}}^{\pm }[j]:=0& {\text{EI}}_1\le j\le {\text{EI}}_K\end{array}$$ E26

The elements of the array $${\underset{\_}{E}}^{\pm }$$ are initialized to unit matrices at the start point of calculation for each line, while the array $\underset{\_}{S}$ should be saved before cleared to start the next line. Note that each line is strictly computed in the sequence of increasing $y$ in order to precisely count attenuation caused by the terrain objects ahead.

Double and triple scatterings between terrain objects and the underlying ground are regarded as shown in Figure 3. First, the corresponding ground patch of multiple scattering can be located by the ray tracing technique. Then the length of propagation path, i.e. the effective range $r_{ef}$ can be given as

$$r_{ef}=\{\begin{array}{ll}\left(r_{\text{Object}}+r_{\text{Ground}}+r_{\text{Path}}\right)/\text{2}\hfill & \text{for double scattering}\hfill \\ r_{\text{Object}}+r_{\text{Path}}\hfill & \text{for triple scattering}\hfill \end{array}$$ E27

In the same way, multiple scattering of a segment of scatterers is mapped into the array $\underset{\_}{S}$ by $r_{ef}$ . The same step of dividing into sub-segments is employed if the mapping interval exceeds one cell.

For double scattering, the attenuation or shadowing suffered through the propagation paths of $r_{\text{Object}}$ , $r_{\text{Ground}}$ are the same with single scattering of the terrain object and ground patch, respectively. It means that the array $${\underset{\_}{E}}^{\pm }$$ can be used as well. The attenuation during the way from the terrain object to ground patch $r_{\text{Path}}$ , is omitted for simplicity. Therefore, the double scattering contribution can be written as

$$\underset{\_}{S}[{\text{SI}}_{ef}]:=\underset{\_}{S}[{\text{SI}}_{ef}]+{\underset{\_}{E}}^{+}[{\text{EI}}_{\text{Ground}}]\cdot {\overline{S}}_{\text{Ground}}^{\text{2+}}\cdot \rho {\overline{S}}_{\text{Object}}^{2-}\cdot {\underset{\_}{E}}^{-}[{\text{EI}}_{\text{Object}}]+{\underset{\_}{E}}^{+}[{\text{EI}}_{\text{Object}}]\cdot \rho {\overline{S}}_{\text{Object}}^{2+}\cdot {\overline{S}}_{\text{Ground}}^{2-}\cdot {\underset{\_}{E}}^{-}[{\text{EI}}_{\text{Ground}}]$$ E28

where ${\text{EI}}_{\text{Ground}}$ , ${\text{EI}}_{\text{Object}}$ are the projection indices of the terrain object and the ground patch respectively. $${\overline{S}}_{\text{Ground}}^{2\pm }$$ , $${\overline{S}}_{\text{Object}}^{2\pm }$$ are the Mueller matrices of the ground patch and the terrain object along the forward and backward directions of double-scattering, respectively. The coefficient $$\rho =v_{ef}$$ or $$a_{ef}$$ is calculated similarly to Eqs. (19,22,25).

For triple scattering (object-ground-object), the expression can be directly written as:

$$\underset{\_}{S}[{\text{SI}}_{ef}]:=\underset{\_}{S}[{\text{SI}}_{ef}]+{\underset{\_}{E}}^{+}[{\text{EI}}_{\text{Object}}]\cdot {\overline{S}}_{\text{Object}}^{2+}\cdot {\overline{S}}_{\text{Ground}}^{\text{3}}\cdot \rho {\overline{S}}_{\text{Object}}^{2-}\cdot {\underset{\_}{E}}^{-}[{\text{EI}}_{\text{Object}}]$$ E29

where $${\overline{S}}_{\text{Ground}}^3$$ is the Mueller matrix of ground patch along the triple scattering path; $$\rho =v_{ef}$$ or $$a_{ef}$$ .

Correct sequence for calculation of each grid unit is first to count single scattering, then multiple scattering, at last its attenuation or shadowing. Notice that the ground patch of multiple scattering could be located at any place around the terrain object. For simplicity, we take the projection cell in the current incidence plane as an approximate substitute. Some other coding techniques are adopted to guarantee the correct computation and make it more efficient.

3.2. Scattering models for terrain objects

For vegetation canopies, VRT model of a layer of random non-spherical particles ( Jin, 1994 ) is employed. Leaves, small twigs or thin stems are modeled as non-spherical dielectric particles under the Rayleigh or Rayleigh-Gans approximation, while branches, trunks or thick stems are modeled as dielectric cylinders. A tree model is composed of crown and trunk. The crown is a cloud with a simple geometrical shape containing randomly oriented non-spherical particles. The trunk is an upright cylinder with the top covered by the crown. Oblate or disk-like particles and elliptic crown are adopted for broad-leaf forest, while prolate or needle-like particles and cone-like crown are for needle-leaf forest. In addition, the crown shapes take a small random perturbation. Similarly, a farm filed is modeled as a layer of randomly oriented non-spherical particles with perturbed layer depth. If the grid units are small enough, the segment of crown or crops cut by one grid unit can be regarded as a segment fully filled with particles. Differences between the trunk and crown are: (a) the trunk is solid within a grid unit i.e. fractional volume is set to 1; (b) double scattering between the trunk and ground is taken into account.

Scattering of the building is seen as the surface scattering from its wall and roof surfaces and multiple interactions with the ground surface. The integral equation method (IEM) is employed to calculate surface scattering. The building is modeled as a parallelepiped with an isoscales triangular cylinder layered upon it. Due to the orientation of wall surface and roof surface, the incident angle must be transferred into the local coordinates of the surface, as well as the polar basis of wave propagation. Geometrical relationships among the wall surfaces, the roof surfaces, as well as the double and triple scatterings between the wall and ground can be found in Franceschetti et al. (2002). Figure 2 illustrates the building model, its image and shadowing in SAR imagery.

Figure 4.

Image and shadow of a building in SAR imagery.

For ground facet, the tangent vectors along $x$ , $y$ axes on the current grid unit are used to construct the facet plane.

The MPA approach is based on the VRT theory for incoherent scattering power, not coherent summation of scattering fields. However, speckle due to coherent interference is one of the most critical characteristics of SAR imagery, especially for studies of SAR image filtering and optimization. Assuming Gaussian probability distribution of the scattering vector $k_L={\left[S_{hh},S_{hv},S_{vh},S_{vv}\right]}^T$ , random scattering vector is then generated as follows:

$$k_L=k^0\cdot {\overline{C}}^{\frac{1}{2}},k^0{}_i=I_i+j{Q}_i,i=\mathrm{1,..,}d$$ E30

where $I_i,Q_i$ are independent Gaussian random numbers with the zero mean and unit variance. Given the positive semidefinite property of $$\overline{C}$$ , its square root can be obtained by Cholesky decomposition.

3.3. Simulation results

The configurations and parameter settings of the radar and platform in simulation are selected follow the AIRSAR sensor of NASA/JPL.

A virtual comprehensive terrain scene is designed, as shown in Figure 5(a), as according to a true DEM of Guangdong province, south China. It contains different types of forests covering on the hill, crops farmland, ordered or random buildings in urban and suburban regions, roads and rivers.

A simulated scattering map (decibel of normalized scattering power, from -50dB to 0dB, pseudo color: R-HH, G-HH, B-HV.) at L band with 12m resolution (pixel spacing is 5m) is displayed in Figure 5(b). The simulated SAR image at L band is shown in Figure 5(c).

Figures 6(a,b) give the simulated scattering map and SAR image at C band, respectively. In this higher band, vegetation scatterings are significantly increased, as well as attenuation is increased and wave penetrable depth is decreased. It causes that trees become thinner while the shadowing becomes darker.

Figure 5.

Simulated SAR image at L band.

Figure 6.

Scattering and Simulated SAR image at C bands.

At the upper-right corner, the blocks images appear like parallel lines, which reveal the dominance of double scattering in the urban area. However, on the other side below the river, the buildings are oriented nearly 45 degrees to the radar flying direction. As a result, the cross-pol scattering (blue) is stronger, while double scattering is reduced. Timberlands can be classified into two classes: broad leaves in yellow color due to balance between HH and VV scattering, and narrow leaves in green color due to weaker scattering of HH. The narrow leaves trees make weaker attenuation of HH and therefore stronger HH scattering of the shadowed ground, which is the reason why most areas of narrow leaves forest become red. Overlay and shadowing effects caused by mountainous topography are particularly perceptible. Croplands are generally uniform and dense, and appear like patches in the image. Its VV scattering (green) is always stronger than HH. They have distinguishable brightness due to different density, sizes, shapes etc. More detail discussion can be seen in ( Xu and Jin, 2006 )

3.4. Bistatic SAR imaging

Bistatic SAR (BISAR) with separated transmitter and receiver flying on different platforms has become of great interests. Most of BISAR research is focused on the engineering realization and signal processing algorithm with few on land-based bistatic experiments or theoretical modeling of bistatic scattering. The MPA provides a fast and efficient tool for monostatic imaging simulation. It involves the physical scattering process of multiple terrain objects, such as vegetation canopy, buildings and rough ground surfaces (Xu and Jin, 2008a).

Figure 7.

Bistatic MPA

Similar to mono-static MAP, the bistatic-MAP steps are described as follows,

(1) The arrays $${\underset{\_}{E}}^{+},{\underset{\_}{E}}^{-}$$ are initialized as unit matrices, $\underset{\_}{S}$ is as zero matrices. Then, visit the grid units, sequentially, along $+y$ direction.

(2) Perform 3D projections of the current grid unit along incidence and scattering directions to the cells $$p_d,p_u$$ of $${\underset{\_}{E}}^{+},{\underset{\_}{E}}^{-}$$ , respectively.

(3) Determine the mapping position of the current grid unit based on the information of its synthetic aperture and Doppler history, which corresponds to the cell $m$ of the array $\underset{\_}{S}$ .

(4) Obtain or calculate the scattering matrix ${\overline{S}}_0$ and the upward/downward attenuation matrix ${\overline{E}}_0^{\pm }$ of the current grid unit.

(5) Refresh the elements of $\underset{\_}{S}$ as

$$\underset{\_}{S}[m]:=\underset{\_}{S}[m]+{\underset{\_}{E}}^{+}[p_u]\cdot {\overline{S}}_0\cdot {\underset{\_}{E}}^{-}[p_d]$$ E31

Refresh the elements of $${\underset{\_}{E}}^{+},{\underset{\_}{E}}^{-}$$ as

$${\underset{\_}{E}}^{-}[p_d]:={\overline{E}}_0^{-}\cdot {\underset{\_}{E}}^{-}[p_d],{\underset{\_}{E}}^{+}[p_u]:={\underset{\_}{E}}^{+}[p_u]\cdot {\overline{E}}_0^{+}$$ E32

(6) Return to Step (2), and continue to visit the next grid at the same coordinate of $y$ or, if all grids at this coordinate have been visited, step forward to a larger coordinate of $y$ till the whole scene is exhausted.

Figures 7, 8 explain the bistatic imaging process of a tree and a building. One of the differences between bistatic and monostatic cases is the shadowing area projected in both incidence and scattering directions. Additionally, due to the split incidence and scattering directions, double scattering terms of object-ground and ground-object are different. One reflects at the ground and diffusely scatters from the object, the other scatters at the object and then reflects from the ground. These two scattering terms have different paths which give rise to separated double-scattering images.

Figure 8.

Composition of a tree and a house in BISAR image.

A comprehensive virtual scene as shown in Figure 9(a) is designed for simulation and analysis. Simulated BISAR images of the configurations, AT (across track) and TI (translational invariant) are shown in Figures 9(b,c), respectively.

Figure 9.

Simulated BISAR images of AT and TI.

3.5. Unified bistatic polar bases and polarimetric analysis

It is found in the BISAR simulation that in AT BISAR image, bistatic scattering of terrain objects preserves polarimetric characteristics analogous to monostatic case. while the major disparity is on total scattering power.

However, the polarimetric parameters behavior is different in the image of general TI mode. First, $H$ is generally higher, probably due to the fact that large bistatic angle is more likely to reveal the randomness and complexity of the object, and makes the scattering energy more uniformly distributed on different polarizations. Second, almost all terrain surfaces show $\alpha {45}^{\circ }$ . It is found that scattering energy is concentrated on the 4-th component of Pauli scattering vector $$k_P$$ all over the image.

It seems that $\alpha$ might lose its capability to represent the polarimetric characteristics under certain bistatic configuration. Generally speaking, all parameters $\alpha ,\beta ,\gamma$ in monostatic SAR image cannot reflect polarimetric characteristics in the general TI case.

Conventional polarimetric parameters such as $$\alpha ,\beta ,\gamma$$ may largely depend on angular settings of the sensors rather than instinct properties of the target. It becomes inconvenient to interpret the physical meaning of polarimetric parameters when they are severely involved with the bistatic configuration.

As proposed in many studies, inverse problems of bistatic scattering are usually discussed in a coordinates system determined by the bisectrix. We believe it is important to first transform the conventional bistatic polar bases to a new one defined by the bisectrix. The unified bistatic polar bases are defined as

$${{\widehat{h}}^{\prime }}_i=\frac{\widehat{b}\times {\widehat{k}}_i}{\left|\widehat{b}\times {\widehat{k}}_i\right|},{{\widehat{v}}^{\prime }}_i={{\widehat{h}}^{\prime }}_i\times {\widehat{k}}_i,\widehat{b}=\frac{{\widehat{k}}_s-{\widehat{k}}_i}{\left|{\widehat{k}}_s-{\widehat{k}}_i\right|}$$ E33

where $\widehat{b}$ denotes the bisectrix, which is defined as the bisector of the incident and scattered wave vectors in the plane, as shown in Figure 10(a).

Figure 10.

Unified bistatic polar bases and bistatic (TI) pseudo color images coded by NPC after unified bistatic polar bases transform.

The relationship between the unified bistatic polar bases and the conventional polar bases can be described by polar basis rotation, i.e. $$E,\overline{S}$$ defined in $$({\widehat{v}}_i,{\widehat{h}}_i,{\widehat{k}}_i)$$ are re-defined as $$E^{\prime },\overline{S^{\prime }}$$ in $$({{\widehat{v}}^{\prime }}_i,{{\widehat{h}}^{\prime }}_i,{\widehat{k}}_i)$$ . It is written as

$${E^{\prime }}_i={\overline{U}}_i\cdot E_i,\overline{S^{\prime }}={\overline{U}}_s\cdot \overline{S}\cdot {\overline{U}}_i^{\text{T}},{\overline{U}}_i=\left[\begin{array}{cc}{\widehat{v}}_i\cdot {{\widehat{v}}^{\prime }}_i& {\widehat{h}}_i\cdot {\widehat{v^{\prime }}}_i\\ {\widehat{v}}_i\cdot {\widehat{h^{\prime }}}_i& {\widehat{h}}_i\cdot {\widehat{h^{\prime }}}_i\end{array}\right]$$ E34

where ${\overline{U}}_i$ is the rotation matrix. Figure 10(b) gives the bistatic TI image after transfom of Eq. (34). The normalized Pauli components (NPC) are defined as $$P_i=\left|k_{P,i}\right|/\left|k_P\right|,i=\mathrm{1,..,4}$$ . The appearance of the NPC color-coded image after transform agrees with the convention of human vision.

We modify the definition of Cloude’s parameters as

$$k_P=\Vert k_P\Vert \cdot \left[\begin{array}{c}\cos \alpha \cos \gamma \exp (j{\varphi }_1)\\ \sin \alpha \cos \beta \exp (j{\varphi }_2)\\ \sin \alpha \sin \beta \exp (j{\varphi }_3)\\ \cos \alpha \sin \gamma \exp (j{\varphi }_4)\end{array}\right]$$ E35

The redefined $\alpha ,\beta ,\gamma$ are calculated and plotted in Figure 11 for bistatic image of Figure 10(b). It can be seen that the redefined $\alpha ,\beta ,\gamma$ are more capable to represent or classify the polarimetric characteristics of different scattering types.

Redefinition of Cloude’s parameters is only a preliminary attempt of bistatic polarimetric interpretation. It remains open for further study in a systematic way and for verification using both simulation and experiments.

Figure 11.

Redefined $\alpha ,\beta ,\gamma$ after bistatic polar bases transform.

The MPA is also applied to simulation of SAR image of undulated lunar surface (Fa and Jin, 2009). Based on the statistics of the lunar cratered terrain, e.g. population, dimension and shape of craters, the terrain feature of cratered lunar surface is numerically generated. According to inhomogeneous distribution of the lunar surface slope, the triangulated irregular network is employed to make the digital elevation of lunar surface model. The Kirchhoff approximation of surface scattering is then applied to simulation of lunar surface scattering. The SAR image for cratered lunar surface is numerically generated. Making use of the digital elevation and Clementine UVVIS data at Apollo 15 landing site as the ground truth, an SAR image at Apollo 15 landing site is simulated.

4. Terrain Surface Classification using De-Orientation Theory

Classification of complex terrain surfaces using polarimetric SAR imagery is one of most important SAR applications. An unsupervised method based on the entropy $$H$$ and target decomposition parameters has been well developed by Cloude et al., which extracts the target decomposition parameter $$\alpha$$ and entropy $$H$$ from eigen-analysis of coherence matrix and construct an unsupervised classification spectrum on the $$\alpha -H$$ plane.

Scattering of the terrain targets functionally depends on the scatter orientation, shape, dielectric property, and scattering mechanism etc. Scatter targets of complex terrain surfaces are often randomly oriented and cause randomly fluctuating echoes. It is difficult to make classification of randomly oriented and randomly distributed scatter targets. Different scatters with different orientations can make the similar scattering, and vice versa, the same scatters with random orientation can make different scattering to make confused classification.

In this section, a transformation of the target scattering vector to rotate the target along the sight line is derived. Deorientation is introduced to transform the target orientation into such fixed state with minimization of cross-polarization (min-x-pol). And meanwhile, the angle is extracted to indicate the angle deviation of the target orientation from this min-x-pol state.

A set of new parameters $$u,v,w$$ from the target scattering vector is defined to indicate the ratio and phase difference of the two co-polarized (co-pol) backscattering terms, and the significance of cross polarized (x-pol) term. All these parameters as well as the entropy $$H$$ are applied to classification of random surface targets.

Numerical simulations of polarimetric scattering of a single small non-spherical particle, and a layer of random non-spherical particles above a rough surface are studied to show the effectiveness of the parameters $$u,v,\psi ,H$$ and the capability of $$u,v,H$$ for classification of complex terrain surfaces (Xu and Jin, 2005).

As examples, the terrain surface classifications for a SIR-C and an AirSAR images are presented.

4.1. De-orientation and parameterization

From Eq. (1), the scattering vector is usually defined as

$$k_P={[S_{hh}+S_{vv},S_{hh}-S_{vv}\mathrm{,2}{S}_x]}^T/\sqrt{2}\text{or}{k}_L={[S_{hh},\sqrt{2}{S}_x,S_{vv}]}^T$$ E36

where $$S_x\equiv S_{hv}=S_{vh}$$ . Cloude et al. defined the parameterization of $$k_P$$ in terms of the parameters $$\alpha ,\beta$$ etc. as

It yields

$$a={\tan }^{-1}(|S_{vv}|/|S_{hh}|),b=\frac{1}{2}\arg (S_{vv}/S_{hh}),c={\cos }^{-1}(\sqrt{2}|S_x|/|k_L|)$$ E38

Rotating the angle $$\psi$$ of the polarization base along with the sight line, the electrical field vector $$E$$ becomes $$E^{\prime }$$ as

$$E^{\prime }=[R]\cdot E,[R]=\left[\begin{array}{cc}\cos \psi & -\sin \psi \\ \sin \psi & \cos \psi \end{array}\right]$$ E39

In backscattering direction, this rotation makes the scattering vector $$k_P$$ to be $${k^{\prime }}_P$$ , which can be expressed as

$${k^{\prime }}_P=[U]\cdot k_P,[U]=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos 2\psi & \sin 2\psi \\ 0& -\sin 2\psi & \cos 2\psi \end{array}\right]$$ E40

Applying minimization of cross polarization (min-x-pol) to $${k^{\prime }}_P$$ , i.e. let

$$\frac{\partial |{k^{\prime }}_{P\mathrm{,3}}{|}^2}{\partial \psi }=\mathrm{0,}\frac{{\partial }^2|{k^{\prime }}_{P\mathrm{,3}}{|}^2}{\partial {\psi }^2}0$$

it yields the deorientation angle $${\psi }_m$$ is obtained as

$${\psi }_m={\left[\mathrm{sgn}\left\{\cos ({\varphi }_2-{\varphi }_3)\right\}\cdot \frac{\beta }{2}\right]}_{\frac{\pi }{2}}$$ E41

It means that such rotation of the angle $${\psi }_m$$ along the sight line makes $$k_P$$ to the min-x-pol status as deoriented $$k_P^d$$ .

New parameters $$u,v,w$$ are then defined from the deoriented scattering vector of Eqs. (38,40) as:

$$u=\sin c\cos 2a,v=\sin c\sin 2a\cos 2b,w=\cos c$$ E42

In the case of non-deterministic targets, we obtain the uncorrelated scattering vectors, i.e. eigenvectors, through eigen-analysis of coherency matrix. Here the most significant eigenvector i.e. principal eigenvector is considered as the representative scattering vector of non-deterministic targets. Thus, the deorientation of non-deterministic targets is merely conducted on the principal eigenvector.

4.2. Numerical simulation

Based on the model of Figure 1, polarimetric scattering is calculated, and the scattering vector and Mueller matrix are obtained. A model of a layer of random non-spherical particles above a rough surface for scattering from terrain surfaces, as shown in Figure 1, is applied to numerical simulation of $$u,v,H$$ at L band in various cases.

Figures 12(a,b) show numerical relationship between the parameters $$u,v$$ and the particle parameters: Eular angle orientation and particle shape (oblate or prolate spheroids). It can be seen that $$\left|u\right|$$ indicates the non-symmetry of particle’s projection along with the sight line, and the dielectric property of the particle, $$\left|{\epsilon }_s\right|$$ , affects $$\left|u\right|$$ , and $${{\epsilon }^{″}}_s/{{\epsilon }^{\prime }}_s$$ affects $$v$$ .

Figure 12(a) presents the distributions of some typical terrain surfaces on the $$\left|u\right|-H$$ plane, which are obtained from the scattering model of Figure 1. It can be seen that $$H$$ indicates the complexity of layer structures of terrain surfaces. Different dielectric property of the soil land and oceanic surface also makes $$u$$ identifiable on the $$\left|u\right|$$ axis.

Figure 12(b) presents the distribution of scattering terms on the plane $$\left|u\right|-v$$ , where $$M_i$$ denotes the i-th order scattering. Concluded from these figures, $$H$$ indicates the canopy randomness, $$u$$ is useful for distinguishing different terrain targets and $$v$$ is helpful to take account of scattering mechanisms of different-orders.

Figure 12.

Distributions of typical surfaces and different scattering mechanism.

4.3. Surface classification

The flow-chart of the classification method of PolSAR Data based on deorientation and new parameters $$u,v,\psi$$ is shown in Figure 13. The classification decision tree of parameters $$u,v,H$$ is displayed on Figure 14.

Figure 13.

Flow-chart of the surface classification.

Figure 14.

Classification decision tree of parameters $$u,v,H$$ .

An AirSAR data at L band of the Boreal area in Canada with rich resources of vegetation is chosen for classification and orientation-analysis, as shown in Figure 15(a). Figure 15(b) is the deorientation classification over this area. Terrain surfaces are divided into 8 classes following the decision-tree: Timber, Urban 2, Urban 1, Canopy 2, Canopy 1, Surface 3, Surface 2 and Surface 1.

Figure 15.

A SIR-C SAR image and surface classification.

Figure 16 shows the data distributions on the planes $$\left|u\right|-H$$ for parameter $$v0.2$$ and $$v-0.2$$ , and corresponding classification.

Figure 16.

Classification by $$u,v,H$$ for 19 classes.

The orientation distribution of several classes are selected to be displayed in Figures 17(a,b). It can be seen that

1. While in the Urban 1 (sparse forest of vertical trunks, instead of artificial constructions), there are uniformly vertical orientations indicating orderly trunks.

2. Random orientation in the Canopy 1 means that the vegetation canopy in this area might be the disordered bush. Note that the region with the roads inside the forest might show randomness confused by the bush vegetation on the roadsides.

Following the above orientation analysis, the terrain surfaces are further classified into the subclasses and are identified by their types.

Figure 17.

Classification of AirSAR data, Boreal area.

4.4. Faraday rotation and Surface classification

To obtain the moisture profiles in subcanopy and subsurface, it needs UHF/VHF bands (435MHz, 135MHz) to penetrate through the dense subcanopy (~20kg/m² and more) with little scattering and reach subsurface. However, when the wave passes through the anisotropic ionosphere and action of geomagnetic field, the polarization vectors of electromagnetic wave are rotated, which is called the Faraday rotation (FR). The FR depends on the wavelength, electronic density, geomagnetic field, the angle between the direction of wave propagation and geomagnetic field, and the wave incidence angle.

Assuming that the propagation direction is not changed passing through homogeneous ionosphere, $ds=\sec {\theta }_{i}dz$ (where $$\widehat{z}$$ is the normal to the surface), and geomagnetic field keeps constant as one at 400 km altitude, FR is simply written as

$${\Omega }_F\approx -2620{\rho }_{e}B(400){\lambda }^2\cos {\Theta }_B\sec {\theta }_i\left(\text{radians}\right)$$ E43

where ${\Theta }_B$ is the angle between the directions of electromagnetic wave propagation and geomagnetic field, ${\rho }_e$ is the total electron content per unit area ( ${10}^{16}$ electron/m²), and $B(400)$ is the intensity of geomagnetic field (T) at 400 km altitude.

The FR may decrease the difference between co-pol backscattering, enhance cross-polarized echoes, and mix different polarized terms. Thus, the satellite-borne SAR data at low frequency becomes distorted due to FR effect. Since the FR is proportional to the square power of the wavelength, it yields especially serious impact on the SAR observation operating at the frequency lower than L band. The FR angle at P band can reach dozens of degrees.

As a polarized electromagnetic wave passes through the ionosphere, the scattering matrix $${\overline{S}}^F$$ with FR (indicated by superscript F ) is written by the scattering matrix $$\overline{S}$$ without FR as follows

$$\begin{array}{l}\left[\begin{array}{cc}{S}_{hh}^F& S_{hv}^F\\ S_{vh}^F& S_{vv}^F\end{array}\right]=\left[\begin{array}{cc}\cos \Omega & \sin \Omega \\ -\sin \Omega & \cos \Omega \end{array}\right]\left[\begin{array}{cc}{S}_{hh}& S_{hv}\\ S_{hv}& S_{vv}\end{array}\right]\left[\begin{array}{cc}\cos \Omega & \sin \Omega \\ -\sin \Omega & \cos \Omega \end{array}\right]\\ S_{hh}^F=S_{hh}{\cos }^2\Omega -S_{vv}{\sin }^2\Omega ,S_{vv}^F=-S_{hh}{\sin }^2\Omega +S_{vv}{\cos }^2\Omega \\ S_{hv}^F=S_{hv}+(S_{hh}+S_{vv})\sin \Omega \cos \Omega ,S_{vh}^F=S_{hv}-(S_{hh}+S_{vv})\sin \Omega \cos \Omega \end{array}$$ E44

The measured polarimetric data with FR, $${\overline{S}}^F$$ or $${\overline{M}}^F$$ , are distorted and cannot be directly applied to terrain surface classification. However, fully polarimetric 4×4-D $$\overline{M}$$ without FR can be inverted from $${\overline{M}}^F$$ . And it shows that the remained $$\pm \pi /2$$ ambiguity error does not affect the classification parameters: $$u,v,H,\alpha ,A$$ . Based on intuitive assumption of gradual change of FR degree with geographical position, a method of 2D phase unwrapping method with random benchmark can be employed to eliminate the $$\pm \pi /2$$ ambiguity error is designed. Some example shows the fully polarimetric SAR without FR at low frequency, such as P band, can be fully inverted (Qi and Jin, 2007).

5. Change Detection of Terrain Surface from Multi-Temporal SAR Images

Multi-temporal observations of SAR remote sensing imagery provide fast and practical technical means for surveying and assessing such vast changes. One direct application is to detect and classify the information on changes in the terrain surfaces. It would be laborious to make an intuitive assessment for a huge amount of multi-temporal image data over a vast area. Such assessment based on qualitative gray-level analysis is not accurate and might lose some important information. How to detect and automatically analyze information on change in the terrain surfaces is a key issue in remote sensing.

In this section, two-thresholds EM and MRF algorithm (2EM-MRF) is developed to detect the change direction of backscattering enhanced, reduced and unchanged regimes from the SAR difference image (Jin and Luo, 2004, Jin and Wang, 2009).

On May 12, 2008, a major earthquake, measuring 8.0 on the Richter scale, jolted southwestern China’s Sichuan Province, Wenchuan area. To evaluate the damages and terrain surface changes caused by the earthquake, multi-temporal ALOS PALSAR image data before and after the earthquake are applied to detection and classification of the terrain surface changes. Using the tools of Google Earth for surface mapping, the surface change situation after the earthquake overlapped the DEM topography can be demonstrated in multi-azimuth views as animated cartoon. The detection and classification are also compared with the optical photographs. It is proposed that multi-thresholds EM and MRF analysis may become traceable when multi-polarization, multi-channels, multi-sensors multi-temporal image data become available.

5.1. Two thresholds expectation maximum algorithm

Consider two co-registered images $$X_1$$ and $$X_2$$ with size $$I\times J$$ at two different times $$(t_1,t_2)$$ . Their difference image is $$X=(X_2-X_1)$$ , and denote $$X_D=\{X(i,j),1\le i\le I\mathrm{,1}\le j\le J\}$$ (in some cases, $$X=X_2/X_1$$ also can be used).

As usually, one-threshold expectation maximum (EM) algorithm has been employed to classify two opposite changes: the unchanged class $${\omega }_n$$ and the changed class $${\omega }_c$$ . The probability density function $$P(X),X\in X_D$$ was modeled as a mixture distribution of two density components associated with $${\omega }_n$$ and $${\omega }_c$$ , i.e.

$$p(X)=p(X|{\omega }_n)P({\omega }_n)+p(X|{\omega }_c)P({\omega }_c)$$ E45

Assume that both the conditional probabilities $$p(X|{\omega }_n)$$ and $$p(X|{\omega }_c)$$ are modeled by Gaussian distributions. The iterative equations for estimating the statistical parameters and the a priori probability for the class $${\omega }_n$$ are the following

$$\begin{array}{l}{P}^{t+1}({\omega }_n)=\frac{{\displaystyle \sum _{X(i,j)\in X_D}\frac{P^t({\omega }_n)p^t(X(i,j)|{\omega }_n)}{p^t(X(i,j))}}}{IJ}\\ \\ \text{mean}{\mu }_n^{t+1}=\frac{{\displaystyle \sum _{X(i,j)\in X_D}\frac{P^t({\omega }_n)p^t(X(i,j)|{\omega }_n)}{p^t(X(i,j))}X(i,j)}}{{\displaystyle \sum _{X(i,j)\in X_D}\frac{P^t({\omega }_n)p^t(X(i,j)|{\omega }_n)}{p^t(X(i,j))}}}\\ \\ \text{variance}=\frac{{\displaystyle \sum _{X(i,j)\in X_D}\frac{P^t({\omega }_n)p^t(X(i,j)|{\omega }_n)}{p^t(X(i,j)}{[X(i,j)-{\mu }_n^t]}^2}}{{\displaystyle \sum _{X(i,j)\in X_D}\frac{P^t({\omega }_n)p^t(X(i,j)|{\omega }_n)}{p^t(X(i,j))}}}\end{array}$$ E46

where the superscripts t and t +1 denote the current and next iterations.

Analogously, these equations can also be used to estimate $$p({\omega }_c)$$ with $${\mu }_c$$ and $${\sigma }_c^2$$ .

The estimates are computed starting from the initial values by iterating the above equations until convergence. The initial value of the estimated parameters can be obtained by the analysis of the histogram of the difference image. A pixel subset $$S_n$$ likely to belong to $${\omega }_n$$ and a pixel subset $$S_c$$ likely to belong to $${\omega }_c$$ can be obtained by preliminarily selecting two-threshold $$T_n$$ and $$T_c$$ on the histogram. This is equivalent to solving the ML boundary $$T_o$$ for the two classes $${\omega }_n$$ and $${\omega }_c$$ on the difference image. The optional threshold $$T_o$$ requires

$$\frac{P({\omega }_c)}{P({\omega }_n)}=\frac{p(X|{\omega }_n)}{p(X|{\omega }_c)}$$ E47

Under the assumption of Gaussian distribution, Eq. (3) yields

$$({\sigma }_n^2-{\sigma }_c^2)T_o^2+2({\mu }_n{\sigma }_c^2-{\mu }_c{\sigma }_n^2)T_o+{\mu }_c^2{\sigma }_n^2-{\mu }_n^2{\sigma }_c^2+2{\sigma }_n^2{\sigma }_c^2\ln \left[\frac{{\sigma }_{c}P({\omega }_n)}{{\sigma }_{n}P({\omega }_c)}\right]=0$$ E48

to obtain optimal $$T_o$$ .

Figures 18 (a,b) are the ALOS PALSAR images (L band, HH polarization) in February 17 and May 19, 2008 before and after earthquake in Beichuan County, Wenchuan area, respectively. Spatial resolution is 4.7m4.5m. Figure 18(c) is the difference image of Figures 18(a,b), which seems very difficult to evaluate the terrain surface changes only using man’s vision, especially to accurately classify the change classes in large scale.

Figure 18.

ALOS PALSAR image in in Beichuan.

The pixels of the difference image (i.e. $$X_D={\sigma }_D^0={\sigma }_2^0-{\sigma }_1^0$$ ) are classified into three classes: $${\omega }_{c1}$$ of $${\sigma }_D^0$$ -enhanced, $${\omega }_n$$ of $${\sigma }_D^0$$ -unchanged and $${\omega }_{c2}$$ of $${\sigma }_D^0$$ -reduced. Thus, the probability density function $$P(X)$$ is modeled as a mixture density distribution consisting of three components:

$$p(X)=p(X|{\omega }_{c1})P({\omega }_{c1})+p(X|{\omega }_n)P({\omega }_n)+p(X|{\omega }_{c2})P({\omega }_{c2})$$ E49

The parameter $$T_{o1}$$ is first obtained by application of the EM algorithm to the enhanced class $${\omega }_{c1}$$ and no-enhanced class $${\omega }_{n1}$$ ( $${\omega }_n$$ and $${\omega }_{c2}$$ ). Then, we obtain $$T_{o2}$$ by applying the EM to the reduced class $${\omega }_{c2}$$ and no-reduced class $${\omega }_{n2}$$ ( $${\omega }_n$$ and $${\omega }_{c1}$$ ), where $${\omega }_{n1}={\omega }_n\cup {\omega }_{c2}$$ , $${\omega }_{n2}={\omega }_n\cup {\omega }_{c1}$$ . Finally, the two results are superposed to get the final classification.

Figure 19.

The histogram of the difference image and change detection using 2EM method.

As an example, Figure 19(a) shows the grayscale histogram of a small part chosen from the difference image, Figure 18(c). Three classes of the probability distributions are outlined in this figure. The histogram is normalized in accordance with the probability density distributions.

The initial statistical parameters of the subsets of

$$\begin{array}{l}{S}_{n1}=\{X(i,j)|0X(i,j)T_{n1}\},S_{c1}=\{X(i,j)|X(i,j)T_{c1}\\ \\ \text{are first derived},\text{ and then derive the parameters of}\\ \\ S_{n2}=\{X(i,j)|T_{n2}X(i,j)0\},S_{c2}=\{X(i,j)|X(i,j)T_{c2}\}\end{array}$$ E50

The initial statistical parameters are related to the classes $${\omega }_{n1}$$ , $${\omega }_{c1}$$ and the classes $${\omega }_{n2}$$ , $${\omega }_{c2}$$ , respectively. Then, Eqs. (46a-c) are sequentially used to perform the iterations on the four classes, i.e. the a priori probability and other statistical parameters [ $$P({\omega }_{n1})$$ , $${\mu }_{n1}$$ , $${\sigma }_{n1}^2$$ ], [ $$P({\omega }_{c1})$$ , $${\mu }_{c1}$$ , $${\sigma }_{c1}^2$$ ], [ $$P({\omega }_{n2})$$ , $${\mu }_{n2}$$ , $${\sigma }_{n2}^2$$ ] and [ $$P({\omega }_{c2})$$ , $${\mu }_{c2}$$ , $${\sigma }_{c2}^2$$ ]. Solving Eq. (48), the thresholds $$T_{o1}=3.1643$$ dB and $$T_{o2}=-2.8381$$ dB are obtained.

Figure 19(b) is the change detection using EM algorithm, where red color indicates the area with scattering enhanced, blue color indicates the area with scattering reduced, and green color denotes no-changed. It makes three classes change of the difference image.

5.2. The change detection using 2EM-MRF

Actually, the pixels of the image are spatially correlated, i.e. a pixel belonging to the class $${\omega }_k$$ is likely to be surrounded by pixels belonging to the same class. To take account of spatial texture may yield more reliable and accurate change detection.

Let the set $$C=\left\{C_{\ell }\mathrm{,1}\le \ell \le L\right\}$$ represent the possible sets of the labels in the difference image:

$$C_{\ell }=\left\{C_{\ell }(i,j)\mathrm{,1}\le i\le I\mathrm{,1}\le j\le J\right\},C_{\ell }(i,j)\in \left\{{\omega }_{c1},{\omega }_{c2},{\omega }_n\right\}$$ E51

where $$L=3^{IJ}$$ . According to Bayes rule for minimum error, the classification result should maximize the posterior conditional probability,

$$C_k=\underset{C_{\ell }\in C}{\arg \max }\left\{P(C_{\ell }|X_D)\right\}=\underset{C_{\ell }\in C}{\arg \max }\left\{P(C_{\ell })p(X_D|C_{\ell })\right\}$$ E52

where $$P(C_{\ell })$$ is the prior model for the class labels, and $$p(X_D|C_{\ell })$$ is the joint density function of the pixel values in the difference image $$X_D$$ given the set of labels $$C_{\ell }$$ .

The MRF algorithm is employed with a spatial neighborhood 55 pixels system to take account of spatial texuture. The generation of the final change-detection map involves the labeling of all the pixels in the difference image so that the posterior probability of Eq. (52) is maximized.

The MRF algorithm is equivalent to the minimization of the Gibbs energy function. The MRF algorithm is iteratively carried out.

Figure 20(a) presents the final result of 2EM-MRF for detection and classification of the terrain surface changes from the difference image, Figure 18(c). The numbers 1~8 indicate some areas with typical changes.

It can be seen that, for example, in the area 1 the river was largely blocked up with landslide; in the area 2, there were landslides causing large scale blocks; the area 3 is Beichuan town, where the terrain surface was significantly undulated or roughed due to landslide and building collapse; in the areas 4 and 5 the river was blocked up due to landslides along river lines; in the area 6 the highway was significantly blocked up with landslide; in the area 7 there was collective landslides almost towards one azimuth direction, and the flat area with reduced scattering along the river might be due to seasonal risen water in May than February; in the area 8 collapse of mountain blocks caused the terrain surface undulating to enhance stronger scattering.

Figure 20(b) gives the topography and contour lines of the same area from the google website. It is useful to locate where and which kind of the changes are happening. It might be seen that the landslides, e.g. in the areas 1 and 7, are correlated with direction and magnitude of slopes. A quantitative analysis to assess the landslides and surface damage can be further developed.

Figure 20.

Automatic detection and classification of terrain surface changes using 2EM-MRF.

As an example, Figure 20(c) gives a comparison of the change at the region 3 with a photograph distributed to public from the website of the Ministry of the State Resources of China. It can be clarified in both figures that A shows lanslides and makes the surface smooth, and B shows the highway blocked. It can be seen that optical shadowing actually does not confuse these classifications.

Using the tool of Google Earth mapping with the 2EM-MRF results, the terrain surface change situation classified by three types overlapped the DEM topography can be showed in multi-azimuth views as a animated cartoons. Since all process of 2EM-MRF are automatic and carried out on real time, it should be helpful, especially for commanding the rescue works in disaster scene. Figure 21 shows four views.

It can be seen that information retrieval from a single amplitude image (scattered power) with one-frequency, mono-polarization at one time is very limited. In fully polarimetric technology, for example, using target decomposition theory, deorientation transform etc., the physical status of the terrain surfaces such as vegetation canopy, road, building, river etc. can be well classified and would be of greatly helpful to change detection. Four Stokes parameters of polarimetric measurement can be feasible to show the surface slope variation and anisotropy, and INSAR has been well applied to retrieval of surface digital elevation. All of these progress show superiority over a single amplitude image analysis and manual vision estimation. It is also helpful to fuse CFAR (constant false alarm rate) for detection of the object, e.g. edge and block, from a SAR image.

As multi-polarized, multi-channels and multi-temporal image data become available, 2EM-MRF can be further extended to utilization of principal characteristic parameters of the difference image, such as VV+HH, VV-HH, entropy, angular $$\alpha ,\beta ,\gamma ,\psi$$ for change detection and classification. More information about terrain surface changes can be retrieved.

Figure 21.

Multi-azimuth views in animated cartoon showing terrain surface changes on Google Earth mapping.

6. DEM Inversion from a Single POL-SAR Image

Co-polarised or cross-polarised backscattering signature is the function of the incidence wave with the ellipticity angle $$\chi$$ and orientation angle $$\psi$$ . Recently, polarimetric INSAR image data has been utilized to generate digital surface elevation and to invert terrain topography. When the terrain surface is flat, polarimetric scattering signature has the maximum largely at the orientation angle $$\psi =0$$ . However, it has been shown that as the surface is tilted, the orienattion angle $$\psi$$ at the maximum of co-polarized (co-pol) or cross-polarized (cross-pol) signature can shift from $$\psi =0$$ . This shift can be applied to convert the surface slopes. Making use of the assumption of real co-pol and zero cross-pol scattering amplitude functions, the $$\psi$$ shift is expressed by the real scattering amplitude functions. Since both the range and azimuth angles are coupled, two- or multi-pass SAR image data are required for solving two unknowns of the surface slopes. This approach has been well demonstrated for inversion of digital elevation mapping (DEM) and terrain topography by using airborne SAR data.

However, scattering signature is an ensemble average of echo power from random scatter media. Measurable Stokes parameters as the polarized scattering intensity should be directly related to the $$\psi$$ shift. In this section, using the Mueller matrix solution, the $$\psi$$ shift is newly derived as a function of three Stokes parameters, $$I_{vs},I_{hs},U_s$$ , which are measurable by the polarimetric SAR imagery.

Using the Euler angles transformation between the principal and local coordinates, the orientation angle $$\psi$$ is related with both the range and azimuth angles, $$\beta$$ and $$\gamma$$ , of the tilted surface pixel and radar viewing geometry. These results are consistent with Lee et al. (1998) and Schuler (2000), but are more general.

It is proposed that the linear texture of tilted surface alignment is used to specify the azimuth angle $$\gamma$$ . The adaptive threshold method and image morphological thinning algorithm are applied to determine the azimuth angle $$\gamma$$ from image linear textures. Thus, the range angle $$\beta$$ is then solved, and both $$\beta$$ and $$\gamma$$ are utilized to obtain the azimuth slope and range slope. Then, the full multi-grid algorithm is employed to solve the Poisson equation of DEM and produce the terrain topography from a single pass Polarimetric SAR image (Jin, 2005; Jin and Luo, 2003 ).

6.1. The Ψ shift as a function of the Stokes parameters

Consider a polarized electromagnetic wave incident upon the terrain surface at $$(\pi -{\theta }_i,{\phi }_i)$$ . Incidence polarization is defined by the ellipticity angle $$\chi$$ and orientation angle $$\psi$$ ( Jin 1994 ). The 22-D complex scattering amplitude functions are obtained from polarimetric measurement as

$${\overline{I}}_i(\chi ,\psi )={[I_{vi},I_{hi},U_i,V_i]}^T={\left[\frac{1}{2}(1-\cos 2\chi \cos 2\psi ),\frac{1}{2}(1+\cos 2\chi \cos 2\psi ),-\cos 2\chi \sin 2\psi ,\sin 2\chi \right]}^T$$ E53

$${\sigma }_c=4\pi \cos {\theta }_i{P}_n$$ E54

where

$$P_n=0.5\left[I_{vs}(1-\cos 2\chi \cos 2\psi )+I_{hs}(1+\cos 2\chi \cos 2\psi )+U_s\cos 2\chi \sin 2\psi +V_s\sin 2\chi \right]$$ E55

When the terrain surface is flat, co-pol backscattering $${\sigma }_c$$ versus the incidence polarization $$(\chi ,\psi )$$ has the maximum at $$\psi =0$$ . However, it has been shown that as the surface is tilted, the orientation angle $$\psi$$ for the $${\sigma }_c$$ maximum is shifted from $$\psi =0$$ .

Let $$\partial P_n/\partial \psi =0$$ at the maximum $${\sigma }_c$$ and $$\chi =0$$ of symmetric case, it yields

$$0=-(I_{hs}-I_{vs})\sin 2\psi +U_s\cos 2\psi +0.5(I_{hs}+I_{vs})\text{'}+0.5{U^{\prime }}_s\sin 2\psi +0.5(I_{hs}-I_{vs})\text{'}\cos 2\psi$$ E56

where the superscript of prime denote $$\partial /\partial \psi$$ . It can be shown that

$$\begin{array}{c}0.5(I_{hs}+I_{vs})\text{'}~(M_{13}+M_{23})\cos 2\psi ~(\mathrm{Re}{S}_{vv}{S}_{vh}^{*}+\mathrm{Re}{S}_{hv}{S}_{hh}^{*})\cos 2\psi \\ 0.5{U^{\prime }}_s~M_{33}\cos 2\psi ~\mathrm{Re}{S}_{vv}{S}_{hh}^{*}+S_{vh}{S}_{hv}^{*}\cos 2\psi \\ (I_{hs}-I_{vs})~0.5(M_{11}+M_{22})\cos 2\psi ~0.5({\left|S_{vv}\right|}^2+{\left|S_{hh}\right|}^2)\cos 2\psi \end{array}$$ E57

$$\begin{array}{c}0.5(I_{hs}-I_{vs})\text{'}~0.5(M_{13}-M_{23})\cos 2\psi ~0.5(\mathrm{Re}{S}_{vv}{S}_{vh}^{*}-\mathrm{Re}{S}_{hv}{S}_{hh}^{*})\cos 2\psi \\ U_s~0.5(M_{31}+{\text{M}}_{\text{32}}\text{)~Re}{S}_{vv}{S}_{hv}^{*}+\mathrm{Re}{S}_{vh}{S}_{hh}^{*}\end{array}$$ E58

It can be seen that (57a) and (57b) are much less than (57c), and (58a) is much less than (58b), so the last three terms on RHS of Eq. (56) are now neglected. Thus, it yields the $$\psi$$ shift at the $${\sigma }_c$$ maximum expressed by the Stokes parameters as follows

$$\tan 2\psi =\frac{U_s}{I_{hs}-I_{vs}}$$ E59

It can be seen that the third Stokes parameter $$U_s\ne 0$$ does cause the $$\psi$$ shift.

By the way, if both $$U_s$$ and $$I_{hs}-I_{vs}$$ approach zero, e.g. scattering from uniformly oriented scatterers or isotropic scatter media such as thick vegetation canopy, $$\psi$$ cannot be well defined by Eq. (59).

6.2. The range and azimuth slopes and DEM inversion

The polarization vectors $${\widehat{h}}_i,{\widehat{v}}_i$$ of the incident wave at $$(\pi -{\theta }_i,{\varphi }_i)$$ in the principal coordinate ( $$\widehat{x},\widehat{y},\widehat{z}$$ ) are defined as ( Jin, 1994 )

$$\begin{array}{c}{\widehat{h}}_i=\frac{\widehat{z}\times {\widehat{k}}_i}{|\widehat{z}\times {\widehat{k}}_i|}=-\sin {\varphi }_i\widehat{x}+\cos {\varphi }_i\widehat{y}\text{and}{\widehat{v}}_i={\widehat{h}}_i\times {\widehat{k}}_i\\ \\ \text{where the incident wave vector}\\ \\ {\widehat{k}}_i=\sin {\theta }_i\cos {\varphi }_i\widehat{x}+\sin {\theta }_i\sin {\varphi }_i\widehat{y}-\cos {\theta }_i\widehat{z}\end{array}$$ E60

As the pixel surface is tilted with local slope, the polarization vectors should be re-defined following the local normal vector $${\widehat{z}}_b$$ as follows

$${\widehat{h}}_b=\frac{{\widehat{z}}_b\times {\widehat{k}}_{ib}}{|{\widehat{z}}_b\times {\widehat{k}}_{ib}|}\text{and}{\widehat{v}}_b={\widehat{h}}_b\times {\widehat{k}}_{ib}$$ E61

By using the transformation of the Euler angles $$(\beta ,\gamma )$$ between two coordinates ( $$\widehat{x},\widehat{y},\widehat{z}$$ ) and $$({\widehat{x}}_b,{\widehat{y}}_b,{\widehat{z}}_b)$$ ( Jin, 1994 ), it has

$$\begin{array}{c}\widehat{x}=\cos \gamma \cos \beta {\widehat{x}}_b+\sin \gamma {\widehat{y}}_b-\cos \gamma \sin \beta {\widehat{z}}_b\\ \widehat{y}=-\sin \gamma \cos \beta {\widehat{x}}_b+\cos \gamma {\widehat{y}}_b+\sin \gamma \sin \beta {\widehat{z}}_b\\ \widehat{z}=\sin \beta {\widehat{x}}_b+\cos \beta {\widehat{z}}_b\end{array}$$ E62

Substituting above relations into Eqs. (60,61), it yields

$${\widehat{h}}_i=-\cos \beta \sin (\gamma +{\varphi }_i){\widehat{x}}_b+\cos (\gamma +{\varphi }_i){\widehat{y}}_b+\sin \beta \sin (\gamma +\varphi ){\widehat{z}}_b$$ E63

$$\begin{array}{l}{\widehat{k}}_{ib}=(\cos \beta \sin {\theta }_i\cos (\gamma +{\varphi }_i)-\sin \beta \cos {\theta }_i){\widehat{x}}_b+\sin {\theta }_i\sin (\gamma +{\varphi }_i){\widehat{y}}_b\\ -[\sin \beta \sin {\theta }_i\cos (\gamma +{\varphi }_i)+\cos \beta \cos {\theta }_i]{\widehat{z}}_b\end{array}$$ E64

Using Eq. (64) to Eq. (61), the polarization vector $${\widehat{h}}_b$$ for local surface pixel is written as

$${\widehat{h}}_b=\frac{a{\widehat{x}}_b+b{\widehat{y}}_b}{\sqrt{a^2+b^2}}$$ E65

where $$a=-\sin {\theta }_i\sin (\gamma +{\varphi }_i)$$ , $$b=\sin {\theta }_i\cos \beta \cos (\gamma +{\varphi }_i)-\sin \beta \cos {\theta }_i$$ . Thus, Eqs. (63,65) yield the orientation angle as

$$\begin{array}{c}\cos \psi ={\widehat{h}}_b\cdot {\widehat{h}}_i=\frac{\cos \beta \sin {\theta }_i-\cos (\gamma +{\varphi }_i)\sin \beta \cos {\theta }_i}{\sqrt{a^2+b^2}}\\ \\ \tan \psi =\frac{\tan \beta \sin (\gamma +{\varphi }_i)}{\sin {\theta }_i-\cos {\theta }_i\tan \beta \cos (\gamma +{\varphi }_i)}\end{array}$$ E66

Thus, the range and azimuth slopes of the pixel surface can be obtained as

$$S_R=\tan \beta \cos \gamma ,S_A=\tan \beta \sin \gamma$$ E67

Since a single $$\psi$$ shift cannot simultaneously determine two unknowns of $$\beta$$ and $$\gamma$$ in Eq. (64), two- or multi-pass SAR image data are usually needed.

If only the single-pass SAR image data are available, one of two unknown angles, $$\beta$$ or $$\gamma$$ , should be first determined. The azimuth alignment of tilted surface pixels can be visualized as a good indicator of the azimuth direction. We apply the adaptive threshold method and image morphological thinning algorithm to specify the azimuth angle in all SAR image pixels. The algorithm contains the following steps:

(1) Make speckle filtering over the entire image.

(2) Apply the adaptive threshold method to produce a binary image. The global threshold value is not adopted because of the heterogeneity of the image pixels.

(3) Apply the image morphological processing for the binary image, remove those isolated pixels and fill small holes. Referring to the part of binary’s “1” as the foreground and the part of binary’s “0” as the background, the edges from the foreground are extracted.

(4) Each pixel on the edge is set as the center of a square 2121 window, and a curve segment through the centered pixel is then obtained. Then, applying the polynomial algorithm for fitting curve segment in the least-squares sense, the tangential slope of the centered pixel is obtained. It yields the azimuth angle of the centered pixel. Make a mark on that pixel so that it won’t be calculated in the next turn.

(5) Removing the edge in Step (4) from the foreground, a new foreground is formed. Repeat Step (4) until the azimuth angle of every pixel in the initial foreground is determined.

(6) Make the complementary binary image, i.e. the initial background now becomes the foreground. Then, the Steps (4) and (5) are repeated to this complementary image until the azimuth angle of every pixel in the initial background is determined.

This approach provides a supplementary information to firstly determine the angle $$\gamma$$ over whole image area if there is no other information available.

Figure 22.

Morphological thinning algorithm to determine the azimuth angle.

As an example for DEM inversion, the L-band polarimetric SAR data is shown in Figure 23(a). As the azimuth angle $$\gamma$$ of each pixel is obtained by the adaptive threshold method and thinning algorithm as described in the above steps, the $$\beta$$ angle of each pixel can be determined by Eq. (66b). Taking $${\varphi }_i=0$$ , it yields

$$\tan \beta =\frac{\tan \psi \sin {\theta }_i}{\sin \gamma +\tan \psi \cos \theta {}_{i}cos\gamma }$$ E68

where the orientation angle $$\psi$$ is calculated using Eq. (59), while the incident angle $${\theta }_i$$ is determined by the SAR view geometry.

Substituting $$\beta$$ and $$\gamma$$ into Eq. (67), the azimuth slope $$S_A$$ and range slope $$S_R$$ are obtained. We utilize the two slopes to invert the terrain topography and DEM. The slopes $$S_A$$ and $$S_R$$ of all pixels in the SAR image are calculated as shown in Figures 23(b,c).

Figure 23.

SAR image and azimuth/range slopes.

The DEM can be generated by solving the Poisson equation for a $$M\times N$$ rectangular grid area. The Poisson equation can be written as

$${\nabla }^2\varphi (x,y)=\rho (x,y)$$ E69

where $${\nabla }^2$$ is the Laplace operator. The source function $$\rho (x,y)$$ consists of the surface curvature calculated by the slopes $$S(x,y)$$ , where the $$\widehat{x},\widehat{y}$$ -directions are used as the range and azimuth directions, respectively. In this section, the full multi-grid (FMG) algorithm is employed to solve the Poisson equation. The benefits of FMG are due to its rapid-convergence, robustness and low computation load. The inverted DEM of SAR image in Figure 23(a) is given in Figures 24(a,b), and its contour map, Figure 24(b) has been validated in comparison to the map of China Integrative Atlas (1990).

Figure 24.

A DEM inversion.

7. Reconstruction of Buildings Objects from Multi-Aspect SAR Images

Reconstruction of 3D-objects from SAR images has become a key issue for information retrieval for SAR monitoring. 3D reconstruction of man-made object is usually based on interferometeric technique or fusion with other data resources, e.g. optical and GIS data. It has been noted that scattering from man-made objects produces bright spots in sub-metric resolution, but present strips/blocks image in metric resolution. This difference is largely attributed to the different imaging ways employed for different resolutions, for example, spotlight for sub-metric resolution, stripmap for metric resolution.

An automatic method for detection and reconstruction of 3D objects from multi-aspect metric-resolution SAR images. The steps are as follows:

The linear profile of the building objects is regarded as the most prominent characteristics. The POL-CFAR detector, Hough transform, and corresponding image processing techniques are applied to detection of parallelogram-like façade-images of building objects. Results of simulated SAR images as well as real 4-aspect Pi-SAR images are given subsequently after each step.

A probabilistic description of the detected façade- images is presented. Maximum-likelihood estimation (ML) of building objects from multi-aspect observed façade-images is given. Eventually, in association with a hybrid priority rule of inversion reliability, an automatic algorithm is designed to match multi-aspect façade- images and reconstruct the building objects. Moreover, taking advantage of the multi-aspect coherence of building-images, a new iterative co-registration method is presented.

Finally, reconstruction results are presented, and good performance is evaluated comparing with ground truth data. It is also concluded that the reconstruction accuracy closely depends on the number of available aspects. In the end, a practicable scheme of the 3D-object reconstruction from satellite-borne SAR image is proposed (Xu and Jin, 2007).

7.1. Object image detection

As shown in Figure 25, the scattering image of a simple building is produced by direct scattering from the wall/façade, roof and edges, double scattering from wall-ground and shadows projected etc. In metric resolution, the scattering produces bright strips/blocks from respective parts of the object. A key step of reconstruction is to identify and extract these strips/blocks.

For the case of a smooth wall/façade, the only double scattering term to be considered must follow a specific propagation path, i.e. wall (reflect) to ground (diffuse). Simple building object is taken into account and is modeled as a cuboid, and the spatial distribution of the building objects is assumed not crowded, i.e. without serious shadowing and superposition. The cuboid object is described by seven geometric parameters: 3D position, 3D size and orientation. Besides, the flat roof is assumed with much less scattering compared with the edges and façades.

Figure 26 shows (a) a model of the cuboid object, (b) its simulated SAR image, (b) using the mapping and projection algorithm, (c) a photo of real rectangular building, and (d) its image of Pi-SAR observation, respectively.

It can be seen from Figure 26 that the longer wall of the cuboid-like building, called as major wall henceforth, plays the dominant role in a SAR image. In this section, main attention is focused on the major wall image. At the first step, the edge detectors of constant false alarm rate (CFAR) such as ratio gradient are used for edge detection.

Figure 25.

SAR imaging of a simple building model and its image composition.

To our experience, the POL-CFAR detector derived from complex Wishart distribution can fulfill the segmentation requirements of medium- and small-scale building-images. Besides, the ridge filter applied after this step can accommodate segmentation error to some extent.

Figure 26.

A cuboid object in m-resolution SAR images.

In order to improve the segmentation precision, the detected edge by window-operator needs to be further thinned. Using the edge intensity produced in the POL-CFAR detection, an edge thinning method of ‘ridge filtering’ is presented. Taking an 8-neighbor 3×3 window as an example, the center pixel is regarded as on the ridge if its intensity is higher than the pixels along two sides.

4-aspect Pi-SAR images acquired over Sendai, Japan by a square-loop flying path (Flight No. 8708, 8709, 8710, 8711, X-band, pixel resolution 1.25m) are taken as an example of real SAR image study. The region of interest (ROI) is the Aobayama campus of Tohoku University.

Figure 27(a) shows an aspect Pi-SAR images as an example.

Figures 27(b-d) show the results processed by the edge detection of Pi-SAR images.

The most distinct feature of a building object in SAR image is parallel lines. The Hough transform is employed to detect straight segments from the thinned edges, and parallelogram outlines of the façade-images can then be extracted. It is carried out in tiling manner in this paper, i.e. the original picture is partitioned into blocks, each of which is detected independently via Hough transform.

Figure 27.

A PI-SAR image and edge detection.

The detection steps of parallel line segments are: i) find bright spots in transform domain with a minimum distance between every two of them so as to avoid re-detection; ii) search the segments consisted of successive points along the corresponding parallel lines in spatial domain, which is longer than a minimum length, and the distance between two successive points is shorter than a maximum gap; iii) only the pairs of points lying on two lines and facing directly are taken into account for segment searching.

Figure 28.

Extraction of strip-like building-images.

Figures 28(a-d) show the detection from 4 aspects Pi-SAR images. Few building-images are not detected and some false images are wrongly produced. In fact, there is a compromise between over-detection and incomplete-detection. We prefer to preserve more detected building-images, so as to control the non-detected rate, while the false images are expected to be eliminated through the subsequent auto-selection of effective images for reconstruction. However, there always remain some undetectable images, attributable to shadowing of tree canopy, overlapping of nearby buildings or with too complicated wall structures.

The detected parallelogram of a homogenous scattering area could be direct scattering from façade, double scattering of wall-ground, combination of direct and double scatterings, projected shadow of building, or even strong scattering of strip-like objects (e.g. metal fence or metal awning), which is not considered in classification.

Shadowing is identified if scattering power of that area is much lower than the vicinity. Specifically, first set up two parallel equal-length strip windows on its two sides and then calculate the median scattering powers of the three regions. If the middle one is weaker than two sides, it is classified as shadow instead of building image.

The wall-ground double scattering can be differentiated from direct scattering based on polarimetric information. The de-orientation parameter $$v$$ indicating the scattering mechanism is used to identify double scattering.

To implement target reconstruction from multi-aspect SAR data, calibrating multi-aspect data requires that at least one aspect is calibrated beforehand, and other aspects are then calibrated with respect to this calibrated aspect. A natural object, such as flat bare ground, is usually chosen as a reference target, which is expected to preserve identical scattering for all aspects. Thereafter, the channel imbalance factors are estimated from distribution of the phase difference and amplitude ratio of co-polarized, hh and vv, echoes of the reference target, and are then used to compensate the whole SAR images (Xu and Jin, 2008b).

7.2. Building reconstruction

Complexity of objects-terrain surfaces and image speckles incorporate certain uncertainty for detection of object images. To well describe multi-aspect object image, the parameterized probability is introduced for further proceeding of automatic reconstruction. For convenience, the detected building-images are parameterized. Generally, the edge pixels detected by CFAR are randomly situated around the real, or say, theoretical boundary of the object. It is reasonable to presume that the deviation distances of the edge points follow a normal distribution.

The original edge can be equivalently treated as a set of independent edge points. The line detection approach is considered as an equivalent linear regression, i.e. line fitting from random sample points. According to linear regression from independent normal samples, the slope of the fitted line follows the normal distribution.

All parameters have normal distributions and their variances are determined by the variance of the edge points deviation. After counting the deviation of the edge points in the vicinity of each lateral of all detected building-images and making its histogram, the variance can be determined through a minimum square error (MSE) fitting of the histogram using normal probability density function (PDF). The 4-aspect Pi-SAR images are counted.

Figure 29.

Geometry of a wall/façade imaging.

Given a group of building-images detected from multi-aspect SAR images, the corresponding maximum- likelihood probability can be further used as an assessment of the coherence among this group of multi-aspect images. A large maximum-likelihood probability indicates a strong coherence among the group of multi-aspect building-images, and vice versa.

Multi-aspect co-registration, as a critical pre- processing step, is necessary when dealing with real SAR data.

Given the specification of a region of interest (ROI), e.g. the longitudes and latitudes, the corresponding area in SAR image of each aspect can be coarsely delimited according to its orbit parameters. It can be regarded as a coarse co-registration step. Manual intervention is necessary if the orbit parameters are not accurate enough.

Only parameters of the building-images are needed to be co-registered to the global coordinates, but rather than the original SAR images. It is regarded as a fine co-registration step. In this study, only linear co-registrations are considered. A simple and effective method should be manually choosing ground control points. In general, a featured terrain object or building is first selected as the reference of zero-elevation, and then its locations of different aspects are pinpointed. Coherence among multi-aspect building-images is the basis of automatic reconstruction.

Finally, an automatic multi-aspect reconstruction algorithm is designed. The building objects are reconstructed from the 4-aspect simulated SAR images. The inversion accuracy is very good. shows the 3D reconstructed objects on the true DEM. It seems the inverted elevations also match well with the true DEM.

Figure 30.

reconstructed objects and comparisons with optical picture.

Due to the difficulty for authors to collect ground truth data, a high-resolution satellite optical picture (0.6m QuickBird data) is used as a reference. Geometric parameters of each building manually measured from the picture are taken as ground truth data to evaluate the accuracy of reconstruction.

There is a trade-off between the false and correct reconstruction rates. If we increase the false alarm rate of edge detector, relax the requirements of building-image detection and/or increase the false alarm rate of judging correctly matched group, it will raise the reconstruction rate, but also boost the false reconstruction rate. Efficiency will be deteriorated if the false reconstruction rate goes too high. On the contrary, the false reconstruction rate can be reduced and the accuracy of reconstructed buildings can be improved, but the reconstruction rate will also decline.

A critical factor confining the reconstruction precision should be the number of effective aspects, hereafter referred as effective aspect number (EAN).T he reconstruction result will become better if more SAR images of different aspects are available.

Main error of reconstruction is caused by the boundary-deviation of detected building-images, which is originated from complicated scattering and interactions of spatially distributed objects and backgrounds. In probabilistic description of detected building-image, the presented boundary in SAR image is seen as the same as reality. However, the real detected building-image might be biased or even partly lost due to the obstacle and overlapping.

In addition, it is noticed that large-scale buildings might not be well reconstructed. The reason is partly attributed to the premise of Wishart distribution of POL-CFAR detector. Since the images of large-scale buildings might reveal more detail information about the texture feature and heterogeneity, it deteriorates the performance of edge detector. To develop a new edge detector based on certain specific speckle model for high-resolution images can improve the results. Another more feasible way is to employ a multi-scale analysis, through which building-images of different scales can present Gaussian properties in their own scale levels. Of course, the expense for multi-scale analysis might be a loss in precision.

Another issue to be addressed is the exploitation of multi-aspect polarimetric scattering information. In heterogeneous urban area, the terrain objects appear distinctively in different aspects, which gives rise to a very low coherence between their multi-aspect scatterings. Hence, polarimetric scattering information may not be a good option for the fusion of multi-aspect SAR images over urban area

After the ROI is coarsely chosen in each aspect image, edge detection and object-image extraction are carried out, subsequently. Then the object-images are parameterized and finely co-registered. As long as multi-aspect object-images are automatically matched, 3D objects are reconstructed at the same time.

The merits of this approach are the process automation with few manual interventions, the fully utilization of all-aspect information, and the high efficiency for computer processing. Making a reconstruction on a 4-aspect Pi-SAR dataset (500×500) takes less than 10 min CPU time (CPU frequency 3GHz). Complexity of this approach is about $O(K{N}^2)$ , where $K$ is the aspect number and $N$ is the size of SAR images (K=4, N=500 in this case).

It is tractable to extend this approach to reconstruction of other kinds of objects. For other types of primary objects, given the a priori knowledge, new image detection methods have to be developed. It is possible to treat more complicated buildings as combinations of different primary objects, which can be reconstructed separately and then combine together.

Considering a space-borne SAR with the functions of left/right side looking (e.g. Radarsat2 SAR) and forward/backward oblique side looking (e.g. PRISM in ALOS), six different aspect settings are available for ascending orbit. There is a 20 angle between ascending and descending flights for sun-synchronous orbit. Therefore, it can observe from 12 different aspects. Suppose the sensor use different aspect in each visit and the repeat period is 15 days, a set of SAR images acquired from 12 aspects can be obtained through a 3-month observation. Then, application of SBT reconstruction with acceptable precision.

Acknowledgments

This work was supported by the National Science Foundation of China with the grant No. NSFC 406370338.