Modeling and Simulation for InSAR

Lee, S. P., Chan, Y. K., Lim, T. S. and Koo, V. C.

Faculty of Engineering and Technology, Multimedia University, Jalan Ayer Keroh Lama, 75450, Melaka
Malaysia, Email: [Link]. Z@[Link], ykchan@[Link], tslim@[Link],
vckoo@[Link]

Abstract
This paper demonstrates models and simulation methods of interferometry synthetic aperture radar (InSAR).
The InSAR models are presented in detail, which includes modelling of radar signal geometry, wrapped phase
estimation, and phase map estimation. Furthermore, an error model based on interferometric scale is derived
by refering to aforementioned models. To consider the real condition of the return signal echo, Gaussian white
noise model is applied into simulation model as noisy phase image generation. Moreover, simulation also
take scence into challenging terrain conditions like fringes and spiral. The noise or extreme terrain structure
might induce interferometry observation, which may lead to undeterminable phase map. To overcome these
issues, we implement the well-known difficult task, which are two-dimensional phase unwrapping, and
de-noising algorithm. Simulation models discussed in the paper are not only the powerful tools for surface
deformation estimation (or digital elevation map generation), but also useful for advance research on InSAR
phase unwrapping.

1. Introduction
Interferometry synthetic aperture radar (InSAR) has
is wrapped within , however, the
been emerged as a popular topography measurement
principal phase value that represent length of echo
method for geophysical applications i.e. digital
range to terrain surface could be more than
elevation map generation and terrain movement
thousands radian. The principal phase value
investigation. Such remote sensing techniques
recovery job is almost impossible due to
promise a consistent approach to provide the high
non-availability of data recovery necessary
resolution and good precision interferometry
information. Such phase value reconstruction task is
observations in all-weather condition. A number of
called as phase unwrapping. Over decades, despite a
excellent overview papers of InSAR concept and its
large number of research efforts have been done,
applications have been published. (Rosen et al.,
such research problem remains open. Such issue
2000) provides a comprehensive InSAR review
become supporting motive of this research work.
covering the working theories, major algorithmic
The main objective of this research is to develop and
steps and its issues . Balmer treats entire SAR
demonstrate an appropriate simulation tool for
processing as inverse problem (Balmer and Hartl,
advanced research of 2D phase unwrapping.
1998). Thus, he emphasizes on InSAR mathematical
Modelling and simulation works of InSAR are
imaging models and its statistical studies. By
essential for researchers to gain a better handle on
relating InSAR measured data to terrain height,
those phase ambiguos conditions. It also enables the
(Richard, 2007) established an InSAR equation
investigator to concoct different imaginative models
based on geometry considerations . He has delivered
(.i.e. phase unwrapping models) and perform
a clear explanation on the relationship of InSAR
advanced analysis based on its simulated outcomes.
measurement and topography elevation. Thus, we
Various research works related to InSAR simulation
have built and simulated an InSAR geometry model
have been conducted by scientists across the world.
based on Richard’s perspective of geometry
(Xu and Cumming, 1997) introduced an InSAR data
considerations. Phase ambiguity is the central
simulation technique by generating SLC image pairs
problem of InSAR techniques. Overall,
with geometry model and noise model. (Li et al.,
interferometric measured phase are ambiguous by
2005) proposed an iterative interpolation method to
noises, under sampling due to extreme terrain
simulate airborne InSAR complex images. (Lin and
structure, and mathematical ill pose of wrapped
Yang, 2009) developed a simulation model for
phase measurement. Among that, the
distributed satellite-borne InSAR SLC images.
mathematical-ill circumstance has long been viewed
However, the simulation methods mentioned above
as a typical difficult task. SAR measurement phase
are only focus on radar signal or geometry model,
further processing step like absolute phase map
estimation is not taken into consideration. Various
commercial and non-commercial InSAR software
like GMTSAR, Gamma and Doris are able to
complete the further processing steps from above Equation 3
methods (Gens, 1999). However, such InSAR
softwares do not cover artificial data simulation. Where f denotes radar frequency, is initial
Therefore, we proposed a simulation method, which
phase of pulse and indicates pulse length. The
cover complete InSAR processing steps. The
proposed simulation method considers at least two received signal , is expressed on following
distortion cases (which lead to measurement error) equation:
simultaneously while other works consider only one
case. In this paper, we have reviewed and compared
existing algorithms. Also, we have simulated and ,
compared several terrain conditions by the proposed
simulation method. Section 2 introduces the
necessary mathematical framework for InSAR
processing, which covers geometry modelling of ,
InSAR measurement, wrapped phase estimation,
phase map estimation, and noise model (additive
Gaussian white noise). Noise reduction methods Equation 4
include multilook and wiener filter are presented.
The statistical analysis of interferometric phase is Where ρ is the complex reflectivity of Po and we
shown by multilook phase distribution based on its
coherence coefficient. A two-dimensional assume ρ is a fixed value. is a complex- value
phase-unwrapping model is introduced into this that covers amplitude and all radar range factors
section. The proper implementation steps of these other than radar cross section σ (since σ is
models are described in section 3. Also, the proportional to ). Speed of light is
simulation parameters and its taken terrain scenes indicated by c. The ideal baseband received signal
are discussed in this section. Section 4 demonstrates after demodulation can be modelled as below:
simulation results and analysis. Comparisons with
other InSAR software are also made. Section 5
concludes the paper.

2. Modelling of Insar
2.1 Radar Signal Model
During InSAR operation, radar antenna propagates
frequency to target area and collect its echo signal
for data observation. The coherent radar receiver
consists of in-phase (I) and quadrature (Q) channels
that lead to I term and Q term of radar measurement.
The amplitude A and phase φ of the received echo
can be calculated by refer to these I and Q data.
Equation 5

By exploiting the phase of received signal, the

measured data sample is defined as equation 6:
Equation 1

Equation 2 Equation 6

The echo signal can be expressed as complex Where arg(.) is the function operating on complex
sinusoids . The transmitted signal number to return angle value.
can be written as:
2.2 InSAR Geometry Models Equation 7
Topography measurement can be done by
space-borne or air-borne InSAR. Such system can Differential of the function are calculated:
be implemented in two ways which are along –track
and across-track. Along track InSAR system require
at least two pass operations since it is operated in
such a way that one antenna record an aperture
image at first pass operation, then following
aperture images are recorded after certain timespan Equation 8
using same viewing geometry at another pass
operations. Such technique allows the detection of By rewriting the differential result in this form, we
spatial changes of land surface for temporal can relate the range variation respect to height
evolution purpose. The following InSAR equations variation.
are derived by referring to Richard’s perspective of
geometry considerations (Richard, 2007). His
geometry figures have been adapted into this paper
as figure 1 and figure 2. Figure 1 shows the Equation 9
geometry of a radar antenna located at ground plane By considering such effect on ground plane as P0
(reference plane) coordinate y = 0 and altitude where z = 0, so δh = 0. It becomes:
coordinate z = h. x coordinate is normal to the page,
therefore it is not shown in the figures.

Equation 10

In this scenario, the measured phase in equation (6)

is can be expressed as:

Figure 1: Geometry of along-track InSAR system. Equation 11

This diagram based on figure 5 found in (Richard,
2007) During InSAR processing, the real three dimension
world is projected into two dimension plane. Range
Po, P1 and P2 indicate positions of scatterer point. value of scatterer P1 is lower than P0 since the time
For example, Po indicates a position on ground delay arrived earlier due to its higher location. As a
plane z = 0, range dimension y = y1 with grazing result, ground coordinate of P1 is located at y3.
angle or depression angle of θ radian. R0 is the
range to Po. In geometry model of figure 1, P1 is
elevated by the height of δh from P0. The variation
of relevant height δh is represented by following
Equation 12
function:
Equation (11) and (12) shows the single sample of
processed phase information with its coordinate in
one dimensional line (y-plane). The InSAR systems
travel in x-plane, record samples line-by-line and
lead to two-dimensional array of phase image
formation. Recall equation (6), G(x,y) indicates
aperture image of phase information from received
signal. =

Equation 15

By referring to the geometry of figure 2, we obtain:

Equation 13

Figure 2 illustrates the geometry model of

across-track InSAR system. Across-track InSAR
system consists of two antennas separated by a
baseline distance B normal to the plane of moving
track. G1(x,y) and G2(x,y) represent two aperture
image information from two antenna of different
position. In this system, Po and P1 indicates the same
target position from apertures of G1(x,y) and
G2(x,y).The different location of the antenna sensors Equation 16
are lead to a height different of δh. Such principle
could be applied for elevation map generation. By Considering δθ as very small, equation 16 can be
refering to figure 2, the different between G1(x,y) written as:
and G2(x,y) in term of range is B cosθ.

Equation 17

By solving equation 15 and equation 17, we get a

topography map of elevation measurement.

Equation 18

From equation 18, we derive an interferometric

Figure 2: Geometry of across-track InSAR system.
This diagram adapted from figure 8 found in scale as = . Since the
(Richard, 2007) maximum phase variation of is 2

Thus, the phase difference between the two , so abrupt change of height among adjacent
apertures becomes: sample beyond ḣ meter will lead to a residue pixel
or whole image corruption. Such circumstance is
model to equation 19.

ḣ= 2π
Equation 19
Equation 14
It is noted that the radar geometry model is not fully
Differentiating equation 14 with respect to grazing represents real world of SAR data acquisition since
angle gives: SAR issues like backscattering and shadowing are
not being taken into consideration. In real
circumstance, each resolution element is also estimates the differences of the phase at each image
constitutes of a lot of arbitrary scatters. Instead of point in either images dimension as an
considering InSAR data as mean value of the approximation to the derivative and integrate the
geometrical point scatterers, we should treat it as result. Once the discontinuity is detected among
imaging process to construct three dimensional adjacent pixel, multiple of 2π will be added or
terrain structure. deducted in order to remove the discontinuity. By
referring to equation 21, we rewrite the absolute
2.3 Wrapped Phase Estimation Model and Its phase determination function where
Problem Formulation
During InSAR processing, estimation of phase is the estimate absolute phase value,
difference between apertures of SAR images denoted wrap phase. The
(wrapped phase estimation) is represented by number of cycle k is determined by summation of
interferogram. Interferogram can vertical adjustment parameter,
be generated by computing and horizontal adjustment parameter,
multiply conjugate of . .∆h and ∆v denoted the phase difference
and among adjacent pixel horizontally and vertically,
while m and n indicate column and row number of
are apertures of SAR images. the interferogram.

Equation 20

The phase values extracted via receiver are wrapped

in the interval of (-π,π]. However, absolute phase
value which represents the range to ground surface
can be many radians (Zebker and Lu, 1998). As
shown in equation 21, Φ denoted absolute phase. k
is number of wavelength and Φp stands for principle
phase value.

Φ(x,y) = 2πk(x,y) + Φp(x,y)

Equation 21
Equation 22
The phase difference estimated from equation 20 is
not the absolute phase difference (it is wrapped
2.5 Noise Model
phase difference) due to 2π modulus issue. In
In real circumstance, such phase measurements are
addition, extreme terrain structures (e.g fringes) and
induced by unknown noise from echo signal. To
noise will further corrupt the interferometry
simulate this scenario, we generate a noisy phase
measurement. This makes phase determination
image by adding the complex-value Gaussian white
ambiguous. Thus, two dimensional phase
noise Ɲ(.) to the interferogram. Equation 23 show
unwrapping (2D- phase unwrapping) to recover
the PDF model of Gaussian white noise where σ
unambiguous phase value becomes the key problem
of InSAR processing. Strictly saying, data recovery stands for variance and indicate mean of input
work is an impossible task due to its mathematical variable (we employed zero- mean Gaussian model,
ill pose. Literature study reveals that absolute phase ).
can be estimated via interferometry spatial or
temporal context (Ghighlia and Pritt, 1998). Thus,
this can be estimated based on assumptions.

2.4 Two Dimensional Phase Unwrapping Model

The following equations shows 2D-phase Equation 23
unwrapping model based on Itoh condition
(Munther and Francise, 2013). Such method 2.6 Noise Reduction Models
In this paper, two- dimensional wiener filter (Matlab Equation 27
help file) is employed as denoising method by
estimating the local mean and variance around each Where
pixel.
and F(.) is a hypergeometric
function. is gamma function. The
coherence coefficient is defined as follow:
Equation 24

where local mean is estimated by refering to nearby

adjacent pixels within N x M Equation 28
window.
Where |C| is the coherence and is the phase of
coherence coefficient.

3. Proposed InSAR Simulation Technique

Figure 4 shows the proposed InSAR simulation
Equation 25 method. In this simulation, various artificial
three-dimensional terrain structures (include fringes,
Another noise reduction method, multilook spiral and hill plane) is generated. We also simulate
techniques is also frequently employed to overcome interferometry observed data from different aperture
phase noise by averaging neighbouring pixels. views (G1 and G2) by referring to aforementioned
radar signal geometry model.
The interferogram is generated via computation of
the interference pattern caused by the phase
difference between two SAR images acquired from
G1 and G2. The interferogram is an estimation map
of wrapped phase difference. Additive
complex-valued Gaussian white noise is generated
Equation 26 to produce an ambiguous wrapped phase image. In
this step, multilook techniques could be
It is noted that and implemented to overcome the phase noise and
reducing phase unwrapping problem. Then, we
are the pair of complex single look SAR
implement a 2D phase unwrapping method, i.e., Itoh
images, n is the number of look. This is a promising
method for absolute phase map estimation. 2D
technique that Lee have derived following pdf for
Wiener filter with 5 x 5 window is employed for
multilook phase distribution (Lee et al., 1994)
de-noising purpose. Finally, terrain height
(Figure 3):
information can be calculated by relating the
unwrapped and de-noised outcome as geophysical
measurement to radar signal geometry. In this paper,
we implemented these simulation steps in two ways,
with (4 look) and without multilook (single look) by
refer simulation parameter of Table 1.
Figure 3: Pdf Plotting of Multilook Phase Distribution (c=0 and =0)

Table 1: Simulation Parameter

Radar Frequency 10 GHZ(X-band)

Wavelength 0.03m
Antenna Height 5000m
Baseline distance 1m
Grazing Angle 30 degree
Noise Variance 0.4
 Figure 4: Simulation flow chart for InSAR

Figure 5: Plotting of interferometric phase vs coherence coefficient

4. Simulation Result and Comparison Table whole structure, fringe structure, and spiral
4.1 Simulation Result structure. As shown in Figure 6, Figure 7 and Figure
Proposed simulation techniques are implemented by 8, each set of terrain structure simulation results are
referring to simulation parameters of Table 1. Such organised as follow: first diagram represents
simulations take into consideration of three terrain relevant artificial terrain structure (or expected
structure types including terrain hill and sinking outcome). Second result indicates a noisy
interferogram after InSAR processing.

Figure 6: InSAR Simulation of terrain hill & sinking holes

 Figure 7: InSAR Simulation of fringe

Figure 8: InSAR Simulation of spiral plane

Then, it follows by reconstruction plots of relevant encounter the abrupt height variation more than
terrain structure after phase unwrapping and 600m. By referring to equation 19, we substituted
denoising. Such final plotting could be compared the simulation parameter and defined our
with expected outcome (first diagram) for interferometric scale at =[Link] geometry
performance evaluation. Figure 6 shows simulation limit of such InSAR system would be ḣ=44.35 x
results of terrain hill and sinking hole. We observed 2π≈ 260m. In this scenario, the height variation of
that the shape of final simulated result is very close fringe is far beyond its geometry limit. Thus, final
to the expected outcome. In this case, noise is processing outcome offringe structure has
well-suppressed and the phase unwrapping method reconstructed to a wrong shape at its effected region.
has well-approximated and reconstructed the terrain Figure 8 shows simulation of spiral plane with a
structure into proper value. This simulation results failure final simulated outcome that both Itoh
shows great performance since all pre-assumptions method and wiener filter failed to reconstruct the
are well-hold. During fringe simulation in Figure 7, terrain map of spiral shape. Such spiral plane
we observe that terrain map is not being considered as extreme structure since more than 80
reconstructed properly at edge area during Itoh percent of its region encounter height variation
unwrapping. We notice that relevant edge area beyond the limit. The proposed phase-unwrapping
method failed to define the spiral edge since its by incorporating multilook technique. Figure 5 plot
working principal define the integral path vertically of interferometric phase versus coherence
and horizontally while relevant edge shaped spirally. coefficient is matching with the theoretical result of
In this scenario, all pre-assumptions are no longer multilook phase distribution as Figure 3. Both
held. Such extreme terrain structures simulations showing that the standard deviation decreases as
show that the Itoh unwrapping approach allows coherence coefficient and number of look increases.
local errors due to abrupt phase variations beyond For comparing the final outcome with the actual
interferometric limit, causing errors across full value of the simulated terrain structure, we could
unwrapped terrain image. This solution is sensitive use root mean square error to assess the
to inconsistencies in wrapped phase, leading to performance of this simulation method. Such metric
significant errors. Since two-dimensional phase are used for comparison of Table 2 on multilook
unwrapping is a well-known unsolved issue, technique. For terrain hill & sinking holes structure,
advance research could be conducted for extreme the multilook technique has reduce 2.04 RMSE.
terrain structures based on proposed simulation Likewise, 1.08 RMSE has been improved in fringe
method. Based on coherence evaluation of Figure 5 structure. Due to high destructive simulated
and RMSE studies, performance of proposed outcome of spiral, its RMSE result is not appropriate
simulation method shows significant improvement to be taken into multilook studies.
Table 2: Root Mean Square Error of InSAR Simulation
Terrain Structure Type Terrain hill & sinking holes Fringe Spiral
Single Look 8.28 5.18 >100
Multi-Look (4 look) 6.24 4.10 >100

Table 3: Comparison of InSAR Software

Features Aforementioned Simulation Processing Proposed
techniques (Xu and Cumming, Software e.g Simulation
1997), (Li et al., 2005), (Li and GMTSAR, Doris method
Yang, 2009)
Artificial input data Yes No Yes
Noise simulation Yes No Yes
Denoising and phase No Yes Yes
unwrapping
Less complexity Yes No Yes

Our simulation results have demonstrated the contributions of aforementioned parties are sincerely
complete InSAR processing. We also reviewed and appreciated and gratefully acknowledged.
compared existing simulation approaches of InSAR References
processing software with our proposed method as
table 3. Bamler, R. and Hartl, P., 1998, Synthetic Aperture
Radar Interferometry, Inverse Problems, 14,
5. Conclusion 1-54.
The modelling and processing of InSAR are Gens, R., 1999, SAR Interferometry: Software, Data
investigated in this paper. Besides, an efficient Format, and Data Quality, Photogrammetric
simulation method is proposed and demonstrated. Engineering and Remote Sensing, 65(12),
The simulation results of artificial terrain structures 1375-1378.
and the process of interferogram, unwrapping and Ghighlia, D. C. and Pritt, M. D., 1998,
denoising terrain maps are presented. Such methods, Two-Dimensional Phase Unwrapping: Theory,
also, can be great valuable for performance analysis Algorithms, and Software, 1st edition, (New
and advance 2D phase-unwrapping research of York: Wiley-Interscience).
InSAR system. Lee, J. S., Dewaele, P., Wambacq, P., Oosterlinck,
A. and Jurkevich, I., 1994, Speckle Filtering of
6. Acknowledgement Synthetic Aperture Radar Images - A Review,
This work are supported by JICA (Japan Remote Sensing Reviews, Vol. 8, 313-340.
International Coorperation Agency) and JST under Li, W. Q., Li, J. W. and Chen, J., 2005, Fast
SATREPS Project, Ministry of Education Malaysia Simulation of Complex Image for
(MOE) under Grant 203/PJJAUH/6711279. The Interferometric SAR, Journal of Beijing
 University of Aeronautics and Astronautics, 31(
1), 31-35.
Lin, L. and Yang, W., 2009, Modeling and
Simulation of Single-Look Complex Images for
Distributed Satellite borne Interferometric
Synthetic Aperture Radar, Image and Signal
Processing, 2009. CISP '09. 2nd International
Congress.
Matlab help file [Online]. Available:
[Link]
[Link]
Munther, G. and Francise, L., 2013, Two
Dimensional Phase Unwrapping Problem,
[Online].Available:
[Link]
t-us/faculties-and-schools/tae/geri/two_dimensio
nal_phase_unwrapping_finalpdf.pdf?la=en
Richards, M. A., 2007, A Beginner’s Guide to
Interferometric SAR Concepts and Signal
Processing, IEEE A and E Systems Magazine,
22(9), 5-29.
Rosen, P. A., Hensley, S., Joughin, I. R., Madsen, S.
N., Rodriguez, E. and Goldstein, R. M., 2000,
Synthetic Aperture Radar Interferometry.
Proceeding of IEEE, 88(3), 333-382.
Xu, W. and Cumming, I., 1997, Simulator for
Repeat-Pass Satellite InSAR Studies, Proceeding
of International Geoscience and Remote Sensing
Symposium 1997, 1704-1706.
Zebker, H. A. and Lu, Y. P., 1998, Phase
Unwrapping Algorithms for Radar
Interferometry—Residue-cut, Least-squares, and
Synthesis Algorithms. Journal of the Optical
Society of America, 15(3), 586–598.