InSAR Time Series Analysis of L-Band Data for Understanding Tropical Peatland Degradation and Restoration

In this study, satellite radar observations are employed to reveal spatiotemporal changes in ground surface height of peatlands that have, and have not, undergone restoration in Central Kalimantan, Indonesia. Our time series analysis of 26 scenes of Advanced Land Observation Satellite-1 (ALOS-1) Phased-Array L-band Synthetic-Aperture Radar (PALSAR) images acquired between 2006 and 2010 suggests that peatland restoration was positively affected by the construction time of dams—the earlier the dam was constructed, the more significant the restoration appears. The results also suggest that the dams resulted in an increase of ground water level, which in turn stopped peat losing height. For peatland areas without restoration, the peatland continuously lost peat height by up to 7.7 cm/yr. InSAR-derived peat height changes allow the investigation of restoration effects over a wide area and can also be used to indirectly assess the relative magnitude and spatial pattern of peatland damage caused by drainage and fires. Such an assessment can provide key information for guiding future restoration activities.


Introduction
As one of the most important natural wetland ecosystems in the world, peatlands cover about 2-3% of the Earth's land surface and store 500 to 700 Gt of carbon [1], nearly equal to the atmospheric carbon pool (750 Gt) [2]. Therefore, peatlands play an important role in the global carbon cycle [3,4]. In the tropics, peatlands cover around 40 million hectares (Mha) and accumulate around 90 Gt of carbon, while more than 54% and 76% of this tropical peatland area and carbon storage occur in South-East Asia [5]. Although tropical peatlands are rich in carbon storage, land use changes, drainage, and recurrent fires have resulted in extensive peatland degradation in this region during the past 20 years [6,7]. From 1990 to 2010, the forest coverage of South-East Asia decreased from 268 million ha to 236 million ha (an average net loss of 1.6 million ha yr −1 , 0.6% yr −1 ) [8], which results in a warmer, drier environment at local scale [9]. Drainage of peatlands causes peat oxidation and releases carbon into the atmosphere as CO 2 , which results in peat subsidence [10]. It is estimated that an additional drainage depth of 10 cm results in 0.9 cm peat subsidence and approximately 9 t CO 2 ha −1 yr −1 emissions [11]. Another severe consequence of drainage is that peat becomes dry and susceptible to fires during dry seasons (typically from June to October). Fires are the most severe in long dry seasons, such as the 1997/98 El Niño Southern Oscillation (ENSO) event. It was reported that 0.12-0.15 Gt of carbon was emitted due to peat combustion from Central Kalimantan, Indonesia during the ENSO dry season of 1997 [12]. The total carbon released to atmosphere due to the above emission from peat combustion in Indonesia in 1997 was between 0.81 and 2.57 Gt, which is equivalent to 13-40% of the mean annual global carbon emission from fossil fuels of 2001 [12]. Therefore, the restoration of peat hydrology and conservation of tropical peatlands is critical in order to prevent peat oxidation, peat fires and mitigate excess CO 2 emission [13][14][15]. For example, from 2004 to 2008, hydrological restoration commenced in Central Kalimantan (of the Mega-Rice project drained landscape) by building large dams [15,16], and understanding the effectiveness of the restoration activity is important for guiding future restoration activities.
Interferometric Synthetic Aperture Radar (InSAR) measures surface displacements in the line of sight to the satellite by differencing the phase measurements between two complex radar images acquired at two epochs but with similar satellite geometries, and the resulting difference phase forms a new kind of image called an interferogram [17]. In the last two decades, InSAR has been proven to be a powerful tool for measuring the Earth's surface movements [18]. It can be accurate to a few millimeters, with high spatial resolution (pixels in the order of sub-meters to tens of meters) over a wide area (e.g., 100 km × 100 km or even wider) [19].
The potential of InSAR to measure peat height changes has been explored [20][21][22][23][24][25][26][27]. For example, in 2011, using four scenes of European Remote Sensing (ERS, C-band) tandem pairs (1 day interval) from October 1997 to January 2000 and the four-pass (four SAR images) InSAR method, a subsidence rate of 2 cm/year was estimated for Central Kalimantan, Indonesia [21]. In another study, using two pairs of Sentinel-1A image, the first pair acquired in June 2015 and June 2016, and second pair in June 2016 and July 2017, subsidence and uplift signals are revealed by two-pass InSAR (D-InSAR) method in peatland of Pelalawan Regency, Riau, Indonesia [27]. The four-pass InSAR method only uses one pair of SAR images to detect the surface deformation; the other pair is used to generate topography information and then applied to remove the topographic phase contained in the first pair. Although C-band tandem pairs promise high coherence in tropic forest, coherence is low for a longer image acquisition interval (such as 35 days) in areas with high dense vegetation [28]. The temporal decorrelation likely makes it difficult to apply C-band SAR images in tropical peatlands area, comparison to L-band SAR images [28]. Because shorter wavelength (C-band) has weaker penetration capability than longer wavelength (L-band). This leads to C-band radar signals are primary reflected back from the tree canopy and only a few signals are returned from the trunk, while the situation is on the contrary for L-band. For both four-pass and two-pass InSAR method, the surface movement rate estimation is likely to be limited as well, because atmospheric effects, orbital ramps, and decorrelation effects are not mitigated by the methods [17]. In addition, no information on the surface movement history can be obtained as well.
To reduce the temporal decorrelation and various errors, InSAR time series algorithms are developed by using multitemporal SAR images since 2001, such as Permanent Scatterer InSAR (PSInSAR TM ) [29], Small BAseline Subset (SBAS) [30], SqueenSAR TM [31]. For more details about the concept and procedure, please refer to Batuhan et al., [18], the paper provides good reviews on InSAR time series method. Different InSAR time series methods and datasets are applied in measuring peatlands surface change. Using Intermittent Small Baseline Subset (ISBAS) DInSAR algorithm, ERS C-band images (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) and Sentinel-1 C-band (2015-2016), Alshammari et al., [25] detected subsidence and uplift change over peatland areas in the northeast of Scotland, UK. Using 161 Sentinel-1 SAR (C-band) images and ISBAS method, Marshall et al., [24] also measured the deformation patterns of tropical peat swamp around Kuala Lumpur International Airport. Fiaschi et al., [32] detected uplift signal of up to 8.9 mm/year and subsidence signal of up to −11.6 mm/year in Irish peatlands by using 124 Sentinel-1A/B images acquired from 4 May 2015 to 1 March 2018 and PS-InSAR method. Although, C-band shows potential in measuring peatlands surface motion, it is likely hard to extract sufficient number of PS points and provide good coverage, particularly over a tropical peatlands area [24]. So far, the application of L-band data and InSAR time series methods for understanding peatlands degradation and restoration in tropical zone is rarely documented in peer-reviewed scientific literature. Therefore, the main goal of this study was to employ ALOS-1 PALSAR L-band SAR images (free of charge) acquired between 2006 and 2010, and InSAR Time Series with Atmospheric Estimation Model (InSAR TS+AEM) method to investigate the spatial and temporal evolution of the peatland surface change in Central Kalimantan, Indonesia. Since latest ALOS-2 PALSAR data are commercial data (~2200 USD for a standard archived Stripmap image), it restricts their use to many end users because of budget constraints. And then assess the effects of hydrological peatland restoration through the use of dams in the drainage systems, which is important knowledge for the planning of future restoration projects.

The Study Area
In 1996, the Indonesian Government initiated the One Million Hectare Mega Rice Project (MRP) in Central Kalimantan (Indonesia). The principal objective of the MRP was food (rice) production in the area between the Sabangau River in the west, the Kahayan and Kapuas Rivers in the south and the Barito River in the east (Figure 1). The area was divided into four blocks (A, B, C, and D) according to the rivers' boundaries (light blue curves) in the MRP project, and the city of Palangkaraya is located in the north of Block C [16]. Large-scale deforestation was initiated in January 1996 to prepare the land for rice paddies, which involved mostly clear-felling the peat swamp forest and evacuating a massive network of in excess of 4400 km of drainage channels across the deep peat domes; some drainage channels were up to 30 m wide and 10 m deep [15,16,33]. Indonesia [34]. The red square in the top right corner is the location of this study area.  [34]. The red square in the top right corner is the location of this study area.
In July 1999, the MRP was closed down because the area had been over-drained by the channels and unable to store water needed for agriculture irrigation. Consequently, the area was largely abandoned and is currently 'unproductive' and susceptible to flooding [35]. The drainage of the peatland caused by the network of extensive channels results in peat oxidation, subsidence, and regular outbreaks of large-scale fires [13,35,36]. Seven types of land cover dominate this area (Figure 1), i.e., shrub swamp, shrub, settlement, secondary swamp forest, rice field, primary swamp forest and clearing [34]. Most of the shrub swamp and rice field areas can be characterized as destroyed, abandoned and/or fire-damaged, and only a small part of the rice field area is still in use [14]. Note that the land cover data shown in Figure 1 [15,16]. The dams in Block A are similar to those in Block C in terms of size and design [15,16,37].

InSAR Interferometric Processing and Time Series Analysis
In this study, two adjacent tracks of L-band (a wavelength of 23.6 cm) ALOS-1 PALSAR images acquired between 2006 and 2010 were processed and analyzed (Table 1), covering a total area of 11,780 km 2 . L-band was used as this has a better penetration capability than shorter wavelength signals (such as C-band), and the ground and the trunks contribute to the signal and to avoid decorrelation caused by changes on the canopy. In this study, the coherence in heavy forest is generally low for interferograms with a perpendicular baseline (satellite separation) greater than 1000 m and/or a temporal baseline (i.e., time interval) longer than 1 year. To minimize the temporal decorrelation effects, only those image pairs with small baselines (i.e., perpendicular baseline < 1000 m and temporal baseline < 1 year) were selected for interferometric processing and further analysis ( Figure 2). All the interferograms were formed from raw radar images using the JPL/Caltech ROI_PAC (Version 3.1b) [38] with multi-look factor of 8 (in range direction) by 16 (in azimuth direction) to reduce the speckle noise, with a final resolution of about 80 m by 80 m. The 3-arcsecond (about 90 m resolution) SRTM digital elevation model (DEM) [39] was used to remove topographic phase contribution. The interferograms were unwrapped using the SNAPHU algorithm [40] to obtain surface displacements in the satellite line of sight (LOS). A phase closure technique [41] was employed to identify major phase unwrapping errors that can then be manually corrected. The basic idea of the phase closure technique is that phase contributions behave in a conservative manner, i.e., ∅ AB − ∅ AC − ∅ BC = 0, where ∅ AB is the phase contribution of interferogram AB constructed from acquisitions A and B. Using the phase closure technique, unwrapping errors were identified by summing around a loop and checking the residuals. To determine whether a dam can effectively block a drainage system, one way is to study data collected during dry seasons when there is less surface water and the surface area would be less affected by flooding. This avoids the study area being flooded, a period when there is no coherence. Therefore, only 12 and 14 PALSAR images acquired during dry seasons between 20 December 2006 and 1 September 2010 were used in this study (Table 1), and 27 and 30 small baseline interferograms were generated from Tracks 421 and 422 respectively ( Figure 2). Finally, these 57 interferograms were geocoded with the SRTM DEM with a pixel spacing of~90 m for further time series analysis.    (Table 1)

176
The unwrapped phase consists of six terms: time invariant constant motion, non-linear surface 177 motion, topographic errors, orbital ramps, atmospheric effects (atmospheric phase screen, APS) and 178 noise. Several time series approaches have been developed to separate deformation signals (i.e., the 179 first two terms) from the other four terms using multiple interferograms [29][30][31]42,43]. In this study, 180 the InSAR TS+AEM [44] was used for time series analysis, which is based on the small baseline subset 181 algorithm (e.g., [30]). In addition to full coherent pixels (i.e., those coherent pixels in every individual 182 interferogram), partially coherent pixels (i.e., pixels with coherence greater than 0.3 in at least 50% 183 independent interferograms) were also selected for time series analysis to increase the coverage of 184 the final mean velocity maps [45,46]. The mean coherence of all the 27 interferograms from Track 185 Figure 2. Small baseline interferograms generated from ALOS PALSAR Tracks 421 (a) and 422 (b) respectively. Each black triangle denotes one SAR image, and black solid lines between two black triangles represent interferograms. The red square marks the reference date for time series analysis.
The unwrapped phase consists of six terms: time invariant constant motion, non-linear surface motion, topographic errors, orbital ramps, atmospheric effects (atmospheric phase screen, APS) and noise. Several time series approaches have been developed to separate deformation signals (i.e., the first two terms) from the other four terms using multiple interferograms [29][30][31]42,43]. In this study, the InSAR TS+AEM [44] was used for time series analysis, which is based on the small baseline subset algorithm (e.g., [30]). In addition to full coherent pixels (i.e., those coherent pixels in every individual interferogram), partially coherent pixels (i.e., pixels with coherence greater than 0.3 in at least 50% independent interferograms) were also selected for time series analysis to increase the coverage of the final mean velocity maps [45,46]. The mean coherence of all the 27 interferograms from Track T421 is presented in Figure 3b. Comparing with the landcover map (Figure 1) visually, we can find that shrub swap area holds high coherence (higher than 0.5), the rice field area has moderate coherence (ranges from 0.3 to 0.5), while the primary and second swamp forest areas present low coherence (lower than 0.3). This suggests that the selected partially coherent points are expected in shrub swap and rice field areas, other than in primary and second swamp forest areas. For a given set of unwrapped interferograms, the topographic errors are correlated with the perpendicular baselines [17] and can be estimated and removed from the unwrapped phase. To remove orbit ramps, a network approach was employed to remove a best-fit plane from every interferogram [17]. The main difference between the remaining three terms is that the nonlinear displacements are correlated both in space and in time, the atmospheric contribution is correlated in space only, while the noise is spatially and temporally decorrelated in the unwrapped interferograms. Taking into account the spatial structure of atmospheric effects (correlated in space only) (e.g., the power-law process [47,48]), APS can be estimated using a temporarily linear velocity (TLV) model [44] and distinguished from the non-linear surface motion. The algorithm is applied iteratively until convergence is achieved [44,49].
To assess the precision of the InSAR time series method used in this study, a correlation analysis between two mean height change rate maps from the two adjacent tracks was performed with data from their overlapping area. A correlation coefficient of 0.67 was obtained with an RMS difference of 0.7 cm/yr (Figure 3c). Based on the assumption that the mean rate uncertainties in both tracks are identical, the RMS of the mean velocity is 0.49 cm/yr, suggesting the InSAR mean rate maps are reliable. Since the difference in incidence angle between the two tracks is small, its influence on the surface deformation projection on the line of sight can be neglected. Thus we combined data from both tracks to perform time series analysis using InSAR TS + AEM, and from this produced the final mean peat height change rate map and time series results. This was done with 26 dry-season SAR images acquired from 20 December 2006 to 1 September 2010, with the mid-timeline image acquired on 9 September 2008 chosen as the reference image, and a stable location at the city of Palangka Raya as the reference site. Within the overlapping area of the two tracks, the number of images included in the time series analysis is almost double, and hence the temporal resolution is increased significantly, providing a better estimate of the mean rate and displacement history. Outside the overlapping area, only one track of SAR images was used, and the InSAR results are almost identical to the previous time series results (i.e., one track included in time series only). The combined mean velocity map (rate of change of surface height) from both tracks is shown in Figure 3a, which will be used in the following analysis and discussions.
In areas with extensive channels and no restoration, four points (cleared (CL), fire scar (FS), rice field (RF) and secondary swamp forest (SF)) distributed within various land use categories were selected for further investigating peatland height change rates and history (Table 2, Figure 3a).

246
The spatial pattern of the surface height changes appears to be associated with land use. No   Figure 3a shows the final mean height change rate (mean velocity) map along line of sight. Note that the negative value means the surface is moving away from the satellite ('subsidence'), while the positive value means the surface is moving towards the satellite ('uplift'). The total number of partially coherent pixels is 859,118 (1,374,913 in total), covering~62% of the study area, and approximately 73 points/km 2 at the case of 90 m resolution. With extensive channels areas in all blocks, subsidence signals can be observed clearly, the largest subsidence rate is more than 5 cm /yr. While the central part of Block C shows uplift signals. Between 2006 and 2010, the mean velocities for Block A, B, C and D are −1.1 ± 1.0 cm/yr, −0.8 ± 1.2 cm/yr, 0.03 ± 1.4 cm/yr and −0.3 ± 1.2 cm/yr respectively. In terms of area, 49% of the study area showed loss in surface height (less than −0.7 cm/yr), 15% retained similar surface height (between −0.7 cm/yr and 0.7 cm/yr), and 37% increased in surface height (greater than 0.7 cm/yr). Figure 3d shows the final root mean square (RMS) of the mean velocity with an average RMS of 0.45 cm/yr. The spatial pattern of the surface height changes appears to be associated with land use. No obvious surface height change can be observed in Palangkaraya city (NW corner of the study area) during the observation period, which is as would be expected. The drained shrub swamp area in Block B has the largest rate of decrease in surface heights. Much of this rice field areas in Blocks A, C and D show decrease patterns in surface heights. The cleared area has relatively small decreases in surface heights ranging from 0 to 2 cm/yr. In the center of Block C, signals in the secondary swamp forest area show increased surface heights during the observation period, with rates of up to~3 cm/yr.

Peatland Height Changes in Area with Restoration
The locations of the dams and their corresponding mean velocities are shown in Figure 4a Table 3. The mean velocities in Block C are shown in Table 3. Note that the closer the mean velocity is to zero, it suggests the bigger impact of the dam construction on reducing the rate of surface height loss along the dams.

Discussions
Although L-band data is used in this study, coherence is still lost in Block C. As show in Figure 3a, the InSAR results have good spatial coverage across the study area, covering various landscapes (Figure 1). By comparing the final valid pixel distribution in Figure 3a, the average coherence map of Track 421 in Figure 3b and land use map in Figure 1 visually, we observe that loss of coherent pixels is often associated with primary or secondary swamp forest, such as at the western part of Block C, the northern part of Block B, and a small area in the northern part of Block A. The loss of coherence is likely because the density of primary/second swamp forest is higher than 100 t/ha [50], where the radar signals are not able to penetrate it canopy of the forest, it returned radar signals are most likely from the canopy of the forest [51,52]. Therefore, temporal decorrelation occurs even over short time periods. While for forest density below 100 t/ha, L-band radar signals can penetrate canopy and are predominantly reflected from the ground with strong double-bouncing contributions from the trunks [53]. Therefore, stable signals are expected from forest density below 100 t/ha in our final results, such as signals from secondary swamp forest in the center of Block C. The eastern part of Block A also shows loss of coherent pixels, but this is because these areas are relatively-low and close to the river, still resulting in frequent summer flooding in this eastern part of Block A even in summer seasons [36].
L-band data would provide much better coverage than that of C-band. Using C-band data and intermittent coherent pixel algorithm in tropical peatlands, Marshall et al., [24] obtained data coverage over 94.35% of the study area by using a coherence threshold of 0.25 and a recurrence of 28% in all interferograms. However, the study also claims that if a coherent pixel is recurrent in each interferogram, only 25.35% of the study area would be covered by coherent pixels. In this study, a coherence threshold of 0.3 and a recurrence of 50% in all interferograms strategy is applied, the final partially coherent pixels covers 62% of the study area. In addition, using C-band data and PS-InSAR method, only 333 stable PS points are found in 6 km 2 boreal peatlands (~56 PS points /km 2 ), where vegetation is dominated by moss [32].Therefore less stable PS points per km 2 are excepted for C-band in tropical peatlands.
The accuracy of the InSAR time series measurements in tropical peatlands is applicable to measure the peatland surface change. Without detailed ground data of this area, cross-validation is achieved by using the overlapping area data from the two independent adjacent tracks. A correlation coefficient of 0.67 was obtained with an RMS of the mean velocity is 0.49 cm/yr implies the InSAR measurement has potential in measuring the tropical peatland height change. The InSAR time series methodology used in this study is also validated by using ground data in other studies, which provide an accuracy of 1.57 cm/yr with a correlation efficient of 0.93 [54]. The accuracy obtained in these studies is also higher than that of ground measurements (e.g., perforated PVC tube inserted vertically in peatlands) [55]. Therefore, InSAR time series and L-band data could be a valuable tool in mapping magnitude and spatial pattern of peatland surface change.
Restoration requires time to take effect [14,15], implying that the earlier the dam is constructed, the more effective its presence will be. This is evidenced by the dams in Blocks A and C (Figure 4b).
In Block A, dams completed in January 2005 appears to have more significant restoration effects compared to dams completed in May 2007 and June 2008 (both p < 0.05, Table 4). Similar effects of successful rewetting were considered observed by investigating the backscatter changes of three PALSAR images acquired in August 2007, August 2008 and August 2009 respectively [15]. Our SAR images were acquired from December 2006 to September 2010, nearly two years after the dam construction in Block C. The mean velocities of C1, C2 and C5 are close to zero, which means the areas near the dams are no longer losing height and thus loss of peat is slower. Increases in surface heights are observed in Points C3 and C4. The most plausible reason is that the ground water level increased in dry seasons and the peatland was rewetted, reducing loss of carbon through oxidation and, or, causing the peatland surface to rise or grow. Since the ALOS-1 PALSAR SAR images are only available after December 2006, we obtain the peat height change rates before the completion of the dam construction from a previous study [15]. By analyzing the backscatter changes of Envisat ASAR images acquired between 2004 and 2009, it was observed that backscatter increased in Block C after the completion of the dam construction in September 2005, implying that the peat soil moisture continuously increased with ground water level, even during the very prolonged dry season of 2006 [15]. Before the completion of the dam construction, a small decrease of the backscatter was observed between July 2004 and September 2005, implying that the peatland was continuously drained with a loss of carbon through oxidation, and resulted in losing peat height.
The swelling and shrinking of peat surface caused by changes in water availability (such as precipitation) had limited effects on the subsidence trend. In Block A (Figure 4b Compared to surface height mapping with GPS or other methods, InSAR provides greater spatial extent in monitoring surface height changes, whilst keeping good precision. However, the InSAR TS + AEM algorithm has its disadvantages. For tropical peatland areas, there is often lack in sufficient number of L-band images within the same path and frame required for InSAR measurement. The lack of sufficient number of SAR images may make it difficult to detect adequate coherent pixels, and separate displacement signals from noises (e.g., the APS noise). For example, there are less than 20 L-band radar images acquired in this study area within the same path and frame by the current ALOS-2 satellite from 2014 to now. By comparison, there are more than 100 Sentinel-1A/B C-band radar images within the same path and frame over the same period. In addition, the use of L-band images could not detect low rate of motion in peatlands, as longer wavelengths tend to be less sensitive to small changes.

394
such security of understanding the community could proceed to using InSAR time series to make 395 temporal and regional assessment of a landscape response to anthropogenic activity and guide future 396 restoration activities. With this capacity InSAR offers a great advantage over field-based point 397 observations by offering regional coverage, and provides a practical way to monitor peatland height  Figure 6a,b, for each box, the central line is the median, the edges of the box are the 25th and 75th percentiles, and whiskers represent the minimum and maximum values. In Figure 6c, the blue dots present the TRMM record for Block C, while the light blue dots present for Black A. Note: (1) LOS surface displacements correspond to the mean LOS surface displacement of all the points within a 270 m by 270 m window, and (2) the reference date is 9 September 2008.

Conclusions
Here we show that InSAR time series analysis of L-band data may be of considerable value in monitoring restoration effects in peatlands. Our analysis suggests that the restoration does have an impact and its effectiveness increases with construction time-the earlier the dam is constructed, the more significant the restoration is. The InSAR time series also provides insight in the likely relative magnitude and spatial pattern of damage to peatland caused by drainage and fires. To confirm the value of using InSAR time series, targeted fieldwork would help corroborate the InSAR understanding and also offer insight into why there are periods and areas of loss of coherence. With such security of understanding the community could proceed to using InSAR time series to make temporal and regional assessment of a landscape response to anthropogenic activity and guide future restoration activities. With this capacity InSAR offers a great advantage over field-based point observations by offering regional coverage, and provides a practical way to monitor peatland height changes continuously and globally, with the archived JERS/ALOS-1 dataset and the operational ALOS-2.