Minimizing the Residual Topography Effect on Interferograms to Improve DInSAR Results: Estimating Land Subsidence in Port-Said City, Egypt

Abstract

The accurate detection of land subsidence rates in urban areas is important to identify damage-prone areas and provide decision-makers with useful information. Meanwhile, no precise measurements of land subsidence have been undertaken within the coastal Port-Said City in Egypt to evaluate its hazard in relationship to sea-level rise. In order to address this shortcoming, this work introduces and evaluates a methodology that substantially improves small subsidence rate estimations in an urban setting. Eight ALOS/PALSAR-1 scenes were used to estimate the land subsidence rates in Port-Said City, using the Small BAse line Subset (SBAS) DInSAR technique. A stereo pair of ALOS/PRISM was used to generate an accurate DEM to minimize the residual topography effect on the generated interferograms. A total of 347 well distributed ground control points (GCP) were collected in Port-Said City using the leveling instrument to calibrate the generated DEM. Moreover, the eight PALSAR scenes were co-registered using 50 well-distributed GCPs and used to generate 22 interferogram pairs. These PALSAR interferograms were subsequently filtered and used together with the coherence data to calculate the phase unwrapping. The phase-unwrapped interferogram-pairs were then evaluated to discard four interferograms that were affected by phase jumps and phase ramps. Results confirmed that using an accurate DEM (ALOS/PRISM) was essential for accurately detecting small deformations. The vertical displacement rate during the investigated period (2007–2010) was estimated to be −28 mm. The results further indicate that the northern area of Port-Said City has been subjected to higher land subsidence rates compared to the southern area. Such land subsidence rates might induce significant environmental changes with respect to sea-level rise.

1. Introduction

1.1. SAR Interferometry

Differential Interferometry Synthetic Aperture Radar (DInSAR) is a technique which can measure the vertical ground movement of an area in millimeter accuracy [1]. In the last decade a number of Synthetic Aperture Radar (SAR) space-borne satellites have been launched and operated using different wavelengths and polarizations. These repeat-pass interferometric SAR (InSAR) systems have repeatedly proven their efficiency due to their high precision and high spatial resolution [1,2,3]. When the SAR system images the ground, both amplitude (strength) and phase (time) of the backscattered signals are recorded by the receiving antenna. By computing the phase difference from two SAR images acquired at different times, it is possible to generate a radar interferogram, which contains information about the (static) topography and any displacement in the slant range direction that may have occurred between the two SAR image acquisition dates.

Radio interferometry was developed after the Second World War. U.S. military performed the first experiments with airborne radar (SAR) interferometry for topographic mapping in 1971 as reported by Graham in 1974 [4] who used two vertically deployed antennas to record the relative phase difference of the backscattered signals. The application of interferometry derived from satellite data using the repeat-pass method, where the satellite revisits an acquisition area after a certain period, was first demonstrated by Li and Goldstein using historic Seasat data [5].

The main focus of radar interferometry in early applications was on topographic estimation, yielding elevation accuracies comparable with optical methods. A problem in those early applications was that for an effective baseline larger than zero, the topographic signal was always mixed with the deformation signal. A suitable solution to this problem was differential interferometry, where the topographic signal obtained from a so-called reference elevation model was scaled to the baseline conditions of the deformation interferogram and subtracted from it, yielding a differential interferogram [6]. The first demonstration of differential radar interferometry for mapping the displacement field of the Landers earthquake was reported by Massonnet et al. [7]. It follows that an accurate reference digital elevation model is crucial to minimize the effect of the residual topography and improve the generated interferograms, and thus, enhance the final DInSAR results.

Most of the different factors limiting interferometry applications are related to ground surface features and the geometry of radar acquisitions [8]. The interferogram shows the phase difference (ϕ), which represents the difference in distance measured in the radar line of sight (LOS) and includes mainly topography, orbital shifting, surface deformation, and atmospheric effects:

$$\Delta \varnothing =\Delta {\varnothing }_{\mathbf{Topo}}+\Delta {\varnothing }_{\mathbf{Defo}}+\Delta {\varnothing }_{\mathbf{Atm}}+\Delta {\varnothing }_{\mathbf{Noise}}$$

(1)

Several sources of uncertainty arise when calculating the deformation of the surface between the two acquisitions, namely (${\varnothing }_{\mathrm{Deformation}})$, (${\varnothing }_{\mathrm{Topography}}$) which is related to errors in the topographic model, (${\varnothing }_{\mathrm{Orbit}}$) represents the orbital errors, (${\varnothing }_{\mathrm{Noise}}$) represents the total noise introduced due to the instrument and (${\varnothing }_{\mathrm{Atmosphere}}$) is related to errors produced by atmospheric delay mainly caused by the differences in tropospheric water vapor content between scenes at low surface elevations [9].

SAR Interferometry (InSAR) calculates the difference of the phases of two SAR images acquired over the same point from two different positions in space and at different times. Differential interferometry (DInSAR) provides an accurate surface deformation by subtracting the topographic phase from SAR interferograms using an independent DEM, as performed in this study or from another InSAR image depending on the method the interferograms were generated. There are different methodologies in differential interferometry (e.g., two-pass, three-pass, and four-pass interferometry) and with relatively new advanced algorithms (e.g., the Permanent Scatterers (PS) and the Small BAseline Subset (SBAS)) [9,10,11]. These techniques enable the mapping of subtle vertical surface movements (e.g., due to subsidence) down to millimeters [12,13,14,15,16,17].

In this study, the differential synthetic aperture radar interferometry technique (DInSAR), using the Small BAseline Subset (SBAS) algorithm [11] was used to monitor the land subsidence rates in Port-Said City of Egypt, in an attempt to obtain precise detection of land deformation of the study area.

1.2. Land Subsidence in the Nile Delta

The Nile Delta of Egypt is among the most densely populated agricultural areas in the world, with 1360 inhabitants per km² [18]. Its coastline is about 400 km long, from Alexandria in the west to Port-Said and the outlet of the Suez Canal in the east. The low elevation and flat relief of the Nile Delta make it a vulnerable area which is easily subjected to different environmental hazards such as shoreline retreat and sea level rise [19,20,21,22,23,24]. The combination of land subsidence and sea level rise in coastal areas may result in serious problems and has a significant and direct impact on the economy, environment and society [22,25,26].

Hereher [22] stated that around one quarter of the Nile Delta’s area would be flooded if the sea-level raises 1 m. Consequently, over two million people would be enforced to abandon their homes and around 214,000 jobs and over US $35 billion in land value, property, and tourism would be lost [27]. Furthermore, Stanley and Clemente [26] add that the construction of new dams in the upper Nile Basin will only exacerbate the problem of subsidence by reducing Nile water and sediment supply to the lower Nile Basin and Egypt’s delta margin. Consequently, updated maps of the Nile Delta and coastal cities of Egypt provide decision makers with useful information for integrated development and sustainable use of natural resources.

Numerous studies have recorded the vertical land subsidence velocities of the Nile Delta by using in situ observations [23,26,28,29,30,31,32] as well as remote sensing techniques [15,19,33,34,35,36,37,38,39,40]. Across the northern Delta, about 87 continuous sediment drill cores from the National Museum of Natural History (NMNH) of the Smithsonian Institution were used to provide estimates of its subsidence [41]. These drill cores are between 10 and 45 m long and were recovered along transects across the Nile Delta at intervals of 5 to 15 km [41]. The average vertical ground movement rate was recorded as over 4 mm/year across the Nile Delta, while the maximum rates of subsidence were estimated at the east of the Damietta promontory (the eastern down-stream branch of the River Nile). Recently, Stanley and Clemente [26] reported that the maximum subsidence rates have been estimated as 8.4 mm/year in the NE Delta. The highest subsidence velocities were attributed to the major eastern Mediterranean fault system with accelerated velocities attributed to the thick belt of Holocene sediments (Figure 1). This unique geological setting has important implications on subsidence rates in Port-Said City and will be discussed in a later section.

Becker and Sultan [19] used a set of 14 ascending ERS-1 and ERS-2 satellite radar images of the Nile Delta covering the period from 1992 to 1999 to estimate the land subsidence rates. They reported the maximum vertical subsidence velocity as ranging from 6 to 8 mm/year. The study concluded that the maximum rates of subsidence correlate with relatively recent deposition and distribution of young sediments [19]. In another study, a total of 46 descending SAR scenes were used for measuring land subsidence rates on the north portion of the Nile Delta [39]. Scene acquisition dates range from 1992 to 2000 covering a time span of approximately 8 years. In this study, the delta was subdivided into three tracks, i.e., northwestern, north-central, and northeastern tracks.

In the northwestern track, the maximum measured rate of uplift was 12.6 mm/year and maximum measured rate of subsidence was 11.0 mm/year with average subsidence velocity of 0.7 mm/year for the entire section. In the north-central track, the maximum measured rate of uplift was 12.7 mm/year and maximum measured rate of subsidence was 16.8 mm/year. In the northeastern track, the maximum calculated rate of uplift was 4.1 mm/year and maximum calculated rate of subsidence was 10.2 mm/year.

Repeated GPS observations for six stations over the Nile Delta covering the period from 2012 to 2015 were compiled, processed, and analyzed to estimate the land subsidence rate by Hoda et al. [32]. Results of the Hoda et al. [32] study show that land subsidence exists in Port Said and Tanta (Figure 2), while the remaining four stations of the Nile Delta show uplift with vertical land movements ranging from +4.9 to −6.5 mm/year, while the horizontal ground movements vary from 4.2 to 6.6 mm/year in the northeast direction [32].

No precise geodetic measurements have been conducted within Port-Said City to estimate the modern land subsidence rate using the Differential Interferometry SAR method. The only detailed work that has been conducted so far in this area was performed by Gaber et al. [43] where they used the DInSAR technique based on a reference DEM derived from freely available SRTM data, which has limited vertical and horizontal accuracy. The other few existing DInSAR studies in the Delta region do not cover Port-Said City, even though it is considered as one of the most important coastal cities in Egypt. The motivation to extend our previous work [43] was driven by the need to improve the generated interferograms by minimizing the effect of the residual topography component in order to obtain more precise subsidence measurements. Thus, in this work, a field calibrated and more accurate reference DEM was derived from the ALOS/PRISM sensor for Port-Said City with its almost flat topography, resulting in much higher accuracies for land subsidence rates than previously possible by applying the differential interferometry SAR method.

2. Materials and Methods

2.1. Study Area and Geological Setting

Port-Said Governorate lies in the northeastern side of the Nile Delta extending about 30 km along the coast of the Mediterranean Sea and north of the Suez Canal between latitudes 31°00′N and 31°20′N, and longitudes 32°00′E and 32°30′E (Figure 3). The Governorate is bounded by the Mediterranean Sea from the north, Lake Manzala from the west, Lake Malaha from the east, and has a total area of about 1351 km² with an approximate population of 517,216 [44].

The geology of Port-Said is composed of thick layers of sands and pebbles at its base (Pleistocene Mit Ghamr Formation) which are overlain by coastal sands and deposits from the Nile floods (Holocene Bilqas Formation) [45]. The plate tectonic development of the eastern Mediterranean plays an important role in shaping the Nile Delta region. It is bounded by the northern margin of the African Plate, which extends from the subduction zone adjacent to the Cretan and Cyprus arcs to the Red Sea where it drifted apart from the Arabian plate.

The Nile Delta has been affected by small-scale earthquakes along existing minor faults. The Pelusium Line Fault (Figure 1 bottom map) is the main structural feature along the study area, trending NE to SW and bounding the eastern part of Port Said City [29,46,47]. The implications of the Pelusiac Fault system on current subsidence rates in Port-Said City are largely unknown although there might be a possible contribution by this fault system.

2.2. Field Work

This work builds on and significantly refines previous work done by Gaber et al. [43] in several aspects. For instance, a significantly higher number of 347 well distributed Ground Control Points (GCPs) were collected from the field (Figure 4). The increased number of GCPs was necessary to improve the accuracy of calculated ground elevations. Leveling instrument (Sokkia B40) together with the GPS (Leica viva GS15) were used define field points as reference GCPs to convert the height values of the ALOS/PRISM derived DEM for Port-Said from relative to absolute values. GPS (Leica viva GS15) was used for calculating two Bench Marks (BM), one of them in Port-Said at 31°15′31.32′′N and 32°16′35.83′′E and another one in Port-Fouad at 31°14′46.59′′N and 32°18′54.66′′E (Figure 4). Precision leveling (Sokkia B40) was conducted to achieve precision measurements in vertical height and horizontal distance resolutions with respect to the BM (Figure 4).

2.3. DEM Generation from ALOS/PRISM Data

In order to obtain precise detection of land deformation of Port-Said area with millimeter accuracy through the DInSAR technique, a very accurate DEM with proper spatial resolution should be used to minimize the effect of the residual topography on the generated interferograms [48]. Previous DInSAR subsidence studies of the Nile Delta have used widely available DEMs such as SRTM and ASTER as reference DEMs with low spatial resolution and low height accuracy [19,39,43]; consequently, accurate elevation values and deformation results were not possible to achieve. As shown in Equation (1), generating an accurate DEM is very important for accurate estimation of the topographic model (${\varnothing }_{\mathrm{Topography}}$) and thus detection of small deformations.

In this work, a stereo pair of ALOS/PRISM (nadir and backward) images acquired on the 9th of September 2014 over Port-Said City with a wavelength range of 0.52 to 0.77 μm, level 1B2R (geo-reference data), and gain mode 2, was used to generate an accurate DEM with 2.5 m spatial resolution. It should be noted that ALOS/PRISM data are available for the entire world at moderate costs. Extracting a DEM from a stereo pair of ALOS/PRISM data requires image processing software [49], and is a multi-step decision making process involving the setting of many parameters. They are briefly summarized here.

During the epipolar generation process, the nadir viewing image was used as the left stereo pair while the backward viewing image was used as the right stereo pair. To define the relationship between the two stereo images and get an absolute DEM, 50 ground control points (GCPs) and 23 tie points were used with 0.825 maximum y parallax (maximum displacement of pixels between two images in y direction). Subsequently the epipolar geometry was calculated and generated automatically by the ENVI software. It is very important to re-project the stereo pair so that the left and right images have a common orientation. The DEM generation was carried out using epipolar images to minimize the effect of noise in the PRISM data. To validate the quality of the generated DEM from ALOS/PRISM stereo-pair of Port-Said City, two different tests were carried out: (1) checking the quality of the generated DEM with 16 field collected GCPs; and (2) comparing it with the widely available SRTM DEM. It was found that the accuracy of the output ALOS/PRISM DEM was significantly higher than the SRTM one with a RMS error of 2.598 m. This finding is further elaborated in Section 3 (Results). Such obtained accuracy is consistent with [50] results that evaluated the accuracy of the generated DEMs from ALOS/PRISM at Kanagawa Prefecture (Japan), and found that the vertical accuracy ranges from 2 m (plain area) to 5 m (mountainous area) with a horizontal accuracy of around 10 m.

2.4. ALOS/PALSAR Data Processing

The differential synthetic aperture radar interferometry technique (DInSAR), using the Small BAseline Subset (SBAS) algorithm [11] was performed to monitor the land subsidence rates in Port-Said City of Egypt. The same set of eight ALOS/PALSAR images that were used in Gaber et al., 2014 (observation number 608 and center frame number 610), covering the period from 12 November 2007 to 4 April 2010 were used from an ascending orbit with right looking mode and 34.3° incident angle. PALSAR data were obtained in FBS (Fine Beam Single polarization observation mode with HH polarization) with high spatial resolution on the ground surface (about 7 m), and were used to estimate the required displacements information (

Table 1. The ALOS/PALSAR data characteristics (note that TB of 2009 images is very large).

Year	Date Day/Month	Perpendicular Baseline (m)	Critical Baseline (m)	Temporal Baseline (Days)
2007	12 November 28 December	501.919	13,066.59	46
2008	14 November 30 December	631.823	13,059.13	46
2009	28 May 17 November	828.926	13,054.03	184
2010	2 January 4 April	620.006	13,071.28	92

). These eight ALOS/PALSAR images were converted from the Committee on Earth Observing System (CEOS) format to Single Look Complex format (SLC). Thereafter the ALOS/PALSAR data were clipped to the area of interest in Port-Said City to focus on the urban areas only which have high coherence pixels values to reduce the computational time. Figure 5 summarizes the processing workflow applied in this study using the SBAS module.

The correlation between the two complex SAR images decreases systematically with increasing baseline length until it completely disappears [51]. The critical baseline is the baseline length for which the two SAR images become completely decorrelated. Thus, calculating perpendicular, temporal and critical baseline were necessary to examine the validity of DInSAR and can be expressed mathematically for flat surfaces using the following equation [52,53]:

$${\mathbf{B}}_{\perp \mathrm{crit}}=\frac{\mathsf{\lambda }{\mathbf{r}}_0\tan \left({\mathsf{\theta }}_{\mathbf{i}}\right)}{\mathbf{m}{\mathsf{\delta }}_{\mathbf{r}}}$$

(2)

where the carrier wavelength is $\mathsf{\lambda }$, ${\mathrm{r}}_0$ is the slant range, the incident angle ${\mathsf{\theta }}_{\mathrm{i}}$ and the effective critical baseline ${\mathrm{B}}_{\perp \mathrm{crit}}$ are all defined in a plane perpendicular to the flight paths. Therefore, ${\mathsf{\delta }}_{\mathrm{r}}$ is the slant-range resolution and $\mathrm{m}$ is the factor equal to two for a repeat-pass, where each radar image is acquired with its own illuminator. Based on the previous equation, the SLC PALSAR image, which had been acquired on 14 November 2008, was chosen automatically as the super master image by generating a connection graph between the eight ALOS PALSAR images using the SARscape software. Figure 6 shows the relative positions of the used slave SAR images with respect to the super master position, and was used to understand and evaluate the connections between the slave and master images. Every acquisition should be connected with others and if a group of acquisitions cannot be connected to the main network, all the acquisitions belonging to this group are discarded. For this reason, it is suggested to collect images at regular time intervals with small baselines (Figure 6). All slant ranges of the eight PALSAR images were co-registered with super master geometry, which was used as reference image.

During the interferometric generation stage, it is important to define the slant range multi-look. This is typically calculated in order to obtain an almost square pixel of the same sensor resolution. It is necessary to work with a multi-look image greater than 1:1 to increase the signal to noise ratio (SNR) of the interferograms and provide a robust coherence estimation and speed up the processing of large areas.

Each master and slave SAR images acquired with a certain temporal separation were combined to generate 22 interferograms. Ideally, a zero-baseline configuration would result in interferograms whose phase information would only be related to the LOS displacement in the scene. The orbit error contribution and atmospheric parameter $\Delta {\varnothing }_{\mathrm{Atm}}$ were corrected during SAR analysis, and the phase noise was reduced by applying an adaptive filter. The orbital correction was carried out after processing the SAR data into displacement time series using the small baseline approach of [11]. In addition, the atmospheric correction was performed by estimating the statistical properties of SAR signals. An initial set of 22 wrapped interferogram-pairs was generated during the interferometric generation stage, and wrapped interferogram-pairs were subsequently filtered and used together with the coherence data to calculate the phase unwrapping. The wrapped interferograms were analyzed in order to discover unwanted behaviors or data problems (Figure 7) such as inaccurate orbits, un-coherent pairs, atmospheric artifacts, phase jumps, phase ramps, etc.

The wrapped interferograms with strong residual phase ramps originated from the orbit inaccuracies, together with large atmospheric artifacts were corrected by removing the residual phase frequency. The very large temporal or perpendicular baseline between the two acquisitions, especially the SAR data that was acquired in 2009 (Table 1) resulted in generating a total of four wrapped interferograms with very low coherence and these pairs were discarded. Small atmospheric artifacts and phase jumps were removed during the refinement and the re-flattening steps. The phase unwrapping was performed and the 18 unwrapped interferogram pairs were evaluated to discard the interferograms that were affected by large phase jumps and ramps. Subsequently, the 18 well-unwrapped interferograms were chosen for further processing for estimating subsidence rate. Unwrapped interferogram phases were refined and re-flattened by the residual phase method to estimate and remove the remaining phase constants and phase ramps, in order to relate the changes in slant range to the deformation only (due to subsidence).

Therefore, about 50 well distributed GCPs with an RMS horizontal error of 0.029 m were selected to refine and re-flatten the unwrapped interferogram images in order to remove the remaining phase and phase ramps (Figure 8). The ALOS/PRISM generated DEM was used as a reference DEM because of its high accuracy. The locations of the selected GCPs were carefully chosen on the geocoded unwrapped phase interferogram so that they would not be located in the residual topography fringes, phase jumps or displacement fringes. Only four reference points from the selected 50 GCPs were used in this study and their displacement are assumed to be 0. The locations of these reference points were selected from the unwrapped phase file based on the following criteria: (a) no phase jump or ramps; (b) no residual topography fringes; and (c) no displacement fringes if known. Thus, the final DInSAR output measurements are relative to these selected reference points [54].

The phase information of such refined unwrapped interferograms was converted into height to estimate the residual topography. As phase information characterizes the displacement of the reflecting surface in the radial (line of sight) direction and land subsidence is most likely in the vertical direction, the LOS displacement was converted into vertical displacement using the following equation [55]:

$${\mathbf{\Delta }}_{\mathit{S}}=\frac{\mathbf{\Delta }\mathit{R}}{\mathit{c}\mathit{o}\mathit{s}{\mathit{\theta }}_{\left(\mathit{i}\mathit{n}\mathit{c}\right)}}$$

(3)

where ${\Delta }_S$ is the surface displacement in the vertical direction, $\Delta R$ is the LOS displacement and ${\theta }_{\left(inc\right)}$ is the incidence angle. Linear mode was used for first estimation of displacement rate, because there is no prior information about the displacement behavior of the study area. The resulting displacement pattern using the first inversion step is affected by the atmospheric artifacts. However, by applying the second estimation inversion, these artifacts were minimized [54].

The next step provides a more refined estimation of the velocity. At this stage, the displacement location and its spatial coverage can be already retrieved. The estimated residual topography may be affected by ramps, atmospheric artifacts or other distortion, thus the wavelet functionality can be used for correct residual topography estimation [54]. Since in this study the SBAS processing chain was used [54], the wrong residual topography estimation will compromise the final SBAS geocoding, causing horizontal shifts in the results. Thus the wavelet number of levels, which refers to the power of two bases, determines what is kept of the estimated residual topography [56]. Berardino et al. [11] suggested setting this value as a function of the reference DEM (which was used for the interferogram flattening) resolution.

Thus, the wavelet number used in this study was set to zero level because of the high fidelity of the reference DEM derived from ALOS/PRISM data. Due to the DEM’s high spatial resolution, its sensitivity is large enough to estimate the average residual topography to almost 0 m with a standard deviation equal to 9.5 m (Figure 9a). Moreover, the accuracy of the SRTM DEM data was checked in order to compare its results with the field calibrated ALOS/PRISM DEM. The 0 wavelet numbers was also used with the SRTM data and showed average residual topography equal to −52.5 m and high standard deviation equal to 35.1 m (Figure 9b). Meanwhile the two wavelet numbers were also used with the SRTM data and showed average residual topography equal to −1.5 m and high standard deviation equal to 25.03 m (Figure 9c). It is obviously clear that even after using large wavelet number with the freely available SRTM data, the PRISM data with 0 wavelet number showed an improved estimation for the average residual topography.

The second estimation of the displacement using the Small Baseline Subset (SBAS) technique [11] was performed to recover the final processed displacement time series from the atmospheric artifacts, taking into account a coherence threshold applied to the input stack of the unwrapped phase (Figure 10). Finally, all the obtained results were geocoded to adopt two constraints: the velocity precision threshold and the height precision threshold. To estimate meaningful thresholds, in terms of coverage and precision, the statistic tool was used (Figure 11). Consequently, the velocity precision threshold value of seven and the height precision threshold value of three were used (Figure 11).

The deformation component can be isolated from the non-deformation component by addressing phase noise due to changing scattering properties over time. This is achieved by using the phase behavior of the radar signals to select pixels with minimal decorrelation [57]. Pixels that retain a high degree of correlation (high coherence) were used to show their subsidence. The first step output deformation results were calculated for all pixels, which show high and low coherence in terms of signal phase behavior. Finally, the deformation of the pixels that have equal to or more than (≥) 0.3 and 0.4 coherence values were calculated and plotted on the WorldView-2 optical satellite images to accurately estimate the vertical displacement rates and avoid the man-made changes.

3. Results

In recent years, a number of researchers have been interested in measuring land subsidence rates in the Nile Delta region by using radar interferometric techniques [15,19,31,38,39,43]. However, most of those studies were based on inaccurate topographic measurements introducing uncertainties in the resulting subsidence rates. Therefore, in this study we focus on modern land subsidence rates in Port-Said City as measured using the SBAS DinSAR technique. DInSAR technique is becoming a very reliable quantitative geodetic tool for monitoring deformation, compared to other more simple quantitative tools (leveling, total station and GPS). Furthermore, SBAS is a radar interferometric technique that relies on an appropriate combination of differential interferograms created by SAR image pairs characterized by a small orbital separation (baseline), which reduces the spatial decorrelation phenomena [11].

In this study, the land subsidence was calculated only within the urban areas of Port-Said City where coherent pixels could be located to ensure that the obtained land subsidence results would be more accurate. Moreover, generating an accurate DEM was considered crucial to detect small deformations accurately.

Figure 12 illustrates the comparison between the extracted DEM from ALOS/PRISM data with respect to the calculated points in the field using the total station and the SRTM DEM. The figure shows around 10 m height differences between the SRTM and PRISM DEM data. This is because the generated DEM using ALOS/PRISM has around 2.5 m vertical height accuracy, while the SRTM has around 16 m absolute vertical height accuracy and around 10 m relative vertical height accuracy [58].

The higher spatial resolution of the ALOS/PRISM DEM provides more surface variation details than the SRTM DEM. Figure 13 shows the differences between the DEMs, i.e. the one derived from SRTM (Figure 13a) and the one generated from ALOS/PRISM (Figure 13b). The estimated average residual topography and the standard deviation are clearly minimized using the field calibrated ALOS/PRISM data. Consequently, the corresponding estimated velocity was significantly improved by the adopted new method.

In this work, the 2009 SAR data pairs were not found appropriate for DInSAR applications due to their high temporal decorrelation and changes in the urban areas. However, the other DInSAR pairs of 2007, 2008 and 2010 show well-defined fringes in filtered differential interferograms.

It can be clearly observed in Figure 14 that most of the subsidence rates are located in the northern part of Port-Said City, where dense buildings and skyscrapers are located, but some are in the southwestern part of the study area. The subsidence was checked in the field by observing many cracks in the buildings and their lower ground level position compared with the street level. Subsequently, the Line of Sight (LOS) deformation and subsidence rate for each of the fringes were calculated, and finally the weighted average of subsidence rate was calculated for continually subsiding areas, i.e., areas where overlapping fringes were seen in multiple interferograms. The subsidence rates were extracted from the different interferograms, which represent different years, and combined together to obtain the total subsidence rate affecting the study area.

Figure 14 shows that the vertical displacement mainly takes place in the northern part of the urban area in Port Said City as well as in Port-Fuad. During the period from 2007 to 2008, the maximum measured uplift rate was 6 mm/year with a maximum rate of 4 mm/year for 99% of all uplifted areas. The maximum measured vertical subsidence rate was −11 mm/year with a maximum rate of −9 mm/year for 99% of all subsidence areas.

From 2007 to 2009, the maximum measured uplift rate was 13 mm/year with a maximum rate of 10 mm/year for 99% of all uplifted areas. The total measured vertical subsidence rate from 2007 to 2009 was −21.0 mm with a maximum rate of −19 mm for 99% of all subsidence areas (Figure 14). The cumulative measurements of uplift rate from 2007 to 2010 were 18 mm with a maximum rate of 16 mm for 99% of all uplifted areas. The total measured vertical subsidence rate from 2007 to 2010 was −28.0 mm with a maximum rate of −25 mm for 99% of all subsidence areas (Figure 14). These measured rates of subsidence show the cumulative trend of the subsidence during the investigated period as well as the subsidence rate is not constant over the different years.

Figure 15 and Figure 16 show the estimated deformations at the high coherence pixels with values ≥0.3 and ≥0.4 and plotted on the WorldView-2 image. The obtained results were checked in the field by locating, photographing and recording the coordinates of some deformed buildings in Port-Said City. Many buildings show their ground level at lower elevation than the street level and have cracks, which suggest that they have suffered from the land subsidence. Buildings affected by subsidence and checked in the fields were plotted on the final land subsidence maps that were derived by using SBAS DInSAR technique. The location of the affected buildings show high consistence with the SBAS DInSAR estimated results (Figure 17).

4. Discussion

Land subsidence in Port-Said is a growing concern because it will induce marked environmental changes particularly with respect to coastline retreat and sea-level rise in the Mediterranean Sea. Port-Said, one of the most important ports in Egypt, has witnessed significant urban development over the past century. Natural gas extraction in the northern part of the Delta and within the urban footprint of Port-Said has dramatically increased in the last decade. Gas pumping, natural sediment compaction and seismic activities as well as the surface infrastructure (buildings density and height) are factors that might cause or accelerate the observed land subsidence.

Furthermore, records from before the construction of Aswan High Dam show that the average sediment load arriving to the Mediterranean coast was 134 Million tons [59,60]. By the time the construction of the Aswan High Dam was finished in 1964, only 2% of the sediment load in the Nile River makes it through the dam, leaving 98% of the sediment entrapped behind the dam in Lake Nasser [28]. Thus, the construction of the Dam and barrages not only control the water flow, but also the sediment load delivered to the Nile Delta. In addition, future dams such as the construction of Ethiopia’s Grand Renaissance Dam (GERD) are expected to diminish Nile flow and sediment supply to Egypt’s delta region even further [26].

In this study the estimated land subsidence rate for Port-Said City was found to be consistent with previous research that have used a variety of techniques [19,23,26,28,29,30,32,39]. The results of this work show that the northern part of Port-Said (coastal area) is suffering the highest subsidence rates compared to other parts of the city. This might be due to the existence of numerous towers and high-rise buildings and road traffic as well as intensive urban expansions and offshore natural gas pumping not far from Port-Said’s coast. In addition, this part of the city is underlain by relatively thick compactable fine-grains sediments compared to the southern part. Port-Said City is located at the northeastern part of the Nile Delta and surrounded by two lagoons indicating in geological terms a clearly depositional environment.

The proximity of Port-Said City to El-Manzala and El-Malaha lagoons implies the existence of shallow groundwater table conditions. The increased pressure applied to the ground by high-rise buildings may alter groundwater flow patterns and plays a very important role during the compaction of the upper sediments. Altered groundwater flow is restricted within compacted sediment areas where the land subsidence occurs. Regions without extensive development have not yet experienced groundwater flow alteration, thus resulting in observed land swelling. Thus, the uplifted areas identified in this work have been most probably caused by exceeding aquifer poroelasticity under the compaction forces of the urban buildings.

Sediment compaction and swelling are the most common subsidence and uplifting mechanisms, respectively, and were both identified in this study. It follows that the compaction of the relatively young Holocene sediments is expected to be the dominant factor affecting subsidence rates. The identified pattern of the land subsidence rates has a strong spatial correlation with the distribution thickness of the Holocene sediments (Figure 1). These Holocene sediments reach a thickness of up to 47 m in the Port-Said area [41]. Stanley [41] stated that large scale tilting of the delta toward the northeast resulted in the deposition of thicker accumulations of Holocene sediments in the northern part of Port-Said City than any other parts of the Nile Delta. Moreover, active tectonic and seismicity events have also been suggested as a mechanism driving increased subsidence rates in Port-Said [26].

The combined effects of groundwater depletion, natural gas pumping and urban expansion may cause pronounced subsidence effects in large populated coastal cities such as Port-Said, especially if insufficient sediment supplies counteract the effect and hidden normal faults enhance the land surface drop. One such hidden fault is the Pelusiac Fault system (Figure 1 bottom) discussed in Section 2.1 (Study Area and Geological Setting). All aforementioned conditions support the importance of land subsidence estimations in Port-Said City to provide decision-makers with critical information for ensuring sustainable development and land-use.

5. Conclusions

Port-Said City is one of the most important coastal cities in Egypt. Located at northeastern side of the Nile Delta, it was chosen as a representative coastal city to estimate the modern land subsidence rate. Until now, no precise geodetic measurements have been conducted within Port-Said City to capture the present-day subsidence using the differential interferometry SAR method. In order to accurately estimate small subsidence rates, a very accurate DEM has to be used, thus a stereo pair of ALOS/PRISM data was used to generate a DEM with 2.5 m spatial resolution and ±2.5 m high vertical accuracy. A total of 347 well distributed Ground Control Points (GCPs) were collected in the field using leveling instrument and differential GPS in order to calibrate the height values of the generated DEM. Furthermore, the PRISM DEM was compared to the widely used SRTM DEM and the results showed a 10 m height differences with the PRISM DEM, which is matching the field GCP measurements. Such accurate DEM was used as a reference elevation model to minimize the effects of the residual topography and to improve the generated interferograms and consequently achieve more accurate DInSAR results.

The total vertical subsidence rate during the investigated period from 2007 to 2010 was estimated to be −28 mm. The obtained results show that the northern part and coastal areas of Port-Said City exhibit the highest land subsidence rates compared to the southern part. This research found that subsidence estimations in Port-Said City as based on DInSAR data analysis after minimizing the effect of the residual topography, resulted in significantly improved rate estimations. Furthermore, such improved subsidence estimations correlated very well with known distribution of subsurface geology and structures as well as surface infrastructure (buildings density and height) providing valuable information on sediment compaction and swelling patterns to decision makers for sustainable development and land-use planning in Port-Said Governorate.