CENTIMETER RANGE MEASUREMENT USING AMPLITUDE DATA OF TERRASAR-X IMAGERY

The SAR (Synthetic Aperture Radar) imagery are largely used for the environmental, structures and infrastructures monitoring. In particular, Differential SAR Interferometry (DInSAR) is a well known technique that allows producing spatially dense displacement maps with centimetre to millimetre accuracy. The SAR signal is characterized by phase and amplitude value and the DInSAR remote sensing technique allows to analyse deformation phenomena affecting both extended natural areas and localized man-made structures, by exploiting the phase difference of SAR image pairs. New SAR satellite sensors such as COSMO-SkyMed, TerraSAR-X and PAZ offer the capability to achieve positioning in a global reference frame accuracies in the meter range and even better, thanks to the higher image resolution (up to 0.20 m pixel resolution in the Staring SpotLight mode for TerraSAR-X and PAZ) and to the use of on board dual frequency GPS receivers, which allows to determine the satellite orbit with an accuracy at few centimetres level. The goal of this work is to exploit the slant-range measurements reaching centimetre accuracies using only the amplitude information of SAR images acquired by TerraSAR-X satellite sensor. The leading idea is to evaluate the positioning accuracy of well identifiable and stable natural and man-made Persistent Scatterers (PS’s) along the SAR line of sight. The preliminary results, obtained on the Berlin area (Germany), shown that it is possible achieve a slant-range positioning accuracy with a bias well below 10 cm and a standard deviation of about 3 cm; the results are encouraging for applications of high resolution SAR imagery amplitude data in land and infrastructures monitoring.


INTRODUCTION
New earth observation SAR (Synthetic Aperture Radar) satellite sensors, as COSMO-SkyMed, TerraSAR-X and PAZ, acquire imagery on any point of the Earth with high resolutions, in terms of phase and amplitude value.Therefore, these data are routinely and conveniently used in order to monitoring deformation phenomena impacting the Earth surface (e.g.landslides, subsidences, volcano deformations and glacier motions) and infrastructures (e.g.buildings, dams, bridges).The main remote sensing technique to extract centimeter information from SAR data is the Differential SAR Interferometry (DInSAR), based on the phase information only (Berardino et al., 2002).Specifically, this technique is based on an appropriate combination of differential interferograms generated by image pairs characterized by a small orbital separation (baseline), in order to limit the spatial decorrelation phenomena.On the other hand, it is well known that DInSAR technique may suffer for lack of coherence among the considered images, and that it is particularly suited for slow deformation phenomena, due to the intrinsic need to unwrap the phase, what can be difficult in presence of displacements (much) higher than the wavelenght of the SAR signal (few centimetres).In addition, the new high resolution SAR satellite sensors offer the capability to achieve positioning accuracies in a global reference frame in the meter range and even better, thanks to the very high image resolution (up to 0.20 m pixel resolution in the Staring SpotLight mode for TerraSAR-X and PAZ) and to the use of on board dual frequency GPS receivers, which allow the determination of the SAR satellite orbit with an accuracy at few centimetres level.This is the why it appeared advantageous in the last years the development of a different approach for land and infrastructures monitoring by SAR, exploiting the signal amplitude and avoid- * Corresponding author ing coherence deficiency and displacement magnitudes restriction (Nascetti et al., 2014).The very first idea was to evaluate the positioning accuracy of well identifiable and stable natural and man-made Persistent Scatterers (PS's) (e.g.corner reflectors) along the SAR line of sight.This information is of course of crucial importance in order to assess the potential of this technique to monitor the possible displacements of the considered PS's.The first experiments, both using suitable corner reflectors positioned by high precision GPS surveys and natural PS's, were carried out with stacks of TerraSAR-X SpotLight imagery, reaching a slant-range measurements accuracy with a bias of about 90 cm and a standard deviation of about 4 cm (Eineder et al., 2011).Then, these experiments were repeated and refined, reaching similar results as regards the standard deviation but with a more and more reduced bias of about 29 cm (Cong et al., 2012) and 10 cm (this result was obtained with corner reflectors) (Balss et al., 2013) respectively.Here we propose a partially different approach, considering only natural PS's, where, at first, the approximate reference position of each PS is determined through a stereo approach.Then the positioning accuracy of a PS along the SAR line of sight is evaluated starting from this reference position and the orbital information supplied in the SAR imagery metadata, and accounting both for signal propagation delays and for geophysical effects causing not negligible PS displacements.In this respect, three are the major problems to face with.At first, it is necessary to recognize and match the same PS across the images of the considered stack with sub-pixel accuracy (millimitres to centimetre level).To this aim, an area where a PS is clearly visible is manually selected and subsequently an automatic matching procedure is adopted.This procedure is based on a template matching applied on small patches, surrounding the PS, magnified by a suitable factor (100x to 200x).Secondly, it is necessary to correct for the signal propagation delays through The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-7, 2014 ISPRS Technical Commission VII Symposium, 29 September -2 October 2014, Istanbul, Turkey the troposphere and ionosphere.Finally, it is necessary to filter out the known geophysical effects inducing periodic and secular ground displacements.In section 2 it is outlined our methodology for the positioning accuracy assessment of a PS along the SAR line of sight; in section 3 the employed stereo model is recalled; in section 4 the corrections applied to the measured slant-ranges for atmospheric delays and geophysical effects are shortly summarized; in section 5 they are presented and discussed the main results of the investigation, carried out on a stack of 20 TerraSAR-X SpotLight images acquired on Berlin urban area (Germany) between February 2008 and March 2009; finally, some conclusions are drawn.

METHODOLOGY
The methodology proposed in this work is devoted to evaluate the positioning accuracy of a PS along the SAR line of sight using the amplitude information only.

PS absolute position estimation
Starting from an interferometric repeat-pass stack (composed by N SAR images) and using one single same side image acquired on a different orbit, it is possible to retrieve the 3D position of a PS in a global reference frame (ECEF -Earth Centered Earth Fixed) using a stereo approach.To this aim, it is at first necessary to measure the image coordinate of the PS both in all the stack images and in the single image and subsequently, solving the radargrammetric equations (see section 3), the PS position for each pair (one stack image with single image) can be computed.In this way, a good redundancy is introduced, computing N PS positions using N different pairs; their median is taken as robust estimation of the PS reference position.

Image coordinates measurement
An automatic matching procedure has been implemented to measure the PS image coordinates in all the stack images.A patch is defined in the first stack image centered on the selected PS and a template matching is performed between this patch and a corresponding area in all other stack images.A two steps approach was used in order to achieve a sub-pixel accuracy.First we used a moving window of 151 pixels x 151 pixels and a searching window of 351 pixels x 351 pixels and we computed the Normalized Cross-Correlation (NCC) on image intensity after oversampling by a factor of two both in range and azimuth directions.Secondly, after the determination of the position of the maximum NCC value, in order to refine the PS image position, we took the two corresponding patches in the master and in the slave images and we oversampled them by a factor of 100 and 10 in the range and azimuth directions respectively.Using this procedure it is possible to achieve a raster precision of 1/200 samples in range corresponding to 2.5 mm and a 1/20 samples in azimuth corresponding to 20 mm.

Range measurement comparison
Starting from the ECEF positions of PS and SAR sensor when imaging the PS itself (it can be computed by orbital information supplied into the metadata), it is possible to easily compute the the SAR sensor-to-PS distance and to compare it with the slantrange distance measured by the SAR sensor.The corresponding differences ∆: As briefly reported in the section 4.2 the range measures are affected both for atmospheric and geophysical effects, which must be considered.In particular the ∆R values have been corrected accordingly with this equation: where: • TD is the Tropospheric Delay (see section 4.1.1) • ID is the Ionospheric Delay (see section 4.1.2) • EM is the Earth Motion (combination of the geophysical effects projected on the SAR Line of Sight -LOS, see 4.2) • δ is the residual displacement of the PS respect to the reference estimated position • is the residual error due to both orbit errors (probably the largest part) and atmospheric and geophysical model errors It is important to underline that the residual error is likely to be spatially correlated over the entire area, otherwise the point displacement is strictly dependent to the point location.Considering these differences, one o more stable points could be conveniently used to eliminate by differences the model error highlighting the point displacements of the points affected by displacement.

SAR STEREO MODEL
As well-known, SAR images are focused in zero-Doppler slant range geometry, that is the target appears as acquired on a heading perpendicular to the flying direction of SAR satellite sensor and the slant range constrain represents the SAR sensor-to-PS distance.Here we recall the stereo (radargrammetric) model used to estimate the reference position of each considered PS, then to compute the SAR sensor-to-PS distance.The implemented stereo model is based on the equation of radar target acquisition and zero Doppler focalization (Leberl, 1990).The radargrammetry technique performs a 3D reconstruction based on the determination of the sensor-object stereo model, in which the position of each point on the object is computed as the intersection of two radar rays coming from different positions and therefore with two different look angles.Actually, these radar rays can be simply modeled as two segments of measured lengths starting from two different positions (along two different satellite orbits), so that the intersection generating each object point is one of the two possible intersections between two circumferences centred in the two different positions and laying into two planes orthogonal to the two different satellite orbits whose radii are equal to the segment measured lengths (Fig. 2).Consequently, the stereo model is  Therefore, the couple of equations in an ECEF (Earth Centered Earth Fixed) system (for example the WGS84) reads: (3) where: • XP , YP , ZP are the coordinates of the generic ground point P in the ECEF coordinate system • XS, YS, ZS are the coordinates of the satellite in the ECEF coordinate system • uX S , uY S , uZ S are the Cartesian components of the satellite velocity in the ECEF coordinate system • Ds is the so-called "near range" • ∆r is the slant range resolution or column spacing • I is the column position of point P on the image Since the satellite angular velocity can be considered constant along the short orbital arc related to the image acquisition, the epoch of line acquisition t can be related to the corresponding line number J through the following linear function: where start time is the time of start of acquisition, P RF is the Pulse Repetition Frequency.The orbit computation consists of estimation of satellite position at the epoch t of each line acquisition according to zero Doppler geometry.In the metadata file of each SAR imagery the ECEF position and velocity of satellite related to the time are supplied through state vectors at regular intervals, whose number N depends on the considered SAR sensor.
The orbit interpolation has been performed by Lagrange polynomials ( 5), whose degree depends on the state vectors number N .
Lagrange polynomial interpolation is enough accurate to model the short orbital segment and its well-known oscillation problems at the edges do not affect the modeling since the images are acquired just in the central part with respect to the supplied orbital segment.Additionally, using a standard divide and conquer algorithm it is possible to find in a fast and accurate way the seeked epoch t when satellite orbit is perpendicular to the line of sight between the sensor and the generic ground point (Capaldo et al., 2014).

RANGE CORRECTIONS
As mentioned, in order to evaluate the positioning accuracy of a PS along the SAR line of sight it is necessary to correct for the signal propagation delays through the troposphere and ionosphere and for known geophysical effects inducing periodic and secular ground displacements (Fig. 1).
The tropospheric delay (TD), which is by far the most important correction to be applied, is close to 2.5 meters in zenith direction at sea level at mid latitude, and the needed correction along the slant path at better than centimetre level can be estimated by GNSS phase observations collected in the area imaged by SAR, using a proper mapping function (Bonafoni et al., 2013).The ionospheric delay (ID) in X-band is in the order of few centimetres and can be computed by a global model (e.g.Klobuchar) (Klobuchar and John, 1987).The geophysical effects consist primarily in the solid Earth tides (SET), polar tides (PT), crustal deformations due to ocean loading (significant close to shorelines) (TOL), which globally can be at the level of half meter and can be corrected using the International Earth Rotation Service Conventions (McCarthy and Petit, 2003), (Petit and Luzum, 2010); moreover, geodynamics (G), crustal deformations due to atmospheric loading (variation of atmospheric pressure) (AL), glacial isostatic adjustment (significant mainly in north Europe and Canada) (GIA) and the effect of the seasonal hydrological loading (HL) have to be accounted for.All the geophysical effects are grouped within the Earth motion (EM) effect.
For each PS, both the atmospheric and the geophysical effects are computed in the geodetic local system (East, North and Up) centred in the PS itself, then they are summed up and projected on the SAR line of sight (LOS).

Tropospheric delay
In order to estimate the tropospheric delay (TD) we considered the ZT D (Zenith Tropospheric Delay) hourly estimated at the public GNSS permanent station closest to the imaged area and included in the European Permanent Network (http://igs.bkg.bund.de/rootftp/EUREF/products/). ZT D is the excess path length due to the signal travel through the troposphere at zenith.Obviously, a more refined solution could consider a spatial interpolation including at least three GNSS permanent stations around the imaged area.Then the ZT D was corrected for the time and for the (mean) height difference between the imaged area and the GNSS permanent station.
The ZT DGP S relative to the time of the SAR image acquisition was obtained through a linear interpolation between two consecutive hourly values.The height correction was computed applying the the NATO Standardisation Agreement -STANAG 4294 (NATO, 1997) tropospheric model.Both for the GNSS station and for the (mean) height of the imaged area the ZT D were computed according to STANAG model, and their difference ∆ZT D h SAR −h GN SS was applied to get ZT DSAR from the ZT DGNSS: The TD is then derived from ZTD by a mapping function; if the image incidence angle θ is lower than 60 degrees, a suitable mapping function is just the simple inverse cosine ( 1 cosθ ).Therefore, in this case the relationship between the ZTD and the SAR image TD reads:

Ionospheric delay
The ionospheric delay has a very small influence, within few centimeters, in the X-band (e.g.9.65 GHz, for TerraSAR-X) and it was just computed using the Klobuchar model and the related coefficient supplied in the GPS broadcast navigation message on a two-hours basis (Klobuchar and John, 1987).

Geophysical effects
The SET, TOL and PT models are well-described in the International Earth Rotation Services (IERS) Conventions 2010 (Petit and Luzum, 2010).The open source library requests for the computation of the TOL a set of coefficients available on line and derive from one among the many ocean loading models.In our case the model used for the ocean tide is the GOT00.2 that is a pure hydrodynamic tide model tuned to fit tide gauges globally using TOPEX/Poseidon data and are given on a 0.5 by 0.5 degree grid.The effect of the TOL is negligible for the site far away from the coast; the largest magnitude has for the areas near to the ocean.On the other hand, in the areas far away the ocean the changes of the weight of the column of atmosphere due to variations of pressure result in crustal deformations called atmospheric pressure loading.These variations on average have the rms of 2.6 mm for the vertical component and 0.6 mm for the horizontal component, but peak to peak variations can reach 40 mm for the vertical component and 7 mm for the horizontal one.The series of the 3-D displacements due to pressure loading with a 6 hour time resolution at 2.5 by 2.5 degrees grid are computed for 751 GPS station, they are available at http://gemini.gsfc.nasa.gov/aplo/.The correct value at the time of the SAR image acquisition is obtained with linear interpolation.

Global Geodynamics
Considering that the time span of the images stack could be several years long, the effect of global geodynamics has to be accounted for.It can be evaluated on the basis of the closest GNSS permanent station included in one of the international networks (e.g.European Permanent Network, International GNSS Services, SIRGAS), for which public data concerning position and velocity are available in the current International Terrestrial Reference Frame (at present ITRF2008).

Glacial Isostatic Adjastement and Hydrology Loading
The glacial isostatic adjustment (GIA) and the present-day rateof-change crustal uplift must be considered specially for areas of Northern Europe and Canada; they are available on line at http://www.psmsl.org/trainand info/geo signals/gia/peltier/.In addition mass loading is known to cause surface deformations.Changes in water mass or ice mass can result in crustal deformations of the surface.Hydrology loading (HL) is usually very seasonal and is most significant in South America, Southern Africa and Asia, where the variation in water mass is large.It is computed on 1 x 1 degree resolution with monthly time resolution.Usually the peak to peak variation of the hydrology loading is between 3-10 mm in the vertical component.The horizontal displacements are much smaller and the peak-to-peak variation is usually not more than a few millimeters (http://lacerta.gsfc.nasa.gov/hydlo/).

EXPERIMENT
The available imagery for the experiments on Berlin (Germany) test site are an repeat pass TerraSAR-X stack of 20 High resolution SpotLight images acquired from 2008 and 2009, with a mean incidence angle close to 30 degrees (Fig. 3), a swath of 10 km x 5 km and a resolution of 1.1 m in azimuth and 0.45 m in range.The sensor TerraSAR-X operates in X-Band, that means a frequency of 9.65 GHz and a wavelength of few centimetres.
In Table 1 the acquisition date and the relative mean incidence angles are displayed.The first image of the Table 1, is the one used in the stereo radargrammetric model in order to determinate the 3D point positions.It is important to note that all the images are provided with the science orbit products.These kind of orbit products are computed with a latency of several days necessary to use high-quality scientific GPS ephemerides provided by the Center for Orbit Determination in Europe (CODE).This has a significant impact on the overall quality of the orbit product that achieve a 3D accuracy of about 4 cm.
As regards the data necessary for the evaluation of the tropospheric delay (   For the experiment we have selected three different stable PS's in one of the available images (see figure 3) and we have supposed that this points are not affected by a remarkable displacement (δ = 0) during the entire acquisition period (about one year from 02/2008 to 01/2009).The PS positions has been estimated using the described stereo approach and the results are reported in Table 2; for each of the three selected PS's the standard deviation is less than 10 cm for the X component and is about 20 cm for the Y and Z components.
The values obtained applying the described methodology are reported in the Table 3.As shown the residual errors for all the PS's have a standard deviation of about 3 cm and a bias well below than 10 cm.Moreover the trends for each PS observed in the acquisition period has been investigated and are reported in the Fig. 4; it is clearly visible the quite good correlation between the trends (min CC of 0.72 and max CC of 0.82) that is probably due both to orbit errors and to the not totally modelled atmospheric and geophysical effects.The PS absolute position has been computed using a stereo approach and the necessary PS image coordinates has been retrieved using an automatic matching procedure.Starting from the positions of PS and SAR satellite when imaging the PS itself (it can be computed by orbital information supplied into the metadata), it is possible to easily compute their distance and to compare it with the slant-range distance measured by the SAR sensor.The signal propagation delays through the troposphere and ionosphere and the known geophysical effects inducing periodic and secular ground displacements have been considered in order to remove these phenomena from the distances measurements.
The available imagery used for the experiments on Berlin (Germany) test site are a repeat pass TerraSAR-X stack of 20 High resolution SpotLight images acquired from 2008 and 2009, with a mean incidence angle close to 30 degrees and one single image acquired with a mean incidence angle close to 50 degrees.Three different and stable PS's have been selected and the results obtained show that it is possible to achieve a slant-range positioning accuracy with a bias well below 10 cm and a standard deviation of about 3 cm.
In the future this methodology could be conveniently adopted to monitoring deformation phenomena affecting both extended natural areas and localized man-made structures.

Figure 1 :
Figure 1: SAR acquisition geometry and slant-range corrections based on two standard equations (3): the former represents zero Doppler projection and the latter the slant range constrains(Capaldo et al., 2014).

Figure 2 :
Figure 2: SAR acquisition system in zero Doppler geometry

Figure
Figure 4: residuals TD, see 4.) the Potsdam GPS station is used; it

Table 1 :
Berlin TSX test site images features has an ellipsoidal elevation of 144 m, around 65 m higher than the Berlin city and it is far around 30 km from Berlin.