Land Subsidence Potential Detection in Yogyakarta International Airport using Sentinel-1 Insar Data

On January 27, 2017, the Indonesian Government started building a new international airport in Yogyakarta Province, named Yogyakarta International Airport (YIA) to replace Adisucipto International Airport. YIA is located near the beach, which means that an awareness of natural disasters, such as coastal flooding, is essential. One of the causes of sea water flooding is land subsidence phenomenon. This land subsidence phenomenon can be monitored by using Sentinel-1 Interferometric Synthetic Aperture Radar (InSAR) data. To monitor the crustal deformation, the data used in this research are from years 2016-2019. The data were processed through LiCSBAS software which is published by the COMET in the UK. In the processing scheme, interferograms with many unwrapping errors are detected and removed via loop closure. Reliable time series and velocities are extracted using several noise indices, with the help of masking. The results show the subsidence phenomenon in the YIA area (up to 25 mm).


Introduction
Airports are critical component of the aviation network. They play an important role, not only within the macro transport environment, but also in enhancing the quality of life and the regional economies, which are directly involved in the production of wealth [1]. Evans [2] stated that airports offer all of the requisite facilities for moving passengers and freight from surface to air transportation systems and for airlines to be able to take off and land their aircraft. Airports are much more than just places for flights, corporate transit events or tax free shopping; they are one of the biggest investments for cities and regions. In the heart of the regions in which they are located, and in certain sectors of activity (e.g. tourism), airport infrastructure is of the upmost importance. Due to new market demands and trends arising in related sectors, such as the aviation industry, their positioning has been changed in recent years.
Yogyakarta is an economic, tourist, and cultural city that has great potential. The city attracts tourists, both local and international. The present international airport, Adisucipto Airport, is too small and unable to serve large aircraft. Thus, the Government has reacted to this issue by building a new, more feasible airport: Yogyakarta International Airport (YIA). The new airport is in the sub-district of Temon, which is a large region in Kulon Progo, near the beach. It is located in the west of the Central Java frontier, Temon Kulon Progo. District Kulon Progo is one of the most vulnerable districts on the Indian Ocean in terms of tsunami disaster impacts [3]. This risk arises because the area between the Indo-Australian Plates ( Figure 1) and the very active Eurasian plate lies on the coast of Kulon Progo Regency. A tectonic earthquake source that may cause tsunamis can occur in the subduction area of the South Java Sea [4][5][6].
Owing to several factors including groundwater depletion, underground mines, rapid land growth, and natural hazards, land subsidence has been a global issue over the past decades. This can happen gradually and becomes evident whether geological or human activities cause the land subsidence [7]. Land subsidence occurs in more than 150 countries worldwide, particularly in densely populated metropolitan areas, including major cities in China [8], Japan [9], Italy [10], Mexico [11], and Indonesia [12].
Satellite observation is a crucial source of Earth data that enables to monitor the effects of natural disasters, in particular. Data collection in inaccessible locations, the broad coverage that allows a full analysis of world phenomena and the availability of long-standing historical data for wide areas allows all of the phenomena to be temporarily studied. This has many advantages over other monitoring techniques which provides more cost-effective and easier data collection. Interferometric Synthetic Aperture Radar (InSAR) has been developed as an efficient geodetic survey technique and is used in the geophysical investigation of earthquakes, landslides and other geological processes [13]. The InSAR technology can also monitor the deformation of the ground.
Spatial-temporal decorrelation and delay in air typically affect InSAR measurement accuracy [14]. It has been suggested that, persistent scatter interferometry (PS-InSAR) and Small Baseline subset InSAR (SBAS-InSAR) can be used to improve accuracy. In addition to the spatial and temporal decorrelation, the SBAS-InSAR system eliminates the phase decontamination and errors in the atmosphere. It could, therefore, have reliable characteristics for a time series deformation [15].
The European Space Agency (ERS)-1/2 satellite radar data were used by Manunta et al. [16] to detect largescale deformations in Rome, Italy and the SBAS method demonstrated its ability to obtain similar information from aerial, low-resolution data, leading to the identification of various major displacement sites. From 1992 to 2010, Zeni et al. [17] used SBAS technology with full resolution for a broad range of ground deformity measurements and observed potential deformation in historical buildings in Rome. Their findings demonstrate the efficient detection and monitoring of historic and artistic monuments by a two-scale SBAS-InSAR multi-sensor system. Xu et al. [18] used an ENVISAT Advanced Synthetic Aperture radar (ASAR) to review the coastal vertical land recycling movements in Shenzhen, to acquire them by ascending and descending orbit methods. They noticed a major coastal shrinkage (25 mm/year) in the recovered property.
Scientists at Soekarno-Hatta International Airport, Singapore Changi Airport, and Kuala Lumpur International Airport have previously carried out a study focussing on the deformation monitoring of airports using InSAR data [19][20]. Ng et al. [19] explained that the area of settlement at Soekarno-Hatta International Airport (and surrounding industrial area, which is located between Kalideres and Cengkareng) forms a subsidence bowl. Field subsidence was also found by Radhi [20], at the International Airport of Kuala Lumpur, and Singapore Changi Airport.
Based on the above, in this research we would like to monitor the crustal deformation observed at YIA in the field with Sentinel-1 InSAR analysis for the period of 2016-2019. From the InSAR data in the YIA area, the output is the time series and velocity field of the deformation.

Data and Methodology
The Yogyakarta International Airport in Temon (Kulonprogo district, province of Yogyakarta, Indonesia) is the area of study for this research. The data used here are InSAR 2016-2019 Sentinel-1 data, which have been downloaded from the COMET-LiCS web portal (https://comet.nerc.ac.uk/comet-lics-portal/). The frame ID used is 076D 09725 121107, which in descending direction and covers the area around Central Java Province. Figure 2 shows the field of study. The authors used 110.12° E, 110.55° E, -8.17°   S and -7.72 ° S as the boundary co-ordinates of the study area. The total number of InSAR images in this study are 121 and the research workflow is shown in Figure 3.
In Figure 3, the dashed lines denote optional steps to incorporate correction of the atmosphere (generic atmospheric online correction service, GACOS), masking and clipping. Before the small baseline inversion (SB), tests were carried out for incorrectly unwrapped interferograms to produce velocities and time series. The final products were further masked and filtered.
In order to process the Sentinel-1 InSAR data, the authors used the LiCSBAS software, published by Morishita et al. [21]. The LiCSBAS workflow of InSAR time series processor was divided into two main sections: preparation of a stack of unwrapped data and time series analysis. LiCSBAS began downloading LiCSAR products for the region of interest and this was followed by a data format conversion. To improve precision and efficiency of the processing, tropospheric noise correction using the InSAR (GACOS) external Generic Atmospheric Correction Online [22] data and unwrapped data masking/clipping were optional measures. Due to the accuracy and coverage of unwrapped data and control of the loop closure, incorrect unwrapped data was detected and disposed of in the time series study. GACOS uses the ITD model 25, which generates high resolution zenith total delay charts, to be used for correcting InSAR calculation and other applications, in order to derive layer and turbulence signals from tropospheric total delays. Due to its precision, and coverage of unwicked data, and by testing the closure of the loop, erroneously unwrapped data was identified and discarded during the time series analysis. With regard to SB interferograms, unwrapped phases are generally reduced in the STD for each interferogram (from an average of 6.7 to 4.2 rad from 6.0 to a median rad 3.9) which shows that the GACOS correction considerably reduced the tropospheric noise [21].
The refined stack of unwrapped data was inverted for time series and speed displacement, followed by calculation of the STD and noisy pixels, based on multiple noise indices. Finally, a spatiotemporal filter was applied to the time series in order to eliminate the residual noise and extract the filtered time series and speed.
In LiCSAR processing chains, Interferograms are automatically created on the predefined LiCSAR frame base (usually 13 explosions at each of the IWS sub-swaths corresponding to a 250 km x 250 km area). A second subsidiary image (close to the latest image for coherence preservation), which was already recorded (if available) with an enhanced method for spectral diversity, was used to copy newly acquired data into one primary image [23]. For each acquisition with three prior acquisitions and three subsequent acquisitions, interferograms would then be used, by default, although the amount may be increased in future. Interferograms are spatially filtered with an alpha-value of 1.0 adaptive spectrum filter from the GAMMA software [21], to eliminate noise. They are multilooked with an interface factor of 20 x 4 in range x azimuth (46 x 56 m spacing).
SBAS-InSAR is commonly used for the identification of subsidence, on the basis of the following equations: (1) Equation (1) is a quantity range of M differential interferograms, produced in the same area by a N+1 SAR image at a specified time ( 0 , 1 , … , ). Furthermore, Equation (2) represents the interferometric composition of the j (to be generated with two images in and ) interferogram phase in pixels (x, r), where x and r are the azimuth and range coordinates, respectively. The change in distance from the target to the radar along the sight line (LOS) causes the phase ∆ . Moreover, the phases of ∆ , ∆ , ∆ , and ∆ are created through the ground and satellite orbit, the atmospheric effect (and other noise), and, in particular, the effect of tropospheric delay. In order to achieve the deformation process ∆ function, SBAS-InSAR removes the residual components from the interferometric phase. where G is a M x (N-1) zero architecture matrix, representing the interferogram network relationship with incremental displacements, given that the unwrapped interferogram (i.e. displacement between two acquisitions) is the sum of the corresponding incremental displacements [24]. Cumulative displacements (i.e. time series for displacement) are determined by summing up the incremental displacement for each acquisition. The mean displacement velocity is then calculated on at least the quadrature of the cumulative displacements.
The NSBAS approach was adopted [25], which imposes a temporary limitation to obtain the more practical time series of the displacement even with a disconnected network.

(4)
If γ is the scaling (weighting) element in the temporal constraints, a linear displacement (d = vt+c) is assumed. The solution within the linked components of the network is minimally influenced by the low time limit (e.g. 0.0001). The time restriction component, therefore, only affects the connection via network gaps. Equation (4) can be employed for pixels with fully connected networks and those with gaps.

Sentinel-1 InSAR Unwrapped Phase
The downsample will be the one-precision floating point formats of the GeoTIFF file process with no header, in order to further analyze the time series and downsamples (multilooks) data, when set by -n option. Figure 4 provides examples of unwrapped phases. The entire Sentinel-1 InSAR unwrapped process is explained in Figure 4a, while the unwrapped process of the clip using the YIA area is shown in Figure 4b.

Loop Closure and Mask Time Series
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. Fortunately, there are no interferograms in this case (which have major unwrapping errors) and they are now removed. This phase generates the mask with the aid of several noise indices, obtained in previous steps for the time series and speed of displacement. If the value for a pixel of any noise indicator is greater or smaller than the threshold specified, the pixel is masked. Figure 5 shows the time series of the mask. In this figure, the average number of coh_avg are constant, n_unw is the number of unwrapped data, vstd is standard deviation in velocity (mm/yr), maxTlen is the maximum time of linked network (year), n gap is the number of gaps in the network, and stc is spatiotemporary coherence (mm). The masked/unmasked velocities could be seen in the top row and the other pictures being noise indices. The number shown in the brackets next to the noise index titles is the limit used.

Time Series and Velocity Field Analysis
The derived time-series also contains a number of noise-related conditions, such as residual tropospheric noise, ionospheric noise, and orbital errors. An isolation of these components from the displacement time series is possible with a space-time   [26]. In order to complete this filtering, the authors added a single, two dimensionals Gaussian kernel in time and space.
As an interactive time-series display, LiCSBAS has two windows (graphical user interfaces). The first image window shows dimensions, cumulative deformation velocity and noise indexes ( Figure 6). From Figure 6, it can be seen that there was uplift deformation in the northern part. However, the southern part (near to the beach) is dominated by land subsidence. The uplift deformation is affected by the Menoreh Mountains, while the land subsidence near to the beach is affected by the Sunda Megathrust activity, which subducted under the Eurasian Plate.
The corresponding time series, with and without a space time filter, is drawn in Figure 7 immediately when a pixel of interest is clicked (black dot which is shown in Figure 6). Figure 7 shows the time series of black dot (see Figure  6) which means land subsidence is ~25 mm/yr. Figure  8 clarifies the distribution of the pattern of deformation at Yogyakarta International Airport. This figure provides a reference map from google maps for understanding the result of deformation, according to the specific location. Deformation range is typically between -48 mm/yr and 36 mm/yr. The standard deviation of the deformations is 0.1 mm/yr. YIA's main deformation processes comprise ground subsidence. This is especially dangerous for the airport, as this area has a high risk of natural disasters. YIA has a high risk of tsunami disaster, based on work carried out by Fauzi et al. [27].

Conclusions
The aim of this study was to monitor the deformations in the Yogyakarta International Airport (YIA) region. Uplift deformation occured in the northern part of YIA. However, the southern part (near to the beach) showed a dominance of land subsidence. The uplift deformation is affected by the Menoreh Mountains, while the land subsidence near to the beach is affected by the Sunda Megathrust activity, which subducts under the Eurasian Plate. The time series shown in Figure 7, resulted in land subsidence of ~25 mm/yr. The time series displays land subsidence. The deformation range observed in YIA is -48 mm/yr to 36 mm/yr, with a normal difference of 0.1 mm/yr. The deformation in the vicinity of the YIA area is important to understand, as this area is located near to the Indian Ocean and the Sunda Megathrust.