LiCSBAS: An Open-Source InSAR Time Series Analysis Package Integrated with the LiCSAR Automated Sentinel-1 InSAR Processor

For the past five years, the 2-satellite Sentinel-1 constellation has provided abundant and useful Synthetic Aperture Radar (SAR) data, which have the potential to reveal global ground surface deformation at high spatial and temporal resolutions. However, for most users, fully exploiting the large amount of associated data is challenging, especially over wide areas. To help address this challenge, we have developed LiCSBAS, an open-source SAR interferometry (InSAR) time series analysis package that integrates with the automated Sentinel-1 InSAR processor (LiCSAR). LiCSBAS utilizes freely available LiCSAR products, and users can save processing time and disk space while obtaining the results of InSAR time series analysis. In the LiCSBAS processing scheme, interferograms with many unwrapping errors are automatically identified by loop closure and removed. Reliable time series and velocities are derived with the aid of masking using several noise indices. The easy implementation of atmospheric corrections to reduce noise is achieved with the Generic Atmospheric Correction Online Service for InSAR (GACOS). Using case studies in southern Tohoku and the Echigo Plain, Japan, we demonstrate that LiCSBAS applied to LiCSAR products can detect both large-scale (>100 km) and localized (~km) relative displacements with an accuracy of <1 cm/epoch and ~2 mm/yr. We detect displacements with different temporal characteristics, including linear, periodic, and episodic, in Niigata, Ojiya, and Sanjo City, respectively. LiCSBAS and LiCSAR products facilitate greater exploitation of globally available and abundant SAR datasets and enhance their applications for scientific research and societal benefit.


LiCSBAS Overview
Here, we provide a brief LiCSBAS overview prior to the detailed description presented in Sections 2.3-2.5. The workflow of our InSAR time series processor is largely divided into two parts: preparation of a stack of unwrapped data (Step 0) and time series analysis (Step 1; Figure 1). LiCSBAS starts with downloading the LiCSAR products covering the area of interest (Step 0-1) and is followed by data format conversion (Step 0-2). Tropospheric noise correction using the external Generic Atmospheric Correction Online Service for InSAR (GACOS) data [38] (Step 0-3) and masking (Step 0-4)/clipping (Step 0-5) the unwrapped data are optional steps that can be carried out to improve accuracy and make processing more efficient. In the time series analysis, incorrectly unwrapped data, which degrade the results, are identified and discarded based on the coherence and coverage of the unwrapped data (Step 1-1) and by checking loop closure (Step 1-2). The refined stack of unwrapped data is inverted to obtain the displacement time series and velocity (Step 1-3), followed by estimation of the velocity standard deviation (STD) (Step 1-4) and masking of noisy pixels based on several noise indices (Step 1-5). Finally, a spatiotemporal filter is applied to the time series to mitigate the residual noise and derive the filtered time series and velocity (Step 1-6). All steps can be executed from the command line or using a batch script with defined parameters. The derived results (i.e., velocity and time series) can be easily displayed by an interactive time series viewer and exported to GeoTIFF, KMZ, or text format.
LiCSBAS is written in Python3 and the Bourne Again Shell (bash) and operates without relying on any commercial software. The source codes are available on GitHub [31]. Although Steps 0-1 and 0-2 are dedicated specifically to manipulating LiCSAR products, interferograms from satellites other than Sentinel-1 and generated using a variety of software packages can be processed from Step 0-3, assuming the required metadata are available and files are converted into a compatible format Remote Sens. 2020, 12, 424 4 of 29 (4-byte, single-precision floating-point format). Future LiCSBAS enhancements may include functions to prepare input data from InSAR archives other than LiCSAR (e.g., ARIA, etc.).

Prepare Stack of Unwrapped Data
2.3.1.
Step 0-1: Download LiCSAR Products GeoTIFF files of unwrapped phase and coherence are published on the COMET-LiCS web portal [30] on a consistent geographic frame basis (typically 250 × 250 km), consisting of both ascending and descending data (which are treated independently). The GeoTIFF files for all or a specified time period (normally from late 2014 at the earliest) for the LiCSAR frame of interest are automatically downloaded in this step. Some available metadata (e.g., LOS unit vectors and perpendicular baselines) are also included in the download.

2.3.2.
Step 0-2: Convert GeoTIFF (and Downsample) The downloaded GeoTIFF files are converted to a single-precision floating-point format without header information for the following processing steps. Downsampling (i.e., multilooking) can also be applied if the full resolution of the LiCSAR products (~100 m) is not needed. Downsampling is achieved by boxcar averaging and returns null if more than half of the pixels in the boxcar window are null. Downsampling results in faster overall processing, a reduction in file size, and is also useful for quickly testing processing parameters and making adjustments prior to running the processing at a full or higher resolution.

2.3.1.
Step 0-1: Download LiCSAR Products GeoTIFF files of unwrapped phase and coherence are published on the COMET-LiCS web portal [30] on a consistent geographic frame basis (typically 250 × 250 km), consisting of both ascending and descending data (which are treated independently). The GeoTIFF files for all or a specified time period (normally from late 2014 at the earliest) for the LiCSAR frame of interest are automatically downloaded in this step. Some available metadata (e.g., LOS unit vectors and perpendicular baselines) are also included in the download.

2.3.2.
Step 0-2: Convert GeoTIFF (and Downsample) The downloaded GeoTIFF files are converted to a single-precision floating-point format without header information for the following processing steps. Downsampling (i.e., multilooking) can also be applied if the full resolution of the LiCSAR products (~100 m) is not needed. Downsampling is achieved by boxcar averaging and returns null if more than half of the pixels in the boxcar window are null. Downsampling results in faster overall processing, a reduction in file size, and is also useful for quickly testing processing parameters and making adjustments prior to running the processing at a full or higher resolution. Interferograms are contaminated by phase noise from path delays due to spatially-and temporally-varying tropospheric water vapor, temperature, and pressure, which can reach 10 cm or more [39,40]. In standard time series analysis algorithms, spatiotemporal filtering (i.e., high-pass in time and low-pass in space) is used to separate spatially correlated noise from deformation signals under the assumption that the noise component is uncorrelated in time [3,5,7]. This option is also available in LiCSBAS (Step 1-6 explained in Section 2.4.6), but the assumption that noise is uncorrelated in time is not always correct. If the tropospheric delay is correlated in time and/or the deformation is uncorrelated in time, the deformation time series may be poorly recovered [41][42][43]. Furthermore, temporal filtering does not improve the accuracy of the average velocities retrieved from the time series. Mitigation of the tropospheric noise before the spatiotemporal filtering using independent data improves the accuracy of the derived displacement time series and the average velocities.
GACOS provides tropospheric delay maps derived from the high-resolution European Centre for Medium-Range Weather Forecasts (ECMWF) data [44] for correction of the tropospheric noise present in interferograms [38]. The GACOS corrections are globally available in near real-time, and users can request the delay maps for specific areas and times through a web application system [45] and download these corrections once they have been produced.
Step 0-3 applies the GACOS corrections to the stack of unwrapped interferograms. Currently, users need to separately prepare the GACOS data for the processing area and acquisition dates beforehand. However, a one-stop service for LiCSAR is being planned in which the GACOS data for all Sentinel-1 acquisitions will be provided along with the interferograms via the COMET-LiCS web portal.
This step is optional but recommended to improve the accuracy of the time series (see Section 3.4). However, since the GACOS corrections may have a negative impact in some cases [46,47], users should check if reasonable noise mitigation has been achieved after applying the corrections. LiCSBAS calculates the STD of unwrapped phases before and after the GACOS corrections and their reduction rates for each interferogram, as well as creates quick-look images of unwrapped interferograms before and after the correction and the correction itself to aid in checking the GACOS effects ( Figure S1). In case the GACOS corrections have a negative impact, other methods (e.g., linear, power-law [40], or spatially varying scaling [48] methods for the topography-correlated component) could be useful and may be integrated in the future.

2.3.4.
Step 0-4: Mask Interferograms (Optional) This step masks specified rectangular areas (e.g., isolated islands) in the stack of unwrapped data ( Figure S2), which is useful if these areas have unwrapping errors. Masking may improve the network refinement at Step 1-2 (see Sections 2.4.2 and 3.3).

2.3.5.
Step 0-5: Clip Interferograms (Optional) This step clips a specified rectangular area from the stack of unwrapped data ( Figure S3). Clipping is recommended if the area of interest does not cover the entire frame, because it decreases the processing time and results in smaller file sizes (see Section 5.1). Note that a stable reference area should be included in the clipped area. Clipping may also improve network refinement in Step 1-2 by removing unnecessary areas where unwrapping errors might be present (see Section 4).

2.4.1.
Step 1-1: Quality Check LiCSAR products may contain not only useful coherent interferograms but also severely decorrelated interferograms due to factors such as vegetation or snow cover that result in a very low percentage of valid unwrapped pixels ( Figure S4a,b). Moreover, Sentinel-1 SLC data distributed by ESA do not always cover an entire LiCSAR frame (i.e., some bursts may be missing in a frame).
Interferograms derived from data with missing bursts may also have low coverage ( Figure S4c), and these incomplete data with low coverage negatively impact LiCSBAS processing and should be flagged as bad data. In Step 1-1, the bad data are identified based on statistics including the average coherence and the percentage of valid unwrapped pixels.

2.4.2.
Step 1-2: Network Refinement by Loop Closure Unwrapped data may include unwrapping errors, which can cause significant errors in the derived time series and should therefore be removed or corrected beforehand. Several approaches have been proposed to identify and correct unwrapping errors, taking advantage of the redundancy of a network of interferograms and loop phase closure [21,49,50]. However, manual identification or correction is too time-consuming for a full stack of LiCSAR interferograms [49]. Automatic correction can work well only if the redundancy of the network of interferograms is sufficiently high and not too many unwrapping errors exist [21,50], which may not always be the case for large-scale, automatic processing. Moreover, automatic unwrapping error correction can produce false corrections, which are difficult to identify and result in worse errors in the derived time series. Therefore, for now, we have decided to take a conservative approach; we automatically detect and completely remove interferograms that have many unwrapping errors by evaluating the loop phase closure on an image-by-image basis, rather than correcting unwrapping errors on a pixel-by-pixel basis.
Supposing that we have three images (φ 1 , φ 2 , and φ 3 ) and three unwrapped interferograms (φ 12 , φ 23 , and φ 13 ), a loop phase (also called closure phase or phase triplet) is calculated by [49] If there are no unwrapping errors in the three interferograms, the loop phase should be close to zero (but not exactly zero due to the effects of multilooking, filtering, soil moisture change, etc. [51]). On the other hand, if one (or more) of the interferograms contains unwrapping errors, the loop phase will usually be close to an integer multiple of 2π ( Figure 2). The root mean square (RMS) of the loop phase image indicates how many pixels with unwrapping errors are included in the loop. All possible loop phases and the RMS are calculated, and problematic loops with larger RMS than a defined threshold (e.g., 1.5 rad) can be identified. If all loops associated with an interferogram are problematic, this is considered a bad interferogram that likely includes many unwrapping errors, and it is removed from further processing ( Figure 2).
Note that this step evaluates unwrapping errors for each interferogram but not for each pixel. Interferograms that are removed during Step 1-2 may contain correctly unwrapped pixels. However, this is not a significant problem if sufficient interferograms still remain after network refinement. Conversely, remaining interferograms may still contain some unwrapping errors, which would degrade the derived time series. These errors generate unclosed loops on a pixel-by-pixel basis. We count the number of unclosed loops for each pixel and can use this for masking pixels with unwrapping errors after the small baseline (SB) inversion at Step 1-5 (Section 2.4.5).
If two interferograms in a loop have opposite unwrapping errors in the same area, the loop phase is falsely closed, and the bad interferograms may not be automatically identified. However, network redundancy assures that other loops associated with the unidentified bad interferograms are not closed and will be counted on a pixel-by-pixel basis and identified at the later masking step (Step 1-5, Section 2.4.5). Users can also manually identify and remove the unidentified bad interferograms by checking the loop closure images ( Figure 2) and then re-run Step 1-2.
After network refinement, a preliminary reference point is selected. At the reference point, all interferograms must have valid (unmasked) and error-free unwrapped data for the following SB network inversion. To identify the reference point, the RMS of all the loop phases for each pixel is calculated, and the pixel with the smallest RMS is selected.  (1), respectively. Black and red lines between them denote which three interferograms are used to calculate each loop phase. Three loop phases with a red frame (Φ 124 , Φ 234 , Φ 245 ) have areas with ±2π phase and a large (>2 rad) root mean square (RMS). Then, the interferogram 24 with a red frame can be identified as a bad interferogram including significant unwrapping errors, because all loops associated with it are problematic, whilst the other interferograms generate at least one nonproblematic loop.

Step 1-3: Small Baseline Network Inversion
In order to perform an estimate of the velocity of a surface pixel through time based upon a series of displacement data, we perform an SB inversion on the network of interferograms. Suppose that we have a stack of M-unwrapped interferograms = [ 1 , … , ] T produced from N images acquired at ( 0 , … , −1 ), N-1 incremental displacement vector = [ 1 , … , −1 ] T (i.e., mi is the incremental displacement between time ti-1 and ti) can be derived by solving

= Gm
(2) where G is a M×(N-1) design matrix of zeros and ones describing the relationship between the network of the interferograms and incremental displacements, considering that the unwrapped interferogram (i.e., displacement between two acquisitions) is the sum of corresponding successive incremental displacements [52]. Cumulative displacements for each acquisition (i.e., the displacement time series) are simply calculated by summing the incremental displacements. The mean displacement velocity is then derived from the cumulative displacements by least-squares. Even though Sentinel-1 data are constantly being acquired with short spatial and temporal baselines in most areas (unlike with predecessors such as Envisat and ERS), there could still be gaps in the SB network (i.e., separated subsets of the network with no interferograms spanning them) due to decorrelation associated with factors such as vegetation, cultivation, or snow cover, and especially after the network refinement Step 1-2. Additionally, some frames may have large acquisition gaps (several months or longer). It is noteworthy that even if the refined network of interferograms is fully connected on an image-by-image basis, network gaps may still exist on a pixel-by-pixel basis. This is because each interferogram has a different coverage of valid unwrapped data due to the coherencebased masking applied to the LiCSAR products (see Section 2.1). In such cases, the design matrix G is rank deficient; equation (2) can nevertheless be solved using singular value decomposition (SVD). However, the minimum norm solution provided by SVD gives incremental displacements of zero in any network gaps, which is not realistic even for very small gaps (e.g., 6 days) where the true displacement would be almost zero, because noise in the interferograms between the gaps would not be close to zero.  (1), respectively. Black and red lines between them denote which three interferograms are used to calculate each loop phase. Three loop phases with a red frame (Φ 124 , Φ 234 , Φ 245 ) have areas with ±2π phase and a large (>2 rad) root mean square (RMS). Then, the interferogram φ 24 with a red frame can be identified as a bad interferogram including significant unwrapping errors, because all loops associated with it are problematic, whilst the other interferograms generate at least one nonproblematic loop.

Step 1-3: Small Baseline Network Inversion
In order to perform an estimate of the velocity of a surface pixel through time based upon a series of displacement data, we perform an SB inversion on the network of interferograms. Suppose that we have a stack of M-unwrapped interferograms d = [d 1 , . . . , d M ] T produced from N images acquired at (t 0 , . . . , t N−1 ), N-1 incremental displacement vector m = [m 1 , . . . , m N−1 ] T (i.e., m i is the incremental displacement between time t i-1 and t i ) can be derived by solving where G is a M×(N-1) design matrix of zeros and ones describing the relationship between the network of the interferograms and incremental displacements, considering that the unwrapped interferogram (i.e., displacement between two acquisitions) is the sum of corresponding successive incremental displacements [52]. Cumulative displacements for each acquisition (i.e., the displacement time series) are simply calculated by summing the incremental displacements. The mean displacement velocity is then derived from the cumulative displacements by least-squares. Even though Sentinel-1 data are constantly being acquired with short spatial and temporal baselines in most areas (unlike with predecessors such as Envisat and ERS), there could still be gaps in the SB network (i.e., separated subsets of the network with no interferograms spanning them) due to decorrelation associated with factors such as vegetation, cultivation, or snow cover, and especially after the network refinement Step 1-2. Additionally, some frames may have large acquisition gaps (several months or longer). It is noteworthy that even if the refined network of interferograms is fully connected on an image-by-image basis, network gaps may still exist on a pixel-by-pixel basis. This is because each interferogram has a different coverage of valid unwrapped data due to the coherence-based masking applied to the LiCSAR products (see Section 2.1). In such cases, the design matrix G is rank deficient; Equation (2) can nevertheless be solved using singular value decomposition (SVD). However, the minimum norm solution provided by SVD gives incremental displacements of zero in any network gaps, which is not realistic even for very small gaps (e.g., 6 days) where the true displacement would be almost zero, because noise in the interferograms between the gaps would not be close to zero. To obtain the more realistic time series of the displacement even with a disconnected network, we follow the NSBAS method [19,50,53] that imposes a temporal constraint, where γ is a scaling (weighting) factor of the temporal constraint, and we assume displacements are linear (d = vt + c). Solutions within the connected parts of the network are minimally affected by the temporal constraint provided γ is small (e.g., 0.0001). Thus, the temporal constraint part has an impact only on the connection across gaps in the network. Therefore, Equation (3) can be used for pixels with fully connected networks, as well as those with gaps. If the actual displacements significantly deviate from a linear function (e.g., a coseismic jump due to an earthquake or acceleration/deceleration in displacement due to volcanic deformation), the solution between the gaps would become highly unreliable. In that case, functions other than linear, which we assume here (e.g., Heaviside, exponential, sinusoidal, etc.), might be better, but a priori information regarding the style of deformation is required to choose the appropriate function. Other functions are not currently implemented. Regardless, we can never retrieve a solution supported by the real observed data between the gaps. We should carefully consider the gaps when we interpret the time series of the displacement. LiCSBAS identifies and records the gaps for all pixels and increments, and this information can be easily displayed in the time series viewer to aid in the interpretation (see Section 2.5).

Step 1-4: Estimate Standard Deviation of the Velocity by Bootstrap
This step estimates the STD of the velocity from the cumulative displacements derived in Step 1-3 using the percentile bootstrap method [22,54]. The bootstrap approach creates a randomly resampled dataset with data replacement. We repeat the bootstrap resampling 100 times, compute the velocities during each iteration, and then calculate the STD for the ensemble of velocity solutions. High STD estimates indicate a noisy or nonlinear displacement time series. Note that the STD could be underestimated if the network is not fully connected, due to the temporal constraint in the SB inversion (Section 2.4.3).

Step 1-5: Mask Noisy Pixels in the Time Series
Time series can be obtained for all pixels with a sufficient number of available unwrapped data. However, the quality of a stack of unwrapped data may differ significantly from pixel to pixel. Some pixels may be very unreliable due to factors including remaining gaps and unwrapping errors.
To address this problem, Step 1-5 generates a mask for the time series to differentiate between "good" and "bad" pixels using several noise indices derived during the previous steps ( Figure 3 and Table 1). The pixel is masked if any of the values of the noise indices for a pixel is worse (larger or smaller) than a specified threshold. The threshold for each noise index can be manually adjusted to obtain a better (i.e., harsher or gentler) mask for the time series and velocity.
Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 30  Table 1). Numbers shown in the parentheses next to the titles of each noise index are the applied threshold value.    Table 1). Numbers shown in the parentheses next to the titles of each noise index are the applied threshold value.  2.4.6.
Step 1-6: Spatiotemporal Filtering of Time Series The derived time series still includes some noise terms, including residual tropospheric noise, ionospheric noise, and orbital errors, which are, in general, spatially correlated and temporally uncorrelated. A spatiotemporal filter (i.e., high-pass in time and low-pass in space) can be applied to attempt to separate these noise components from the displacement time series [3,6]. To achieve this filtering, we apply a one-and two-dimensional Gaussian kernel in time and space, respectively. Moreover, if the spatial scale of the target deformation signal is sufficiently small compared to the processing area (e.g., subsidence, landslides, or volcanic deformation), removing a ramp (linear, bilinear, or quadratic) from the entire area may help mitigate the noise.
It is noteworthy that too-strong a spatiotemporal filter (i.e., with long and short filter widths in time and space, respectively) can remove not only the noise components but also the actual displacement signals of interest, especially for cases where the displacement time series is nonlinear and the temporal width of the filter is too long compared to the time scale of the nonlinear displacement [6]. On the other hand, the temporal filter width needs to be sufficiently long, compared to the sampling interval, or no effect is produced. These issues have made separating noise from nonlinear displacements using infrequent data provided by earlier InSAR satellites (e.g., Envisat and ERS) challenging. However, with Sentinel-1, we can apply a relatively short temporal filter width and potentially more easily detect real nonlinear displacement. In LiCSBAS Step 1-6, we use 3dt as the default temporal filter width (1σ of Gaussian), where dt is the average temporal sampling interval. For example, if the average temporal baseline of interferograms is 12 days, a temporal filter width of 36 days is applied to the displacement time series, which is sufficiently short to highlight seasonal displacements associated with the hydrological cycle [56,57].

Visualization of the Results
Interpretation of the derived time series is as important as the accurate derivation of the time series. However, the huge amount of associated data (e.g.,~4 GB for 100 images of an entire frame) often makes smooth visualization and intuitive interpretation problematic. Selection of an appropriate reference area is also a critical step towards obtaining meaningful and accurate displacements for interpretation, as all measurements are relative rather than in an absolute reference frame [23]. However, it is not possible to interactively examine and find out the appropriate reference by trial and error in most software.
LiCSBAS is equipped with an interactive time series viewer composed of two windows (graphical user interfaces). The first image window displays velocity, cumulative displacements, and noise indices (Figure 4a). When clicking on a pixel of interest, the associated time series with and without the spatiotemporal filter is promptly plotted in the second time series window (Figure 4b). Any network gaps are clearly displayed, which alerts the user to carefully interpret the time series connected by the temporal constraint (see Section 2.4.3). The reference area is shown by a black rectangle, and it can be interactively changed and examined in the image window.
The viewer can smoothly display the results of the time series analysis, even with a huge file size, without consuming a large amount of computer memory, owing to the HDF5 file format, from which the requested part of a large dataset can be quickly read.

Southern Tohoku, Japan
We processed an entire LiCSAR frame covering Southern Tohoku and Northern Kanto, Japan in a descending orbit (LiCSAR frame ID 046D_05292_131313; acquisition time of day 20:42 UTC and 05:42 JST; Figure 5a). The observation period is from 25 November 2014 to 14 July 2019 (~4.6 years), and the associated network consists of 104 images and 306 interferograms. The acquisition interval is primarily 24 days prior to 18 February 2017 and 12 days after this date due to the increased observational capacity from the availability of Sentinel-1B, although all the data were acquired by Sentinel-1A (Figure 5b).
Inland areas within this frame are mostly mountainous and covered with dense vegetation. Furthermore, it snows during winter months. Therefore, it is quite difficult to obtain the full time series using C-band SAR data due to severe decorrelation. There are several plains and urban areas that are suitable for C-band data: the Kanto Plain in the south and the Echigo Plain in the northwest (Figure 5a).
A portion of the megathrust fault that slipped during the 2011 Mw 9.0 Tohoku earthquake is located offshore, just beyond the frame (Figure 5a) [58], although the frame encompasses a region that experienced a few meters of the eastward coseismic motion based on continuous Global

Southern Tohoku, Japan
We processed an entire LiCSAR frame covering Southern Tohoku and Northern Kanto, Japan in a descending orbit (LiCSAR frame ID 046D_05292_131313; acquisition time of day 20:42 UTC and 05:42 JST; Figure 5a). The observation period is from 25 November 2014 to 14 July 2019 (~4.6 years), and the associated network consists of 104 images and 306 interferograms. The acquisition interval is primarily 24 days prior to 18 February 2017 and 12 days after this date due to the increased observational capacity from the availability of Sentinel-1B, although all the data were acquired by Sentinel-1A (Figure 5b).
Inland areas within this frame are mostly mountainous and covered with dense vegetation. Furthermore, it snows during winter months. Therefore, it is quite difficult to obtain the full time series using C-band SAR data due to severe decorrelation. There are several plains and urban areas that are suitable for C-band data: the Kanto Plain in the south and the Echigo Plain in the northwest (Figure 5a).
Step 1-2, we masked the northwestern portion of the frame including the Echigo Plain and the adjacent mountainous regions in Step 0-4 ( Figure S2; its impact is discussed in Section 3.3). As a result, no interferograms were removed at Step 1-1, and only two interferograms (20141125-20141219 and 20160224-20160506) were removed at Step 1-2; there were no gaps in the network. We also applied the GACOS correction at Step 0-3 (its impact is discussed in Section 3.4). The filter widths used in Step 1-6 are 2 km in space and 49 days (3 ��� ) in time.  A portion of the megathrust fault that slipped during the 2011 M w 9.0 Tohoku earthquake is located offshore, just beyond the frame (Figure 5a) [58], although the frame encompasses a region that experienced a few meters of the eastward coseismic motion based on continuous Global Navigation Satellite System (GNSS) observations [58]. The region has also experienced significant, ongoing postseismic deformation [59]. The rate of postseismic deformation has roughly converged to a constant velocity since around 2014 [60]. GNSS-derived velocities across the northern portion of the frame reveal significant eastward displacements at a rate of~50 mm/yr from 2014-2019, compared to the south (Figure 5a). Vertical displacements are generally smaller than horizontal, but an east-to-west transition from uplift to subsidence with rates on the order of >10 mm/yr occurs across the region covered by the northern portion of the frame.
As the inland mountainous regions tend to be decorrelated by vegetation and snow, the Echigo Plain is often isolated from the Kanto Plain in the unwrapped interferograms, which results in many unwrapping errors. To reduce the impact of these errors and improve the network refinement during Step 1-2, we masked the northwestern portion of the frame including the Echigo Plain and the adjacent mountainous regions in Step 0-4 ( Figure S2; its impact is discussed in Section 3.3). As a result, no interferograms were removed at Step 1-1, and only two interferograms (20141125-20141219 and 20160224-20160506) were removed at Step 1-2; there were no gaps in the network. We also applied the GACOS correction at Step 0-3 (its impact is discussed in Section 3.4). The filter widths used in Step 1-6 are 2 km in space and 49 days (3dt) in time.

Results
Despite having masked almost all pixels in the mountainous areas at Step 1-5 due to decorrelation, more than a million valid pixels, mainly distributed on the Kanto Plain and along northern basins and plains, still remain ( Figure 6a). As a measure of InSAR data quality, we compare the LOS velocities with those stemming from the GNSS measurements with the same reference point of the 950217 GNSS station. The InSAR-derived velocities successfully image the broad tectonic deformation (i.e., LOS decrease in the northeast due to the combination of eastward and vertical uplift motions) and are largely consistent with the GNSS velocities ( Figure 6a). The average and STD of the difference between InSAR and GNSS (n = 56) is 0.1 mm/yr and 1.9 mm/yr, respectively. We detect no significant systematic errors caused by factors such as ionospheric noise and orbital errors, despite no ionospheric correction having being applied in the LiCSAR products (Figure 6a,b). We also tried estimating the best-fitting quadratic polynomial surface through the differences between InSAR and GNSS velocities, at all GNSS stations, and subtracting it from the InSAR-derived velocities, which slightly decreased the STD to 1.6 mm/yr (Figure 6b). In this case, the impact of ionospheric noise and orbital errors is quite small, although, in some cases, these error sources could become significant, particularly at low latitudes on ascending tracks [61].
We also compare the LOS time series data and find that the InSAR-derived results are quite consistent with the GNSS measurements, even at the farthest GNSS site (950179) located 150 km away from the reference (950217; Figure 6c, Figure S5 and S6). The STD of the difference of the time series at each GNSS station shows no clear correlation with distance from the reference point (Figure 6d) once the GACOS correction has been applied, and the average of the STD of the difference at all GNSS stations is 9.0 mm, which represents the root-sum-square of the uncertainties of the InSAR-and GNSS-derived time series for each epoch. We also calculated the STD of the GNSS measurements for each epoch from the differences from those of 30-days moving average, and the average of the STD among all GNSS stations is 5.0 mm. Considering these, the uncertainty of the InSAR-derived time series for each epoch is 7.5 mm.
The results for Southern Tohoku demonstrate the potential of LiCSAR products and LiCSBAS to detect broad tectonic deformation independently without the help of GNSS and to contribute to dense and detailed tectonic strain mapping and global earthquake hazard assessment [62].

Impact of Masking and a Network Gap
Here, we investigate the impact of masking during Step 0-4. Without masking the northwestern area at Step 0-4, many unwrapping errors remain in the interferograms, and 47 interferograms are excluded during the loop closure at Step 1-2, generating a gap in the network between 18 and 30 June, 2017 ( Figure S7). This means the masking has improved the network refinement and avoided the gap at the cost of a much-reduced extent of spatial coverage.
We compare the time series and velocities derived from the network, with and without the gap (i.e., with and without masking at Step 0-4). The time series at 950179 with the gap is consistent with the no-gap time series and the GNSS data, but slightly (~2 mm) biased across the gap (Figure 6c and Figure S6). This bias is caused by the interpolation associated with SB inversion linear temporal

Impact of Masking and a Network Gap
Here, we investigate the impact of masking during Step 0-4. Without masking the northwestern area at Step 0-4, many unwrapping errors remain in the interferograms, and 47 interferograms are excluded during the loop closure at Step 1-2, generating a gap in the network between 18 and 30 June, 2017 ( Figure S7). This means the masking has improved the network refinement and avoided the gap at the cost of a much-reduced extent of spatial coverage.
We compare the time series and velocities derived from the network, with and without the gap (i.e., with and without masking at Step 0-4). The time series at 950179 with the gap is consistent with the no-gap time series and the GNSS data, but slightly (~2 mm) biased across the gap (Figure 6c and Figure S6). This bias is caused by the interpolation associated with SB inversion linear temporal constraint (see Section 2.4.3). In this case, the bias is very small, because the overall trend is almost linear (i.e., constant velocity). However, if the actual displacement is highly nonlinear, the bias would become large and need to be carefully considered in the interpretation of the derived time series. The impact differs from point to point (Figure 6d and Figure S5).
The overall velocity, with and without the network gap, is also broadly similar, although significant differences up to~5 mm/yr exist in some areas ( Figure S8). The STD of the difference of the LOS velocity between InSAR with the gap and GNSS is 2.9 mm/yr, which is significantly larger than that without the gap (1.9 mm/yr). This is because the time length of the connected network, which has an impact on the precision of the velocity, is shortened from 4.6 years to 2.6 years by the gap (see Section 3.5).

Impact of the GACOS Correction
We investigate the impact of the GACOS correction at Step 0-3. In terms of SB interferograms, the STD of unwrapped phases for each entire interferogram is generally reduced (from 6.7 rad to 4.2 rad on average, from 6.0 rad to 3.9 rad on median), which indicates the GACOS correction significantly mitigated the tropospheric noise (Figure 7 and Figure S1).
The overall velocity, with and without the network gap, is also broadly similar, although significant differences up to ~5 mm/yr exist in some areas ( Figure S8). The STD of the difference of the LOS velocity between InSAR with the gap and GNSS is 2.9 mm/yr, which is significantly larger than that without the gap (1.9 mm/yr). This is because the time length of the connected network, which has an impact on the precision of the velocity, is shortened from 4.6 years to 2.6 years by the gap (see Section 3.5).

Impact of the GACOS Correction
We investigate the impact of the GACOS correction at Step 0-3. In terms of SB interferograms, the STD of unwrapped phases for each entire interferogram is generally reduced (from 6.7 rad to 4.2 rad on average, from 6.0 rad to 3.9 rad on median), which indicates the GACOS correction significantly mitigated the tropospheric noise (Figure 7 and Figure S1).
Regarding the result of the LiCSBAS processing, without the GACOS correction, the overall velocity is quite similar, and the STD of the difference in LOS velocity between InSAR and GNSS is 2.4 mm/yr, which is slightly larger than that with the GACOS correction (1.9 mm/yr; Figure 8).
The time series without GACOS at the distal position 950179 shows clear false seasonal displacements, presumably due to tropospheric noise, which has a seasonal correlation and cannot be removed by the spatiotemporal filter at Step 1-6 ( Figure 6c, Figure S5, and Figure S6). The impact of the residual tropospheric noise depends on the distance from the reference point; far points tend to have larger time series STDs (Figure 6d). These results suggest that the GACOS correction plays an important role in more accurately retrieving the time series, especially at distances far from the reference point.  Regarding the result of the LiCSBAS processing, without the GACOS correction, the overall velocity is quite similar, and the STD of the difference in LOS velocity between InSAR and GNSS is 2.4 mm/yr, which is slightly larger than that with the GACOS correction (1.9 mm/yr; Figure 8).
The time series without GACOS at the distal position 950179 shows clear false seasonal displacements, presumably due to tropospheric noise, which has a seasonal correlation and cannot be removed by the spatiotemporal filter at Step 1-6 ( Figure 6c, Figure S5, and Figure S6). The impact of the residual tropospheric noise depends on the distance from the reference point; far points tend to have larger time series STDs (Figure 6d). These results suggest that the GACOS correction plays an important role in more accurately retrieving the time series, especially at distances far from the reference point.

Evolution of Velocity Uncertainty
The velocity uncertainty can be estimated by

= 2√3
0.5 (4) where , N, and T are the uncertainty at each epoch, the number of image data, and the observation time period length, respectively, under the assumption that the measurement errors are uncorrelated [63]. It is noteworthy that, if network gaps are present, T is not equivalent to the time length of the observation period but, rather, the timespan of the connected network. Therefore, the presence of a gap can significantly degrade the accuracy of the velocity estimation (see Section 3.3). We calculate InSAR LOS velocities starting at the first Sentinel-1 acquisition (i.e., 25 November, 2014) to each subsequent acquisition epoch; GNSS LOS velocities for the corresponding periods; and the STD of the difference between these estimates. We also calculate the STD of the InSAR LOS velocity at each subsequent acquisition epoch using the percentile bootstrap method [22,54].
The temporal evolution of the InSAR-GNSS and InSAR-bootstrap STDs are largely consistent with the theoretical velocity uncertainty model calculated using Equation (4) with of 9.0 mm or 7.5 mm (calculated in Section 3.2 by comparison with GNSS), respectively, and a sampling interval ( = /( − 1)) of 24 days for ~2.3 years (i.e., until February 2017; Figure 9). An excursion in InSAR-GNSS STDs around 1.5-2.0 years (i.e., May to November 2016) might be caused by residual noise in the InSAR time series (Figure 6c, Figure S5 and S6). The STD of InSAR-GNSS converges to ~2 mm/yr, not 0 mm/yr, presumably because of a contribution from colored noise components in InSAR and/or GNSS [63] or discrepancies in the location of the measurement points between InSAR (i.e., coherent sum of the scatters in each resolution cell) and pointwise GNSS.
Our prediction of the velocity uncertainty is useful for estimating how small displacements can be detected given a certain amount of data or how long a time period is required to achieve a certain velocity precision. For example, if the sampling interval is 70 (corresponding to past Envisat  Figures 6a, 7a and 8a). Note that the color scale is different from Figures 6a and 8a.

Evolution of Velocity Uncertainty
The velocity uncertainty σ v can be estimated by where σ e , N, and T are the uncertainty at each epoch, the number of image data, and the observation time period length, respectively, under the assumption that the measurement errors are uncorrelated [63]. It is noteworthy that, if network gaps are present, T is not equivalent to the time length of the observation period but, rather, the timespan of the connected network. Therefore, the presence of a gap can significantly degrade the accuracy of the velocity estimation (see Section 3.3). We calculate InSAR LOS velocities starting at the first Sentinel-1 acquisition (i.e., 25 November, 2014) to each subsequent acquisition epoch; GNSS LOS velocities for the corresponding periods; and the STD of the difference between these estimates. We also calculate the STD of the InSAR LOS velocity at each subsequent acquisition epoch using the percentile bootstrap method [22,54].
The temporal evolution of the InSAR-GNSS and InSAR-bootstrap STDs are largely consistent with the theoretical velocity uncertainty model calculated using Equation (4) with σ e of 9.0 mm or 7.5 mm (calculated in Section 3.2 by comparison with GNSS), respectively, and a sampling interval (dt = T/(N − 1)) of 24 days for~2.3 years (i.e., until February 2017; Figure 9). An excursion in InSAR-GNSS STDs around 1.5-2.0 years (i.e., May to November 2016) might be caused by residual noise in the InSAR time series (Figure 6c, Figures S5 and S6). The STD of InSAR-GNSS converges to~2 mm/yr, not 0 mm/yr, presumably because of a contribution from colored noise components in InSAR and/or GNSS [63] or discrepancies in the location of the measurement points between InSAR (i.e., coherent sum of the scatters in each resolution cell) and pointwise GNSS. acquisitions every two cycles), 24, 12, or 6 (current best-case with Sentinel-1A/B) days, it would take 3.1, 2.2, 1.8, or 1.4 years to achieve a velocity STD of 2 mm/yr, respectively.

Echigo Plain
The Echigo (Niigata) Plain is located in the northwest of the frame ID 046D_05292_131313 (Figure 5a and Figure 10). It is an alluvial plain formed by the Shinano and Agano Rivers, surrounded by sea to the north and higher topography elsewhere ( Figure 10). Snow cover is present in winter months in both the mountainous areas and the lowlands.
In addition to processing the whole frame (Section 3), we performed LiCSBAS processing after Step 0-3 (GACOS correction) for a rectangular subset covering the Echigo Plain region by clipping this area from the entire frame (Step 0-5). In addition, we masked (Step 0-4) the mountainous areas to the southeast, where most of the interferograms are decorrelated by vegetation and snow. Only one interferogram (20171203-20171227) was then excluded, due to unwrapping errors at Step 1-2, leaving 305 unwrapped interferograms for time series and velocity estimation. We selected GNSS station 960566 as a reference point, because the observed GNSS displacements are stable compared to other stations in this region and the associated InSAR-derived time series are reliable. As a result, clear displacement signals are detected at Niigata, Ojiya, and Sanjo cities (Figure 10b). In the following sections, we describe the details for each area as examples to test the extraction of both long-term subsidence signals and seasonally varying recoverable signals from the InSAR time series corroborated by independent datasets. Our prediction of the velocity uncertainty is useful for estimating how small displacements can be detected given a certain amount of data or how long a time period is required to achieve a certain velocity precision. For example, if the sampling interval is 70 (corresponding to past Envisat acquisitions every two cycles), 24, 12, or 6 (current best-case with Sentinel-1A/B) days, it would take 3.1, 2.2, 1.8, or 1.4 years to achieve a velocity STD of 2 mm/yr, respectively.

Echigo Plain
The Echigo (Niigata) Plain is located in the northwest of the frame ID 046D_05292_131313 (Figures 5a and 10). It is an alluvial plain formed by the Shinano and Agano Rivers, surrounded by sea to the north and higher topography elsewhere ( Figure 10). Snow cover is present in winter months in both the mountainous areas and the lowlands.
In addition to processing the whole frame (Section 3), we performed LiCSBAS processing after Step 0-3 (GACOS correction) for a rectangular subset covering the Echigo Plain region by clipping this area from the entire frame (Step 0-5). In addition, we masked (Step 0-4) the mountainous areas to the southeast, where most of the interferograms are decorrelated by vegetation and snow. Only one interferogram (20171203-20171227) was then excluded, due to unwrapping errors at Step 1-2, leaving 305 unwrapped interferograms for time series and velocity estimation. We selected GNSS station 960566 as a reference point, because the observed GNSS displacements are stable compared to other stations in this region and the associated InSAR-derived time series are reliable. As a result, clear displacement signals are detected at Niigata, Ojiya, and Sanjo cities (Figure 10b). In the following sections, we describe the details for each area as examples to test the extraction of both long-term subsidence signals and seasonally varying recoverable signals from the InSAR time series corroborated by independent datasets.

Niigata City
Niigata City was previously characterized by rapid subsidence (~50 cm/yr) associated with groundwater pumping and natural gas extraction around 1960 [64]. Since then, prevention measures have been taken to prevent subsidence, including regulation of extraction and restoration of groundwater. As a result, the subsidence rates have decreased, although slow subsidence (~2 cm/yr or less) is still detected around the mouth of the Agano River using conducted dense leveling surveys (Figure 11a) [65].
Here, we compare the vertical velocity between InSAR and leveling. We calculate the levelingderived velocity at the benchmarks from five leveling surveys spanning September 2014 to September 2018. We project the InSAR-derived LOS velocity to vertical by taking into account the incidence angles under the assumption of no horizontal displacement. The reference point for both InSAR and leveling is set at benchmark II2025 (Figure 11a).
The derived vertical velocity shows an obvious correlation between InSAR and leveling with a degree of scattering (Figure 11b). The average and STD of the differences (n = 185) are 0.7 and 2.5 mm/yr, respectively. The time series of InSAR at benchmark A, where the most rapid subsidence has been observed, shows almost constant subsidence at a rate of 15.0 mm/yr, which is also consistent with leveling (13.8 mm/yr; Figure 11c). Note that the leveling results consist of only five data points and that the associated velocity uncertainty could be quite large compared to the GNSS and InSAR results. Therefore, the uncertainty of the InSAR-derived velocity estimate could be smaller than 2.5 mm/yr.

Niigata City
Niigata City was previously characterized by rapid subsidence (~50 cm/yr) associated with groundwater pumping and natural gas extraction around 1960 [64]. Since then, prevention measures have been taken to prevent subsidence, including regulation of extraction and restoration of groundwater. As a result, the subsidence rates have decreased, although slow subsidence (~2 cm/yr or less) is still detected around the mouth of the Agano River using conducted dense leveling surveys (Figure 11a) [65].
Here, we compare the vertical velocity between InSAR and leveling. We calculate the leveling-derived velocity at the benchmarks from five leveling surveys spanning September 2014 to September 2018. We project the InSAR-derived LOS velocity to vertical by taking into account the incidence angles under the assumption of no horizontal displacement. The reference point for both InSAR and leveling is set at benchmark II2025 (Figure 11a).
The derived vertical velocity shows an obvious correlation between InSAR and leveling with a degree of scattering (Figure 11b). The average and STD of the differences (n = 185) are 0.7 and 2.5 mm/yr, respectively. The time series of InSAR at benchmark A, where the most rapid subsidence has been observed, shows almost constant subsidence at a rate of 15.0 mm/yr, which is also consistent with leveling (13.8 mm/yr; Figure 11c). Note that the leveling results consist of only five data points and that the associated velocity uncertainty could be quite large compared to the GNSS and InSAR results. Therefore, the uncertainty of the InSAR-derived velocity estimate could be smaller than 2.5 mm/yr.

Ojiya City
Ojiya city is located on a terrace of the Shinano River (Figures 10 and 12). The Quaternary layer is composed of alluvium (thickness of 10-20 m), Holocene terrace deposits (5-20 m), lower terrace deposits (10-25 m), and Uonuma formation (> 450 m) from the top [66]. The deposits are mainly composed of gravel, sand, silt, and clay. The largest volumes of yearly groundwater usage in this area are devoted to snow-melting (38%) and industrial purposes (27%) [66]. The snow-melting system using the groundwater is laid along roads throughout the city (Figure 12a). The groundwater usage for snow-melting is concentrated in winter, which results in a sudden drop of the groundwater level every winter (Figure 12d and Figure S9).

Ojiya City
Ojiya city is located on a terrace of the Shinano River (Figures 10 and 12). The Quaternary layer is composed of alluvium (thickness of 10-20 m), Holocene terrace deposits (5-20 m), lower terrace deposits (10-25 m), and Uonuma formation (> 450 m) from the top [66]. The deposits are mainly composed of gravel, sand, silt, and clay. The largest volumes of yearly groundwater usage in this area are devoted to snow-melting (38%) and industrial purposes (27%) [66]. The snow-melting system using the groundwater is laid along roads throughout the city (Figure 12a). The groundwater usage Remote Sens. 2020, 12, 424 20 of 29 for snow-melting is concentrated in winter, which results in a sudden drop of the groundwater level every winter (Figure 12d and Figure S9). A continuous GNSS station (950240) has recorded large seasonal vertical displacements since 1996 (Figure 12d and Figure S9) with subsidence starting around December, when winter snowfall begins [67]. When the accumulated snow melts away around March, recovery of the subsidence begins, exponentially decays, and then converges in autumn. This temporal evolution of vertical displacements is consistent with the groundwater level in Well I, which is located 300 m northeast of the GNSS station (Figure 12d, Figures S9 and S10) [67]. The periodic subsidence is almost completely recovered every year, and interannual subsidence hardly remains. The depth of Well I is 200 m, and the screens are located at the depth of 134-189 m, corresponding to the Uonuma formation. These facts suggest the predominant mechanism of this vertical fluctuation is elastic compaction and expansion of the aquifer in the Uonuma formation, rather than inelastic [67,68].
Here, we reveal the spatial extent of the vertical fluctuations, as well as the detailed temporal evolution. We computed the InSAR-derived LOS time series and velocity with reference to the GNSS station 960566 (~26 km north-northeast from 950240). Although the velocity shows no significant long-term subsidence (Figure 12a), the time series clearly captures the periodic subsidence and retrieval temporally and spatially in detail (Figure 12c,d and Figure S11). The distribution of the fluctuating areas is almost the same every year.
The InSAR-derived time series at 950240 closely agrees with the time series measured by GNSS (Figure 12d). The STD of the difference between the two datasets is 9.5 mm, which is consistent with the result for the entire frame (Section 3.2), even with the large nonlinear periodic displacement. In February 2018, the subsidence was larger than the other years, presumably because the amount of snowfall was the largest that winter of 2018, and this is also clearly captured in the InSAR results (Figure 12c,d and Figure S11; see also Section 4.4).
In terms of the spatial distribution, there are two peaks of the subsidence either side of the river (Figure 12b,c). The eastern one has a smaller amount of the subsidence than the western one, which is consistent with the fact that the groundwater level change at the eastern Well II is smaller than that at the western Well I ( Figure S9). The shape of the western subsidence bowl is not asymmetric; the northeast and southwest edges of the subsidence area have sharp boundaries that may be controlled by subsurface geological boundaries.

Sanjo City
Large subsidence rates measured in Sanjo City (Figures 10b and 13a) correspond to time series (Figure 13b) that do not exhibit linear displacements like Niigata City (Figure 11c) but, rather, are characterized by episodic subsidence (~4 cm) in the winter of 2017/2018, which is not recovered thereafter. This suggests a mechanism that is not elastic but inelastic. In addition to the episodic displacement, periodic displacements, as in Ojiya City (i.e., subsidence in winter and recovery until autumn), are also seen in the time series (Figure 13b). This is expected, because Sanjo City employs the same type of snow-melting system that we previously described. The height difference of the periodic component is~3 cm, which is smaller than the Ojiya City signal (~6 cm or more).
One possible explanation for the episodic subsidence in the winter of 2017/2018 relates to overpumping of the groundwater for snow-melting. An unusually heavy snowfall occurred in this region on 30 January 2018, and the snow remained uncharacteristically deep (>100 cm) until the end of February, whilst the normal depth of snow is~60 cm (Figure 13b), recorded at the Nagaoka snow gauge,~25 km SSW from Sanjo City (Figure 10b). The episodic subsidence observed from 1 to 25 February is temporally consistent with the heavy snowfall.  Figure 10b).

Annual Displacement in Echigo Plain
Here, we briefly present the annual displacements over the Echigo Plain, which we estimate from the time series using least-squares. Although the actual periodic displacement is not strictly sinusoidal (Figure 12d and 13b), for simplicity, we use a simple sinusoidal function, along with the linear function, to estimate the amplitude and time offset (e.g., [57,69]).
Large amplitudes are detected in Ojiya City (~22 mm) and Sanjo City (~13 mm), as was previously shown in Sections 4.3 and 4.4 (Figure 14a). Some other areas also exhibit relatively large amplitudes, including Nagaoka City (~17 mm) and Shibata City (~11 mm). Elsewhere, only moderate amplitude signals ranging from 4-7 mm/yr are detected. The spatial pattern revealed by the amplitude map might reflect variations in geological properties. We should note that if the reference point had a periodic displacement, an apparent inverse periodic displacement would erroneously be generated for all the other areas. However, it is unlikely, because some of the valley plains (e.g.,  Figure 10b).

Annual Displacement in Echigo Plain
Here, we briefly present the annual displacements over the Echigo Plain, which we estimate from the time series using least-squares. Although the actual periodic displacement is not strictly sinusoidal (Figures 12d and 13b), for simplicity, we use a simple sinusoidal function, along with the linear function, to estimate the amplitude and time offset (e.g., [57,69]).
Large amplitudes are detected in Ojiya City (~22 mm) and Sanjo City (~13 mm), as was previously shown in Sections 4.3 and 4.4 (Figure 14a). Some other areas also exhibit relatively large amplitudes, including Nagaoka City (~17 mm) and Shibata City (~11 mm). Elsewhere, only moderate amplitude signals ranging from 4-7 mm/yr are detected. The spatial pattern revealed by the amplitude map might reflect variations in geological properties. We should note that if the reference point had a periodic displacement, an apparent inverse periodic displacement would erroneously be generated for all the other areas. However, it is unlikely, because some of the valley plains (e.g., southeast of Sanjo City), which should have different geological properties from the Echigo Plain, show very low amplitudes (< 2 mm).
Remote Sens. 2019, 11, x FOR PEER REVIEW 24 of 30 southeast of Sanjo City), which should have different geological properties from the Echigo Plain, show very low amplitudes (< 2 mm). The estimated time offset is almost constant over the Echigo Plain (Figure 14b). The lower peak of the annual displacements occurs around February and March, implying that periodic displacements occur simultaneously, and the common cause is groundwater pumping for snowmelting.

Processing Time and Disk Usage
For most existing InSAR time series analysis software packages (e.g., StaMPS or GIAnT), users generally need to produce the interferograms over the area of interest themselves. If we start by downloading a stack of Sentinel-1 SLC data and then produce a whole set of the interferograms, the process takes a large amount of time and disk space. For example, to produce ~100 images and ~300 interferograms used in the case study, ~1 TB of Sentinel-1 SLC products must be downloaded, and ~600 GB of coregistered SLC data must be created, along with ~300 GB of interferograms and other intermediate files, corresponding to ~2 TB in total. In addition to the computational resource requirements, extensive knowledge of InSAR processing and Sentinel-1 data is also needed to prepare the data for time series analysis.
One of the advantages of LiCSBAS is that users do not need to produce interferograms if the LiCSAR products are available in their area of interest, saving significant amounts of time and disk space. The size of GeoTIFF files for ~300 interferograms is only ~5 GB ( Table 2). The size of the files produced by LiCSBAS for the Japan test frame is only ~30 GB. LiCSBAS processing times are also reasonable, even for an entire frame (~5 hours; Table 2). The estimates presented here are based on The estimated time offset is almost constant over the Echigo Plain (Figure 14b). The lower peak of the annual displacements occurs around February and March, implying that periodic displacements occur simultaneously, and the common cause is groundwater pumping for snow-melting.

Processing Time and Disk Usage
For most existing InSAR time series analysis software packages (e.g., StaMPS or GIAnT), users generally need to produce the interferograms over the area of interest themselves. If we start by downloading a stack of Sentinel-1 SLC data and then produce a whole set of the interferograms, the process takes a large amount of time and disk space. For example, to produce~100 images and 300 interferograms used in the case study,~1 TB of Sentinel-1 SLC products must be downloaded, and~600 GB of coregistered SLC data must be created, along with~300 GB of interferograms and other intermediate files, corresponding to~2 TB in total. In addition to the computational resource requirements, extensive knowledge of InSAR processing and Sentinel-1 data is also needed to prepare the data for time series analysis.
One of the advantages of LiCSBAS is that users do not need to produce interferograms if the LiCSAR products are available in their area of interest, saving significant amounts of time and disk space. The size of GeoTIFF files for~300 interferograms is only~5 GB ( Table 2). The size of the files produced by LiCSBAS for the Japan test frame is only~30 GB. LiCSBAS processing times are also reasonable, even for an entire frame (~5 hours; Table 2). The estimates presented here are based on an Intel Xeon E5-2690 CPU (8 cores, 2.90 GHz),~5 GB maximum of used RAM, and a download speed of~30 MB/sec. Further, if the area of interest is a specific part within the frame, clipping at Step 0-5 greatly saves processing time and disk space but also could improve the result because of reducing the impact of potential unwrapping errors (Section 4). Step 0-1 10 min --Step 0-2 10 min -3 min Step 0-3 30 min -5 min Step 0-4 5 min 2 min 1 min Step 0-5 -5 min - Step

Downsampling at
Step 0-2 makes further processing much faster and file sizes much smaller at the cost of resolution (Table 2), whilst broadly, the same result can be derived ( Figure S12). If the full resolution (~100 m) is not required (e.g., for large-scale strain estimation), downsampling is worth incorporating. Even when full resolution is preferred as the final result, quick processing with the downsampled dataset is useful to check if a good result can be obtained and to adjust parameters/thresholds prior to processing at full or higher resolution.

Limitations
The interferograms provided by LiCSAR are multilooked to~100 m resolution and strongly spatially filtered for the main purpose of the LiCS project (i.e., large-scale deformation monitoring and tectonic strain mapping), as mentioned in Section 2.1. Although km-scale deformations can be accurately detected, as demonstrated in Section 4, the LiCSAR products are not optimized for more localized deformation studies such as infrastructure monitoring.
LiCSBAS deals with gaps in the SB network by imposing a temporal constraint at the inversion step (Step 1-3; Section 2.4.3). We note that the interpolated displacement between the gaps might not be realistic if the true deformation is not linear. It is difficult to quantify the uncertainty of the derived Remote Sens. 2020, 12, 424 25 of 29 time series with gap, and the potential bias should be carefully considered. Gaps are easily identified in the time series viewer to aid in interpretation.
Densely vegetated areas, including large portions of Japan, suffer from incoherence in C-band data, even with the short temporal baselines of Sentinel-1, compared to earlier radar satellites (see Section 3). In this study, the minimum acquisition interval is 12 days versus 6 days for Sentinel-1 across Europe [12]. Quite a few of the 12-day interferograms retain sufficient coherence over the frame, whereas most of the 24-day and longer interferograms lose coherence in the vegetated areas ( Figure S13). In general, the coherence decays exponentially with time [70,71]. More frequent acquisition intervals (i.e., 6 days) increase the likelihood of retaining coherence at C-band and detecting deformation over vegetated areas [72].
Another possible solution for the vegetated areas is exploitation of L-band data. Although abundant L-band data like Sentinel-1 are currently not available, several new planned L-band satellites missions, such as ALOS-4 and NISAR, are capable of frequent, wide-swath observations. LiCSBAS can also be used to process data from other satellites (see Section 2.2).

Conclusions
We have developed LiCSBAS as an open-source InSAR time series analysis package integrated with the automated LiCSAR processing chain outputs. LiCSBAS utilizes published LiCSAR products, and users do not need to produce interferograms from SLC data, thus avoiding high data processing time and disk space consumption. LiCSBAS automatically identifies and removes interferograms with many unwrapping errors by the loop closure test and estimates reliable time series and velocities with the help of masking based upon several noise indices. All processing steps are easy to run using a batch script, and users can adjust processing parameters to obtain improved results. The interactive time series viewer aids in the interpretation of the derived displacements.
Case studies presented here over Japan demonstrate that LiCSBAS can detect both large-scale and localized deformation with an accuracy of < 1 cm at each epoch and~2 mm/yr in velocity. These results are validated by comparisons with GNSS and leveling data. Displacements with different temporal characteristics, such as linear, periodic, and episodic, in the Echigo Plain are also successfully and accurately detected.
LiCSBAS enables users to easily perform InSAR time series analysis without producing interferograms, thus facilitating the use of globally available and abundant Sentinel-1 data for a wide range of scientific investigations.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/3/424/s1: Figure S1. Quick-look image example of the GACOS correction. Figure S2. Example of masking at Step 0-4. Figure S3. Example of clipping at Step 0-5. Figure S4. Example of bad unwrapped interferograms. Figure S5. Displacement time series at all GNSS stations. Figure S6. Displacement time series at 950179 after removing the constant velocity. Figure S7. Perpendicular baseline configuration and network of the SB interferograms without masking of the northwestern area. Figure S8. LOS velocity with the gap in the network. Figure S9. Time series of groundwater level at Well I and II. Figure S10. Time series of GNSS displacement at 950240, with reference to 960566. Figure S11. Cumulative LOS displacements in Ojiya City, with reference to the first acquisition of 25 November, 2014 for all acquisitions. Figure S12. LOS velocity derived from 10 × 10 downsampled dataset. Figure S13. Examples of 12-and 24-days interferograms.