Exploring PAZ co-polarimetric SAR data for surface movement mapping and scattering characterization

In this contribution, we investigate PAZ co-polarimetric SAR data applicability for surface movement mapping and scattering characterization. PAZ simultaneously collects SAR imagery in both VV and HH channels. Using a small stack of PAZ data, we apply the real-valued impulse response function correlation to identify constantly coherent scatterers (CCS), separately in VV and HH, in the course of time series InSAR (Interferometric SAR) processing. The proposed method has an advantage to selecting the CCS with minimal incoherent scatterer inclusion and exact radar location, which can eventually lead to the precise deformation time series estimations of all CCS, and a high-precision surface movement map. Moreover, we apply the co-polarimetric phase difference (CPD) method to classify the CCS in terms of scattering mechanisms which provides a new attribute to every individual CCS. We recognize the sibling pairs by both thresholding the spatial distance between any two CCS observed separately in VV and HH, and using common scattering characteristic as a new criterion. The deformation estimates of sibling pairs are used to reduce the biases in the deformation estimates of every ground target. The proposed methods are demonstrated in a test site, in the northern part of the Netherlands, using 10 co-polarimetric SAR data acquired between September 2019 and April 2020. The results show that 83.5% sibling pairs behave a linear deformation trend over time, and that the other pairs show a correlation between their deformation and temperature, and the sibling pairs with the surface, dihedral, volume scattering mechanisms account for 62%, 12% and 26%, respectively. We conclude that by combining data from VV and HH polarization as siblings, PAZ co-polarimetric SAR data are highly suited to map surface changes and characterize surface


Introduction
Radar satellite PAZ ('peace' in Spanish), launched in February 2018, is a promising X-band SAR (Synthetic Aperture Radar) mission of Spain, and by definition is able to regularly deliver radar images with high spatio-temporal resolution and millimeter precision (Palma et al., 2010). It is operated in a constellation with TerraSAR-X (Werninghaus and Buckreuss, 2009) and TanDEM-X (Krieger et al., 2007) on the same orbit. A successful cross-sensor interferometry generation between PAZ and TerraSAR-X SAR data at Burgan, Kuwait, which was reported on 28 February, 2019 (eoportal, 2019), demonstrates a giant leap of satellite InSAR (Interferometric SAR) techniques and boosts the further PAZ SAR data exploration for mapping surface elevation and movement. Moreover, as PAZ possesses co-polarimetric channels, i.e. VV and HH, the scattering mechanism(s) of any observed ground target can be disclosed. Yet, as of now, few works in the literature related to the PAZ-based applications to surface movement mapping and scattering characterization are published, probably due to the limited PAZ data availability to the public or data collection protocol. As availability of PAZ data is increasing, it is crucial to investigate the potential of PAZ data, and expand the scarce experimental material available.
In general, given a stack of SAR data, surface movement can be retrieved using multi-temporal InSAR (Interferometric SAR) techniques, MTInSAR for short. Examples are persistent scatterer interferometric (PSI) (Ferretti et al., 2001;Osmanoglu et al., 2011), small baseline permanent scatterers (SBAS) (Berardino et al., 2002;Mora et al., 2002;Lanari, 2003;Tizzani et al., 2007;Lanari et al., 2010;Bateson et al., 2015) SqueeSAR (Ferretti et al., 2011), and CAESAR (Fornaro et al., 2014). In essence, MTInSAR primarily identifies all coherent scatterers by using an appropriate and ad-hoc SAR image combination and later on extracts the movement (deformation in other words) of those scatterers. An MTInSAR technique applies distinct method to select coherent scatterers. For example, the normalized amplitude dispersion (NAD) (Ferretti et al., 2001) is widely used for selecting constantly coherent scatterers (CCS), i.e. persistent scatterers, during PSI processing. The signal-to-clutter ratio (SCR)  and amplitude thresholding  and supervised selection (Hanssen and van Leijen, 2008), can be considered as alternative selection methods. Furthermore, the selection methods using phase characteristics (Hooper et al., 2004), maximum likelihood estimation (Shanker and Zebker, 2007), and top eigenvalue of coherence matrix (Navneet et al., 2017), were proposed for selecting constantly high-coherent scatterers. However, given a small stack of SAR images (e.g. ⩽12 SAR images), the selected CCS using some of these methods may contain a number of incoherent scatterers and may not be sensitive to identify sidelobes against mainlobes. The incoherent scatterers and the scatterers at sidelobes in the entries do not enhance the stableness and robustness of the unwrapping network in time and space. On the contrary, these scatterers weaken the stableness and robustness of the network and degrade the CCS deformation estimates. To circumvent this problem, a CCS selection method using an impulse response function correlation will be introduced and demonstrated in this study.
The presence of VV and HH polarization channels, offers an opportunity to attribute the scattering mechanisms to individual selected coherent scatterers. We distinguish three main scattering mechanisms, surface scattering, dihedral scattering and volume scattering mechanisms (Freeman and Durden, 1998;Lee and Pottier, 2017). A plethora of works have shown the methods of scattering characterization, mostly for land-use classification and land-use change identification, e.g. coherent and incoherent polarimetric decomposition and machine/deep learning, (Van Zyl, 1989;Cloude and Pottier, 1996;Cloude and Pottier, 1997;Pottier and Lee, 2000;Lee et al., 2004;Yamaguchi et al., 2011;Leinss et al., 2014;Shimoni et al., 2009;Mullissa et al., 2017). For instance, Leinss et al. (2014) proposed to straightforwardly employ the copolarimetric phase difference (CPD) between the VV and HH channels to categorize the coherent scatterers in terms of scattering mechanisms (i.e. surface, dihedral and volume scattering), and quantitatively detect snowfall according to the dependence of the CPD on fresh snow in the sequel. The use of the CPD is not limited to land use/change identification. As the outcome of using the CPD, the identified scattering mechanisms of the CCS can be used to assist for improving the CCS deformation estimates by integrating the MTInSAR and CPD measurements. However, the research in this regard is not yet fully explored.
The objective of study is to explore the usability of PAZ copolarimetric SAR data to map surface changes and characterize surface features. The study is illustrated with an application in the northern Netherlands. This paper is organized as follows. Section 2 is dedicated to the methods of constantly coherent scatterer selection (especially for a small stack of SAR data), and scattering characterization and deformation estimate improvement. Section 3 presents all data used over a test site. The results of the test site are presented in Section 4, followed by discussion and conclusions in Section 5.

Methods
This section first reviews the constantly coherent scatterer selection by using a real-valued impulse response function correlation in Section 2.1, and presents a quality indicator of kinematic deformation time series estimates, i.e. ensemble temporal coherence, in Section 2.2. Next it briefly introduces the basis of co-polarization phase difference (CPD) in Section 2.3, and then elucidates how to precisely classify scatterers in terms of scattering mechanisms in Section 2.4. Finally, it describes how to improve deformation estimates and select sibling pairs in Section 2.5 and time series modelling in Section 2.6.

Constantly coherent scatterer selection
An MTInSAR, such as PSI, employs a stack of SAR images to identify all constantly coherent scatterers (CCS) and generates the associated kinematic deformation time series. A CCS has relatively stable phase measurements in time irrespective of temporal and spatial decorrelation. Buildings, bridges and bare rocks are all good examples of CCS in SAR imagery. Every resolution cell in a SAR image is recorded as a single complex number (SLC), and can be expressed both in polar form and in modulus-argument form In polar form, the complex number Z is composed of two parts, the real part X and imaginary part Y. j denotes the imaginary unit. In modulusargument form, Z is shaped by the amplitude |Z| = (X 2 + Y 2 ) 1 2 , and the SLC phase ψ = atan Y X , (ψ ∈ (− π, + π]), depicted in Fig. 1a). Since ψ is uniformly distributed between − π and +π (Bamler and Hartl, 1998), it is unrealistic to determine the potential CCS using ψ information. Instead, |Z| information is often employed to recognize potential CCS. For instance, the normalized amplitude dispersion (NAD) (Ferretti et al., 2001), which is a ratio of the standard deviation σ a of SAR amplitude time series over its mean μ a , expressed as NAD = σ a /μ a ≈σ ψ , has been widely used to successfully detect the potential CCS in a relatively long time series (e.g. more than twenty SAR images) (Cuenca et al., 2013;Chang et al., 2014;Wang et al., 2018). However, it is not always optimal and reliable to extract all CCS in a relatively short time series (e.g. ⩽12 SAR images), as the selected CCS may contain a number of incoherent scatterers. Moreover, it is insensitive to detect sidelobes and determine the exact CCS location in radar coordinates. Therefore, we review an impulse response function (IRF) correlation method to select the CCS with minimal incoherent scatterers inclusion and an exact radar location, especially in a short time series.
In a one dimensional domain, the complex impulse response function (IRF) correlation at position x can be formed as, c.f. (van der Torren, 2011) where s(x) is the received signal in an SLC form at position x, f(κ) is the IRF, and κ is a dummy distance variable. (*) is the conjugation operator, and w(κ) is the weight function that is used to form a weighted window. The weighted window merely includes the local neighbourhood and mainlobes and excludes sidelobes. The upper and lower limits c and b are the maximum and minimum distance w.r.t. the position x, which can be customized by the user. The term ∫ c b s(x +κ)⋅f * (κ)⋅w(κ)dκ in the numerator represents the correlation of the received signal and IRF, given a weight function w(κ), later denoted as ρ(x) (Cumming and Wong, 2005).
in the denominator is applied to normalize ρ(x). When the weight function w(κ) is defined as the amplitude of f(κ), i.e. |f(κ)|, in a weighted window matrix, the position with a zero weight has no influence on the scatterer's signal at position x, whereas the position with a higher weight has more influence. The real part of Eq.
(2) is recommended by (van der Torren, 2011) to further compute the correlation of the received signal and the ideal signal with an equal phase value. The real valued IRF correlation is obtained by If the signal at position x is considerably influenced by the surrounding clutter, ρ IRFR (x) decreases, as the phase of the received signal largely diverts from the average phase within the weighted window matrix. The highest real valued IRF correlation occurs in the center of a point-wise scatterer. Hence the exact scatterer location in radar coordinates can be determined by computing the local maxima (Perissin, 2006).

Ensemble temporal coherence
The temporal model of a kinematic deformation time series per CCS is often not known completely. To facilitate the interferometric phase unwrapping processing, a linear function of time is often used as a temporal model for all ρ IRFR (x)-based selected CCS, c.f. (Chang and Hanssen, 2015).
To assess the phase time series model misspecification and temporal coherence, the ensemble temporal coherence (ETC) is proposed by (van Leijen, 2014), where m is the number of SAR images, m − 1 is the number of the interferometric or slave images. The unwrapped interferometric phase estimate and model value (e.g. using a linear function of time) between the master SAR image and ith slave image, at position x, are denoted by φ i x and ϕ i x,model . The ETC ranges between 0 and 1. A higher ETC value implies a better match between the estimate and model.

Co-polarized phase difference (CPD)
Given two SLC in the co-polarised VV and HH channels, denoted by Z VV and Z HH respectively, which are acquired simultaneously, their interferogram, Z c , is a complex multiplication of Z VV and the conjugate of Z HH , where |Z c | indicates the amplitude of the interferogram, while ϕ c decipts the phase difference between VV and HH channels, and ranges between − π and + π, named as the co-polarized phase difference (CPD). |Z VV | and ψ VV are the amplitude and phase of the complex number Z VV , while |Z HH | and ψ HH are the amplitude and phase of the complex numer Z HH .
Theoretically, abide by the classical theory of polarimetry in (Lee and Pottier, 2017), three scattering mechanisms can be categorized by ϕ c : 1) surface scattering ϕ c = 0; 2) dihedral scattering ϕ c = π; and 3) volume scattering ϕ c ∼ unif(− π, +π) (unif() represents the notation for a uniform distribution). Description of these three scattering mechanisms is presented in Appendix A.

Scattering characterization
In practice, a SAR resolution cell (where a single CCS may be located) is described by a superposition of these three scattering mechanisms, due to the meter-level spatial resolution of SAR imagery, and SAR imagery contains speckle noise. Hence, it is not straightforward to characterize pure surface, volume and dihedral scatterers and associate an appropriate scattering mechanism with every resolution cell in a single CDP ϕ c image.
As temporal averaging is a way to mitigate the noise in ϕ c measurements, and the CPD coherence γ c is indicative of the correlation level of VV and HH channels, we propose to use the temporally-averaged CPD, ϕ c , which is weighted by γ c , to characterize scatterers. Suppose there are m CPD, ϕ i c (∀i, i ∈ [1, m]), yielded from m SAR image acquisitions in the co-polarimetric channels, ϕ c is formulated as where Here γ i c ranges between 0 and 1. The large value of γ i c , corresponds with a high correlation between VV and HH SLC images. Z i c,q represents the CPD value for the qth resolution cell in the ith acquisition date. Z i VV,q and Z i HH,q present the complex number of the qth resolution cell in VV and HH in the ith acquisition date, respectively. N is the total amount of resolution cells that are used to compute γ i c . For instance, when γ i c = 1, Assuming an equal uncertainty of individual ϕ i c measurement, the estimate uncertainty of ϕ c is reduced by 1/ ̅̅̅̅ m √ . The precision of scattering characterization result can be improved by the increase of the number of acquisitions m, when all CPD coherence can retain an allowable range (e.g. > 0.7) over time.

Deformation estimate improvement and sibling pair selection
In response to natural and human activities, every single scatterer on Earth may retain the same, move towards or away from a satellite SAR sensor which captures this scatterer in radar imagery, namely deformation. When a SAR sensor has multiple polarization channels such as VV and HH, multi-measurements per scatterer are obtained. As such, the kinematic deformation time series estimates per scatterer can be improved. To be more specific, assuming a CCS is observed by both VV and HH channel, its associated deformation estimate in the ith acquisition is denoted as ŷ i VV for VV, and ŷ i HH for HH, we improve the deformation estimate in the ith acquisition (i ∈ [1, m], m is the total number of the acquisitions) by computing their weighted average, where ω i VV and ω i HH are the allocated weight of the deformation in VV and HH, respectively. ω i VV and ω i HH can be shaped by the CCS coherence in the ith interferogram. The large weight is assigned for the CCS when the coherence is high. The quality of the deformation estimates, ŷ i VV and ŷ i HH , is highly dependent on the uncertainty of SAR raw observation and InSAR time series processing errors (Teunissen et al., 2005). If ω i VV and ω i HH are assigned as a unit weight, the variance of the deformation estimate is σ 2 yy for both VV and HH, the variance of y i achieves a 50% decrease and reaches to σ 2 yy /2 according to the error propagation law. Identifying a pair of CCS, which separately appear in VV and HH, bear the deformation, ŷ i VV and ŷ i HH , and both represent a common ground target, is crucial to the deformation estimate improvement by using Eq. (7). We label such a pair of CCS as a sibling pair, as this pair represents a common ground target and is much alike in scattering mechanism. If two CCS representing two distinct ground targets are mistakenly grouped as a sibling pair, we may bias their deformation estimate by using Eq. (7). Therefore, to exclude the occurrence of this unacceptable case, we propose to first detect the potential sibling pairs by thresholding the euclidean distance between them, and then screen more qualified sibling pairs by thresholding the CPD standard deviation . If the scattering mechanism of a CCS does not strongly vary over time, ϕ i c (i, ∈ [1, m]) would not divert much from ϕ c , resulting in a small σ ϕ c . Moreover, with the increase of the total amount of SAR image acquisitions in both VV and HH, σ ϕ c drops. Fig. 2 illustrates the CPD standard deviation σ ϕ c changes of 1000 simulated individual CCS in response to the number of SAR image acquisitions (m ∈ [10, 300]), given 0.2 radians noise level. Each colored line shows the fluctuation of σ ϕ c values as the number of SAR acquisitions increases, for a simulated CCS. The simulation results show that having more SAR acquisitions, the range of σ ϕ c shrinks, which signifies the CPD standard deviation (σ ϕ c ) values of more simulated CCS are close to their average (or expectation), given a larger amount of acquisitions. For instance, the inset shows the CPD standard deviation σ ϕ c histograms of the 1000 simulated CCS, when the number of SAR acquisitions m equals to 10, 100, 200, and 300 (indicated by the dashed line in the main figure). The spread of the histogram distribution is wider when m gets smaller, and the most probable value of σ ϕ c is ∼ 0.12 radians, ∀m.

Sibling pair deformation time series modelling
In practice, the apriori knowledge on the sibling pair temporal behavior is often incomplete. Then we empirically parametrize the sibling pair deformation time series by a linear function of time. Because a linear function of time, which facilitates the phase unwrapping but may not always be the most suited option to all sibling pairs, we evaluate the default linear model and determine the suitable model per sibling pair by hypothesis testing (Teunissen et al., 2005;Chang and Hanssen, 2015). We define an alternative (physically-realistic) deformation model including temperature-related parameter as an alternative hypothesis H 1 , against the linear model as a null hypothesis H 0 . As such, H 1 considers the thermal dilation of ground scatterers such as concrete buildings and steel-made structures.
The functional and stochastic models for H 0 and H 1 can be separately expressed as where E{} denotes the expectation operator, and D{} denotes the dispersion operator. B T represents the temporal baselines between all slave images and master, v represents the unknown deformation constant velocity, B temp is the temperature baseline which is the temperature difference between the slaves and master, and η is the temperaturedependent parameter. y is the m × 1 deformation measurement matrix of a sibling pair, [y 1 , y 2 , ⋯, y m ] T . The measurement noise of y is described by the variance-covariance matrix Q y y (= σ y y 2 R y y where σ y y 2 is the variance of unit weight and R y y is the cofactor matrix). m is the number of observations. To determine the suitable deformation model from H 0 and H 1 , we follow a three-step procedure: 1) we first apply the overall model test to the null hypothesis H 0 , using the test statistic T =ê T 0 Q y y − 1ê 0 ∼ χ 2 (m − 1,0). Here ê 0 is the column matrix of the error estimations between y and its estimate ŷ in H 0 . We define the level of significance α by calculating α = 1/2m, and compute the corresponding critical value K α (Chang and Hanssen, 2015). If T < K α , H 0 is determined to be the suitable deformation model, and if T > K α , H 0 is rejected and H 1 is further assessed. 2) Next, we apply a least-squares estimation to estimate the unknown parameters in H 1 . 3) Finally, we use the quality metrics such as the posterior variance, to determine whether H 1 has a better agreement with the sibling pair measurements, compared to H 0 , and can be considered as the suitable deformation model.

Data description and test site
Various data are employed in this study. The basic information on these data is presented in this section. The test site was chosen as it is a complex region and covers highly dynamic and decorrelated areas.
SAR data: 10 Spanish PAZ SAR image acquisitions which were acquired between September 2019 and April 2020 are employed. It covers 15 × 50 km 2 , including some areas in the province Friesland, the islands Terschelling and Ameland and Waddenzee, which is highlighted in red in Fig. 3. Each acquisition has two polarization channels, VV and HH, and both in stripmap mode with 3 × 3 m resolution in range and azimuth direction. The repeat cycle of PAZ mission is 11 days, however, the data used have 22 days time difference between adjacent acquisitions, as shown in Table 1 with the temporal baseline B T information.
Tide gauge data: The tide gauge data are collected at two local stations, Nes and Holwerd, which are closest to the PAZ scene (Fig. 3). The tide height per station is measured every ten minutes, w.r.t. the NAP (Normaal Amsterdams Peil, the Dutch vertical datum). Fig. 4 depicts the tide height observations at Nes and Holwerd on the PAZ SAR acquisition dates. Fig. 4 shows that the tide has a semi diurnal pattern, with two high tides and two low tides almost each day. The tide cycle is approximately 24 h 20 min (2 × 12 h 10 min). The extra 20 min compared to the 24 h per solar day, leads to the tide pattern slips day by day. As every two adjacent SAR acquisitions have 528 h (22 days × 24 hours) time difference, the tidal dynamics in time behaves oppositely, e.g. on 28 September in red and 20 October 2019 in black in Fig. 4. Moreover, the acquisitions with (a multiple of) 44 days have the tidal dynamics in concert, to some extent. An example is that the tidal dynamics on 11 November and 25 December 2019, in blue and yellow respectively, reach the high and low tides at roughly same moments.
Air temperature data: The daily temperature information on the SAR acquisition dates is recorded by the Leeuwarden station (nr. 06270, 53.22 • N, 5.75 • E), see Table 1. The maximum temperature difference is 12 • in the past eight months. The Leeuwarden location is situated in the suburb of the city of Leeuwarden, indicated by the white star in Fig. 3.
AHN3 data: The airborne and terrestrial laser scanning devices are used to produce the DEM (digital elevation model) and DSM (digital surface model). Particularly, the AHN3, referred as to Actual Height Netherlands, version 3, is a laser-scanning-measured height product, for the entire Dutch territory, which has 5 cm height uncertainty and 50 cm spatial grid size (AHN3, 2020; van der Zon, 2013). The AHN3 measurement campaign took place between 2014 and 2019. An example is shown in Fig. 7a).
Test site: The soil map, as the backdrop map in Fig. 3), indicates that roughly 30% area is covered by water bodies, e.g. the Wadden Sea, for a full scene of PAZ SAR image shown in red. Most of the inland areas are comprised of light clay (in dark green), heavy clay (in cyan) and heavy gravel (in light green), while the two islands (Terschelling and Ameland) are mainly covered by sand in yellow. Some reclaimed salt-marsh areas along the Wadden sea shoreline (Ven, 2004), are indicated as water, but they fall dry twice during low tide and are inundated during high tide.

Results
The results of the test site are presented in a step-wise manner as per SAR data processing procedure, in this section.

Interferogram generation
We used Doris (v4.10 beta version, (DEOS, 2012)) to generate two stacks of 9 interferograms in VV and HH from the 10 PAZ SAR image acquisitions (see Table 1). As mentioned before, the PAZ SAR data has an equivalent data format w.r.t. TerrSAR-X data, merely readinput.cc, slcimage.cc and rpocessor.cc in this Doris version are modified to accommodate the new syntax of PAZ data.
To inhibit the high spatio-temporal decorrelation and considering the large spatial extent of water bodies in the test site, a SAR image with relative short temporal and spatial perpendicular baseline to all other SAR images and a low tide height, is preferred to be selected as the reference image. Knowing the spatial and temporal baseline listed in Table 1, and all SAR images were illuminated the whole area at 17:18 h indicated by the dashed black line in Fig. 4, and the tide height on 11 November 2019 at Nes and Holwerd retains low, we then deliberately selected the SAR image acquired on 11 November 2019 as the reference image for VV and HH modes. There were 9 interferograms in VV and 9 interferograms in HH separately generated and registered to the reference SAR acquisition. The AHN3-derived DEM, as shown in Fig. 7a), was used to estimate topographic phase components from interferograms. According to the AHN3 data, the PAZ-data-covered area is flat and does not have a great difference in elevation.
Note that we defined a 13 × 22 km 2 as the area of interest (AOI), marked by the red rectangle in Fig. 3), for further InSAR time series analysis. As this AOI embraces large water bodies and agricultural regions and reclaimed areas (i.e. salt marsh, indicated by the black rectangle in Fig. 7 small part of Ameland island, it is challenging to detect the CCS and estimate the kinematic deformation time series. In this endeavour, we applied a redundant network instead of a Delaunay network (Kampes, 2005) to obtain more redundant spatial arcs for spatial unwrapping, we loosed the maximum arc length in the network to avoid isolated network clusters, and tightened the threshold for the potential CCS selection.

CCS selection
Since for both VV and HH channels, only a small number of in-  terferograms is available, we employed the real-valued IRF correlation method (using Eq. (3)) to extract the CCS. In total 14944 CCS in VV and 17935 CCS in HH with ρ IRFR ⩾0.60 are extracted, and none of them are located in the Wadden sea area. To demonstrate that using the realvalued IRF correlation method for selecting high-quality CCS is highly suitable for short InSAR time series as compared with using a classical and standard method (i.e. the NAD (normalized amplitude dispersion)), over a complex and decorrelated region, we discuss the quality of the real-valued IRF correlation-and NAD-derived kinematic deformation time series results of the CCS, by means of the ensemble temporal coherence γ t (using Eq. (4)  and NAD (normalized amplitude dispersion) in VV and HH. Using the NAD method, with the threshold NAD ⩽0.23, there are 12877 CCS in VV and 20571 CCS in HH. A handful of falsely selected CCS, accounting for ∼ 0.27% of 12877 CCS in VV and ∼ 0.28% of 20571 CCS in HH, are located in the Wadden sea area, as they are coincidentally keeping constantly stable and high amplitude over time. Fig. 7b) illustrates a multi-reflectivity map in VV channel in radar coordinates, which shows that some scatterers in the Wadden sea area have comparably high amplitudes w.r.t. the scatterers in inland and islands. Those scatterers in the Wadden sea area are probably mainly reflected by the offshore drilling rig (TES-01) (Duin et al., 2006) and sandbanks. For VV mode, using the real-valued IRF correlation method, 77% CCS have ρ IRFR ⩾0.90 and 34% CCS have γ t ⩾0.90 as shown in Fig. 5a).
The rightmost and bottommost insets of Fig. 5a) present the scatterer population density w.r.t. ρ IRFR and γ t separately, in which monotonously increasing trend is manifested. It implies that γ t improves as ρ IRFR rises.
Using the NAD method, 49% CCS have NAD ⩽0.20 and 5% CCS have ρ IRFR ⩾0.90, as shown in Fig. 5b). Fig. 5b) also shows that albeit all NAD values are already small enough to assure the scatterers possessing stable amplitude and phase time series, the result quality of the corresponding kinematic deformation time series cannot be all guaranteed, as γ t presents a climax at 0.75 and considerably population between 0.30 and 0.80 (see the bottommost inset). Likewise, for HH mode, 78% CCS have ρ IRFR ⩾0.90 and 47% CCS have γ t ⩾0.90 as shown in Fig. 5c). In contrast, 46% CCS have NAD ⩽0.20 and 5% CCS have γ t ⩾0.90, and γ t mainly spans between 0.60 and 0.75. With this evidence, using the realvalued IRF correlation method yields better-quality results of the selected CCS, in comparison with using the NAD method.

CPD-based scattering characterization of CCS
In order to attribute the scattering mechanism to every single CCS, we used 10 CPD measurements of each CCS, ϕ i c (here i is the PAZ SAR acquisition date), to compute the temporally-averaged CPD ϕ c (using Eq. (6)). Figs. 6a) and 6b) show the temporally-averaged CPD ϕ c histogram in polar coordinates for 14944 CCS in VV and 17935 CCS in HH. By accepting 2σ n (2 × 0.2 radians) noise level in ϕ c at the confidence level 95%, we categorized the CCS into three classes, those are surface scattering (|ϕ c |⩽2σ n ), dihedral scattering (|ϕ c |⩾π − 2σ n ), and volume scattering (|ϕ c | > 2σ n &|ϕ c | < π − 2σ n ). As a result, the CCS with surface, dihedral and volume scattering characteristics, account for 30%, 24% and 46% in VV, 36%, 17% and 47% in HH. It implies that the vertically polarized radar signal is more sensitive to detect dihedral scatterers, while horizontally polarized radar signal is a little bit more sensitive to detect surface and volume scatterers.

Sibling pair identification and improvement on deformation estimates of CCS
Having this additional information on the scattering characterization of individual CCS in VV and HH, we can identify sibling pairs by means of the proposed strategy in Section 2.5. To be more specific, we first extracted the potential sibling pairs if any pair of two CCS obtained from VV and HH respectively have ⩽1 m spatial distance, and then we computed σ ϕ c of the potential sibling pairs and defined the more qualified sibling pairs that bear σ ϕ c ⩽0.3 radians. Eventually, 4150 sibling pairs were selected in the AOI. The associated deformation estimates of the sibling pairs were improved by using Eq. (7). Here we applied an equal weight, i.e. 1, to ω i VV and ω i HH (i ∈ [1, 10]), by accepting an allowable estimation bias of y i . This bias is induced by the difference in coherence in VV and HH in the ith acquisition, whose absolute values range between 0.01 and 0.14. Fig. 7c) gives the line-of-sight deformation velocity map of the 4150 sibling pairs. Most of the sibling pairs behave relatively stable between September 2019 and April 2020, w.r.t. the spatial reference point denoted by the magenta cross (53.2804 • N, 5.6602 • E). Most of them have [− 5, 5] mm yr − 1 deformation velocity rates. Note that all CCS deformations in VV and HH were aligned to this spatial reference point before the sibling pair selection. The spatial reference point noise was mitigated by applying the Shenzhen algorithm (Chang and Hanssen, 2015). There are 217 sibling pairs (with 0.6⩽ρ IRFR < 0.7) in the salt marsh region (shown as the black rectangle), where the earth dikes are built. The signal reflectance strength of the earth dikes seems fairly good, see Fig. 7b). Among these 217 sibling pairs, 47 sibling pairs (out of 127 sibling pairs over the entire AOI) show an uplift trend with > 5 mm yr − 1 , while the rest have the linear deformation velocity between − 5 and 5 mm yr − 1 . Fig. 7d) shows the classification map of the sibling pairs in terms of scattering mechanisms. 2577, 475, and 1098 sibling pairs are characterized as surface, dihedral and volume scattering respectively. The dihedral scatterers are primarily detected in the civil infrastructure, such as buildings and lamp poles. In the salt marsh region, in the black rectangle, for instance, there are 144 volume scatterers, out of 217. Along the south coastline of the Ameland island, the sibling pairs are classified as volume and surface scatterers, see Fig. 7d).

Examples of sibling pairs deformation time series modelling and interpretation
We used the hypothesis testing (see Eq. (8)), to decide the suitable deformation model per sibling pair. The temporal baseline B T information can be found in Table 1, and the temperature baseline B temp can be obtained by computing the temperature difference of the slave T slave temp and the master acquisition on 11 November 2019, T master temp − T slave temp . We initially defined the stochastic model Q y y ( = σ y y 2 R y y) of deformation time series of sibling pairs, as a diagonal matrix, in which no temporal correlation between deformations is assumed, i.e. R y y = diag(1, 1, ⋯, 1), and the variance of unit weight, σ y y 2 , was assigned to 2 2 mm 2 . We defined the level of significance α, by α = 1/2m = 0.05, then the corresponding critical value K α=0.05 = 17. The results show that 83.5% remains in H 0 , and the others accept H 1 in which the estimated temperature-related parameter η ranges between − 0.8 and + 1 [mm/ • C]. Figs. 8b), c) and d) separately illustrate the line-of-sight deformation time series of 3 dihedral scattering sibling pairs on a same single building, in the town of Sint Annaparochie, 4 surface scattering sibling pairs in front of a farmhouse, in Stiens, and 8 volume scattering sibling pairs on the salt marsh region. The location of the 3 dihedral scattering, 4 surface scattering and 8 volume scattering sibling pairs is depicted on the left zoomed-in images of Figs. 8b), c) and d), as well as indicated in the red, green and blue squares respectively in Fig. 7b). These three zoomed-in images also separately show the dihedral, surface, and volume scattering sibling pairs in red, green and blue, according to the classification results c.f. Fig. 7d). The suitable deformation model per sibling pair (SP) is shown as the colored dashed line. For instance, the dihedral scatterers SP1, SP2, SP3, show a temperature-related pattern in time. The temporal deformation fluctuation of these three sibling pairs positively resonates the temperature variation in time (see the local air temperature of all acquisitions in Fig. 8a)). The surface scattering SP4, SP5, SP6 and SP7 are nearly stable, and have small-scale deformation dynamics. The volume scattering sibling pairs SP8-SP14 behave unstable over time with a 23.4 mm yr − 1 maximum upward movement. Table 2 lists the unknown parameter estimates (deformation constant velocity v and temperature-dependent parameter η) and the quality indicators of the results (the posterior variance σ 2 of the suitable deformation model and the posterior variance ratio of H 1 and H 0 , σ 2 1 /σ 2 0 ), of the suitable deformation model, for these 14 sibling pairs. It shows that SP1, SP2, SP3, SP8, SP10, SP11, SP12, and SP14 follow the alternative hypothesis H 1 . Their posterior variance ratios between H 1 and H 0 (σ 2 1 /σ 2 0 ) are all smaller than 1, which implies that the quality of their results in H 1 is improved compared to H 0 . SP4, SP5, SP6, SP7, SP9, SP13 remain the null hypothesis H 0 , and have σ 2 1 /σ 2 0 > 1 which signifies that the additional temperature-dependent parameter is unnecessarily included in the functional model as it degrades the quality of the results, for those sibling pairs.

Discussion and conclusions
In this study we investigated and demonstrated the usability of PAZ co-polarimetric SAR data to map surface movement and categorize scattering characteristics. An IRF correlation method is recommended to be used to select constantly coherent scatterers (CCS) in InSAR time series processing, as it yields the CCS with exact radar location and maximum exclusion of incoherent scatterers and scatterers at the sidelobes. The results showed that 77% and 78% CCS in VV and HH respectively hold ⩾0.9 real-valued IRF correlation values, possibly benefiting from good coherence of all interferograms in both VV and HH channels. The results also showed that more CCS that were selected by the real-valued IRF correlation method have large ensemble temporal coherence, than those that were selected by the NAD method. Hence, the real-valued IRF correlation method leads to an increase in the quality of The co-polarimetric phase difference (CPD) method is proposed for categorizing the CCS in terms of scattering mechanisms. The results showed that by predefining an allowable tolerance/noise range, i.e. 0.4 radians for the test site, and using the temporal averaged CPD, we can achieve a reliable CPD-based classification. In addition, the results showed that a higher percentage of CCS in VV are classified as dihedral scatterers (24%), while a highly higher percentage of CCS in HH are classified as surface scatterers (36%) and volume scatterers (47%). This means that the radar signals have distinct sensitivity to various scatterer classes, in VV and HH channels. This clearly enriches the scatterer characteristic knowledge of CCS. The study further showed that the CPD method is able to categorize surface, dihedral and volume scatterers, but as yet it is unable to fine-classify e.g. to distinguish dipole, complex structure, random surface, and Bragg surface scatterers. It will require further research to obtain such a fine-classification, to which the coherent/incoherent polarimetric decomposition can be resorted to. Possibly, machine learning methods can be useful here, although the CPD-derived classification results may be affected by noise, leading to biased results.
The identified scattering mechanisms of all CCS in VV and HH are, for the first time, being used to assist in selecting sibling pairs. As such, only the pair that contains a constantly coherent scatterer in VV and its counterpart in HH, that also have a small distance in space and holds a common scattering characteristic, is identified as a sibling pair. Such a pair offers a geometric and geophysical CCS linking between the different polarimetric channels. We note that in our study area, we identified 4150 sibling pairs among the 14944 CCS in VV and the 17935 CCS in HH, with every sibling pair having deformation time series estimations in both VV and HH. As the result further showed that the weighted average deformations of every sibling pair yielded the improved estimations of deformations, these have a smaller posterior variance of the deformations as compared to the posterior variance of deformations in either VV or HH. The improved estimations in turn lead to an improvement of the kinematic time series, its interpretation and modelling.
As concerns the deformation time series modelling, we tested two deformation models, both including linear deformation velocity, and  Fig. 7b). The classification result is also superimposed on these three zoomed-in images, where the dihedral, surface, and volume scattering SP are respectively shown in red, green and blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 2 Unknown parameter estimations, posterior variance σ 2 , of the suitable deformation model, for the 14 sibling pairs. σ 2 1 /σ 2 0 indicates the posterior variance ratio between the alternative model H 1 and the default model H 0 . SP nr.
Model either including or not a temperature-dependent parameter. The results for this testing showed that 83.5% of the sibling pairs behaved as a linear deformation trend over time, and that the other sibling pairs had some degree of correlation between the deformation time series and temperature variations over time. This correlation is manifested by the temperature-dependent parameter ranging [− 0.8, +1] mm/ • C. These findings can guide us to further hypothesize the driving mechanisms and early detect spatio-temporal anomalies in deformations.
In conclusion, this paper has shown that PAZ co-polarimetric SAR data are highly suited to map surface movement and characterize ground targets in complex areas, in particular by taking siblings from VV and HH polarization for data combination. From the case study we observed that when using the co-polarimetric phase difference (CPD) method for categorizing the CCS in terms of scattering mechanisms, a reliable classification is obtained.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.