Bridging Loss-of-Lock in InSAR Time Series of Distributed Scatterers

We introduce the term loss-of-lock to describe a specific form of coherence loss that results in the breakage of a synthetic aperture radar interferometric (InSAR) time series. Loss-of-lock creates a specific pattern in the coherence matrix of a multilooked distributed scatterer (DS) by which it may be detected. Along with identification, we introduce a new DS processing methodology that is designed to mitigate the effects of loss-of-lock by introducing contextual data to assist in the time-series processing. This methodology is of particular relevance to regions that suffer from severe temporal decorrelation, such as northern peatlands. We apply our new method to two subsiding cultivated peatland regions in The Netherlands which previously proved impossible to monitor using DS InSAR techniques. Our results show a very good agreement with in situ validation data as well as spatial correlation between regions and the natural terrain.

Bridging Loss-of-Lock in InSAR Time Series of Distributed Scatterers Philip Conroy , Simon A. N. van Diepen , Freek J. van Leijen , Member, IEEE, and Ramon F. Hanssen , Senior Member, IEEE Abstract-We introduce the term loss-of-lock to describe a specific form of coherence loss that results in the breakage of a synthetic aperture radar interferometric (InSAR) time series.Loss-of-lock creates a specific pattern in the coherence matrix of a multilooked distributed scatterer (DS) by which it may be detected.Along with identification, we introduce a new DS processing methodology that is designed to mitigate the effects of loss-of-lock by introducing contextual data to assist in the timeseries processing.This methodology is of particular relevance to regions that suffer from severe temporal decorrelation, such as northern peatlands.We apply our new method to two subsiding cultivated peatland regions in The Netherlands which previously proved impossible to monitor using DS InSAR techniques.Our results show a very good agreement with in situ validation data as well as spatial correlation between regions and the natural terrain.

I. INTRODUCTION
L AND subsidence in The Netherlands is becoming an increasingly critical issue as it is closely linked with sealevel rise, flooding risks, and greenhouse gas emissions due to peat oxidation [1], [2], which is abundant in the region.Despite the importance of this issue, it is very difficult to accurately monitor subsidence rates across the country, and currently, no land surface time-series data exist with the required levels of accuracy, length, and spatial extent.Synthetic aperture radar interferometric (InSAR) is a very promising technique for monitoring land surface motion at large spatial scales with frequent temporal sampling.While InSAR techniques employing stable point scatterers (PS) have been successfully used to monitor subsidence in The Netherlands [3], [4], [5], these PS's are usually found at greater depths and the movement of the surrounding landscape has had to be indirectly inferred.
So far, it has been impossible to directly observe land surface motion using distributed scatterer (DS) techniques in The Netherlands because rapid soil movement, seasonal landuse changes, and high noise levels result in sudden losses The authors are with the Department of Geoscience and Remote Sensing, Delft University of Technology, 2628 CN Delft, The Netherlands (e-mail: p.n.conroy@tudelft.nl;s.a.n.vandiepen@tudelft.nl;f.j.vanleijen@tudelft.nl;r.f.hanssen@tudelft.nl).
In this article, we present a novel methodology for dealing with irreparable losses of interferometric coherence, which we term loss-of-lock.These events are almost always observed when attempting to monitor the motion of peatland surfaces with DS InSAR, however, the term is more general and can be applied to other circumstances as well.In Section II, we provide a definition for loss-of-lock as it relates to InSAR, as well as examples with real data.Section III describes the methodology we have developed to combat this problem and enable InSAR monitoring of these challenging regions.Section IV provides the results of several test areas and their validation against in situ measurements.Section V provides discussion, and finally Section VI concludes the article.

II. LOSS-OF-LOCK A. Definition
Interferometric coherence locks subsequent SAR acquisitions together.While losses of coherence are a common phenomenon, we differentiate between intermittent losses of coherence, where an event results in the loss of one or more epochs, but the overall time series is unaffected [Fig.1(a)], and loss-of-lock, a more serious loss of coherence which results in an irreparable discontinuity of the time series [Fig. 1

(b)].
A loss-of-lock is defined at a given epoch t L such that a sustained loss of coherence is observed, and no coherent interferometric combination exists which connects observations across the epoch in question.At low sample coherences (<0.1), the distribution of the interferometric phase approaches a uniform distribution [10].This means that all useful interferometric information (i.e., the displacement component of the interferometric phase) at that epoch is lost.Thus, the time series is effectively cut at t L , and the information content in the SAR image stack alone is not sufficient to estimate a connected set of interferometric phases spanning the entire observation period without additional information or assumptions.
Loss-of-lock is a diagnostic term in that it is defined based on conditions in an observed coherence matrix, rather than the occurrence of a particular scattering phenomenon, although the coherence losses are ultimately related to physical changes in the scattering object(s).For instance, a short-lived snowfall on an otherwise undisturbed and stationary grassland will result in an intermittent loss of coherence, while a loss-oflock may be caused by agricultural activities such as plowing, or changes in vegetation such as harvesting of crops, resulting in drastic reconfiguration of the scattering geometry of the ground, without implying any subsidence.
It is important to note that the presence of a loss-of-lock event may not be readily apparent from inspecting a displacement (or phase) time series.If one considers an event in which the mean surface level of the region under observation remains constant, but the scattering geometry changes drastically (for instance, by plowing), due to the wrapped nature of phase observations, the wrapped phase observation following the event may quite likely fall close to that of the previous epoch, and both phase unwrapping algorithms and manual inspection will overlook the change [as shown in Fig. 2(a) and (b)].Alternatively, it is also possible for large phase differences to be observed, due to changes in the scattering surface which are then misattributed as displacements [as shown in Fig. 2(c) and (d)].
Different sensors will be sensitive to different phenomena occurring on the ground and in the atmosphere, i.e., a lossof-lock observed at C-band may not be observed at L-band.A more practical description of the phenomena affecting our study areas is provided in Section II-B.

B. Observed Loss-of-Lock Events
A practical example of loss-of-lock is shown in Fig. 3, showing sample coherences for the period between January 2020 and May 2021.Two coherent periods (identified with the red dashed lines) are separated by a substantial amount of time, but more importantly, it can be seen that there is no significantly coherent interferometric combination linking them.This means that the two periods are effectively disconnected, and additional information will be required to estimate a time series spanning the entire observation period.Depending on the threshold used, an additional loss-of-lock event can be observed, resulting in three coherent periods, as indicated by  the magenta lines.This shows that: 1) detecting a loss-of-lock is dependent on the choice of allowed level of noise versus the amount of data used and 2) a loss-of-lock can be both a sudden and/or a sustained condition.This choice of threshold is discussed further in Section III-E.

III. DS PROCESSING METHODOLOGY A. Overview
A high-level end-to-end process flow diagram is given in Fig. 4. The system makes use of two well-established InSAR software packages, DORIS [11] and DePSI [12] (blue and red sections of Fig. 4, respectively).DORIS is used to align, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.resample, and geolocate the level-1 single-look complex (SLC) SAR image stack.DePSI is used to create a network of PS and estimate an atmospheric phase screen (APS), which can be applied to both the PS and DS phase observations.The remaining parts of this section are dedicated to describing the "Contextual Data" and "DS Processing Workflow" sections of the diagram (yellow and blue sections of Fig. 4, respectively), which are the novel aspects of this methodology.

B. Spatial Contextual Data and DS Pixel Identification
Three publicly available spatial datasets are combined via a spatial join operation: cadastral land use (parcel) polygons, soil maps, and groundwater management zoning (Dutch: peilgebied) [13].The datasets are provided as vector geometries in geopackage format, allowing for a straightforward combination of the data.This is accomplished by taking the land parcel delineations in the cadastral dataset as the base layer for the spatial join and performing a one-to-one attribution with the features with the largest overlap.While the cadastral and groundwater management zones follow similar geographic boundaries, the soil map has a different spatial structure, and this one-to-one attribution results in some information loss.For instance, a parcel mostly composed of peat with a smaller vein of sand running through it will simply be labeled with the peat soilcode.This is done to constrain the problem variables and have only one value per region.This method could be extended for regions with larger or more heterogeneous parcels by subdividing them either geographically or based on other relevant contextual information.Finally, the nearest weather station is found by Voronoi polygonization, and its corresponding ID is added to enable the attribution of (temporal) meteorological data.An example of such a combined dataset is shown in Fig. 5.
DORIS provides a coregistered SLC stack along with the geolocation of each pixel in a grid.Each pixel in the stack can, therefore, be assigned an ID corresponding to the polygon it belongs to (Fig. 5).Each polygon is assigned a coordinate according to its centroid.

C. Coherence Matrix and ESM Phase Estimation
Multilooking is performed on a per-polygon basis.As can be seen in Fig. 5, the Dutch peatlands are divided into rectangular parcels surrounded by drainage ditches, which provides us with a natural set of multilooking boundaries.While parcel sizes vary in shape and size, in general, the groundwater level and land cover within a parcel are consistent, atmospheric delay variability will be negligible (at the sub-mm level) [14], [15], and a parcel will typically contain 100 pixels.A minimum number of 50 pixels per polygon is enforced for noise suppression.Thus, we can ensure ergodicity and representativity while maximizing the number of equivalent looks.Following this parcel selection, we also optionally apply a statistically homogeneous pixel (SHP) test as outlined in [16].This can filter out misattributed pixels due to geolocation errors in the radar and contextual data, as well as the effects of unwanted scatterers within the region, such as electrical masts, light posts, trees, and so on.
The complex sample coherence matrix of a multilooked region, Ĉ, consisting of elements ĉi j is given by where S i, j n contains the complex value of the nth pixel in SAR images acquired at epochs i and j, the asterisk denotes the complex conjugate, and is the set of all selected pixels within the multilooked region.We differentiate between the complex coherence matrix C, and the coherence matrix, , which is the matrix of the magnitudes of the elements of C.An example of ˆ of a multilooked parcel is shown in Fig. 3.An equivalent single main (ESM) [17] set of phases is estimated using the "Eigendecomposition-based Maximumlikelihood-estimator of Interferometric phase" (EMI) method, as described in [18].This procedure reduces the full set of all interferometric combinations to a single set of consistent phases, φ, as estimated by the phase of the minimum eigenvector of the Hadamard product of the inverse of the sample coherence matrix with the complex sample coherence matrix, as given by where • denotes the Hadamard product, λ is the eigenvalue, and ξ is the eigenvector.The estimated interferometric phases are given by Strong decorrelation can hinder the effectiveness of the ESM phase estimation.In cases in which coherence is completely lost during a loss-of-lock event, it may be advantageous to only perform the estimate within the identified coherent blocks.This can reduce the amount of noise at the input to the estimator, however, one risks losing useful long-term coherent information.The decision to perform block-wise estimation could also be driven by contextual data, that is, a priori information about the land use/cover which indicates that a loss-of-lock has occurred, such as knowledge of plowing or harvesting events.
The estimated ESM phases φ per polygon are then imported into DePSI as virtual points into the secondary network of scatterers to apply APS filtering.The locations of the virtual representative points are given by the centroids of the given parcel polygon.The APS estimation is based on an initial primary network of PSs.The filtered phases φAPS are read back out of DePSI following the APS filtering stage.

D. Contextual Enrichment and Grouping
We have now obtained a set of wrapped, multilooked, and filtered DS phases which are each characterized by the set of attributes shown in Fig. 5, along with a point coordinate given by the polygon centroid.The parcels and their estimated phases are grouped together according to their shared attributes, establishing contextual groups.We contend that parcels that share the same land use, soil classification, and belong to the same groundwater management regime should be expected to behave in a similar fashion.That is, although we expect to see variations in phase according to differing noise and clutter conditions, local variations in soil stratigraphy, and variations in optical depth due to land cover, we expect that the parcels in a contextual group can be expected to move according to the same displacement model in the mean sense.This grouping becomes critical in the context of bridging loss-of-lock in the parcel time series, which is described in Sections III-E-III-G.The contextual group corresponding to the red highlighted parcel of Fig. 5 is indicated in yellow.The identified contextual groups are then filtered by their number of members: we have found a minimum of 30 members is needed to ensure sufficient coverage throughout the year; however, this value will change depending on the coherence behavior of the area under investigation.

E. Segment Identification
Due to the loss-of-lock phenomenon, attempting to interpret the entire ESM time series of phases at once is not possible and will result in several types of error, such as interpreting a noise-dominated signal as real deformation, or phase unwrapping errors when transitioning from incoherent to coherent interferograms [6], [8], [9], [19], [20].Thus, a different approach is required.
We begin by identifying which parts of a time series are of sufficient quality that they contain physically interpretable information that can be unwrapped with an acceptable degree of error.Despite using a full-rank method to estimate the ESM phases (Section III-C), we find that the best quality indicator we have available is the so-called daisy-chain coherence, γ DC , which is the magnitude of the first off-diagonal of the coherence matrix [corresponding to the indices j = i − 1 in (1)].These are the coherence magnitudes of the interferograms with the shortest temporal baseline in the dataset, which for Sentinel-1 data is six days.In general, we expect these to be the most coherent interferograms in the dataset, as less time has passed for decorrelation effects to occur [21], while orbital baselines for Sentinel-1 are always small, resulting in negligible baseline decorrelation [14].We threshold the daisy-chain coherence to identify sufficiently coherent subsections of the full time series, which we term (temporal) segments.A segment is a contiguous subset of a time series in which the coherence is sufficiently high to estimate a consistent set of interferometric phases.Thus, a segment is defined by two thresholds: the minimum coherence, and the minimum number of consecutive coherent epochs, which can be determined experimentally.In our case, a minimum of five consecutive epochs with γ DC > 0.12 is used as a threshold.
Each contextual group, therefore, contains many coherent segments: one for each contiguously coherent period of a parcel, times the number of parcels, times the number of tracks.By considering such a large number of segments, we can span loss-of-lock events in one parcel time series with coherent observations from a neighboring one from the same contextual group.

F. Temporal Ambiguity Resolution
The identified segments are initially treated as independent time-series.Temporal phase unwrapping (or ambiguity resolution) is performed independently on each segment using a method aided by a machine learning model, as described by [22].The ground surface level of peatlands is extremely unstable and prone to rapid fluctuations depending on temperature and precipitation levels, so we use a recurrent neural network (RNN) to aid in making predictions about which Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
ambiguity level is correct.This RNN model uses temperature, precipitation, and day of the year as inputs, which is publicly available daily weather data.Detailed information about the implementation and testing of the methodology is provided in [22].

G. Displacement Model
We have now obtained a collection of temporally unwrapped segments, which are internally consistent but disconnected from one another by an unknown vertical shift, z, which represents the unknown displacement history of the DS during the loss-of-lock period.Thus, to recombine the coherent segments, this unknown shift must be estimated.This can be accomplished with the aid of a displacement model, which can be used to align all the segments of a contextual group.
A parametric model that relates precipitation and evapotranspiration to soil surface displacement at a particular location has been developed in [23].This model has been shown to accurately model shallow soft soil movement in several locations in The Netherlands with different Holocene lithologies.While we prefer this model because it is accurate and requires very little input data, in principle any model can be used with this method.The model, M, is a function of precipitation, evapotranspiration, and the Holocene stratigraphy at the modeled location.It is a combination of reversible processes, such as shrinkage and swell, and irreversible processes, such as soil oxidation M x, P(t), E(t) = R x, P(t), E(t) + I x, P(t), E(t) where the model is parameterized by the lithology-dependent unknowns in x, t is the time, P is the daily mean precipitation [mm], and E is the daily mean reference evapotranspiration [mm] [24].R represents the reversible component and I is the irreversible component of the relative soil surface position.Daily values for P and E are provided at every weather station in The Netherlands [25].The reversible component is estimated by the scaled cumulative difference between precipitation and evapotranspiration R x, P(t), E(t) = where x P and x E are empirical scaling factors and τ is the integration time.The irreversible component is approximated as a linear rate, which is only considered active when R is negative, indicating drying soil conditions where x I is a constant, and f x, P(t), E(t) = 0, for R x, P(t), E(t) > 0 1, for R x, P(t), E(t) ≤ 0.
Thus, the model is parameterized by the four unknowns These parameters depend on the depth and stratigraphy of the Holocene sequence at a given location, that is, the lithology of that location.Details on the validation of the model are provided in [23].For the test locations shown in Section IV, the RMSE of the model with validation data is 6.9 mm in Zegveld, and 4.1 mm in Rouveen.Now we will show how to accurately estimate these model parameters, given the sparse unwrapped measurements we have available.This result can then be used to align the unwrapped segments of a given contextual group and estimate a continuous displacement time series.The relationship between the unwrapped phases of the mth segment of a given DS polygon, φ m , and the group displacement model, M, is given by where t is the time, θ is the incidence angle, λ is the wavelength, z m is the unknown vertical shift (constant for a given segment), and ϵ is a combination of noise, phase unwrapping errors, and model residuals.Equation ( 9) cannot be solved in its current form, as the model parameters x must be known a priori to evaluate the correct z.While they can theoretically be estimated simultaneously, the high degree of correlation between these unknowns can result in a very poor estimation.Instead, we note that z is common for all phases within a given segment.Thus, by taking the difference in time between phases, the z term drops out and the model parameters x can be estimated directly by solving Now that the model parameters have been estimated, z for each coherent segment can subsequently be estimated by taking the average difference between the model and the phase time series over the coherent period where x are the estimated model parameters, and ⟨•⟩ denotes averaging.This process is repeated for each contextual group described in Section III-D, so there is one model for every identified contextual group.This method can also be used to align the phase observations of multiple satellite tracks together, provided there is no significant horizontal motion, or else that the vertical component of the displacement phase can be accurately estimated, and that care is taken to ensure that the same object is used as a reference point across all tracks.

H. Spatial Ambiguity Estimation
The typical approach to DS InSAR processing involves applying a minimum cost flow spatial unwrapping algorithm to the data, such as the well-known SNAPHU algorithm [26].However, this approach is not well-suited to peatland observations due to rapid soil movements and the high degree of multilooking required [8], [20], [22].Heterogeneity in both Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the type and depth of the soft soil layer of the Holocene will result in different responses to the seasonal weather conditions to which the ground is exposed, leading to spatial differences in the seasonal amplitude of the reversible displacement.When combined with high degrees of multilooking, this can create sharp discontinuities in the downsampled interferogram, which will essentially lead to aliasing if strong spatial continuity constraints, such as those in SNAPHU, are applied [27].
For these reasons, direct spatial comparison of phase changes between adjacent parcels is an error-fraught process and could result in introducing additional phase unwrapping errors instead of improving the result.We therefore take advantage of the mean displacement of the contextual group, which is the best estimate of how the contextual group of parcels should behave on average.Since a time series for the expected mean behavior of the contextual group has already been estimated, it is now a straightforward process to apply integer least-squares (ILSs) estimation [28], [29] to refine the estimated ambiguities of each DS polygon belonging to the contextual group.The ambiguities are estimated by first obtaining a float solution, given by where â are the real-valued float ambiguity corrections, A is 2π times the n × n identity matrix, φDS is the vector of unwrapped interferometric phases of the given DS polygon, and φgroup is the vector of mean unwrapped phase of the entire contextual group.Q y is the variance-covariance matrix of the phase observations and can be approximated by the Cramer-Rao Bound (CRB) [19].Thus, the covariance between two interferometric phases φ i j and φ kl is estimated by where L is the effective number of looks [14], and γ is the magnitude of the sample coherence [as determined by ( 1)] of the given interferometric combination, as indicated by the epoch subscripts i, j, k, and l.Next, integer bootstrapping [30] can be applied, which provides the most likely integer ambiguities as where [•] is the rounding operator, ǎ ∈ Z n is the vector of estimated integer ambiguities, and l are the entries of a lower triangular matrix L obtained by decomposing the matrix Q −1 y into L and a diagonal matrix D, such that LDL T = Q −1 y .

I. Overall Model Test
Finally, a quality check is performed on the estimated contextual group results to ensure reliability.Overall statistics of the estimated contextual group parameters are generated for the entire AOI, and groups are flagged whose parameters deviate significantly (i.e., greater than 2σ ).The unwrapped parcel phases are compared to the estimated group model in flagged groups in which it is suspected that the contextual group model has been poorly estimated by means of an overall model test (OMT).
The OMT is performed by comparing the model residuals ê to the estimated precision of the observations Q y to generate the test statistic T for each DS polygon where N is total number of epochs and the nth element of ê is given by The operator refers to the fact that we use the differential daisy-chain phase as defined in (10) to remove the estimated vertical displacement shifts (the displacement occurring during the loss-of-lock periods) from the equation.
The test statistic T follows a central chi-squared distribution with four degrees of freedom, corresponding to the four unknown model parameters [see ( 8)], and is compared to a critical value which follows from a chosen significance level α.If T exceeds the critical value, then the model does not follow the observations to within the estimated precision of the observations at that significance level.In our case, the precision estimation comes from the CRB, which is the theoretical lower bound on the best achievable uncertainty.Thus, while it is correlated with the true uncertainty, estimating the CRB based on the sample coherence [see (13)] will systematically overestimate the uncertainty of the phase observations.Therefore, the significance level is chosen more strictly to compensate for this.
The OMT is performed recursively on flagged groups by choosing an initial α and removing points that are rejected by the test.The model parameters of the contextual group are then re-estimated with the rejected points removed.If the new model parameters fall below the acceptable threshold then the group is sustained.If the parameters still deviate, α is slowly increased and the procedure is performed again.If after several iterations (ex.5) the estimated model parameters still fall outside the accepted bounds, it is concluded that the model is not suitable for the terrain in question, and the group result is discarded.In a multiple-hypothesis testing context, this procedure could be reiterated with an alternative model.

A. Description of Satellite Data Used
Sentinel-1 imagery of two 10 × 10 km regions around Zegveld and Rouveen, The Netherlands, are used as test sites for the time period spanning January 2017-December 2022.In Zegveld, four tracks are used: ascending 088 and 161, and descending 037 and 110.In Rouveen, three tracks are used: ascending 015 and 088, and descending 037.The unwrapped segments of all available tracks are combined (as discussed in Section III-G) by projecting them onto the vertical after ensuring that all phases are referenced to the same object.The common reference point is found by identifying the PS with Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

B. Multilooking Based on Contextual Data Versus SHP Test Only
A comparison between a standard multilooking approach that employs 300 × 300 m regions and our parcel-based multilooking approach is shown in Fig. 6(a) and (b), respectively.As can be seen in the standard approach, despite the use of an SHP test, pixels from a number of objects that we do not expect to behave the same way are still averaged together.This is particularly apparent in the NE and NW corners of the image, where agricultural fields, residential yards, and greenery along roadways are all grouped together.
By including parcel cadastral information, we can help ensure that we are indeed averaging pixels that belong to the same objects or regions.An SHP test can also still be applied to remove unwanted pixels from within the parcel boundaries.

C. Coherent Segment Identification and Commonalities
Fig. 7 illustrates the advantage of grouping similar parcels together into contextual groups.While almost all regions provide sufficiently coherent (i.e., γ DC > 0.12) data over the winter period, from approximately October to April, the coherence of most regions drops significantly in the spring and is only intermittently present throughout the summer period until the following October.However, by combining the observations of enough similar parcels, we can have a year-round set of data with which to estimate the parameters of the displacement model as per (10).The coherence threshold of γ DC > 0.12 was experimentally found to be the highest value which still ensured sufficient data coverage year-round.
It is interesting to note that there is both a systematic and a random aspect to the coherence behavior of these regions.A systematic loss of coherence from April to October is clearly visible in the majority of parcels, however, the exact timing of this loss, as well as the intermittent recovery of coherence during the summer, seems to be a random event.This distribution is visualized by the shading of the background of Fig. 7.While it is clear that losses of coherence in these regions are caused by agricultural activities such as mowing and grazing, as well as changes in the scattering properties of the medium [21] caused by the drying of the soil and vegetation over the summer periods, it is unclear why some parcels seem to show higher coherence levels than others from the same contextual group at the same moments in time.This may be caused by some fields being used more intensively for agriculture than others, for instance, differences in the level of grazing between various fields.

D. Time-Series Estimation
An example group time-series result is given in Fig. 8.This result demonstrates how the displacement estimates of several temporarily coherent regions can be combined together to produce an unbroken time series of the overall region.The result matches very well with the available in situ validation data.Note that the validation data is not available for the entire span of the time series due to their installation dates.The difference between the contextual group median result and the validation data is quantified by the root mean square difference (RMSD) in Table I.However, it should also be noted that we do not expect an exact match between the InSAR and ground-based results, because the InSAR result shows the average behavior over a large spatial extent, whereas the ground-based measurement is of a single point.Moreover, the ground-based results do not capture the influence of the top five centimeters of soil, due to the position of the extensometer anchors.Nevertheless, as the major factors driving the motion are the same for both cases, we see that the agreement between them is very close, particularly in the observed short-term dynamics.

E. Effective Number of Looks Over Time
An important factor governing the accuracy of the result is the effective number of looks [14], shown in Fig. 9.This number fluctuates throughout the year due to the availability of coherent segments in the contextual group, as discussed in Section IV-C.It is important to ensure that there remain  where OSR is the oversampling rate given by where PRF is the pulse repetition frequency, f s,R is the range sampling rate, and BW az and BW R are the azimuth and range bandwidths, respectively.

F. Estimated Linear Rates
Approximate linear subsidence rates are shown in Fig. 10.These rates are estimated by linear regression of the contextual group mean time series results shown in Section IV-D, however, it should be noted that the total length of the observation period (five years) is too short to establish a robust estimate of the rate.Thus, these results provide an order of magnitude estimate and can be used to assess the spatial distribution of subsidence in the area.

A. On the Absence of Contextual Data
It is quite likely that in some cases, additional contextual data may not be available for the region under investigation, for instance, in peatland regions in remote locations.In such a case, additional remote-sensing data may be integrated into the processing workflow to identify and group common pixels together, such as the SAR backscatter data, as is done in the established squeeSAR [16] method, or through the use Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. of semantic segmentation techniques on colocated optical imagery.A posteriori technique such as T-SNE [31], [32] can potentially be used to group similarly behaving scatterers together into contextual groups.

B. Model Reliability and Goodness of Fit
Steps are taken to ensure goodness of fit and reliability of the estimated model in the OMT (Section III-I).The T -score of a given parcel [see (15)] is shown in Fig. 11, however, whether or not that parcel is used in the final rate estimation depends on the procedure outlined in Section III-I.When comparing Figs. 10 and 11, it can be seen that some parcels with a high T -score, and therefore a poor agreement with the contextual group model, are discarded from the final result.These are the ones detected by the iterative testing procedure.Other parcels with high T -scores are flagged for further re-evaluation but are not discarded immediately because their corresponding groups fall within expected bounds.There are several main causes of error that make a parcel deviate from the estimated contextual group model: • Misattribution within the contextual dataset: for example, errors in the soil map, or incorrect land use classifications.
• Model parameter estimation errors.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Often these causes are correlated; a misattributed parcel may be grouped with a set of other parcels to which it should not belong and introduce error into the contextual group model estimation.One region where this is evident is in a group of central-northern parcels (approximate coordinates: 52.62 • N, 6.13 • E) in the Rouveen area.Although classified as grassland, it is in fact a large rewilded "Natura-2000" region.Some of the parcels in this region are more densely covered with vegetation as opposed to being simple grasslands.This means that the phase behavior in these parcels is possibly different from the surrounding areas, and in some cases, the estimated displacement model may not be valid there.The OMT procedure can identify this and re-estimate a valid model with the remaining parcels not discarded by the test.A similar situation is visible in the SW corner of the Zegveld region (approximately coordinates: 52.095 • N, 4.75 • E).

C. Mean Displacement Model Versus Mean Phase Change
The use of a mean displacement model is a choice that in theory could be omitted.One could, for example, simply take the mean of the daisy-chain differential phase φ(t) of the entire contextual group and integrate it to obtain a relative position time series.However, in that case, one becomes very dependent upon a select few sets of measurements during the low-coherence times, as shown in Figs.7 and 9. Any biases, noise, or phase unwrapping errors in these observations would then be directly propagated into the mean contextual group time series.Therefore, using the set of all φ(t) observations to estimate a set of global model parameters is a safer option, provided the model is valid for the contextual group.

VI. CONCLUSION
Loss-of-lock is a permanent loss of coherence between two or more parts of a time series which is impossible to repair using the data in the SAR image stack alone.While decorrelation is a topic that has been discussed at length in the past, the specific implications of a loss-of-lock event are not well understood nor has a name been given to the phenomenon despite its very common occurrence in certain regions around the world, such as northern peatlands.
We introduced a new DS processing methodology that makes use of contextual data to reconnect coherent observations separated by loss-of-lock.With this methodology, we perform multilooking based on polygons which mark physically existing divisions in the terrain and assign a set of attributes and multilooked phases to each polygon.As is observable from their coherence matrices, most of these phase histories suffer from loss-of-lock.We combine the observations of different polygons which we expect to behave in a similar manner to parameterize a common functional model.This model is used to align the disparate observations to estimate a single unbroken time series for the contextual group.
Using this methodology, we have successfully been able to estimate accurate InSAR displacement time series in several subsiding peatland regions in The Netherlands which was previously not possible with InSAR.To our knowledge, this is the first time that an accurate and validated time series has been estimated based on direct observation of the peatland pixels using DS techniques.

Manuscript received 27
June 2023; revised 9 September 2023 and 9 October 2023; accepted 1 November 2023.Date of publication 3 November 2023; date of current version 16 November 2023.This work was supported in part by the Living on Soft Soils (LOSS): Subsidence and Society Project; and in part by the Dutch Research Council (NWO-NWA-ORC), URL: nwa-loss.nlunder Grant NWA.1160.18.259.(Corresponding author: Philip Conroy.)

Fig. 1 .
Fig. 1.Canonical coherence matrices showing different types of coherence losses for a set of five subsequent SAR acquisitions.Shaded cells: interferometric combinations that are sufficiently coherent to produce a useful phase estimation.Empty cells: insufficiently coherent combinations.(a) Intermittent loss of coherence at epoch t L that does not produce a loss-of-lock, because coherent interferometric combinations exist that connect epochs preceding and following t L .Archetype: intermittent snow cover.(b) Loss of coherence results in a loss-of-lock.There are no sufficiently coherent interferometric combinations connecting the epochs preceding and following t L .Archetype: plowing, harvesting.

Fig. 2 .
Fig. 2. Two examples (a)-(d) of a loss-of-lock which result in different interpretations, resulting in a perceived change in the displacement phase at epoch 50.Top (a) and (c): interferometric complex phasors, with arrows indicating the mean value pre-and postloss-of-lock (blue and red, respectively).Bottom (b) and (d): resulting unwrapped phases.Blue dots: true (noisy) displacement phase of the ground level.Red dots: observed (noisy) displacement phase in the presence of loss-of-lock.Green dashed line: estimated linear velocity.

Fig. 3 .
Fig. 3. Observed coherence matrix showing loss-of-lock in Sentinel-1 ascending track 88 of a multilooked region near Zegveld, The Netherlands.The red and magenta dashed lines are added to indicate the disconnected coherent periods.When or if a loss-of-lock occurs, it depends on the minimum allowable coherence threshold (discussed in Section II-B).

Fig. 4 .
Fig. 4. Simplified processing flow diagram showing the major steps taken to create DS time-series estimates from Sentinel-1 Level 1 SLC SAR data with the aid of spatial and temporal contextual data.

Fig. 5 .
Fig. 5. Graphical visualization of spatial contextual data in QGIS, based on land parcel polygons of the region surrounding Zegveld, The Netherlands.The attributed polygons are shown in green over a background optical satellite image of the region.A parcel of interest is highlighted with a red border, and the corresponding contextual group is highlighted in yellow.

Fig. 6 .
Fig. 6.Comparison between two multilooking strategies of the same area, color-coded by the region.Solid lines: boundaries of a given multilooking region.Circles: pixels included in the multilooking.(a) Standard square multilooking procedure using 300 × 300 m areas.(b) Parcel-based multilooking.

Fig. 7 .
Fig. 7. Chart showing the availability of coherent data over time for a period of one year for a selection of parcels belonging to the same contextual group.The y-axis indicates the ID number of a given parcel, and the presence of a solid line indicates the presence of sufficiently coherent data.The background is shaded to indicate the relative degree of availability (i.e., the number of coherent parcels divided by the total number of parcels) such that a white background indicates complete availability with darker shading as availability decreases.

Fig. 8 .
Fig. 8. Surface-level time-series results plotted against in situ extensometer measurements for the period January 1, 2017-December 31, 2022, at the a) Zegveld and b) Rouveen regions.Gray lines: all segments of all parcels belonging to the contextual group.Blue line: contextual group median time series.Red line: Mean of time-series segments of a selected parcel in the contextual group.Black line: in situ measurement by extensometer of the same parcel.All the individual coherent segments belonging to the contextual group are shown in gray for readability.

Fig. 9 .
Fig. 9. Chart showing the effective number of looks over time for the period January 1, 2017-December 31, 2022, for the same contextual group shown in Fig. 8(a).

Fig. 10 .
Fig. 10.Color-coded map of estimated linear subsidence rates for the period January 1, 2017-December 31, 2022, at (a) Zegveld and (b) Rouveen regions.Parcels with no fill indicate that no estimation has been made at that location.

Fig. 11 .
Fig. 11.Color-coded map of estimated goodness-of-fit according to the T -score value at (a) Zegveld and (b) Rouveen regions.Parcels with no fill indicate that no estimation has been made at that location.

TABLE I DIFFERENCE
BETWEEN INSAR AND EXTENSOMETER ESTIMATES