Remote Sensing Mapping Land Subsidence Related to Underground Coal Fires in the Wuda Coalfield (northern China) Using a Small Stack of Alos Palsar Differential Interferograms

Coal fires have been found to be a serious problem worldwide in coal mining reserves. Coal fires burn valuable coal reserves and lead to severe environmental degradation of the region. Moreover, coal fires can result in massive surface displacements due to the reduction in volume of the burning coal and can cause thermal effects in the adjacent rock mass particularly cracks and fissures. The Wuda coalfield in Northern China is known for being an exclusive storehouse of prime coking coal as well as for being the site of occurrence of the maximum number of known coal fires among all the coalfields in China and worldwide, and is chosen as our study area. In this study, we have investigated the capabilities and limitations of ALOS PALSAR data for monitoring the land subsidence that accompanies coal fires by means of satellite differential interferometric synthetic aperture radar (DInSAR) observations. An approach to map the large and highly non-linear subsidence based on a small number of SAR images was applied to the Wuda coalfield to reveal the spatial and temporal signals of land subsidence in areas affected by coal fires. The DInSAR results agree well with coal fire data obtained from field investigations and thermal anomaly information, which demonstrates that the capability of ALOS PALSAR data and the proposed approach have remarkable potential to detect this land subsidence of interest. In 1153 addition, our results also provide a spatial extent and temporal evolution of the land subsidence behavior accompanying the coal fires, which indicated that several coal fire zones suffer accelerated ongoing land subsidence, whilst other coal fire zones are newly subsiding areas arising from coal fires in the period of development.


Introduction
Underground coal spontaneous combustion means that underground coal layers are burned by natural or artificial factors.Developing with time, the coal fires spread and grow, and coalfield fire disasters (hereafter called coal fires) gradually emerge.Underground coal spontaneous combustion is a kind of natural disaster, which is a very serious problem in many countries in the world, such as the United States, Australia, India and Indonesia.Very serious coal fires have occurred in coalfields in northern China, and have been ongoing for a long period.The coal fires result not only in the loss of millions of tons of coal resources annually, but also damage the ecological environment, cause air pollution, and influence quality of life in the coal fire areas [1].
Satellite remote sensing methods are quick and effective for the earth's surface and are extensively used to survey the vicinity surface coal fires.The methods mostly used optical mutispectral data to extract surface features related to coal fires and employed thermal infrared data to analyze coal fires according to the direct relationship between coal fires and thermal surface anomalies [1][2][3][4][5][6][7][8].However, few studies focus on the land subsidence accompanying uncontrolled coal fires.Space-borne differential interferometric synthetic aperture radar (DInSAR) has proven to be an effective technique for land subsidence measurement due to its precision, spatial coverage and resolution.The capability of DInSAR for land deformation mapping has been demonstrated in many applications, such as glacier movement [9], volcanic activity [10], crustal movement [11,12], slope instability [13][14][15], groundwater overexploitation [16,17] and underground mining activity [18][19][20].However, few research efforts using the DInSAR techniques have been carried out for coalfields to investigate the land subsidence related to underground coal fires.[21][22][23] tried to apply DInSAR techniques to detect land subsidence related to subsurface coal fires in ERS interferograms, covering the Ke-er Jian coalfield, the Wuda coalfield and the Ruqigou coalfield in Northern China, respectively.[24] tried to use a time series of 34 ENVISAT ASAR datasets over the Wuda coalfield to reveal the potential of detecting the land subsidence accompanying coal fires by means of satellite DInSAR observations.Ground-based geodetic measurements (e.g., spirit leveling or GPS surveys) have not been applied widely to monitor land subsidence caused by coal fires because most of them around the world are located in the places that are hardly able to access, such as Wuda coalfield in Northern China.Compared with ground-based geodetic measurements, the DInSAR technique can determine a two-dimensional (2D) deformation field on a large scale and at a reasonable cost.However, DInSAR observations [21][22][23][24] in coalfield areas were limited, due to the interferometric decorrelation of C-band radar induced by the rapidly changing topography of the coalfield.L-band radar is more effective for observation in such coalfield areas, as it maintains much higher coherence in cases of high deformation gradient due to its longer radar wavelength.The Phased Array type L-band Synthetic Aperture Radar (PALSAR) instrument on the Japanese Space Agency's Advanced Land Observing Satellite (ALOS) provided global coverage with longer wavelength L-band data [25].In [26], PALSAR data has been used to study the subsidence phenomena in the Wuda coalfield using 3-pass DInSAR observations.However, since the deformation detecting ability of 3-pass DInSAR is strongly influenced by the effect of inaccurate compensation of the topographic phase and atmospheric, and hence greatly reducing the precision of ground deformation detected.In addition, 3-pass DInSAR can encounter difficulties in mapping slow-evolving, long lasting deformations and cannot obtain deformations time series.To cope with this, advanced InSAR algorithms, namely interferometric Permanent Scatterers InSAR (PSI) [27,28], has been proposed.This PSI implementation exploits the temporal and spatial characteristics of interferometric signals collected from point-likely targets or so-called Persistent Scatterers (PSs), and allows a precise estimation of the different phase contributions in terms of deformation, atmospheric, topographic and orbital phase.However, due to its insufficient number of archived PALSAR data available and significant non-linear nature of the land subsidence induced by coal fires, it is often impossible to apply PSI techniques using the L-band SAR satellites in coalfields.
The main objective of this study is to evaluate the capabilities and limitations of ALOS PALSAR data for monitoring the land subsidence accompanying coal fires.This study also describes an approach to map the large and highly non-linear subsidence induced by underground coal fires based on a small number of SAR images compared with the PSI technique.We selected the Wuda coalfield as our study area.A time series of seven ALOS PALSAR datasets over the study area, acquired between October 2006 and January 2008, were used to retrieve average Line-of-Sight (LOS) deformation rates and the deformation history of the time series.The DInSAR derived results validated by field survey data of coal fires and thermal anomaly information demonstrate the capability of ALOS PALSAR data and the proposed approach in our study to reveal the spatial and temporal behavior of land subsidence in the Wuda coalfields.
This paper is organized as follows.The study area and data acquisition are presented in Section 2. Section 3 describes the DInSAR method to estimate the LOS deformation rates and the deformation history of the time series.The results obtained from the DInSAR method over the Wuda coalfield using ALOS PALSAR data is then presented in Section 4. The results are then compared with field survey data of coal fires and thermal anomaly information.In Section 5, an analysis of the capabilities and limitations of the ALOS PALSAR data used for monitoring the land subsidence accompanying the coal fires is firstly provided.Subsequently, a comparison between the approach in our study and the PSI technique is discussed.Section 6 presents the conclusions and discusses further development work to be undertaken.

Study Area Location
The study area is located in the 5 km west of Wuda city, which is in the south-central part of the Inner Mongolia Autonomous Region, North China.The area extends from 106°34′E to 106°38′E and from 39°27′N to 39°27′N (see Figure 1).The area of nearly north-south direction, which is 13 km wide from east to west and 16 km long from north to south, has a spatial extent of approximately 208 km 2 with elevation varying from 1,150 m to 1,300 m above sea level and covers the whole Wuda coalfield.The Badain Jaran Desert and the Helan Shan Mountains surround the coalfield in the northern, southern and eastern direction, while the Yellow River in the east supplies enough water for Wuda city.The Wuda coalfield area consists of three underground coal-mining fields: Wuhushan operating in the south, Huangbaici in the east and Suhaitu in the north-west [24].The study area is a low hills landform with well-developed valleys.The high annual average evaporation and low annual average precipitation results in very dry climatic conditions.As a result, the land surface in Wuda coalfield is covered mainly by bare rocks and soils and has sparse dry desert shrubs, i.e., almost no vegetation [29], which means that the coal seam is in a dry state for a long time favoring spontaneous combustion events.

Land Subsidence Caused by Coal Fires
Land subsidence in the Wuda coalfield is caused by the volume loss arising from underground coal-seam burning, which leads to the change in surface morphology.The land subsidence can cause significant surface cracks and fissures [22].Furthermore, cracks and fissures that develop within the overburdened rocks provide ventilation paths for the air; these exhaust pathways for the byproducts of coal combustion and hence accelerate the burning process [30].According to a previous survey, there is a particular distribution of different amounts and types of surface deformation in underground coal fire combustion areas [1].From the perspective of time evolution of the coal fire area, the underground coal fire combustion experiences the four evolutionary processes of combustion (1) beginning, (2) combustion center formation, (3) combustion system development and (4) combustion extinguished.Once the underground coal fire is formed, it will spread continuously along the underground coal seam, and can last for several years or even hundreds of years.During the underground coal fire combustion center formation and combustion system development stages, since coal-seam loss is small, land subsidence is slow leading to land subsidence per year values in the order of centimeters, as illustrated in Figure 2(a).In the underground coal fire combustion extinguished stage, the residue of extensive ash-layers after a coal-seam has been burnt out results in a volume loss underground and thus often leads to land subsidence and cap rock collapses, especially in the Wuda area.This results in sudden surficial cracks and land subsidence per year values in the order of ten centimeters, as presented in Figure 2(b).From the perspective of the spatial distribution of the entire coal fire area, sudden surficial cracks and slow trough-like land subsidence exist at the same time.Land subsidence areas extend from several to hundreds of square meters, with spatial distributions that are generally roughly close to the location of the underground coal fire.In addition, underground coal mining is another deformation factor in Wuda coal field, which results in the order of ten centimeters of surface changes in a few months, as presented in Figure 2(c).

Data Acquisition
For this study, we used a total of 7 L-band ALOS PALSAR images acquired from 11 October 2006 to listed in Table 1.All data were single-look complex images (CEOS format and ALOS-specific L1.1 products) and in HH polarization with the 34.3° look angle.

Method and Data Processing
SAR interferometry utilizes the phase information from SAR images to compute the geometry between the satellite position and the target on the ground.The differential phase in interferogram m, in the pixel of range and azimuth coordinates (x,y), from the SAR acquisitions at times t A and t B after the removal of the topographic phase, is given by [31,32]: where  W  is the wrapping operator, which indicates that the differential interferometric phase is modulo-2π, and λ is the wavelength of the radar signal., is the Doppler centroid phase term, which is related to the azimuth sub-pixel position with respect to the center of the resolution cell, az.It depends on the difference in Doppler centroid frequency between image acquisitions at times t A and t B and the instantaneous velocity of the satellite at the pixel in the earth-fixed coordinate system., m xy  denotes the azimuth sub-pixel position coefficient [33].In this study, ALOS PALSAR data are focused on a common Doppler centroid, therefore the azimuth sub-pixel position component was not considered here.The forth term, Considering its insufficient number of archived PALSAR data available in this study and significant non-linear nature of the land subsidence induced by coal fires, it is often impossible to apply PSI techniques using the L-band SAR satellites in Wuda coalfield.Therefore, this paper describes an approach for monitoring the land subsidence related to underground coal fires using a small stack of DInSAR interferograms.The proposed technique is based on utilizing multiple SAR images, with the shortest revisiting times in order to mitigate the phase saturation problem and also to minimise temporal decorrelation.The main function of the proposed technique is to estimate the DEM error, especially for relatively long baselines pairs, using 1-D regression analysis.The deformation from each differential interferogram can be derived using standard DInSAR processing after the removal of the DEM error phase.Subsidence time series maps and an average subsidence velocity map can then be calculated.The overall processing procedure is described by the block diagram in Figure 3.

Interferogram Formation
To construct interferograms, and because of the coherence loss and high subsidence rates, we limit image pairs to those having temporal baselines lower than 200 days and perpendicular baselines lower than 2,000 m.We build a total of 11 interferograms from 7 images that provide links between all 7 acquisition dates (see Figure 4).
Because the perpendicular baselines of some interferograms constructed by ALOS PALSAR data are relatively large, which results in changes of the scattering energy at the same position on the ground, it is very easy to generate a bias using a general image registration method.Therefore, in this study, we adopt a multi-scale image registration method based on wavelet multi-resolution decomposition [34] to ensure that each ground target contributes to the same (range, azimuth) pixel in both the master and the slave image.Multilooking of a factor 3 was applied to the interferograms along the azimuth, leading to a ground pixel resolution of 10 m × 10 m.During the processing to generate the interferograms, bandpass filtering is used to eliminate the non-overlapping portions of the Doppler signal.After the interferograms are generated, the topographic contribution is removed using the 30-m resolution Aster digital elevation model (GDEM).These DEM data were also applied for geocoding the resultant InSAR products from range-Doppler coordinates into map geometry corresponding to the Universal Transverse Mercator (UTM) coordinate system.An orbit refinement was performed that attempted to remove orbital ramps from each differential interferogram by means of calculating the interferogram fringe rate, which can correct linear and quadratic orbital ramps successfully, especially in the case of large perpendicular baselines.Nonlinear adaptive spatial filtering [35] was applied to each interferogram to increase the signal to noise ratio.Once the filtered differential interferograms are generated, the next step is to exclude the pixels that are strongly affected by noise and decorrelation.Coherence stability [36] was used in this study to exclude the pixels strongly affected by noise that carry out no significant phase information and may have a negative impact on the performance of the unwrapping procedure, which only involves those pixels that exhibit an estimated coherence value higher than an assumed threshold.The next step of the procedure involves the retrieval of the original (unwrapped) phase signals from the modulo-2π (wrapped) phases directly computed from the interferograms.This operation, referred to as phase unwrapping, is based on the Minimum Cost Flow algorithm [37], which is integrated with a region growing procedure to improve the performances in areas with low signal to noise ratio.

DEM Error Correction
In Wuda coalfield, due to the dumping of massive overburden waste materials near the mining areas, which is a fairly common practice when coal mining, the topography of the coalfield changes significantly [24].It is common to find that spatial unwrapping is unsuccessful for pairs with relatively long normal baselines due to inaccurate compensation of the topographic phase.We therefore need to correct the DEM error based on pairs with relatively short normal baselines that have been spatial unwrapped successfully.Once the interferograms with short normal baselines are unwrapped successfully, 1-D regression on this unwrapped phase is used to determine DEM corrections.The (unwrapped) phase value of each selected pixel is relative to the same reference pixel, which can be seen as a star network.The atmospheric signal, which is not modeled, is expected to increase for points further away from the reference point [38].Phase difference observations with long arcs contain more atmospheric signal, and hence it is more difficult to separate the parameters (e.g., DEM error) from the atmospheric signal.The phase differences have to be considered between nearby points and therefore a triangular irregular network is constructed based on the Delaunay triangulation network.Using this network the DEM error can be estimated from the double phase differences between two points for each arc in the interferogram stack via one-dimensional regression analysis of the following system of equations: where M is the number of input unwrapped interferograms, (x 1 ,y 1 ) and (x 2 ,y 2 ) are the pixel coordinates (range and azimuth) at each end of the arc, and is the relative DEM error between the two points.
In order to examine the reliability of the arc (pairs of points), the baseline dependences of the spatial difference in the unwrapped interferograms are analyzed for the arc (pairs of points).For a stack of M records, this means that M spatial differences between differential interferometric phases are considered.The phase model indicates a linear dependence of the topographic phase on the perpendicular baseline component, with the slope of the regression indicating the relative height correction.This height correction is the height that needs to be added to the second point, so that its phase becomes consistent with that of the reference point.The standard deviation of the phase from the regression is used as a quality measure, allowing for the detection and rejection of points that are not suited for regression analysis.The DEM error is estimated from a sparse grid of selected points.It is then interpolated onto the original grid of the DEM using Kriging interpolation.Based on these, a DEM correction can be estimated.Using the improved DEM, it is then usually also possible to unwrap the interferometric phases of pairs with longer normal baselines.The corrected DEM phases for each interferogram can then be calculated.Following this operation, the phase contribution from the corrected DEM is subtracted from each input original differential interferogram.The results of DEM error corrections are shown on Figure 5(a,b) before and after applying the topographic correction.The residual topographic noise is clearly visible on the first image and significantly reduced on the second.Accordingly, the residual phases are believed to be due to atmospheric disturbances.

Correction of the Atmospheric Phase Delay
It is possible to classify the atmospheric contributions into turbulent mixing and vertical stratification contributions [31].The former is the result of turbulent processes in the atmosphere, the latter results from different vertical refractivity profiles during the two SAR acquisitions.Because of their nonlinear nature, turbulent processes affect the interferometric signal on a wide range of scales and are difficult to model.Currently, no methods are able to quantify them accurately and routinely [39], so most studies consider it as a random signal or noise affecting interferograms.The vertical stratification contribution is correlated with elevation [40][41][42], as the delay in radar microwave propagation from the satellite to the ground depends on the integrated atmospheric water vapor content, which is dependent upon the scene elevation.
From one perspective, since the Wuda coalfield is located in an arid and semiarid area of North China, the range of the atmospheric fluctuations is relatively small at the different acquisition times of each SAR image, with the result that turbulent contributions can not affect the interferometric signal seriously.Therefore, turbulent contributions are identified as a random signal that is associated with one particular image acquisition and cancels out in summed "stacks" of interferograms, because turbulent contributions are generally uncorrelated between satellite image acquisitions.On the other hand, the stratified atmospheric artifacts should be related to the range of height difference.Thereby, for the Wuda coalfield, the elevation of which only varies from 1,150 to 1,300 m above sea level, this component should be negligible.To confirm this idea, we performed a linear regression between the interferometric phases and the elevation [43] to estimate tropostatic phase delay from all interferograms.For each interferogram, this method estimates the parameters of a linear model of the residual phase with respect to the altitude values in the original DEM.The output vertical stratification phase has the height dependent unwrapped phase as determined from the regression: The model consists of a phase constant a 0 in radians and phase slope a 1 in units of radians per meter.Estimating, for each interferogram, we found that the statistics of the stratified displacement estimations only less than 0.003 m.Therefore, the stratified atmospheric artifacts were not considered in this study.

Time Series Analysis
After the removal of DEM error from the differential interferograms, the standard DInSAR technique can be applied.The phase noise can be reduced by applying an adaptive filter [35].The phases in the refined interferograms are unwrapped using the Minimum Cost Flow (MCF) method, which is applied independently [37].Once the refined interferograms are unwrapped, the cumulative deformation in the LOS direction can be obtained by solving the linear system: is the accumulated deformation along the LOS for the pixel of range and azimuth coordinates (x,y) at each acquisition time of the SAR image.To solving this linear system, we used singular value decomposition (SVD) to obtain the least squares solution for the deformation time series.A similar approach was proposed in [36].SAR displacement measurements are not 3D; the displacement vector along the LOS is a composite of vertical, easting and northing displacement components.At least two independent LOS measurements are needed in order to solve for the vertical and horizontal displacement vector.However, there is only a single satellite direction available in this study.Sudden surficial cracks and slow trough-like land subsidence in the Wuda coalfield are always perpendicular to the spreading direction of the underground coal fire [4].Therefore, the LOS displacement is directly converted into vertical displacement under the assumption that the subsidence from underground coal fire burning is larger in the vertical direction than in the horizontal direction.Thereby, we obtained the accumulated vertical displacements of each pixel: , , ,     N T T T .Relying on the relationship between the accumulated vertical displacements and the temporal interval, the average vertical displacement rates can be modeled as follows by the least square method [44]:

Comparison of Our Results with Coalfield Fire on-Site Surveys
The average subsidence velocity map generated with the ALOS PALSAR data is shown in Figure 6.It should be noted that no highly coherent points can be found in some regions of the study area, which results in some black holes in some subsidence bowls.The areas where coal mining activities are frequent, which become complex deformation regimes, are easily influenced by temporal decorrelation, which seriously hinders the detection of subsidence in the areas of interest for this study.In spite of this, the DInSAR driven result can still reveal characteristics of land subsidence caused by underground coal fires in the Wuda coalfield.It can be seen that the remarkable subsidence mainly occurred inside the Wuda coalfield, where the land was deforming at rates ranging from −30 cm/yr (subsiding) to 30 cm/yr (uplifting), whilst there was a relatively stable pattern (−5 to 5 mm/yr) in areas near the urban center of Wuda city.Overlaying the coal fire distributions that were mapped in field works provided by Beijing Remote Sensing Corporation of Shenhua Group in 2008 [24] on the displacement velocity maps (Figure 6), reveals that the subsidence signals derived by DInSAR analysis agree well with the coal fire data that were derived via on-site surveying.There are several noticeable subsidence bowls located in the coal fire areas, particularly in coal fire zones 5, 7, 10, 11 and 12, where subsidence rates were up to −30 cm/yr.Several other subsidence bowls outside the coal fire areas have also been observed in the other parts of the Wuda coalfield that are inferred to be coal mining areas.The observed subsidence rates in the coal mining areas are also as high as −30 cm/yr.
The correlation between coal fires and subsidence observations is further demonstrated by a quantitative statistical analysis, in terms of the mean and the standard deviation of the deformation rate corresponding to each of the coal fire zones, as listed in Table 2.The statistical results show relatively evident subsidence signals within the coal fire zones 5, 7, 10, 11 and 12, with mean deformation rate of −8.52 cm/yr, −3.92 cm/yr, −4.6 cm/yr, −8.93 cm/yr and −6.07 cm/yr, respectively.
In order to further validate the difference between the subsidence caused by coal fire burning and by coal mining, we selected two authentication areas from the subsidence rate map; one is the zone illustrated by rectangle A in Figure 6, and the other is the area illustrated by rectangle B in Figure 6, both of which experienced significant subsidence.After the field survey, we found that area A is a new subsiding area caused by coal fire burning, which is identified as coal fire zone 18 by the Wuda Mine Bureau.We also found that the reason for the subsidence in area B is the intensive coal extraction activity with mechanized longwall methods.Since the deformation amount induced by coal mining activity is greater than the critical value utilizing DInSAR technique, some black holes can be seen in the area of the subsidence bowl on the subsidence rate map.The field survey pictures of the authentication area are presented in Figure 7.In addition to the subsidence rate map, we also obtained the subsidence history of the time series from October 2006 to January 2008 (Figure 8) to monitor the displacement evolution in the Wuda area.It clearly reveals an acceleration of the surface deformation within the Wuda coalfields, which can be confirmed by the fact that the coal fire evolution in this area.Moreover, the coal fire distributions surveyed in 2008 [24] are overlaid on the six accumulated deformation maps from the six time intervals to compare the observed displacements with the on-site surveying.From the comparison, it can be observed that the magnitude of the land subsidence within the coal fire zones increased rapidly.The land over some obvious subsidence bowls, such as coal fire zones 5, 10, 11 and 12, up to −40 cm between 18 October 2006 and 21 January 2008 (about 460 days) shows accelerated subsidence.It is also worth noting that small-scale land subsidence occurs in some coal fire zones (e.g., coal fire zones 16, 17 and 18) during the period of coal fire development, which is further evidence that such areas are new subsiding areas.Displacement history for two typical points within the coal fire zones 5 and 18 are illustrated in Figure 9.We have to conclude that the DInSAR technique enables us to reveal the major displacement trend caused by underground coal fire burning in the Wuda area.However, taking into account that underground coal fire burning is not the only factor in the resulting surface deformation, only using the DInSAR technique to distinguish land subsidence caused by coal fire burning from the other sources of surface deformation is difficult, particularly the other mining-related activities.Therefore, some other data are required, such as thermal infrared and high resolution optical multispectral data, to aid the extraction of the subsidence signal related to underground coal fires.6) within the coal fire zones 5 and 18.The circles signs represent the points in coal fire zone 5, which reveals an accelerated land subsidence phenomenon in this zone.The plus signs represent the points in coal fire zone 18, which indicates that this zone are newly subsiding areas during the period of the coal fire development.

Comparison with Thermal Anomaly Information Based on Thermal Infrared Data
It is common knowledge that there exists a direct relationship between coal fires and thermal surface anomalies.Therefore, thermal satellite imagery has been used to delineate and analyze coal fires [21,45,46].ETM, one of the Moderate-resolution Imaging Spectroradiometers, has one thermal band (10.4-12.5 Am) with a 60 m spatial resolution.ETM data was applied in this study for its high spectral resolution in the infrared part of the electromagnetic spectrum.12 ETM scenes from October 2006 to September 2009, covering the Wuda coalfield, were used in this study.Because most of the coal fires showed up in one ETM thermal infrared band, an effective and simple thermal anomaly detection algorithm for ETM had to make use of a single band concept.Here, a density slice method was used to form the thermal anomaly field using the average digital numbers (DN) map generated by 12 thermal infrared bands of ETM data.Thresholds of density slice to distinguish coal fire and non coal fire pixels were derived on a pixel-by-pixel basis using a moving window technique, thus allowing coal fire detection in large, thermally inhomogeneous areas.The thermal anomalies revealed show that even very subtle thermal anomalies are detectable on the thermal channels of the ETM sensor.Of course, these detections will be more reliable the hotter or larger an anomaly is.The thermal anomaly contours were then generated from the 2-D thermal anomaly field and plotted on the DInSAR derived LOS deformation rates map that is overlaid by the coal fire distributions from field investigation, as presented in Figure 10.A qualitative visual inspection of Figure 10 indicates that the thermal anomalies agree well with the coal fire data that were derived via on-site surveying especially in coal fire zones 11 and 12. Furthermore, a clear correlation is quantitatively highlighted between the displacement velocities of pixels and their corresponding average digital numbers (DN) higher than thresholds that were regarded as thermal anomalies in coal fire zones 11 and 12 (see Figure 11).
However, we can observe from Figure 10 that there are some thermal anomalies outside coal fire zones.This is probably because the daytime thermal imagery from ETM (taken at 10:30 am local time) was strongly influenced by solar heating and surface characteristic effects.In addition to this, thermal emissions of large-scale industrial complexes such as cement or coking plants, as well as steel factories in the regions, may be mistaken as false alarms in coal fire detection.It is worth remarking that the area illustrated by rectangle B in Figure 6 did not exhibit thermal anomalies in Figure 10, which demonstrated the reason for land subsidence in this area is not related to underground coal fires, but rather to coal mining activities.Therefore, a synergistic analysis of the combination of the thermal surface anomalies with DInSAR observations was proposed to distinguish land subsidence effects related to coal fires from those related to other subsidence sources.

Figure 11. (a)
A cross-plot between the displacement velocities of pixels and their corresponding average digital numbers (DN) higher than the threshold of 163 that were regarded as thermal anomalies in coal fire zone 11.(b) A cross-plot between the displacement velocities of pixels and their corresponding average digital numbers (DN) higher than the threshold of 160 that were regarded as thermal anomalies in coal fire zone 12.

Analysis of ALOS PALSAR Interferometry
In our study, we showed that ALOS PALSAR interferometry produces excellent results in the Wuda coalfield, relying on the longer radar wavelength of L-band with a relatively high spatial resolution of 10 m after multi-looked.L-band InSAR images generally maintain much higher coherence in cases of high deformation gradient due to their longer radar wavelength [47,48].In addition, the PALSAR data has a larger critical baseline, i.e., 13.1 km in FBS mode, resulting in high coherence interferograms in coalfield areas [49].Therefore, interferometric decorrelation due to rapid subsidence and surface breaking caused by coal fires or coal mining could be relaxed by using L-band InSAR analysis.ALOS PALSAR interferometry is thus less impacted by temporal decorrelation and phase saturation.In addition to this, due to the higher maximum deformation that is detectable, PALSAR data enables the detection of land subsidence signals that cannot be monitored by C-band ENVISAT and RADARSAT-2 data in large and highly non-linear subsidence areas, such as the Wuda coalfield.Thanks to the relatively high spatial resolution compared with ENVISAT data, small-scale subsidence signals can be detected and much more detailed spatial behavior of land subsidence in the Wuda coalfield can be revealed.
It should be noted that there are also some limits of ALOS PALSAR data for mapping subsidence due to coal fires.Firstly, the longer revisit time of 46 days has hampered interpretation of the more detailed temporal resolution of the land subsidence behavior accompanying the coal fires.Secondly, the L-band radar signal is less sensitive to deformation than the C-band and X-band signals because of its comparatively longer wavelength.More precise results can not be obtained with L-band data that have a longer wavelength.It is possible that the limits above could be relaxed by means of a joint analysis of other InSAR systems such as TerraSAR-X.This new-generation satellite, which was launched in 2007, can provide data suited to interferometry at high spatial resolution.Furthermore, TerraSAR-X is operated with an 11-day repeat interval, which is significantly shorter than the 46-day intervals of the ALOS PALSAR satellite, and this relatively shorter revisiting time also favors a fast build-up of interferometric data stacks.In spite of its shorter wavelength, which is influenced more by temporal decorrelation and phase saturation, the temporal and spatial sampling of deformation patterns is significantly better than C-band data and potentially useful for resolving the faster deformation as it is occurring in the case of coal fire burning.
Another limit of ALOS PALSAR interferometry is residual topographic error that results from the significant topographic changes in Wuda coalfield, which is believed to be the single largest source of error limiting the interpretation of results in this study.In [50], it was shown that the residual topographic noise could not be considered random in the case of ALOS PALSAR utilizing the DInSAR technique, because the baseline lengths within the set of interferograms are large and variable.This paper describes an approach for estimating the DEM error adopting a 1-D regression method, which is particularly applicable to L-band ALOS PALSAR data because of the smaller number of images presently available and the observed temporal variation of spatial baselines.
Due to only a single satellite direction being available for this study, the LOS subsidence is directly converted into vertical subsidence based on the assumption that the subsidence from underground coal fires is larger in the vertical direction than in the horizontal direction.However, it should still be noted that underground coal fires can result not only in significant vertical motion but also in lateral movement due to volume loss underground.If SAR data from both ascending and descending satellite orbit directions were available, the subsidence along both the vertical and east-west directions could be derived with significantly higher accuracy, which could dramatically improve the knowledge of land subsidence behavior related to coal fires.

Comparison with PSI Techniques
PSI techniques are capable of measuring ground deformation to sub-centimeter accuracy by using multiple interferometric SAR data, especially in areas with slow terrain motion [51].However, there are several factors that have limited the applicability of PSI techniques to map the land subsidence related to coal fires in Wuda coalfield: (1) Coalfields are usually located in rural areas where only a small number of man-made features are available for imaging.The limited number of such features, the large deformation rate, and the rapidly changing deformation greatly reduce the density of stable scatters identified by PSI techniques, which has limited the ability of PSI techniques to monitor the deformation due to underground coal fires.
(2) The problem of resolving the ambiguity of phase measurements prevents the PSI technique from being able to monitor rapidly changing deformation [52].The coal fires induced subsidence exceeds the absolute limit imposed by the sampling theorem if there is no a priori information available.Many PSI approaches make use of a linear deformation model in their estimation procedures.The linear model assumption can have a negative impact on the PSI estimates for all deformation phenomena characterized by non-linear temporal deformation behavior, i.e., where the linear assumption is not valid.In coalfield areas where the deformation shows "significantly non-linear motion" the PSI products may lack PSs, due to the fact that the PSI phases do not fit the (incorrect) linear model.This may represent a critical limitation because PSI may be unable to provide deformation measurements over the most interesting deformation area [53].(3) The amplitude dispersion criterion [27,28] is used to select PS pixels.In general, at least 25 images are needed to identify PS pixels [54], though this number can be reduced to about 12 [55].Unfortunately, the amount of ALOS PALSAR data available is limited (seven images for this study) due to the long revisit cycle and the insufficient amount of archived data available, which leads to identifying the wrong PS pixels.In addition, it is complex to resolve the temporal phase ambiguity in time with such small datasets, especially when the deformation is very large and its gradient is steep.(4) PSI techniques measure deformation time series on a point-by-point basis, whereas deformation measured on a pixel-by-pixel basis could provide more information to aid understanding of the mine subsidence phenomenon.As a consequence, we have to remark that using PSI approaches is not possible in our study area with its small stack data.
The results in Figure 6 and Figure 8 show that the proposed approach in our study is successful for monitoring underground coal fires induced land subsidence in Wuda coalfield.Even under the circumstances that the amount of data available over the study area is relatively small and deformation due to underground coal fires is very large and highly non-linear, the approach described in our study provided information to improve detection of the land subsidence related to coal fires in the study area.In addition, the deformation time series and deformation rates are measured on a pixel-by-pixel basis in our approach.Compared with the sparse pixels identified by PSI techniques, our approach enhances the understanding of underground coal fires induced subsidence phenomena.However, as mentioned earlier, only using our approach to distinguish the land subsidence caused by coal fires from the other sources causing surface deformation is still a challenge, particularly the various mining-related activities.In our study, a joint analysis of DInSAR subsidence signals and thermal anomalies mapped from thermal infrared images was proposed, which will help us to distinguish between different subsidence causes.Furthermore, some studies of coal fires have suggested that underground subsidence caused by the coal fires can result in systematic orthogonal patterns of opening fissures at the surface, which is different from the patterns of fissures caused by underground mining [24].As a result, high-resolution remote sensing data can be used to extract the surface fissures related to the different subsidence causes, taking account into the different geometrical features of opening fissures.However, normally remote sensing can only acquire information of coal fires on the earth surface, whilst geophysical exploration methods are used to probe for information on the combustion of underground bodies.This is thus a relatively comprehensive and effective method to integrate geophysical and remote sensing methods for coal fire detection.The combined method can overcome the shortcomings of a single probe method [1].

Conclusion
In this paper, the capability of DInSAR interferometry using ALOS PALSAR data for detecting and monitoring long-term ground displacement due to underground coal fire burning was demonstrated.Based on our work we concluded that the ALOS PALSAR interferometry produces more excellent results in the Wuda coalfield compared with the C-band ENVISAT InSAR images, relying on the longer radar wavelength of the L-band with a relatively high spatial resolution.This paper also described an approach to map the large and highly non-linear subsidence induced by underground coal fires based on a small number of SAR images in the Wuda coalfield, one of the largest and most well-known coal fire areas in China and worldwide.The DEM error induced by significant topography changes in the Wuda coalfield was estimated.This effect is particularly clearly observed in the results of deformation signals before and after applying the topographic correction.We also found that the stratified atmospheric component should be negligible, due to the small range of height difference in the Wuda coalfield.The DInSAR results between October 2006 and January 2008 were overlaid with the coal fire distributions that were mapped in field work in 2008.A good agreement between the LOS displacement velocity value and the on-site surveying was found.In addition, a quantitative comparison between coal fire zones and subsidence observations has been presented.The statistical results show relatively evident subsidence signals within coal fire zones 5, 7, 10, 11 and 12, with mean deformation rates of −8.52 cm/yr, −3.92 cm/yr, −4.6 cm/yr, −8.93 cm/yr and −6.07 cm/yr, respectively.These comparisons demonstrate that the capability of ALOS PALSAR data and the proposed approach have remarkable potential to reveal the major land subsidence trend of interest.However, it is worth noting that apart from the land subsidence areas inside the coal fire zones we observed, we also found some areas outside the coal fire zones with similar subsidence signals.This makes it difficult to separate land subsidence effects related to coal fires from those related to other subsidence factors when relying only on the DInSAR technique.In order to validate the observations above, we chose two authentication areas, one from an area inside the coal fire zone, and the other from outside.Both of them had a similar subsidence signal.From the field survey, we found that the former one is induced by coal fires and identified as coal fire zone 18, while the latter one is caused by intensive coal mining activities.We also found that the subsidence history time series could provide a temporal evolution in addition to the spatial extent the land subsidence behavior accompanying the coal fires, which indicated that several coal fire zones suffered accelerated ongoing land subsidence whilst other coal fire zones were newly subsiding areas during the period of the coal fire development.Furthermore, we tried to use thermal surface anomaly information extracted from the thermal infrared images based on ETM data to compare this with DInSAR observations.The comparison not only reveals a good correlation between the thermal surface anomalies and DInSAR results, which proved again that the DInSAR technique is capable of monitoring the land subsidence related to underground coal fires, but also demonstrated a synergistic analysis of combining thermal surface anomalies with DInSAR results to distinguish the land subsidence effects related to coal fires from those related to other subsidence sources.
Note that because of the long revisit time of 46 days, it is not possible to obtain a higher temporal resolution subsidence signal in the area of interest.In addition to this, more precise results cannot be obtained using ALOS PALSAR data, since the L-band radar signal is less sensitive to deformation than C-band and X-band signals.In future work, a joint analysis of other new-generation SAR sensors, such as TerraSAR-X and COSMO-SkyMed, which have higher spatial resolution and relatively shorter revisiting time and more efficient orbital control, will to be carried out to reveal the potential for a more thorough and comprehensive understanding of the spatial and temporal evolution of land subsidence behavior induced by underground coal fires.It is also worth noting that it is still a challenge to separate land subsidence related to coal fires from that related to coal mining activities when using only the DInSAR technique.Although we have made a attempt to identify coal fire induced land deformation using thermal infrared data in this paper, another subject of our future work is a synergistic analysis of remote sensing methods and geophysical exploration methods to improve the separation of subsidence effects due to coal mining from those due to coal fires burning, which will enhance the understanding of underground coal fires induced subsidence phenomena.

Figure 2 .
Figure 2. Field pictures of different types of land subsidence in the Wuda coalfield.(a) A slow trough-like land subsidence related to underground coal fires happened in coal fire areas.(b) A sudden surficial cracks related to underground coal fires happened in coal fire areas.(c) Land subsidence caused by underground mining activities in coal mine areas.
signal due to atmospheric inhomogeneities at the two radar acquisition times, t A and t B .The fifth term, ,, m Orbit x y  , accounts for the orbit error.The last term, ,, m Noise x y  , is caused by many factors, including thermal noise, image miscoregistration, wrong focusing parameter and decorrelation of the interferometric signal.

Figure 3 .
Figure 3. Block diagram of the approach implemented in our study.

Figure 4 .
Figure 4. Baseline-time plots relevant to the PALSAR images and interferograms used in this study.

Figure 5 .
Figure 5. Example of the DEM error correction (here applied to the interferogram covering 18 October 2006 to 18 January 2007).(a) The sensitivity to topography was apparent for the original interferogram because of B⊥ being large, about 1812 m.The spatial unwrapping of the interferogram is unsuccessful due to inaccurate compensation of the topographic phase, which results in wrong deformation pattern.(b) DEM correction was estimated and subtracted from the original interferogram.Based on these, it is then typically possible to unwrap this interferogram with longer baselines.The residual topographic noise is significantly reduced on the deformation map after DEM error correction.


is the unwrapped phase in refined interferogram M for the pixel of range and azimuth coordinates (x,y), and ,

Figure 6 .Table 2 .
Figure 6.DInSAR derived subsidence rate map in Wuda coalfield from October 2006 to January 2008.The colour tones from light yellow to red correspond to subsiding areas; the tones from dark blue to sky blue indicates slight uplift; the light green correspond to stable areas.This result is superimposed on a Google Earth image and the magenta cross represents the location of the reference point, which is loated in Wuda urban area.The red polygons illustrate the coal fire boundaries mapped in the field survey in 2008 and the rectangles labeled A and B indicate field survey spots.

Figure 7 .
Figure 7. Field survey pictures of the authentication area showing (a) a new subsiding area (area A marked in Figure6) caused by coal fire burning, which was identified as coal fire zone 18 by the Wuda Mine Bureau, and (b) a coal mining area (area B marked in Figure6) where intensive coal extraction activity was carrying out with longwall methods.

Figure 8 .
Figure 8. DInSAR derived subsidence time series in Wuda coalfield with respect to 18 October 2006 for 6 dates; 3 December 2006, 18 January 2007, 5 March 2007, 5 June 2007, 6 December 2007, and 21 January 2008.The colour tones from light yellow to red correspond to subsiding areas; the tones from dark blue to sky blue indicates slight uplift; the light green correspond to stable areas.This result is superimposed on a Google Earth image and the magenta cross represents the location of the reference point, which is loated in Wuda urban area.The red polygons illustrate the coal fire boundaries mapped in the field survey in 2008.

Figure 9 .
Figure 9. Displacement history for the two typical points (marked in Figure 6) within the coal fire zones 5 and 18.The circles signs represent the points in coal fire zone 5, which reveals an accelerated land subsidence phenomenon in this zone.The plus signs represent the points in coal fire zone 18, which indicates that this zone are newly subsiding areas during the period of the coal fire development.

Figure 10 .
Figure 10.Examples of detected thermal anomalies, which can be related to coal fires, in the Wuda Coal Field.Thermal anomalies (purple and light purple contours) are plotted on a DInSAR derived subsidence rate map and overlaid by coal fire distributions from field investigation.

Table 1 .
The PALSAR data used in this study ,,

Table 2
are only results of statistical computing, and do not representing the measurement accuracy of the displacements.