TIME SERIES ANALYSIS OF SUBSIDENCE IN DHAKA CITY, BANGLADESH USING INSAR

Despite existing literature suggesting that Dhaka City, Bangladesh is undergoing subsidence, few researches have been carried out to actually measure the subsidence rate. Previously conducted studies either do not provide sufficiently accurate subsidence results, or the study period is not long enough. In this research, we have tried to address that gap by performing time series subsidence analysis of Dhaka City utilizing Interferometric Synthetic Aperture Radar (InSAR) technique for a study period of 20 years. Synthetic Aperture Radar (SAR) C-band images from ERS, ENVISAT and Sentinel-1A were used to obtain the results. We had to use C-band SAR data from multiple sensors considering data availability issue of the period of investigation (i.e. 1992 to 1999(using ERS); 2003 to 2010(using ENVISAT); 2014 to 2017(using Sentinel 1A)). Most parts of the city is found to be subsiding. Mirpur and Uttara have subsided by over 221mm and 232mm respectively over the 20 years. Ramna and Cantonment subsided around 205mm compared to their level in 1992, whereas both Gulshan and Tejgaon have subsided by about 200mm. Demra and Lalbagh show similar subsidence to the Ramna area, whereas Dhanmondi and Mohammadpur show subsidence rates similar to Tejgaon. We have also assessed the parameter sensitivity to perform this time series subsidence analysis. The parameter selection of coregistration, filtering and unwrapping was found to greatly influence the results. The result is being calibrated with the available GPS observation.

structures to even increased waterlogging as evident in Jakarta, Indonesia and in San Joaquin Valley, California [4][5][6][7]. Therefore, it is deemed crucial to know the pattern of subsidence, in order to anticipate and mitigate those impacts. Moreover, subsidence rates estimation will certainly contribute to a better understanding of the process of subsidence and its contributing factors.
Although it is well established that Dhaka is subsiding, only a handful of research has been conducted to actually measure the subsidence rates. Previously, researches only attempted to infer the rates. Such methodologies mainly utilized in situ techniques such as magnetostratigraphic dating, carbon dating, bore/well-logs and geomorphic surveys [8][9][10]. The authors note that the findings have high uncertainties and limited spatial extent. More recently, a scholar aimed at unraveling the subsidence rates through a time series analysis from 2003 to 2008 with GPS velocity measurements [11]. However, only 1 GPS station observation was available for Dhaka city. The results, therefore, do not provide any idea of the spatial variability of subsidence rates that may be present within the city. Subsidence can also be measured by the earth observation technique of Interferometric Synthetic Aperture Radar (InSAR). InSAR permits remote detection of ground deformation to a high degree of precision [12][13][14]. Moreover, it eliminates the need for expensive and time-consuming field measurements while at the same time enabling large spatial coverage (unlike the rest of the in situ techniques). A study by another scholar used L-band InSAR to investigate subsidence rates in Bangladesh and unlike the previous studies, this was able to show any spatial variability that was present. However, the study period was shortened to only 4 years (2008 to 2011). The average subsidence rate measured by GPS at Dhaka was 12.4 mm/year but the InSAR derived average found it to be only 3.8 mm/year. Consequently, the authors recommended that time series analysis for a larger study period will reduce the anomalies.
Fortunately, the European Remote Sensing (ERS) satellite's and Envisat's (C-band Synthetic Aperture Radar) data was recently made freely available for the period 1992 to 2010 [15]. In addition, Sentinel-1 SAR data is also available for free from 2014 for research purpose. In this study, we have performed the time series subsidence analysis for Dhaka city utilizing these SAR data from different sensors using conventional InSAR technique. Synthetic Aperture Radar (SAR) C-band images from ERS, ENVISAT and Sentinel-1A were used to obtain the results. Although Persistent Scatter Interferometry and Small Baseline have certain advantages over Conventional InSAR, i.e., they are less affected by phase noise but conventional InSAR technique has been chosen in this study due to the lack of a large number of images (with small temporal gaps) [16][17][18].
Specifically, we have addressed two research questions through this study in the context of subsidence rate estimation of Dhaka city using InSAR time series analysis. These are (1) what are the results of time series analysis of subsidence in Dhaka city? (2) Which parameters are found to be highly sensitive for this displacement mapping.

STUDY AREA
It is well established that the city of Dhaka has been subsiding for quite some time, driven by the combined influence of factors like expansive seasonal flooding, tectonics, and compaction of sedimentary layers. Given that the individual contribution of those factors still remains unexplored, it denotes that more reliable, long-term data is missing. InSAR is expected to solve that problem. In order for C-band InSAR to work, it requires ground targets that can produce strong backscatter of the radio waves and for them to remain coherent (i.e. relatively unchanged). The presence of urban clusters within the city, many of which have changed little (between the acquisition times of satellite images), and, the area is mostly devoid of large-scale vegetation-cover, ensure that the technique can be applied within the Dhaka City Corporation area. Furthermore, to calibrate the results, few surface displacement data must be known and the green star (in Figure 1) marks the location of the GPS station which was used for the calibration purpose. These formed the basis of our justification for selecting Dhaka City Corporation (DCC) as our study area.

Materials
The process of InSAR requires SAR (Synthetic Aperture Radar) images. SAR Level1 images from ERS-1 & -2 and ENVISAT-1 & -2 and Sentinel-1A have been used for the analysis. When selecting the images (listed in Table  1) several factors that reduce the signal to noise ratio, i.e., induce decorrelation, were taken into account. Amongst them are atmospheric noises, perpendicular baseline, and seasonality. To reduce noise from atmospheric sources (mainly water vapor), scenes from dry seasons have been preferred [19]. When images from the dry period are not available in the archive only then images from the rest of the year are considered.
Furthermore, in order to ensure that atmospheric noises are not too high, scenes that were imaged when the water vapor content was high has been disregarded. Moderate Resolution Imaging Spectroradiometer (MODIS) Level 2 Water Vapour Product gives a good idea of this phenomenon. This product is derived from the ratio of reflectance value of different bands (ranging from 11 to 12 µm) that are highly sensitive to moisture [20]. However, it is necessary to point out that this derived product only exists for years after 2000.
The other factor, Spatial Baseline, refers to the distance between satellite positions when the images are acquired in repeat-pass orbits. These orbits are quite close together but do not overlap (resulting in slightly different viewing geometries). The perpendicular component of that distance is called the perpendicular baseline. The smaller this component is the lesser is the geometric decorrelation. All perpendicular baseline (for InSAR pairs) have been kept below the critical baseline limit of about 400m [21]. Seasonality means the time gap between image acquisitions, has been attempted to be kept constant. This means (when possible) the images were taken from the same month of each year. When images could not be taken from the same month, it is to be understood that the image was not available in the archive.
The images (listed in Table 1) were processed in SNAP (Sentinel Application Platform by European Space Agency) and unwrapping was done in SNAPHU [22]. Both software is freely available over the internet. SRTM 1arcsecond Digital Elevation Models (DEM) was used for topographic corrections [23].

Methods
The phenomenon of extracting the displacement values from the phase difference information (i.e. an interferogram) is the core concept of Interferometric SAR. Each SAR image consists of an amplitude band and a phase band. InSAR only utilizes the phase band. This band contains the number of complete cycles that the radio waves underwent during its travel time (i.e., it represents a sinusoidal function). This phase information, in turn, signifies how much distance the wave has traveled. When another SAR image's phase band is subtracted from the preceding image (i.e. the older master image is subtracted by the recent slave image) it gives the phase difference, which can then be used to find the ground displacement between the duration of the two image acquisitions (after compensating for Geometric decorrelation, Atmospheric Noise etc.). [21] The difference in phase for a given pixel, in an interferogram, is the sum of-land deformation, topography, atmospheric noise, and satellite geometry-induced errors: ΦTotal Phase Difference = ΦDeformation + ΦGeometry + ΦAtmosphere + ΦTopography [21] The aim is to find Deformation. Therefore, in order to make the phase difference value representative of deformation only, the contribution of all other factors must be compensated. Reduction of Geometric decorrelation is done during the Coregistration process and by choosing images with low baseline values.
Firstly, one image pair is being subset and subsequently coregistered (to ensure that the SAR images overlie the same area). For coarse coregistration 2500 to 3000 GCPs were selected, whereas the fine coregistration window was kept as 8 pixels in order to ensure that only corresponding pixels' phase is subtracted in the interferogram formation stage.
To account for Topographic contributions a DEM has to be subtracted which results in a Differential Interferogram. Furthermore, to account for changes in the surface that is not due to subsidence (like changes in water level, building construction etc.) a coherence map is created. It is used in later steps to mask out the areas of low coherence.
Phase filtering is also performed for the reduction of remaining noise. This step utilizes an adaptive filtering algorithm, namely Goldstein [24]. Such an algorithm is self-learning, meaning it will change the filtering coefficients automatically based on the values of neighboring pixels in order to select the most suitable coefficients for a given pixel [25].
The phase difference value of each pixel in an Interferogram ranges from -pi to +pi (i.e. within -180 to +180 degrees). Phase Unwrapping converts the phase difference to the absolute phase, by adding the appropriate integer number of cycles (i.e. multiples of 360 degrees) [26].
These steps ensure that deformation (due to subsidence) is the only contributing factor remaining in the phase difference of the Interferograms and that they are ready for Phase Unwrapping.
The result of SNAPHU is not an image but tiles of unwrapped phase. Those tiles are stitched together to form a continuous image in SNAP. As the absolute phase is proportional to the travel distance, thus the displacement can be obtained from this unwrapped image (in the Line-ofsight direction).
The subsequent generation of Vertical Displacement maps has been done through the "Phase to Displacement" option in SNAP. These maps illustrate the (spatially varying) deformation values of the terrain between the times of acquisitions of two images used in the Interferogram.
However, there will be the inevitable presence of noise in the interferogram and these will, therefore, be transferred into the displacement maps. Hence, the low coherence areas (i.e. areas with high noise) are masked out using a threshold of the coherence band. The presence of ambiguity in displacement readings necessitate that they are calibrated. Fortunately, there is the presence of GPS data from the "DHAK" station (marked with a green star in Figure 1) is situated in the Ramna area. The long-term average subsidence value adjusted for the line-ofsight component (12.4.004 mm/ year) has been used for calibration ( Figure 1).
Time Series Analysis helped determine the spatially and temporally varying nature of subsidence for the entire study period. When all Calibrated and Geocoded maps had been prepared, the time series analysis was performed. All the maps were stacked in a chronological order. Then "Pixel Info" of SNAP window reveals the displacement values for each pair of images at the location where the cursor is pointed at.
Since the coherence threshold for each image is slightly different, the same locations across all the images are not displayed in the displacement maps. Meaning, the same location in one image may have passed the coherence threshold but in another image, it did not. Not passing the threshold results in that location being masked out and so no subsidence values are found at that location. When this is the case, subsidence value from the nearest coherent location has been taken.

Result Interpretation
Through visual inspection of the subsidence maps, 2 profiles were demarcated; one is from south to north, from Ramna (via Tejgaon and Cantonment) to Uttara. The second profile is from Gulshan towards Mirpur and Pallabi. The points of interest (from the two profiles) were marked ( Figure 7) and their displacement values from each map were collected and tabulated to show the cumulative subsidence (Table 2). In general, increasing subsidence trends along both the profiles were noticed. Mirpur and Uttara had subsidence rates of 10mm/year and 10.1mm/year respectively in 1992. By the end of the year 2017, the rates had increased to 14.5 and 16.6 mm/year respectively. Consequently, these areas were found to have subsided the most, by around 221mm and 232mm in 20 years (i.e. 1992 to 1999(using ERS); 2003 to 2010(using ENVISAT); 2014 to 2017(using Sentinel 1A)).  Finally, the profile ends in Uttara which has subsided 232.6mm during the study period of 20 years. The second profile starts at Gulshan and moves west, towards Mirpur. These areas have subsided 200mm and 221mm respectively. It is obvious that the cumulative subsidence has a rising trend along the profiles. Not only is cumulative subsidence rising along the profile but also with the passage of time. This notion is supported by Figure 3 which shows that the cumulative subsidence in the early years of the study period (1992 to 1996) was more or less similar. The drastic rise was only noticed after 1996. Cumulative values only provide an overall scenario of subsidence along the profiles. However, it is also necessary to study how the marked locations have subsided with time. In order to investigate that, the individual subsidence rates for those locations were noted ( Table 3). The general trend shows a rise in subsidence over time. For areas like Uttara and Mirpur subsidence rates have increased by half their original values during the 20 years study period, and rates for Tejgaon, Cantonment and Gulshan have increased by around a quarter times their rates measured in 1992. Demra and Lalbagh were found to have similar subsidence pattern to the Ramna area, whereas Dhanmondi and Mohammadpur displayed rates similar to Tejgaon and that is why their subsidence rates were not explicitly stated. Since the subsidence rates are changing by a fraction of millimeters, the data will be better visualized through a graph (Figure 4). Upon closer inspection of the graph, it was noticed that the subsidence rates were not constantly rising. Starting from 1992 there is a gradual rise up till 1996. Although the identification of the exact cause of the changes in subsidence is well beyond the scope of this study, some changes may be explained logically. For example, the time period (from 2005 to 2007) that is marked by reduced subsidence is actually derived from a very noisy interferogram (A.4). This resulted in a poor quality displacement map (Figure 7 Bottom Right) which in turn manifested as low subsidence rates in the graph (Figure 4). The change in rates from 2015 may be attributed to the change in sensor characteristics.

Figure 4:
The graph shows the yearly subsidence rates at each marked location.
The results discussed above were obtained by using the value of the parameters mentioned in Section 3.2. Amongst them, Coregistration has proven to be one of the most sensitive steps in this study. Using a greater number of Ground Control Points (GCPs) (i.e. increasing them to 3000, from 2000) proved to be effective. In order to understand why it was effective, it is necessary to first understand how GCPs are used in the Coregistration process. A GCP is a distinguishable feature in an image that is identified automatically by the software and has a known coordinate [27]. Initially, GCPs are identified in the master image by automated algorithms. Using their geographical position information, the algorithm searches for those GCPs in the slave image (within the user-defined window size).
Once GCPs in both images are identified, the information is used to find coefficients for transformation equations. These equations, in turn, are used to establish the relationship between master and the slave image. So the improvement of results with a greater number of GCPs is logical [28]. Moreover, a good number of GCPs in urban areas results in better registration [29]. Upon inspecting the image footprints, it was clear that they did not have a spatial shift of more than 600m in either azimuth or slant range direction. So, a window size of 128 pixels was deemed to be sufficient (ERS and ENVISAT SAR pixels had a resolution of 9.5 x 5m). However, keeping larger window size will not render the results inaccurate, but will take a greater amount of time to produce more or less the same results.
This step is followed by fine coregistration. The process attempts to align the images to sub-pixel accuracy which is not possible through coarse coregistration. This is a crucial step since the results directly determine the accuracy of displacement readings. The algorithm is based on a crosscorrelation technique that creates an alignment between the master and the slave by automatically matching similar pixels through the use of a distributed correlation optimization window [30]. The aim is to improve the coherence between the two images. Multiple sub-pixel shifts of the slave GCPs is computed and the one that produces the highest coherence is selected and applied. Having high coherence means that there is little noise in the images. The noise source can be from satellite viewing geometry (influenced by the size of the perpendicular baseline), atmosphere (water vapor) or even a change in the reflective properties of ground targets. This information is helpful in recognizing that InSAR pairs, in which the same parameters have been applied, fail to produce similar results.
Since the results of coregistration cannot be directly visualized, the coherence maps ( Figure A.1) are intended to illustrate them. The coherence map on the left is for InSAR pair #4 (see the InSAR pair column in Table 1). This pair has one of the lowest perpendicular baselines of all the pairs that have been used. The time difference is restricted to one year meaning there is relatively less change in the land cover and the atmospheric condition during acquisitions was also good. The map shows that using a fine correlation window of 2 pixels, have resulted in good correlation in the majority of the pixels representing the urban areas. However, for pairs with larger baselines, such as for pair #7, larger fine coregistration windows (like 4, 8 pixels) needed to be used. This was necessary since the large baseline adds mathematical constraints to the fine coregistration process. Moreover, the large temporal gaps between the image acquisitions and the significant presence of moisture in the atmosphere also add to the noise in the image. The extent of noise is so much that using 2 by 2-pixel windows cause coregistration to fail completely in these specific cases. Even using larger windows does not produce good coherence ( Figure A.2).
Completion of coregistration is followed by generation of a flattened interferogram. The subtraction of slave phase images from their respective masters is a relatively simple procedure. This is done in conjunction with the topographic phase removal. Since the contribution of Satellite viewing geometry was compensated in the coregistration phase, now the topographic contribution is removed from the interferogram using SRTM DEMs. However, as already mentioned the interferogram still consists of atmospheric noise. Therefore, Goldstein filtering is necessary. This ensures that the fine changes (due to deformation) are preserved but the larger changes (due to noise) gradually becomes less significant. This is illustrated in Figure A.3. However, the presence of too much noise cannot be effectively removed by filtering as is shown in Figure A

Limitations
The availability of SAR images was a major limitation. Although ERS and ENVISAT are supposed to have a revisit time of 35 days, the image archive for Dhaka's location tells a different tale. The lack of images in the dry season and those with small baselines has caused some Interferograms to be very noisy.
Since 3 different Satellite data were used to cover the non-continuous 20 year study period, it cannot be said with certainty whether the change in subsidence shown is representative to reality or simply due to the different sensor characteristics.
Whilst, the results were indeed calibrated, there was no way to validate the results since a second GPS station was not available in or near Dhaka City. Moreover, the GPS station used for calibration only had data from 2003 to 2015.
C band SAR is known to cause decorrelation in vegetated areas. This can be seen quite clearly in Figure A.1 So areas with dense vegetation cover were always masked out and no idea of displacements in those areas could be found. As a matter of fact, we applied the same process to map mininginduced subsidence of Boropukuria Mine, Dinajpur, which is a heavily vegetated area ( Figure A.5). The results were not acceptable owing to the high presence of noise in the images.

CONCLUSION
The overall objective of the study was to perform time series analysis in order to find the subsidence pattern in Dhaka city for the last 20 years. In the calculation of those results, minute change in the parameters of coregistration, filter and unwrapping were found to greatly influence the results. Coregistration proved to be one of the most sensitive steps in the analysis process. For images with small baselines and low atmospheric noise, good subpixel coregistration (using the fine window of 2 pixels) was achieved, which resulted in good coherence ( Figure A.1). Fortunately, all except for 2 image pairs (#7 and 9 from Table 1) exhibit low noise. Using the Goldstein filtering also aided in the removal of most of the residual noise ( Figure A.3). As for, images with high noise levels the problem was overcome by using larger coregistration window (4 or 8 pixels) and a smaller (masking) threshold of coherence. However, the noise content in pair #7 was so high that even after filtering the result was not satisfactory (Figure 7 (bottom right)).
The presence of noise can never be completely eliminated, and this means when the displacement maps are created there is still an extent of ambiguity in the displacement values. Thus, those maps are required to be calibrated with the readings from the GPS station. The results of calibration are shown as Geocoded and Calibrated Displacement Maps in Figure 5 through 10. These maps display the subsidence values in different places of the city during specific portions of the study period. It is also worth to note that we have used C-band SAR imagery from multiple sensors to cover non-continuous 20 year study period which has made the analysis quite challenging. Still, the subsidence result is realistic considering the image acquisition temporality, data quality and sensor characteristics.
The study was possible due to the availability of freely available Copernicus Sentinel 1A, ERS & ENVISAT SAR data; and processing software (SNAP -ESA Sentinel Application Platform v2.0.2, http://step.esa.int) which were provided by the European Space Agency. Special thanks go to the STEP Forum, SAREDU.com and to the patient researchers in this field whose feedback and research work served as the motivation for doing this research.

APPENDIX Figure A.1:
The map on the left shows pixels with good correlation in white. Comparing these white pixels with the Google Earth image on the right, it is clearly visible that the urban areas have good correlation and little noise, whereas the vegetated areas have very high noise (and so less coherence).

Figure A.2:
When the coherence map for pair #4 is compared to Figure A1(left) it becomes clear that the presence of noise sources has taken a toll on the fine coregistration process. There are significantly less number of coherent (white) pixels.  Shows GPS data from Dhaka station. The long-term average subsidence rate was found to be 12.4 mm/year from the trendline of this graph [3]