Sinkhole Scanner: A New Method to Detect Sinkhole-Related Spatio-Temporal Patterns in InSAR Deformation Time Series

: Sinkholes are sudden disasters that are usually small in size and occur at unexpected locations. They may cause serious damage to life and property. Sinkhole-prone areas can be monitored using Interferometric Synthetic Aperture Radar (InSAR) time series. Deﬁning a pattern using InSAR-derived spatio-temporal deformations, this study presents a sinkhole pattern detector, called the Sinkhole Scanner. The Sinkhole Scanner includes a spatio-temporal mathematical model such as a 2-dimensional time evolving Gaussian function as a kernel, which moves over the study area using a sliding window approach. The scanner attempts to ﬁt the model over deformation time series of Constantly Coherent Scatterers (CCS) intersected by the window and returns the posterior variance as a measure of goodness of ﬁt. In this way, the scanner searches for subsiding regions resembling sinkhole shapes over a sinkhole prone area. It is designed to detect large sinkholes with a high efﬁciency, and small sinkholes with a lower efﬁciency. It is tested at four different spatial scales, and on a simulated and real set of deformation data. Real data were obtained from Sentinel-1A SLC data in IW mode, over Ireland where a large sinkhole occurred on 24 September 2018. The Sinkhole Scanner was able to identify a pattern of low posterior variance zones consistent with the simulated set. In case of the real data, it is able to identify signiﬁcantly low posterior variance zones near the sinkhole area with the lowest value being 51.1% of the maximum value. The results from Sinkhole Scanner over the real sinkhole site were compared with Multiple Hypothesis Testing (MHT), which identiﬁes Breakpoint and Heaviside temporal anomalies in the deformation time series of CCS. MHT was able to identify high likelihood for Heaviside anomalies in deformation time series of CCS near the sinkhole site about 10 epochs before the sinkhole occurrence. We show that the Sinkhole Scanner is efﬁcient in monitoring a large area and search for sinkholes and that MHT can be used successively to identify temporal anomalies in the vicinity of areas detected by the Sinkhole Scanner. Future research may address other Sinkhole shapes whereas the underlying stochastic model may be adjusted. We conclude that the Sinkhole Scanner is important to be applied at different levels of scale to converge on potential sinkhole centers.


Introduction
The surface of the Earth is in a state of dynamic equilibrium due to a range of geophysical and geochemical phenomena. A particular phenomenon that we study in this paper concerns cavities that are formed due to the dissolution of subsurface bedrock, see Figure 1a. They result in unbalanced gravitational forces due to the Earth mass present above the cavity [1]. The burden of mass gradually sinks into the cavity and thereby leads to subsidence. Regularly, this gradual subsidence manifests itself on the ground surface [2]. This phenomenon is defined as a sinkhole. Sinkholes commonly but not exclusively occur in karst terrains [3]. Subsidence of the top surface varies in shape and size depending upon the soil type and nature of cavity formed underneath [4].   Sinkholes may exhibit precursory subsidence, indicating their formation. Early detection of precursory deformation patterns is important for early warning and sinkhole monitoring [2]. Sinkholes can be monitored using both in situ and remote sensing methods. For instance, leveling [5], extensometers [6] are helpful to monitor micrometer scale deformations in mining galleries, and micro-seismics are used to monitor activity in salt mining caverns [7]. Passive optical [8] and thermal remote sensing using UAVs [9] is used to monitor sinkhole prone regions. However, these methods have limited coverage and are expensive when monitoring large areas. Radar remote sensing using instruments such as Synthetic Aperture Radar (SAR) [10], can be useful to do so at a relatively lower cost. SAR images have a large spatial coverage, a high spatio-temporal resolution, and data from few sensors can be freely accessed like for Sentinel-1 SAR [11]. In the past three decades, conventional InSAR (Interferometric SAR) and time series InSAR (TSInSAR) techniques have been successfully applied to land subsidence monitoring [12][13][14][15][16][17]. TSInSAR is particularly well established in monitoring post-hazard surface deformation patterns [14], but it is only partially developed for observing small spatial scale hazards such as sinkholes [2,18,19]. Moreover, there is a lack of generic methods to detect sinkholes, although attempts have been made to model TSInSAR derived deformations as Gaussian shaped sinkholes [20].
Detection of sinkhole-related precursory deformation patterns can be conducted by either modeling the deformation patterns using known deformation models, e.g., Gaussian subsidence bowl, or by identifying anomalies in InSAR deformation time series. Modeling deformation patterns works well in detecting subsidence in a region of a high Constantly Coherent Scatterer (CCS) density, while identifying anomalies in InSAR deformation time series is recommended if the CCS spatial density is relatively low. Application and efficacy of the methods depend upon the CCS density in the sinkhole area of interest. The CCS density is commonly limited by the spatial and temporal decorrelation and atmosphere heterogeneity, especially in rural areas. Simulation of synthetic sinkholes can thus help in further developing models to detect sinkholes, based upon contextual information. Sinkhole simulations have been conducted by simulating subsurface bedrock and interaction between bedrock particles [1] during the occurrence of a cavity. When the cavity migrates upwards, surface deformation manifests itself on the surface, which could be modeled as a spherical Gaussian function. Modeling of surface deformation due to sinkholes has been done as well using viscoelastic modeling [21]. Currently, SAR-based sinkhole monitoring methods vary with respect to the spatial and temporal resolution of the images. Various SAR satellites provide data at different spatio-temporal resolutions and in different modes. The spatial resolution of currently available SAR products ranges between high (e.g., TerraSAR-X Stripmap mode) to medium (e.g., Sentinel-1 IW mode). High-resolution images allow us to identify the formation of sinkhole-related patterns using interferograms [19]. They are based on Gaussian-shaped patterns that, however, are unlikely to appear on medium resolution images. Sentinel-1A&B [11] offer medium spatial resolution SAR data with resolution equal to 20 m × 5 m in azimuth and slant range directions and a relatively high temporal resolution of up to six days. To study sinkholes, medium resolution Sentinel-1 data [22], unwrapped interferograms, and time series InSAR methods [18,23] have been used. Sinkholes have mainly been studied by analyzing deformation trends over constantly coherent scatterers (CCS) such as persistent scatterers (PS), detecting anomalous behaviour in InSAR deformation time series. Sinkholerelated subsidence can also be studied by searching for instantaneous acceleration or nonsteady-state deformation behavior. Such behavior can manifest as a Heaviside anomaly [24] in a deformation time series. Moreover, deformation behaviour can also manifest itself as a change in the linear deformation velocity model, where the deformation velocity increases a few epochs before the sinkhole occurrence [18]. This is termed as a Breakpoint anomaly. Such anomalies have been modeled using the Multiple Hypothesis Testing (MHT) [25].
The main objective of this study is to design and test a sinkhole detection method called the Sinkhole Scanner. It is based upon modeling the deformation time series, derived from InSAR, as spatio-temporal deformation patterns, using the inverted Gaussian function. The Sinkhole Scanner is designed on the basis of a spatio-temporal deformation model and using a moving window, inside which the deformation time series is modeled. The Sinkhole Scanner is mainly tested by detecting sinkhole related spatio-temporal deformation patterns over simulated and real sinkholes. Furthermore, the test results from this Sinkhole Scanner method are compared with the results from the Multiple Hypothesis Testing.
The remainder of the paper is arranged as follows. Section 2 presents a brief background on various sinkhole shapes and relates them to bedrock properties. Section 3 introduces the method for TSInSAR, sinkhole simulation, Sinkhole Scanner, moving window operation, and MHT. Section 4 presents the details of the study area and Sentinel-1 data used. Section 5 states the results followed by a discussion in Section 6. Finally, Section 7 draws conclusions and also gives future recommendations.

Background on Sinkhole Types and Shapes
The diameter of sinkholes commonly ranges between 1 and 300 m [4], while their evolution depends upon the nature of the terrain and the mechanism of material loss in the region. Usually, sinkholes have a low depth to diameter ratio if they occur over materials which have weaker cohesion, while they have a high depth to diameter ratio if the material has high cohesion [1]. Following the authors of [4], we distinguish four major categories of sinkholes or dolines: Dissolution, Collapse, Dropout, and Suffosion dolines. Dissolution dolines occur by water penetration through a large number of small fissures, followed by dissolving the rocks solutes. This is followed by a collapse caused by weathering of the overburden material [26]. Their shape is often Gaussian, see an example shown in Figure 1b, and the subsidence on the surface is gradual and evident. Collapse dolines occur when the above ground suddenly collapses into the cavity created below. This generally results in a cylindrical or steep sided conical shape. These types of sinkholes occur generally in karst and volcanic regions [27]. An example is shown in Figure 1c. Dropout or subsidence dolines occur majorly in urban areas and are relatively small in size, such as 10 m in diameter. These sinkholes are quite dangerous because they cause sudden failure of the ground. Suffosion dolines are caused through constant loss of underground material through a channel. This results in a funnel-shaped structure. An example is shown in Figure 1d. From the above, we note that the most common shapes of sinkholes are inverted Gaussian, cylinders, and cones. Figure 1e shows an example of a random distribution of 30 CCS which can be used to reconstruct the formation of a sinkhole shapes, as shown in Figure 1f-h for cylindrical, cylindrical, and conical shapes, respectively. In this study, we will use this way of recon-structing sinkhole shapes from TSInSAR-derived CCS location and their corresponding deformation time series.

Methods
This part first reviews the method for CCS location and associated deformation time series extraction using multi-burst SAR data in Section 3.1. It next describes the sinkhole simulation method and its use in extracting simulated sinkhole-related deformation time series in Section 3.2. It then presents the sinkhole scanner method in Section 3.3, and the method for moving window operation in Section 3.4. The MHT method [25] used for detecting temporal anomalies in individual deformation time series of CCS points is discussed in Section 3.5.

Multi-Burst TSInSAR Processing
An overview of the used methods can be seen in Figure 2. Given a stack of SAR images (usually ≥ 20), a standard TSInSAR processing can be applied to obtain the location of the CCS and their corresponding deformation time series in the line of sight (LOS) direction. Every Sentinel-1 SAR image in Interferometric Wide (IW) swath mode contains three sub-swath images, and each sub-swath image consists of nine bursts with sufficient overlap between adjacent bursts [28]. When the study area covers multiple bursts/swaths, multi-burst TSInSAR processing is required. Here, we take two-burst TSInSAR processing as an example and employ the Delft implementation of Persistent Scatterer InSAR (DePSI) as one of the tool packages for PSI processing [29][30][31]. For single-master interferogram stack generation, we first deramp all Sentinel-1A SAR images for each burst to improve the coregistration [28] and next coregister all slave images with the master image using Sentinel-1 precise orbit data and Digital Elevation Model (DEM)-Assisted Coregistration (DAC). Then, we reramp all images and create the interferogram stack by multiplying the complex conjugates of all slave images with the master image [32]. To concatenate bursts in a same swath, we use the deburst method presented in [33]. In this method, the bursts are concatenated using the average burst scan time in azimuth direction and scan time per azimuth line [34]. Specifically, the merging time t merg of two adjacent bursts can be calculated by t merg = (t last zd,be + t first zd,bl )/2. Here, t last zd,be and t first zd,bl represent the zero Doppler time (zd) for the last line of the earlier burst (be), and first line of the latter burst (bl), respectively. After debursting, PSI processing with, e.g., DePSI can be conducted as DePSI tool is customized to spatially divide debursted interferograms into small patches for efficient PSI processing in a large area.

Sinkhole Simulation
As discussed in Section 2, three sinkhole models are used in sinkhole simulation and later in modeling the data. The definitions of the sinkholes shapes are described in Table 1. The deformation function F is defined as a function of the two-dimensional coordinates (x, y) over the ground surface, and time t. Therefore, each function has a spatial and a temporal component. In time, the sinkhole shape can be assumed to evolve linearly, and is defined as where v is the linear deformation velocity of the center of sinkhole model located at (x 0 , y 0 ), and c is the constant term. In space, the sinkholes are defined on the 2-D horizontal surface using the functions described in Table  1. The deformation functions F for Gaussian and Conical models are dependent upon spatial coordinates (x, y), whereas the deformation function for the Cylindrical model is independent of the spatial coordinates, and models the subsidence of points using the temporal model only. For the Cylindrical and Conical functions, there is a constraint placed on the points: ∀P i | x i 2 + y i 2 ≤ r 2 , where P i represents CCS points with i ∈ [1, b], (x i , y i ) represents the CCS coordinates relative to the center of the sinkhole, r is the radius of extent for the cylindrical and conical function on the top surface, and b is the total number of CCS intersected in the simulation window. Therefore, this means that only the points located within the circle with radius r and sinkhole center (x 0 , y 0 ) are influenced by the simulated deformation function F, and later, participate in the inverse modeling of sinkholes, see Section 3.3. Table 1. Sinkhole inventory: mathematical model of shapes used for sinkhole modeling.

Inverted Gaussian Model Cylindrical Model Conical Model
Model equation r: radius of extent for the cylindrical and conical function on the top surface, P i : i th CCS with coordinates (x i , y i ), relative to the center of the sinkhole (x 0 , y 0 ), I(t) = v · t + c, t: time step, v: linear velocity rate, c: constant term, and ζ signifies the extent of the Gaussian function.
The assumptions behind these simulations are that the sinkhole evolves in the shape of one of the sinkhole models defined in the sinkhole inventory, and that the radar signal is able to capture the movement of surface at the location of CCS, albeit in the LOS direction. Synthetic deformation time series is simulated over the estimated CCS locations. The deformation time series of the CCS locations, modulated with artificial Gaussian noise is obtained by , which is extracted from the simulated deformation function F; b is the total number of CCS; and N (0, σ 2 ) signifies a normally distributed noise with zero expectation and σ standard deviation introduced to simulate the error d e in the deformation. The sinkhole is then positioned over a location in the study area using the moving window operation.

Sinkhole Scanner
The Sinkhole Scanner is defined as a two-dimensional testing space shaped by its window size in the x − y plane and the scanning operator. The scanner moves in either directions over the observation space with a predefined stride in the x and y directions. The testing window searches for sinkhole shapes defined in the sinkhole inventory that show an evolving pattern. For instance, for the Gaussian shape this is obtained by solving the mathematical model of the sinkhole shape, e.g., Equation (4) for Inverted Gaussian shape. The following text describes the mathematical models for the three sinkhole models derived using the deformation functions F, defined in Table 1.

Inverted Gaussian shape solution:
As the inverted Gaussian model is described as a nonlinear equation in Table 1, the least squares estimation cannot be employed in a straightforward way. Therefore, we linearize this nonlinear equation to have its linear form, by using the logarithmic transformation, Upon this, we can define the corresponding functional and stochastic model, expressed as [35] E{d} where E{} and D{} are the expectation operator and dispersion operator. The matrix A is the design matrix and x is the vector of unknowns. The stochastic model Q dd is designed to describe the observation errors. By relating the functional model in Equation (3) with Equation (2), we specify the functional model of the inverted Gaussian model as The vector d  To simplify the simulation, one can use a scaled identity matrix, i.e., Q dd = σ 2 · R bm×bm = σ 2 · I bm×bm , which implies deformation of the b points are uncorrelated in time and space and the deformation variance is invariant over time. Here, σ 2 is the variance of unit weight, R bm×bm is the cofactor matrix, and I bm×bm is the identity matrix.

Cylinder shape solution:
The cylindrical and conical function are dependent upon the user defined parameter for the extent of the sinkhole r. The default value for this variable is chosen to be half of the scanner window size. The solution for the cylindrical sinkholes is performed by firstly applying the constraint ∀P i | x i 2 + y i 2 < r 2 . The CCS inside the circle centered at the center of the scanning window and with radius r are modeled using a common linear deformation velocity, expressed as (see in Table 1 as well) where v and c represent the linear deformation velocity and the constant term in the model fitting, respectively. The CCS outside the circle, but inside the scanning window, are assumed to be in a steady state. Assuming b points are intersected due to the constraint, the functional model and stochastic model shown in Equation (3), with b replaced by b , can be used to solve for the cylinder shape. Specifically, in this case the functional model can be expressed as where d b m×1 represents the deformation observation vector. The matrix A, with size b m × 2 is the design matrix and x, with size 2 × 1 is the vector of unknowns.

Conical shape definition:
Similar to the solution for cylindrical shape, the points satisfying ∀P i | are used for solving the model equation. The functional model is derived using Figure 3, Figure 3. Therefore, the observation function of any point P i can be written as Here, the unknown parameters are the velocity of sinkhole center is denoted by v and the constant term associated to the model is denoted by c. The corresponding functional model is described as Here, d has the size of b m × 1, where b signifies the number of points intersected by the constraint circle. The symbols have the same meaning as in the case of Cylindrical model. Moreover, the shapes of A and x are the same as in the previous Cylinder shape case as well. According to the least squares estimation method, the unknown parameters of the mathematical models shown in Equation (3), are estimated bŷ The error after the parameter estimation is obtained asê = d − A ·x. The precision of the least-squares estimations can be estimated by The detectability of the model is assessed using root mean square error RMSE = E[ê 2 ], and posterior varianceê T (Q dd /σ 2 ) −1ê /(bm − n) for inverted Gaussian shape,ê T (Q dd /σ 2 ) −1ê /(b m − n) for conical and cylindrical shapes, after fitting the model to the data. Here, n is the number of unknown parameters, i.e., 2, and σ 2 is the variance of unit weight.

Moving Window Operation
Instead of moving the scanner over the whole study area, and intersecting points with the moving window, we suggest to efficiently identify sinkhole scanning locations based upon the locations of CCS points. We identify the lower-left coordinate of the box encompassing the CCS points in the area by where Box ax,0 represents the initial position of the window or box on axis ax, ax ∈ {x, y}, and the superscript 0 identifies the bottom left index of the window. The parameter P ax i represents the component of the positional vector of P i of the i th CCS point with i ∈ [1, b] on the axis ax, and win ax represents the window size on the corresponding axis. The notation . represents the greatest integer function and min( ) is the minimum operator. Unique windows are identified using Equation (10). The corresponding CCS point locations and deformation estimates are grouped with the boxes, as attributes of the box. The Scanner then performs an iteration over the boxes and their grouped attributes by performing a least squares estimation. This one-dimensional iteration method largely enhances the performance of the Scanner as compared to two-dimensional iterations with successive intersection of CCS points with the moving windows.

Hypothesis Testing Method
The pointwise hypothesis testing is used to identify anomalies in time. The Breakpoint [18] and Heaviside anomalies [24] are chosen to be two temporal anomalies which can resemble sinkhole-related deformation behaviour. Hypothesis testing is used to identify the best epochs in the deformation time series where either a Heaviside model, or a Breakpoint model, with one jump, or break, would fit. The functions for Heaviside (M H ()) and Breakpoint anomalies M Br () are defined as [25] Here, H(t − τ i (P k )), for i ∈ [1, m − 1] is the Heaviside step function, centered at τ i (P k ), ∆ i (P k ) is the Heaviside step size, v k is the deformation velocity at time t k . For instance, for a single Breakpoint model, the deformation time series is divided into two parts with two different linear velocity parameters, i.e., v 1 and v 2 , at the breakpoint epoch l. The corresponding single Breakpoint model can be written as M i . This anomaly detection is done by calculating the test statistic and test statistic ratios for a null hypothesis and other alternative hypotheses (details in Section II in [25]). The null hypothesis H 0 and alternative hypotheses H j , for instance, for the inverted Gaussian shaped sinkhole case, can be expressed as [25] where C j ∇ j in the j th alternate hypothesis is used to capture the incidental shifts in the observations. C j is the specification matrix and ∇ j is the vector of additional unknowns with dimension q, and q ∈ [1, bm − n], where n is the total number of unknowns in the model. For cylinder and conical shaped sinkhole case, b is replaced by b for Equations (11)- (13). The test statistic ratio results from MHT are then compared with the results from the Sinkhole Scanner step. The performance, effectiveness and use of the results is also compared.

SAR Data and Test Site Description
A total of 75 Sentinel-1A SAR images with track number 1 and frame numbers varying between 170 and 175, acquired between 15 April 2015 and 31 December 2018 with a 12-day repeat cycle, were used in this study. All SAR images were obtained in ascending orbital direction, IW mode, and VV polarization channel, see Table 2. Figure 4a shows the study area located in North Eastern Ireland, while Figure 4b shows the zoomed-in area. It is in this area that a sinkhole occurred, see Figure 4c. Figure 4a also depicts the Sentinel-1A image boundary at the first acquisition date, as well as the two adjacent bursts in the black dashed line. Figure 4d shows the cropped region from the merged bursts used for this study using a brown dashed line, with the earlier and latter bursts labeled as Burst-1 and Burst-2 respectively. This sinkhole, shown as the red dot in Figure 4a-c, occurred on 25 September 2018 and its size was reported to be in the order of 100 m [36]. Figure 4c shows an areal picture taken immediately after the sinkhole had occurred. This sinkhole had been attributed to underground activities in a Gypsum mine which can be observed at the bottom of the sinkhole spot Figure 4b [37]. Historically, this is not a sinkhole-prone area. However, the area is situated in a karst region and is rich in gypsum and limestone [38]. For these reasons, the type of sinkhole resembles a combination of a dropout and collapse doline.

Results
This section presents the results of the study, in response to the subsections in Section 3.

TSInSAR Results
A total number of 294,519 CCS (i.e., PS) were identified in the cropped study area with size 63.73 km × 43.84 km in east and north direction, respectively, using standard PSI approach with DePSI. Shuttle Radar Topography Mission (SRTM) 1 arc second resolution DEM data [39] was used to remove the reference DEM phase from the interferograms. The CCS locations were obtained in radar coordinate system, see Figure 5a, and in World Geodetic System WGS-84 coordinate reference system (CRS) see Figure 5c. The coordinates in WGS84 CRS were projected onto the local Universal Transverse Mercator (UTM) coordinate system (EPSG code: 2157) and thus used in sinkhole simulation and scanning steps. However, all the results are presented in the WGS84 CRS to ensure consistency. The pink triangle in Figure 5a,c shows the location of the reference point relative to which the deformations are estimated.
During the postprocessing of the PSI results, additional filtering was performed by applying thresholds on ensemble coherence (>0.2), deformation velocity (<5 mm yr −1 and >−25 mm yr −1 ), and spatio-temporal consistency (STC > 0.008) [40]. The threshold on the ensemble coherence was applied in order to reduce noisy CCS, i.e., CCS with noise in the deformation time series. Moreover, the thresholds on deformation velocities removed CCS with unrealistic deformation velocities, e.g., unusual uplift of points by velocity >5 mm yr −1 . This filtering step resulted in a large reduction of noisy CCS, as the total number of CCS reduced to 46,611. A general stable trend was visualized over the cropped study area, see Figure 5a,c. Interestingly, there were a few points in the center of Figure 5b,d with high deformation velocities (∼ −20 mm yr −1 ). These points are located in the mine area that is ∼ 100 m away from the rim of the sinkhole.

Sinkhole Scanner Results
We now describe the results from the Sinkhole Scanner step.

Case: Simulation
The two-dimensional horizontal space was divided into grids of 10 × 10 m, which became the minimum mapping unit for the simulation and for scanning of sinkholes. The common deformation time series of sinkhole centers was defined between −0.5 and −100 mm varying linearly over the temporal baseline between 0 and 3.65 years, with a discretization of 10 intervals or epochs. This temporal model was used along with spatial simulation models described in Table 1. The spatial models were defined over the grid defined above. These simulation parameters are summarized in Table 3. The experiments were designed to simulate and scan sinkholes with sizes ranging from large ∼200 m, to small <∼10 m. Therefore, simulations were conducted with two different sets of values for window size and stride parameters, i.e., a common size of 2000 m and 100 m in both directions. These simulated sinkholes as well as real deformation data were scanned using the Sinkhole Scanner with common window-size and stride values of 2 km, 1 km, 500 m, and 100 m. Window-sizes and stride parameters were chosen to be ∼10 times the size of the sinkhole which is being looked for. This factor was chosen to ensure the intersection of enough CCS to solve the mathematical model. The least number of CCS required to attempt the mathematical model solution was 3, because the number of unknown parameters were 2 in all sinkhole shapes. With this setup, the aim is to find the optimal parameters to model the simulated sinkholes, and sinkhole-shaped patterns in the real data, using the Sinkhole Scanner. The stochastic model was designed as a scaled identity matrix in which σ 2 = 5 mm 2 based on experimental test. We assume that the spatial correlation of CCS which are further away from each other was assumed to be marginal for this study. Table 4 illustrates the choices for the window sizes and stride parameters in the simulation and scanning space. Case: Scanning Simulated Data Figure 6 shows the simulation and scanning results using the Gaussian model. Figure 6a shows the contour plot of the simulated sinkholes which were simulated using Scenario-L, draped over CCS locations shown as black dots. Figure 6b shows the corresponding simulated deformation velocity map over CCS locations. High deformation velocity circular zones can be observed at regular intervals over CCS locations. Finally, Figure 6c shows the application of Scanner 3 (window-size= 500 m) on the simulated data. The map shows goodness of fit of the Gaussian mathematical model for each 500 m × 500 m scanning window. The quality statistic for goodness of fit is the posterior variance. The scattered white patches in the map show regions with insufficient number of CCS to solve the mathematical model. A regular pattern of low posterior variance zones was observed, as expected over low subsidence zones. This can be seen better in the zoomed in image in Figure 6d, corresponding to the black box in Figure 6c. However, low posterior variance zones were observed in windows where the CCS density was sufficient to recreate the sinkhole shape. This regular pattern shows that this scanning method can highlight low subsidence zones with respect to relatively stable regions. We show the results from Scanner 3, which provided the best results of all the four scanners. Note that the best performing window size in this case was ∼10 times the ζ parameter chosen for the simulated sinkholes. However, limited success was observed while scanning for sinkhole simulation Scenario-S using any of the scanners. This may have been due to either too many sinkholes in case of a large window size (e.g., Scanner 1) or too few CCS intersected in case of a small window size (Scanner 4). In addition, we observed that the model was able to correctly estimate the v parameter with an error of under 5%, whereas the parameter ζ was overestimated which may be due to the finite value of the parameter δ.    Figure 7 shows the result of scanning real deformation data using all four Sinkhole Scanners, i.e., Scanner 1-4, using the Gaussian shape in Figure 7a-d, respectively. In the case of Scanner 1 (Figure 7a), a trend of better fitting grids was observed when moving from west to east. This is also slightly evident in the case of smaller window sizes as well, see Figure 7b,c. As the window size decreases, this trend is observed to converge near the actual sinkhole location.

Case: Scanning Real Data
We also observed that the range for posterior variance values widens as we decrease the window size. This is because the lower bound of the values decreases, but the upper bound remains almost the same. The absolute values of posterior variance are high. This is likely due to the finite value of the δ parameter used in the functional model to elevate the deformation values above zero. This causes a Gaussian to be modeled considerably above the zero deformation value, and causes a bias between the modeled and original observations. However, as the δ is constant and chosen to the smallest deformation value in the observation space, its effect in space is consistent. This, however, leads to a higher dynamic range in the posterior variance values and therefore a higher contrast in the goodness of fit map, which helps in identifying vulnerable areas much better. This offset can be reduced by employing more flexible models, or by using Taylor expansions. Furthermore, we observed that reducing the window size of the scanner resulted in a finer resolution sinkhole scanning map, however at the cost of a reduced scanned area due to the nonuniform CCS distribution. Table 5 shows the variation of scanned area with reducing window size.    Contrasting results were observed for the Cylindrical and Conical sinkholes where the trend of fitting was opposite to the Gaussian model. High posterior variance values were observed for sinkhole zones as compared to stable areas. This suggests that the Gaussian model is best suited for modeling sinkholes in this area. The values of posterior variance were much lower than in the Gaussian case, most likely due to the absence of the offset parameter δ in the observation vector. However, the contrast between fitting zones was less pronounced than in the Gaussian case.
The estimated velocities for the windows also showed high deformation rate near the sinkhole area, see Figure 8b. On comparison with the window-wise (Scanner 4) spatially averaged InSAR derived deformation velocities of CCS, see Figure 8a, we notice that the model overestimates the deformation velocities Figure 8c. This is likely due to spatial averaging of the deformation velocities.    Figure 9 shows the deformation velocity map and deformation time series of CCS located near the sinkhole site. The results show a large deformation velocity near the sinkhole site, as shown in Figure 9a. The points are annotated and their corresponding deformation time series over 75 epochs is shown in Figure 9b. We observed that a few epochs before the sinkhole event, there was a sudden change in the behaviour of the deformation time series of CCS lying over the mining region, close to the sinkhole site. The change is in the form of sudden deformation velocity change approximately 10 epochs before the sinkhole.

Hypothesis Testing Results
In order to validate any such anomalies, we use multiple hypothesis testing (MHT), to automatically detect any such behaviour. The level of significance was initially defined as α = 1/(2m) = 1/(2 × 74) ≈ 0.68% and the power of the test γ was chosen to be 50%. The horizontal axis of Figure 10 shows the temporal baseline, and each colored time series shows the test ratio statistic for the Heaviside and Breakpoint model in Figure 10a,b, respectively, for each CCS. For instance, we can notice a sudden drop in the deformation for PS2 (time series in orange) at the 24th epoch, (temporal baseline = 1.085 yr), in Figure 9b. Corresponding to the 24th epoch, we can observe spikes in the corresponding test ratio statistic time series for both Heaviside and Breakpoint anomaly (also shown in orange). This way, MHT is able to detect the best deformation model from a set of deformation functions. Figure 10 shows that the highest value of test ratio statistic was for PS1, which indicates that these anomalies are more apparent in the deformation time series of PS1 as compared to other CCS. Furthermore, Figure 10a shows that there are high indications of Heaviside jumps for CCS above PS8 and PS9, at the epoch corresponding to the sinkhole event. Moreover, Breakpoint test ratio increases for PS6 just before the sinkhole occurrence. PS6, PS8, and PS9 are situated directly above the sinkhole site.
The results from the Sinkhole Scanner were compared with the results from the Hypothesis testing method. Both the methods give relative information about only temporal anomalies with MHT, and both spatial and temporal anomalies with Sinkhole Scanner. The Sinkhole Scanner gives a map of sinkhole similarity for a large area because of its gridwise operation, and is not limited to give information only over the CCS point locations. The MHT method gives results only in time, which can be projected onto the CCS point location in space, as done in [41], but cannot extend to give information of a spatial trend in the result. However, MHT is not limited to the CCS number threshold required to deliver a result as in the case of Sinkhole Scanner.

Discussion
One of the most interesting observations in this work was the identification of large deformations near the sinkhole site in Ireland. This was observed in the deformation velocity map in Ireland. These observations were valuable and fortunate because of the orientation of mining terraces to the radar line of sight which acted as double bounce scatterers. This would otherwise not have been possible because of the large number of vegetation scatterers. The advantage of using CCS was a higher accuracy in geolocation and deformation estimates, although at the expense of lesser number of measurement locations. The subsidence over the sinkhole area was attributed to a mine collapse [42]. The collapse of the mine was approximately 50 m away from the periphery of the sinkhole. This information, along with the InSAR subsidence data, were sufficiently convincing that the deformation was related to the unstable ground over the mining area, and that its collapse led to the sinkhole. Validation using in situ reference data could not be performed because of their unavailability near the sinkhole spot. As soon as reference data, e.g., from GPS or extensometers are available, the InSAR signals can be validated [43] and thus used to verify actual anomalous activities in space and time. Further ancillary information, e.g., groundwater table can help in identifying signals due to sinkholes from other geo-processes.
The major contribution of this study was the Sinkhole Scanner which was defined to model the LOS deformation time series over CCS. It was assumed that that the LOS direction was sensitive to the sinkhole-related deformation in the vertical direction. The scanner was defined using a 4-dimensional moving window-based method in space and time, where a 3-dimensional sinkhole model, like an inverted Gaussian function evolving in time, is fitted to the CCS deformation time series overlapped by the moving window using least-squares estimation. The Gaussian function was simplified such that there were only two unknown parameters: the extent of the Gaussian function ζ of the circular shaped Gaussian, and the linear deformation velocity v of the center of the Gaussian function. The stochastic model used for estimating the unknown vector was an identity matrix scaled to a common variance.
The results from the Sinkhole Scanner show its effectiveness in detecting sinkhole-like evolving spatio-temporal deformation patterns. The four Sinkhole Scanners, designed to detect sinkholes at different resolutions were all able to detect Gaussian shaped sinkhole subsidence patterns. This was validated using the sinkhole simulation experiment and was especially evident in the area around the real sinkhole location. However, it is not easily observed in the deformation velocity map. Therefore, the Sinkhole Scanner is well applicable when using spatio-temporal deformation data in modeling such subsidence patterns. Furthermore, Scanners 1-3 highlight the grid containing the sinkhole area to be the most sinkhole-like deforming region. This establishes that the Sinkhole Scanner is useful in detecting sinkhole-like features. The Gaussian shape model, in contrast to Cylinder and Conical shape model, showed the best results in terms of fitting better to sinkhole zones as compared to stable areas. This fits well with the Karst geology of the terrain where Gaussian shaped sinkholes are more likely to occur [4].
From Scanners 1-4, the lower-bound of the posterior variance range reduces which may indicate improved Gaussian shaped patterns at a fine spatial scale. Unfortunately, the coarse CCS density does not allow us to investigate this further, but this can be done with higher spatial resolution datasets, e.g., Spanish PAZ SAR data [44].
The mathematical model for the Gaussian shape was defined to identify sinkholeshaped patterns with the center of the sinkhole lying at the center of the search window. It may happen that the sinkhole is located at the peripheries of the search window. In that case, the stride parameter can be reduced to ensure that the center of the search window aligns with the potential sinkhole center. It may also happen that there may be more than one sinkholes in the same search window. In that case, the window size of the scanner can be further reduced. The scanning results were dependent upon the CCS point density. When the scanning window size decreased, the number of CCS points overlapped by the scanning window also decreased and thus reducing the scanned area. For instance, the scanned area reduced by almost 100 times by reducing the window size from 2 km to 100 m. However, the fine-resolution map made with Scanner 4 helped in highlighting the location of sinkhole and also identified the best fitting locations with lower posterior variance values than the other three scanners. Furthermore, the spatial resolution of the radar is important when choosing a window size. A finer resolution image would give a higher density of CCS, which can be worked upon by even smaller Sinkhole Scanner window sizes.
The results also showed that the MHT is able to identify Heaviside anomalies in two CCS a few epochs before the sinkhole occurrence. MHT method can be used as another tool to identify temporal anomalies which can indicate vulnerable areas. The Sinkhole Scanner, in conjunction with MHT can be useful in detecting sinkhole related subsidence. The Sinkhole Scanner can thus be used for performing a global search over a sinkhole affected area, followed by Hypothesis testing to further investigate sinkholerelated anomalies.
In this research, we extended the research in [20]. We used more robust models, better model fitting methods, worked on a coarser resolution data and worked on a larger study area in a non-sinkhole prone region. However, several limitations still exist to the work. The major one is that we have used a general stochastic model for the Sinkhole Scanner. A more detailed stochastic model can be used to take the uncertainty in the deformation estimation into account [45]. Furthermore, a simplified Gaussian function was used to detect sinkholes, whereas the deformation behaviour can be more complex.

Conclusions
This study presents the Sinkhole Scanner, which includes methods to identify spatiotemporal changes in CCS behavior. The method is robust in the sense that different sinkhole shaped models can be used in place of the Gaussian function, e.g., a cylinder, a cone, or even a more flexible Gaussian function. Four different variants of the Sinkhole Scanners with different window sizes were used to identify sinkholes in simulated and real deformation time series. The scanner showed ∼30% lower posterior variance values over the sinkhole areas for both simulated and real datasets. Moreover, it showed a high contrast in the posterior variance values for windows near the sinkhole area. This was not visible from the deformation velocity map. Therefore, the Sinkhole Scanner provides additional information about sinkhole shaped spatio-temporal patterns, as compared to the deformation velocity map. Furthermore, Gaussian shaped sinkhole performed best in highlighting correct contrast between sinkholes and stable zones.
A trade-off was identified between the model's ability to detect smaller sinkholes and the number of windows with no solution to the least squares model. In this study area, there was 100 times reduction in the area returning successful model fitting from a larger Scanner 1 to a finer Scanner 4. The deformation velocity map over Ireland showed high deformation velocity near the mining area. This deformation trend suggested precursory deformation patterns in both space and time around the mining area. Sinkhole Scanners 1-3 were able to identify Gaussian-shaped patterns over the real sinkhole area. In the case study, Scanner 3 performs best role of providing high coverage area, and also highlighting the sinkhole location. In addition to real deformation data, deformation data were extracted from a simulated sinkhole field. The simulated data gave sufficient evidence of providing contrast between sinkhole and non-sinkhole areas. In the simulation case, the window size of 500× 500 m also provided the best results. The study suggests that multiple window sizes should be tested to get an idea of the spatial trend of posterior variance in the area. Smaller window size can be further tested in areas with relatively lower values of posterior variance. This would also considerably reduce the scanning time.
Future studies may consider learning-based algorithms [46] to identify anomaly types, while convolution filters [47] with successive lowering of the window size can be explored to encode information at various scales and then decode it to identify sinkhole shaped deformation patterns. Furthermore, a more complex stochastic model in InSAR processing can be used, see, e.g., in [45]. The addition of DS-InSAR would also be useful in increasing the InSAR measurement point density. Distributed Scatterers [48], and Temporary Coherent Scatterers [49] would help in increasing the number of measurement locations. Their precise role should be explored in future studies. A priori contextual information about the size and shape of the brewing sinkholes can be used for choosing the window size of the scanner and sinkhole shape, respectively. Finally, in the MHT step, more and physically realistic number of the breaks or jumps could be added to the models.
Funding: This research received no external funding.

Data Availability Statement:
The datasets used for this research are available in [50].

Acknowledgments:
The authors would like to thank ESA for providing Sentinel-1 images used for this study. We also thank Delft radar group for sharing DePSI toolbox.

Conflicts of Interest:
The authors declare no conflict of interest.