Natural Hazards and Earth System Sciences Monitoring Subsidence Effects in the Urban Area of Zonguldak Hardcoal Basin of Turkey by Insar-gis Integration

Zonguldak Hardcoal Basin is the largest bituminous coal region in Turkey where extensive underground mining activity exists. Because of this activity subsidence effects have been experienced in different locations of the city. In this study, surface deformations caused by the subsidence have been observed by D-InSAR technique using C-Band RADARSAT data. InSAR data process of 16 RADARSAT images acquired between 24 July 2005–23 October 2006 has resulted in significant deformations in the order of about 6 cm in the most populated region of the city. The deformation map obtained has been integrated with digitized mine production maps and Quickbird Orthoimage into GIS. According to GIS analysis, there are three mine seams at different levels driven below the deformed zone. Many governmental and private buildings located in this area have a high potential risk of subsidence damage. Also, this area covers approximately 12 km of transportation routes.


Introduction
Zonguldak is the capital city of the Zonguldak province, located at the Western Black Sea coast of Turkey, 360 km east of Istanbul and 270 km north of Ankara (Fig. 1).The topography around the city is very undulating and steep; 56% of the city's land is surrounded by mountains, and only 29% with the slope with less than 20% suitable for settlement and agriculture.Also, very heavy forests cover the land surface in the immediate vicinity of the city centre (Zonguldak Province Environment and Forestry Directorate, 2006).
Forestry and fishing are among the important sources of income in the city with a population of 200 thousand while the major industry is the hard coal mining, so that Zongul-Correspondence to: H. S. Kutoglu (kutogluh@hotmail.com)dak is the most famous mining region of Turkey.Underground coal mine extraction in the basin was initiated in 1848 when the Ottoman Empire, before Republic of Turkey, ruled the land.Now, Turkish Hardcoal Enterprises (TTK) is the company authorized for mine production in the basin.According to the official records of TTK, the hard coal production is about 2.5 million tons per year, and has totally reached 400 million tons since 1848.There are widespread coal seams located between the levels of +155 m and −550 m under the city (Turkish Hardcoal Enterprise, 2008).
Due to the above-mentioned mining activities, major subsidence events occurred at some locations of the city in the past (Fig. 2).Some locations are still suffering from subsidence phenomena, and some locations are under the risk of subsidence damage.For the safety of human life and property, monitoring the temporal evolution of the subsidence effects for the city of Zonguldak is very crucial.Also, such a monitoring program can facilitate defining the locations at risk, early warning before major occurrences, and proper urban planning for the future.
Regarding the deformation monitoring, GPS is the most powerful geodetic technique producing the most precise, reliable and exact results to detect pointwise surface deformations.However, to keep wide grounds under control, Differential InSAR is today's most useful geodetic technique.GPS may need thousands of site points to monitor an area of interest which can be controlled only through a pair of InSAR images.
Therefore, D-InSAR was decided as the main technique for this study, since subsidence coming out in the basin is formerly not known.GPS was preferred for validation of the results obtained from InSAR.The surface deformation maps obtained from D-InSAR analysis were integrated with QuickBird Ortho-image into Geographical Information System (GIS).Finally, the deformation areas in different levels were determined, and buildings, roads etc. located in those areas were documented by a GIS analysis.
60 km east of Istanbul and 270 km north of Ankara (Fig. 1).The topography around ery undulating and steep; 56% of the city lands are surrounded by mountains, and ith the slope of less than %20 is suitable for settlement and agriculture.Also, very cover the land surface in the immediate vicinity of the city center [1].

Detecting surface deformations using InSAR technique
InSAR technique has found a wide use since the late 1990s for monitoring ground deformations caused by landslide, tectonic motions, natural and mining-induced subsidence events, and etc. InSAR applications for landslides can be exemplified by Colesanti et al. (2003), Colesanti and Wasowski (2006), Hilley et al. (2004), Kimura and Yamaguchi (2000).This technique has been mostly applied for monitoring tectonic motions (Kimura and Yamaguchi, 2000;Bürgmann et al., 2006;Delouis et al., 2002;Johanson and Bürgmann, 2005;Wright et al., 2004;Bürgmann et al., 2002).It has also been used successfully for surface effects of subsidence induced by natural causes or man made activities (Poland et al., 2006;Amelung et al., 1999;Tomas et al., 2005;Dixon et al., 2006;Fernandez et al., 2009).However, InSAR technique, for deformation monitoring, depends a lot on local conditions such as topography, vegetation, atmosphere as well as resolution of data.Therefore, each study can be regarded unique.For instance, in Turkey, many successful studies were carried out for monitoring tectonic movements, however, no successful study to detect mining-induced deformations has been carried out yet.
locations of the city in the past (Fig. 2).Some locations are still suffering from subsidence phenomena, and some locations are under the risk of subsidence damage.For the safety of human life and property, monitoring the temporal evolution of the subsidence effects for the city of Zonguldak is very crucial.Also, such a monitoring program can facilitate definition of the locations under the risk, early warning before major occurrences, and proper urban planning for the future.As regards deformation monitoring, GPS is the most powerful geodetic technique producing the most precise, reliable and exact results to detect pointwise surface deformations.However, to keep wide grounds under control Differential InSAR is today's most useful geodetic technique.GPS may need thousands of site points to monitor an area of interest which can be controlled through only a pair of InSAR images.
Therefore, D-InSAR was decided as the main technique for this study since where subsidence comes out in the basin is formerly not known.GPS was preferred for validation of the results obtained from InSAR.The surface deformation maps obtained from D-InSAR analysis were Fig. 2. Subsidence scenes from different places of the city.

Data choice
Surface deformations in heavily forested areas can be detected best through L-band SAR data, which can penetrate vegetation (Gens, 2006;European Space Agency, 2009).For that reason, L-band ALOS Palsar data was considered first to detect the surface deformations in the Zonguldak basin.However, sufficient data pairs having a baseline less than 500 m could not be found.Therefore, it was decided to utilize the data from C-band RADARSAT, because orbit stability is quite good.Figure 4 shows the acquisition dates of 16 RADARSAT images used for the study.The image dated 24 July 2005 was chosen as the master image and others were considered slave images (Fig. 3).Therefore, 15 InSAR data pairs were composed to detect surface deformations in relation to the master image.Each pair has the baseline less than 500 m.

InSAR processing for deformation monitoring
InSAR uses phase information of SAR data.The distance change between a sensor and the ground can be measured from phase difference of two observations using phase property in slant-range length: where φ orbit is the orbit fringe caused by baseline distance obtained by two observations while φ topo is the topographic fringe with respect to terrain.These are described by the following equations Surface deformations in heavily forested areas can be detected best through L-band SAR data, which can penetrate vegetation [17,18].For that reason, L-band ALOS Palsar data was considered first to detect the surface deformations in the Zonguldak basin.However, sufficient data pairs having baseline less than 500 m could not be found.Therefore, it was decided to utilize the data from C-band RADARSAT, of which orbit stability is quite good.Figure 4 shows the acquisition dates of 16 RADARSAT images used for the study.

InSAR processing for deformation monitoring
InSAR uses phase information of SAR data.The change of distance between a sensor and the ground can be measured from phase difference of two observations using phase property in slant-range length: where orbit φ is the orbit fringe caused by baseline distance obtained by two observations while topo φ is the topographic fringe with respect to terrain.education campus located in the most densely populated area of the city.Also, the university campus is located in the blue zone.The maximum rate of the deformation is obtained as 5.6 cm per 15 months.A detailed temporal development of the deformation, detected from the each pair, is given in Fig. 6 for the foregoing area.In these equations, B para and B perp represent parallel and perpendicular components of baseline, respectively.Also h, λ, ρ and α represent elevation, wavelength, slant-range length and incidence angle, respectively.Finally, φ atm is a phase delay caused by the reflection of microwaves in the water vapor layer, φ noise is the error component caused by thermal noise, or temporal and spatial decorrelation associated with baseline distance or scattering characteristic change, and φ def represents the amount of surface deformation during the period between two observations (see Massonnet and Feigl, 1998;Franceschetti and Lanari, 1999;Hanssen, 2001).Radar Interferometry -Data Interpretation and Error Analysis, Springer Verlag, New York.The whole process in which the deformation interferograms are obtained is explained by a scheme in Deguchi et al. (2006), illustrated in Fig. 4.
In order to measure the time series of deformation, we applied the smoothness-constrained least-squares method to the pixels satisfying (a) and (b) below: The amount of deformation since 24 July 2005 was solved as an unknown parameter under the condition that U in the Eq. ( 4) was minimized.
where φ  using i and j .The amount of deformation of k-th periodic observation estimated by smoothness-constrained condition is described in t (k) i,j .Finally, α 2 is a parameter to control the smoothness of t (k) i,j and the optimum value is decided by minimizing Akaike's Bayesian Information Criterion (ABIC) given below.
where n is the number of unknown parameters and σ 2 is defined by U/n.Also, I is the identity matrix and D is the matrix of first or second difference of t (k) i,j (Goldstein and Werner, 1998;Akaike and Kitagawa, 1999).

Results
Applying the above-mentioned procedure for the RADARSAT data pairs have resulted in the deformation map given by Fig. 5. Since RADARSAT is a C-band microwave sensor, it is difficult to provide high coherence between master and slave images in densely, vegetated areas.Consequently, fringe with spatial continuity was observed only in the urban district.As seen from the figure, there are significant surface deformations near Kozlu district (Fig. 5).The core (maximum) area of the deformation in red, coincides with the secondary and high school education campus located in the most densely populated area of the city.Also, the university campus is located in the blue zone.The maximum rate of the deformation is obtained as 5.6 cm In this study, a geodetic network with fifteen site points was established in the study area (Fig. 10).This network was observed 4 times through GPS technique, which is the most effective surveying method for monitoring pointwise deformations.Fig. 7 shows the time schedule for the GPS observations.per 15 months.A detailed temporal development of the deformation, detected from the each pair, is given in Fig. 6 for the foregoing area.

GPS validation
D-InSAR can reveal areal deformations over resolution limits of data depending on the factors stated above.InSAR technique can catch an existing deformation, but can not verify non-existent deformation.Moreover, deformation rates can be slightly different point to point in equi-deformation area detected by InSAR.Last, but not least deformations from InSAR are in the direction of slant range.Therefore, the InSAR results must be validated through a more exact method such as GPS.
In this study, a geodetic network with fifteen site points was established in the study area (Fig. 10).This network was observed 4 times through GPS technique, which is the most effective surveying method for monitoring pointwise deformations.Figure 7 shows the time schedule for the GPS observations.
The GPS observations were conducted in static mode with one hour sessions.For datum definition, the network points were connected to the four points of Turkish National Fundamental GPS Network (TUTGA).The observations in The GPS observations were conducted in static mode with one hour sessions.For datum definition, the network points were connected to the four points of Turkish National Fundamental GPS Network (TUTGA).The observations in all the periods were processed into a batch adjustment with respect to the principle of free network adjustment to obtain the temporal coordinates of the site points.In the first process, two TUTGA points were found unstable.Then, these points were removed from the datum definition, and the final results were obtained regarding the other two TUTGA points.Figure 8. Temporal deformations at GPS site points in the significant deformation area.The point 112 is not listed here because it was destroyed during a building construction after the second GPS observation period.
Fig. 8 shows the surface deformations at the site points located in the significant deformation area determined through InSAR technique.The GPS results validate the existence of a vertical deformation in this area.The total deformations are -8.2 cm/ 13 months, -7.3 cm/ 13 months and -4.4 cm/13 months at the points 113, 114 and 115, respectively.In this region, the deformation rate was obtained 5.6 cm/15 months by InSAR.As stated above, the deformations obtained from InSAR are areal while those from GPS are pointwise.For a comparison of both methods, an average deformation rate is obtained 6.6 cm/13 months for the three GPS points.If two months time difference is ignored coherence of the average GPS rate with InSAR result is %85.In this respect, the results from both techniques can be regarded as comparable.

GIS integration and analysis
After InSAR processing, the underground mine production maps for the related area were digitized and incorporated to GIS to establish the relationship between the obtained surface deformations and the mining activity.Also, the final deformation map was integrated to the system, and equivalent deformation zones with 1 cm intervals have been vectorized for defining different levels of deformation (Fig. 9).Then, they were both displayed on QuickBird Orthoimage for the region (Fig. 10).GPS observations were conducted in static mode with one hour sessions.For datum nition, the network points were connected to the four points of Turkish National damental GPS Network (TUTGA).The observations in all the periods were processed into a h adjustment with respect to the principle of free network adjustment to obtain the temporal rdinates of the site points.In the first process, two TUTGA points were found unstable.Then, e points were removed from the datum definition, and the final results were obtained rding the other two TUTGA points.8 shows the surface deformations at the site points located in the significant deformation determined through InSAR technique.The GPS results validate the existence of a vertical rmation in this area.The total deformations are -8.2 cm/ 13 months, -7.3 cm/ 13 months and cm/13 months at the points 113, 114 and 115, respectively.In this region, the deformation was obtained 5.6 cm/15 months by InSAR.As stated above, the deformations obtained from AR are areal while those from GPS are pointwise.For a comparison of both methods, an rage deformation rate is obtained 6.6 cm/13 months for the three GPS points.If two months e difference is ignored coherence of the average GPS rate with InSAR result is %85.In this ect, the results from both techniques can be regarded as comparable.

IS integration and analysis
er InSAR processing, the underground mine production maps for the related area were tized and incorporated to GIS to establish the relationship between the obtained surface rmations and the mining activity.Also, the final deformation map was integrated to the em, and equivalent deformation zones with 1 cm intervals have been vectorized for defining erent levels of deformation (Fig. 9).Then, they were both displayed on QuickBird hoimage for the region (Fig. 10).all the periods were processed into a batch adjustment with respect to the principle of free network adjustment to obtain the temporal coordinates of the site points.In the first process, two TUTGA points were found unstable.Then, these points were removed from the datum definition, and the final results were obtained regarding the other two TUTGA points.
Figure 8 shows the surface deformations at the site points located in the significant deformation area determined through InSAR technique.The GPS results validate the existence of a vertical deformation in this area.The total deformations are −8.2 cm/13 months, −7.3 cm/13 months and −4.4 cm/13 months at the points 113, 114, and 115, respectively.In this region, the deformation rate was obtained 5.6 cm/15 months by InSAR.As stated above, the deformations obtained from InSAR are areal while those from GPS are pointwise.For a comparison of both methods, an average deformation rate is obtained 6.6 cm/13 months for the three GPS points.If two months time difference is ignored coherence of the average GPS rate with InSAR result is 85%.In this respect, the results from both techniques can be regarded as comparable.

GIS integration and analysis
After InSAR processing, the underground mine production maps for the related area were digitized and incorporated to GIS to establish the relationship between the obtained surface deformations and the mining activity.Also, the final deformation map was integrated to the system, and equivalent As seen from Fig. 10, there are mining activities at three different mining levels (-425, -485 a 560 meters) under the deformed zone.This might be the main reason of the subsidence ef monitored on the ground.As seen from Fig. 10, there are mining activities at three different mining levels (-425, -485 and -560 meters) under the deformed zone.This might be the main reason of the subsidence effects monitored on the ground.deformation zones with 1 cm intervals have been vectorized for defining different levels of deformation (Fig. 9).Then, they were both displayed on QuickBird Orthoimage for the region (Fig. 10).
As seen from Fig. 10, there are mining activities at three different mining levels (−425, −485, and −560 m) under the deformed zone.This might be the main reason for the subsidence effects monitored on the ground.
For the risk assessment, superstructure over the region such as buildings, roads etc. was extracted from the Quickbird image, thus, the final product was obtained for further GIS analysis (Fig. 11).All the GIS integration processes applied for the study is illustrated in Fig. 12.Using the final product, statistics of the structures affected by the subsidence were studied by the buffer zone analysis.Thus, buildings and roads located in the different degrees of the deformation zones were determined (Fig. 13).
Based on the analysis, the size of the deformed area detected by InSAR is totally 1.83 km 2 .A total of 51 government and 1040 private buildings (residence) are located in the detected deformation zones.Also total length of the roads in these zones is 11.8 km.A detailed distribution of the structures into the deformation zones can be seen in Table 1.For the risk assessment, superstructure over the region such as buildings, roads etc. was extracted from the Quickbird image; thus the final product was obtained for the further GIS analysis (Fig. 11).All the GIS integration process applied for the study is illustrated in Fig. 12.
Using the final product, statistics of the structures affected by the subsidence were studied by the buffer zone analysis.Thus, buildings and roads located in the different degrees of the deformation zones were determined (Fig. 13).
Based on the analysis, the size of the deformed area detected by InSAR is totally 1.83 km 2 .A total of 51 government and 1040 private buildings (residence) are located in the detected deformation zones.Also total length of the roads in these zones is 11.8 km.A detailed distribution of the structures into the deformation zones can be seen in Table 1.Finally, the buildings were investigated on site; the damaged ones were documented, and the documents were relocated into the GIS (Fig. 14).In Fig 14, a high school damaged due to the subsidence is illustrated from the designed GIS.

Conclusions
As stated previously, Zonguldak lands have very steep topography with heavy forests cover.In this difficult geography, not only a heavy mining activity is carried out underground but also a dense urban life survives on the ground.In this context, Zonguldak city is a rare place in the world.This is the first study that the subsidence effects around the city were successfully observed through the InSAR technique and interpreted by GIS.At the same time, this is also the first successful study in Turkey in which D-InSAR is applied for monitoring the mining-induced surface deformations.Furthermore, an extensive mining activity carried out under a densely populated city makes InSAR-GIS integration necessary.In this regard, the study is one of the first studies in which urban deformations monitored by InSAR results are analyzed within GIS.
For instance, a similar investigation was performed for the city of Naples, Italy [25].Finally, the buildings were investigated on site; the damaged ones were documented, and the documents were relocated into the GIS (Fig. 14).In Fig. 14, a high school damaged due to the subsidence is illustrated by the designed GIS.

Conclusions
As stated previously, Zonguldak lands have very steep topography with heavy forest cover.In this difficult geography, not only is a heavy mining activity carried out underground but also a dense urban life survives on the ground.In this context, Zonguldak city is a rare place in the world.This is the first study that the subsidence effects around the city were successfully observed through the InSAR technique and interpreted by GIS.At the same time, this is also the first successful study in Turkey in which D-InSAR is applied for monitoring the mining-induced surface deformations.Furthermore, an extensive mining activity carried out under a densely populated city makes InSAR-GIS integration necessary.In this regard, the study is one of the first studies in which urban deformations monitored by InSAR results are analysed within GIS.For instance, a similar investigation was performed for the city of Naples, Italy (Lanari, et al., 2004).
In this study, processing of the 15 RADARSAT data pairs resulted in subsidence of about 6 cm/15 months between 24 July 2005-5 September 2006.This is quite a rapid deformation rate.The deformation area monitored surrounds the most populated area in the city, where many buildings and roads are located.The mine production maps of TTK show that there are mining activities at the three different mining levels (−425, −485, and −560 m) under the deformed zone.These activities might be the main reason of the subsidence effects monitored on the ground.In addition, the intense building stock in the region might be another factor triggering the subsidence.According to the GIS analysis, it was determined that 51 government, 1040 private buildings and 11.8 km long roads are in the 1.8 km 2 risky area.Those statistics reveal the dimension of the problem experienced in Zonguldak city.If the deformation rate progresses at the same rate the structures in the zone can be exposed to destructive effects.

Figure 2 .
Figure 2. Subsidence scenes from different places of the city

Figure 3 .
Figure 3. Radarsat data used for the study The image dated 2005/7/24 was chosen as the master image and others were considered slave images (Fig 3).Therefore, 15 InSAR data pairs were composed to detect surface deformations in relation to the master image.Each pair has the baseline less than 500m.
(a) capable of phase unwrapping, (b) coherence values greater than 0.1.

Figure 6 .
Figure 6.Temporal deformation developments from the D-InSAR pairs3.GPS ValidationD-InSAR can reveal areal deformations over resolution limits of data depending on the factors stated above.InSAR technique can catch an existing deformation, but can not verify non-existent deformation.Moreover, deformation rates can be slightly different point to point in equideformation area detected by InSAR.Last, but not least deformations from InSAR are in direction of slant range.Therefore, the InSAR results must be validated through a more exact method such as GPS.

Figure 7 .
Figure 7. Time schedule for GPS campaigns

Figure 7 .
Figure 7. Time schedule for GPS campaigns Temporal deformations at GPS site points in the significant deformation area.The point 112 is not listed here because it was destroyed during a building construction after the second GPS observation period.

Fig. 8 .
Fig.8.Temporal deformations at GPS site points in the significant deformation area.The point 112 is not listed here because it was destroyed during a building construction after the second GPS observation period.

Figure 9 .
Figure 9. Different deformation levels vectorized from InSAR data

Figure 9 .
Figure 9. Different deformation levels vectorized from InSAR data

Figure
Figure 11.Final GIS product

Figure 13 .
Figure 13.Structures located in the deformation zones

Fig. 14 .
Fig. 14.An example of the damaged buildings documented on site (from a high school).

Table 1 . Distribution of the structures into the deformation zones Deformation zone Area (m 2 ) Government building Private building Road length (m)
Figure 14.An example of the damaged buildings documented on site (from a high school)

Table 1 .
Distribution of the structures into the deformation zones.