INTERFEROMETRIC SAR DEFORMATION MONITORING USING PASSIVE REFLECTORS AND ASCENDING AND DESCENDING PASSES

: The article concerns the estimation of land deformation in an area of copper mining using the Persistent Scatterer Interferometry (PSI) technique and C-band data coming from the Sentinel-1 satellite. The aim of the work is a comparison of two different PSI approaches exploiting passive corner reflectors deployed pairwise, oriented towards the satellites in ascending and descending orbits. The comparison was done in the line of sight geometry. In addition, a validation was carried out using GNSS measurements. The results achieved using the two PSI approaches are consistent and mutually coherent. They are also congruent with observations from GNSS stations. These results proves that the method proposed to identify pixels on multi-looked images is correct. In addition, the intermediary results on the corner reflectors shows that their responses, in terms of dispersion of amplitude and the coherence, is slightly different. This needs to be further studied.


INRODUCTION
This work concerns the estimation of land deformation using the Persistent Scatterer Interferometry (PSI) technique and Cband data coming from the Copernicus Sentinel satellites 1A and 1B. It belongs to a research project entitled "Automatic system for monitoring the impact of high-energy para-seismic shocks on the surface of the area using GNSS/PSI satellite observations and seismic measurements (ASMOW)", co-financed by the European Regional Development Fund (XI 2018-XII 2021. The project aimed at the creation of the local service for geohazard (mining induced para-seismic events) monitoring. A dedicated system has to integrate seismic observations, GNSS measurements as well as satellite interferometry to map ground displacements related to underground exploitation of copper, and to present ASAP the consequences of para-seismic incidents. A dedicated software has been developed by the Centre Tecnològic de Telecomunicacions de Catalunya (CTTC) as part of the ASMOW project, analyzing both time series of Sentinel-1 ascending and descending scenes as well as differential interferograms (3-pass or 4-pass) generated using the pre-and post-seismic event SAR acquisitions.
The PSI and artificial Corner Reflectors (CR) have been used since many years as a geodetic measurement technique for ground and infrastructure displacement monitoring. E.g. an array of 40 trihedral corner reflectors-triangle co-located with survey marks, was installed by the end of 2014 in the northern Queensland in Australia (Garthwaite et al., 2015). By performing a series of interferograms formed from many SAR acquisitions the velocity and time-series maps can be generated. Small multi-directional radar reflectors have been used for stereo SAR measurements by the German Space Agency DLR (Gisinger et al., 2016), (Balss et al., 2018). For landslides monitoring in Dubrava city low-cost corner reflectors have been used by (Lazecky and Kacmarik, 2011), for urban deformation monitoring in Hong Kong by (Qin et al., 2013), in Italian Alps corner reflectors have been used by Darvishi et al. (Darvishi et al., 2018). Corner reflectors have been integrated for monitoring of mountain glacier areas with Sentinel-1 time series by Jauvin et al. (Jauvin et al., 2019). In some recent studies (Crosetto et al., 2020) for deformation monitoring both active and passive reflectors have been deployed.
In this work two independent chains for radar scenes processing using the PSI method were compared. One from a commercial source distributed by HARRIS Geospatial Solutions (SARscape) and the other one developed by the CTTC. SARScape is available on the ENVI platform and it includes a toolbox for the processing of radar images and module for PSI to detecting/monitoring terrain displacements. The processing chain created by the CTTC for Differential Interferometric Synthetic Aperture Radar (DInSAR) was designed for the study of displacements on LGOM (Lubin-Glogow Copper Belt) area described in the next chapter. The comparison refers to slow ground displacements determined on a small network (set) of passive corner reflectors deployed for the purposes of the AS-MOW project and manufactured to be effective for SAR interferometry based both on Sentinel-1 and TerraSAR-X data.
The analyses performed have the purpose to justify the necessity of CRs deployment on the areas of strong mining influences but not having sufficient quantity of strong natural scatterers. The comparison focuses on PS "candidates" determination on two different ways, using the same SLC dataset. The comparison was done on "pixel-to-pixel" basis by identifying in the CTTC output the "PS pixels" extracted using SARScape package. The main stage of the comparison has been done on the dataset without unwanted errors of further geoprocessing: geocoding errors, interpolation grid setup, etc. The comparisons have been made on CRs as the reference objects having stable radiometric characteristics and known 3D orientation, mounted pairwise on a common pedestal for ASC and DESC orbits. The results of the comparison in LoS geometry have been next validated using precise permanent GNSS measurements. For the validation the displacements calculated on the CRs were decomposed into UpDown (vertical) and EastWest (horizontal) components.
The content of the paper is: i) the used corner reflectors description and deployment procedure, ii) the description of the PSI processing using both software packages, iii) the description of the comparison methodology, iv) the presentation of achieved results and their validation and v) discussion of the results.

AOI description
The area of interest is a rural area located in South Western Poland in a region of intensive copper mining identified as a source of important geohazard. The Legnica-Glogow Copper Belt (LGOM) is an urban-industrial area and one of the largest copper mining centres in the world and covers the area of Lubin, Polkowice, Sieroszowice and Rudna with total geological resources of 1740 million tonnes (about 75% of national resources). With the expansion of industry in this area, technical infrastructure was developed (including roads, power lines and railways), and especially housing facilities (Glogow, Legnica, Lubin, Polkowice) in Figure 1. The economic structure of the agricultural region has changed to a large extent, which quickly advanced to the rank of the most important in the country. Currently, LGOM resources are estimated at 2.3 billion tons of ore with 1.6% metal content. Deposits of considerable thickness (1-20 m) are located at a depth of several hundred meters to about 1.5 km. Approximately 27 million tonnes of copper ore are recovered annually. In addition, there are deposits of brown coal and salt. Large portions of such an area are covered by vegetation and agricultural fields, where PSI provides low density of deformation measurements.

Procedure of CR deployment
In order to guarantee a sufficient sample of the area of interest, pairs of passive corner reflectors have been installed. The artificial corner reflectors should be designed in such a way as to concentrate the energy incoming from the satellite and send it back towards the receiver. To locate the radar reflector in the SAR image, it is necessary to select the specific dimensions of the reflector, adapted to wavelength of the SAR system. The CR are passive targets with a simple geometrical shape, designed to perform a high radar reflectivity. There are many different types of reflector being used including (Garthwaite et al., 2015): spheres, cylinders, dihedrals, trihedrals, flat plates, top hats and bruderhedrals. A trihedral radar reflector is often known as CR because the reflection is facilitated by the bouncing of the incident radar wave from three mutually orthogonal plates. The orthogonality of all three plates ensures that the CR resembles the corner of a cube (and hence the name "corner reflector") with one baseplate and two 'vertical' plates. Trihedral CR have been used for many years for calibration of SAR images. They are usually constructed with metal plates, with a size that is large with respect to wavelength, and with faces oriented to maximize the energy reflected towards the radar.Two types of trihedral reflectors are used: a triangle trihedral and a square trihedral. Trihedral corner reflectors-triangle (TRC-T) and trihedral corner reflector-square (TRC-S) have different Radar Cross Section (RCS) for the same length. The RCS is a measure of the size of a target as seen by the imaging radar.
The corner reflector to be used should maximize the RCS or at least the SCR of the target to allow accurate phase measurements, while keeping the weight and the size as small as possible. The small reflectors, their radiometric features and usability for ground motion monitoring have been analyzed in details in (Dheenathayalan et al., 2017). Size of corner reflector should be chosen such that the SAR output should not get saturated due to the large amount of reflected energy. As shown by Gisinger in (Gisinger et al., 2016) SCR can by expressed by following formula: The magnitude of clutter depends on terrain type, vegetation density, soil moisture, radar wavelength, incidence angle, polarization and SAR resolution. With reference to typical clutter levels for different land cover types and dependent on incidence angle the used CRs will not encounter clutter levels at C-band greater than -10 --12 dB. In our project we used TRC-T reflectors with dimension of 100 cm. They were installed pairwise in order to provide accurate phase measurements using both ascending and descending SAR orbits. Much useful information on design, deployment and localization of bi-directional corner reflectors can be found in (Saeed and Hellwich, 2017). Following the formula of Freeman for C-band (λ = 5.6cm) the RCS for used 1m length (inner leg) TRC-T reflectors The same shape and dimensions reflectors have been mounted on Fogo volcano (Cape Verde) in 2011 by Bignami et al. (Bignami et al., 2013) taking into account the requirements to be compliant with C-band SAR: Radarsat and Sentinel-1 and X-band systems: COSMO-Sky-Med and TerraSAR-X. Corresponding SCR can be calculated as follows: Where 0.063 equals a linear value of -12dB supposed for the σC for the clutter, 4m and 14m are the range and ground resolution of Sentinel-1 SLC images correspondingly. And expressed in dB The estimation of spatial and temporal SCR for InSAR corner reflectors from SAR time series is largely discussed in (Czikhardt et al., 2022). As considered in (Qin et al., 2013), (Ferretti et al., 2007) following (Ketelaar et al., 2004) SCR is directly related to the dispersion of the phase values. For high SCR values, the relationship between SCR, phase noise, and the dispersion of the displacement measurements along the satellite line of sight (LOS) can be approximated as follows: Taking into account SCR = 238.31 we receive the value of σLOS = 0.20mm which is largely enough to lead this study.
To improve the SCR it is essential to choose the place for CR deployment free from buildings, fences, metal or dome shaped objects. Site selection should be carried out based on the scattering properties of the ground. To reduce the ground reflectivity interfering with CR, it is preferable to choose dry, soil covered field or uniform grass field.
In 2019, a network of 14 corner reflectors was deployed in the LGOM area Figure 2. Suitable locations were selected in areas far from the huge obstacles and the edges of forests if possible. Strong reflectors such as railroad tracks or light poles have been also avoided. The areas of low backscatter and low coherence have been selected.
The elevation and azimuth angle for a CR to given location can be predicted from the Sentinel-1A/B acquisition scenarios over that location. The vertical and horizontal positioning of the CR were assured by precise geodetic measurements.

SAR data used: Sentinel-1 time series
The stack of 27 SLC scenes from the ascending (02.09.2019-05.02.2020) and 26 from the descending orbit (5.09.2019-08.02.2020), in Interferometric Wide Swath mode (IW) at VV polarization have been processed. The Incidence Angle at the scene centre for both orbits 39 degrees.

METHODS
One of the most significant steps of the PSI method is the detection and selection of stable scatterers. Using multiple SAR images, we can analyse the time series of the amplitude values for each pixel in the region of interest, looking for such scatterers (Ferretti et al., 2001). Two independent processing chains have been compared and the methodology of our work is summarized with the flowchart in fig.6.
The processing parameters were chosen in such a way that the entire PSI calculation process could be as similar as possible in both cases (CTTC and SARScape). The same master image was also selected from the first registration day. In both cases the same reference point was chosen, and the same Digital Elevation Model (SRTM 90 m) was used as well as multilooking 10x2 pixels in range and azimuth (pixel footprint 40 by 28m).

PSI processing -CTTC
The CTTC's solution starts with an automatic burst extraction and precise co-registration is performed on the entire bursts stack that covers the area of interest. This is basically based on the information provided by the precise orbits associated with the images. The interferogram generation includes the removal of the topographic component. In addition to the interferograms, this software generates the coherence files, containing values of the interferometric coherence, associated to the The amplitude dispersion index has been introduced by (Ferretti et al., 2001) as the ratio between the standard deviation σA and the mean value of the image intensities in the considered time interval.
The mean amplitude is very useful for visualization of the frame of images in radar geometry, to recognize the localization of points of interest. The next step is interferogram phase unwrapping, unlike in case of SARScape processing chain. The AOI is then chosen on the generated interferogram. Phase unwrapping is a key processing step that proceeds the int2ima transformation. This is a critical processing step: the phase unwrapping errors can affect the quality of the deformation time series. This script performs the selection of the pixels to be unwrapped, which is based on a threshold on either the DA (in the case of full resolution data) or the interferometric coherence (in the case of multi-looked data). The processing approach utilized in this work employs multi-looked interferograms, thus a coherence threshold included between 0.35 and 0.4 has been used for the pixel selection. This script performs the 2D phase unwrapping (Crosetto et al., 2019), (Devanthéry et al., 2018) of each interferogram of the interferograms' stack. The input data are images, in raster format, where the sample values are included between −π and +π. The interferometric phases (not temporally ordered) are transformed in image phases, which are temporally ordered. The image phases are later transformed in deformation time series. In the next step it is employed a spatio-temporal filter to estimate and remove the atmospheric phase screen (APS) component for each of the N image processed phases. The atmospheric signal is characterized by a high spatial correlation and a low temporal correlation thus it is estimated using a spatial low-pass (SLP) filter followed by a temporal high-pass (THP) filter. Geocoding PS results are also possible. This is a key step, which is used to transform the data from the radar image geometry in a geographic or cartographic geometry for further comparison with GNSS measurements.

PSI processing -SARScape
In SARScape software all images are connected to create a network of master and slave pairs and co-registered onto the same master acquisition which has been selected manually (first registration). Then the interferograms are generated for each slave image always using the same master image. Then the Interferogram Flattening is performed using an input reference Digital Elevation Model and calculation of the amplitude dispersion index is done. The amplitude dispersion index (in SAR-Scape it is referred to as DA or MuSigma) is calculated inversely than in CTTC solution as the ratio between the mean value and the standard deviation σA of the image intensities in the considered time interval.
M uSigma = mA σA (7) A low temporal variability in backscattering intensity (i.e., high MuSigma) indicates the presence of a PS. The calculation of the amplitude dispersion index allows at the very beginning, without phase coherence analysis, to select candidates for PSC (Permanent Scatterers Candidates) points. Pixels that have a high dispersion index, and therefore have similar over time and relatively high amplitude values, are potentially good candidates for fixed scatterers. The typing of PS candidates consists in thresholding the value of the amplitude dispersion index. The amplitude dispersion index is calculated on radiometrically calibrated images. Calibration is performed automatically by the program. In the next step called PS First Inversion the first displacement (date by date value), velocity and height (correction values) related products are generated without removing any phase component due to the atmosphere. Whereas in PS Second Inversion the atmospheric corrections, related to spatial and temporal variations, are performed. Then this component is estimated and finally subtracted from the interferogram files in order to generate the final displacements. Atmospheric filtration removes most of the variation in signal propagation delay, which is mainly caused by changes in the troposphere. Atmospheric filtration includes spatial (Low Pass Filter) and temporal (Hi Pass Filter) filtration. Additionally, results can be geocoded when all PS related products (e.g. displacement velocities, residual heights, displacement time series) are projected into the cartographic system of the input "DEM file". No phase unwrapping is performed during interferometry, as Persistent Scatterers works with wrapped interferograms only. The PS inversion looks for the best linear model which best fits the interferometric measures, date by date, for each pixel (Riccardi, 2020). This means that, a velocity is associated to each pixel, which is computed based on the linear model. The time series reconstruction is then made using a sort of "temporal unwrapping" date by date. The accuracy of the time series depends on how good the velocity has been estimated during the inversion.

Comparison methodology
The main part of the methodology applied for the comparison consists in the identification the same pixels in both outputs (CTTC and SARScape) on the ASC and DESC scenes in LoS slant geometry, representing the TRC-T reflectors. The first step is finding the shifts in rows and columns between two datasets coming from two independent processing chains on the master images: DESC and ASC. The main supposition is that the one dual/twin TRC-T reflector is representing by 2 adjacent pixels, due to azimuth resolution of the SLC image (14m) Figure 5. The MuSigma and coherence values for the pixels representing the corner reflectors and their closest surroundings in LoS slant geometry issued from SARScape routine both for ASC and DESC orbits and the applied multilooking (10x2). The reflectors have been identified on Mean Amplitude image from SARScape software and corresponding pixels have been marked with the numbers accordingly to their location on the ASC and DESC scenes as indicated in the figure 6.
In every "neighbourhood window" the Vprecision parameter for pixels 1-12 has been calculated in SARScape following the formula: Where -λ is the wavelength of the imaging sensor and γ -the coherence of given pixel. Vprecision corresponds to the estimate of the velocity measurement average precision in [mm/year]. It is computed considering the acquisition temporal baseline and the multitemporal coherence (Riccardi, 2020). If for given pixel from such a window the Vprecision is marked as "NaN" value in the output from SARScape routine, there is no further calculation of MuSigma neither the coherence value in our approach. Otherwise, both parameters are calculated for the remaining pixels in the window. The values of DA, referred to as MuSigma parameter in SARScape, and coherence levels for the pixels representing the corner reflectors and their closest surroundings in LoS slant geometry, both for ASC and DESC orbits, are shown in the Figure 5 for 9 reflectors deployed on studied area. Further analyses consist in comparison of time series of the LoS displacements both for ASC and DESC orbits only for two pixels marked as 5 and 8 correspondingly for The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France CTTC and SARScape outputs.
How the deployed corner reflectors look like from point of view of their signal stability? It can be estimated from the graphs in the Figure 5 The threshold values of the MuSigma and the coherence adopted for SLC data processing were 3.2 and 0.4 correspondingly. They correspond with the values of DA and coherence adopted by CTTC. The graphs show that generally the values of MuSigma > 10 (strong and stable) and the coherence > 0.6, sometimes exceeding 0.8 (very stable) for the pixels 5 and 8 designating CRs. There are some exceptions: CR9/8 ASC -not visible, CR10/5 -small coherence, weak signal and partially CR11 and CR14. Perhaps, unfortunately the corners 9,10 and 11 are located to close to trees in forested areas? CR 14 is located on the village area where some objects can influence the backscatter. These cases must be analysed on the spot in order to identify potential external factors influencing reflectors response (intentional damage or orientation changing).
The comparison of displacements and velocities has been done on the "matched" pixels (SARScape versus CTTC) in LoS geometry which were additionally georeferenced in the output files of both packages. GNSS measurements were used here as a validation tool after the decomposition of LoS displacement estimates determined from ASC and DESC orbits into UpDown and EastWest components according to the equations: DispEW = cos θDDispA − cos θADispD cos θA sin θD cos αD − cos θD sin θA cos αA (9) DispUD = sin θD cos αDDispA − sin θA cos αADispD cos θA sin θD cos αD − cos θD sin θA cos αA where αA and αD are the azimuth angles of ascending and descending orbits, respectively. θA and θD are the look angles in ascending and descending orbits, respectively. DispA is the displacement on the ascending LOS, DispD is the displacement on the descending LOS.

RESULTS AND DISCUSSION
The results of the comparison are shown on the graphs in the Figures 7 8, 9 and 10 first as the LoS displacements on the CR5/5 and CR145. The trends for both orbits are very similar. The average difference of the displacements for ASC orbit in time between CTTC and SARScape for CR5/5 is -0.98 mm (small bias) and the standard deviation is equal 0.73mm. The average difference of the displacements for DESC orbit is 0.68 mm and the standard deviation is equal 0.79 mm. It is worth to notice that the LoS displacement for ASC orbit (about -6mm) is 3 times smaller than for DESC orbit (about -16mm). The absolute discrepancies between both estimates are smaller in case of more serious displacements. The average difference of the displacements for ASC orbit in time between CTTC and SARScape for CR145 is -4.08 mm (higher bias) and the standard deviation is equal 4.59 mm (high dispersion). The average difference of the displacements for DESC orbit for CR145 is 2.45 mm (smaller bias than ASC) and the standard deviation is equal 1.40mm (smaller dispersion than ASC). It is worth to mention that CR14 is far away from the reference point located in Polkowice, contrary to CR5 -close to it. The dispersion of values for ASC orbit is noticeable and will be  The agreement achieved between PSI, both CTTC and SAR-Scape, is very good for vertical components (UpDown) for two CRs -about 3/15mm and 5/35mm. East-West components are also convergent for the CR5, but not for CR14. It is the result of LoS displacement discrepancy for that case. The best term and the measure for PSI technique performance assessment is the displacement velocity. The absolute values of the displacements can be charged by negatively influencing factors like atmosphere. The velocity is estimated with redundant information, so is more reliable. The velocities achieved on the analysed corners are the following shown in the

CONCLUSIONS
The results achieved using two different approaches in PSI processing are consistent and mutually coherent. The deployed corner reflectors permitted to extract proper pixels in dataset to compare the LoS displacement estimates on "pixel-to-pixel" basis between CTTC and SARScape. The convergence of the values did proof that the method proposed to identify pixels on multilooked images is correct. They are also very congruent with permanent observations on close GNSS stations. Nevertheless, the intermediary results showed that the responses (in terms of dispersion of amplitude and the coherence) of some reflectors are slightly different. In further research the authors will check the physical status of the reflectors deployed on the AOI. After the assessment of reflectors' condition (orientation, presence of the objects increasing the clutter response) the second attempt of the comparison of both software packages will be made on another dataset. The special attention will be paid to the corners 9,10 and 11 exhibiting outlying values of coherence and DA. The sources of observed discrepancies will be further studied.