Land subsidence prediction using recurrent neural networks

In an environment, one of the natural geological hazards is land surface subsidence. Underground mining and subsurface coal fires are primarily responsible for subsidence of land. Activities, such as, over-exploitation of coal, minerals, groundwater and petroleum resources, depillaring of the existing galleries and water logging of the relinquished galleries are the major factors resulting in subsidence. The deformation is primarily measured in terms of change in ground elevation values (Z-dimension) at different time intervals at identified ground locations. All the conventional and exiting techniques have certain limitations in monitoring and predicting land surface subsidence. In this work, we predict the land subsidence in Jharia Coalfield, Dhanbad, India for one year in the interval of twelve days on the datasets collected through a monitoring technique called Modified PSInSAR. The sample datasets contains 14 locations and 67 previous land subsidence value calculated from each location. We train and test predictive models and perform the prediction of the land subsidence using Vanilla and Stacked long short-term memories. Finally, we demonstrate the predicted deformation values of the 14 locations for one year. The prediction model shows the subsidence rate in Nai-dunia basti near Jharia, Dhanbad is alarming as 93.8 mm/year where as Digwadih and Godhar showed the critical rate as 82 mm/year.


Introduction
In an environment, one of the natural geological hazards is land surface subsidence. Underground mining and subsurface coal fires are primarily responsible for subsidence of land. Activities, such as, over-exploitation of coal, minerals, groundwater and petroleum resources, depillaring of the existing galleries and water logging of the relinquished galleries are the major factors resulting in subsidence (Kumar et al. 2020;Pandey et al. 2016;Chatterjee 2006;Jianjun et al. 2012). Subsidence vulnerability becomes more in those areas where large underground voids have been created by extracting coals, ores, etc. (Ishwar and Kumar 2017;Qin and Perissin 2015;Engelbrecht and Inggs 2013;Paradella et al. 2015;Gupta et al. 2014;Guang et al. 2009;Miao et al. 2008). It is very terrible as it involves human loss and great loss of national properties. It also affects the surface and subsurface water resources and the ultimate result is the degradation of the environment. Mine subsidence can take a shape of disaster in inhabited areas if preventive measures are not taken in time. Unfortunately, the increasing demand of energy and mineral resources worldwide has brought mechanization and rapid expansion of mining activities. With the increase of mining activities, there will be a corresponding increase in mine subsidence problems causing more damage unless proper subsidence control measures are taken. Control measures are directly dependent on detection, monitoring and prediction pattern of subsidence area. Spatial-temporal monitoring and need for precise calculation of land subsidence for mapping in zonal management and corresponding control of surface deformation caused by both underground mining and subsurface fires. However, the effectiveness of the preventative and protective measures of subsidence greatly depends upon the accuracy of subsidence monitoring and associated prediction parameters.
Most of the existing techniques of subsidence monitoring are based on ground-level that rely on the instruments as Precise Level (PL), Auto Level (AL), Digital Level (DL), Total Station (TS), etc. These instruments with associated field survey techniques provide a highly relevant measurement with a millimeter accuracy, but very cumbersome in comparison to modern geo-spatial techniques. Ground-based subsidence monitoring methods are also not safe because measurements are required to be taken along subsidence prone areas. Global Navigation Satellite System (GNSS) based techniques have made the measurement quite easier in terms of portability of the instruments and satellite dependencies for the data acquisition (Chatterjee et al. 2015;Wang et al. 2011;Jing-Xiang and Hong 2009;Lü et al. 2008). However, GNSS suffers from some of the same vexing problems as physical movement required in subsidence prone areas with very costly instruments. The limitations of GNSS techniques have been overcome by spaceborne imaging techniques. Spaceborne subsidence monitoring has emerged as a better technique after the development of Synthetic Aperture Radar (SAR) Interferometry. The radar satellites can observe almost anywhere on the surface with ease, even during darkness and cloudy conditions, which makes it invaluable for subsidence monitoring. Interferometric Synthetic Aperture Radar (InSAR) technique uses two or more SAR images to generate maps of surface deformation or digital elevation, using differences in the phase of the waves returning to the satellites. Its limitation is accuracy as atmospheric errors introduce several errors. When images are acquired at different times (temporal baseline), utilizing the Differential SAR Interferometry (DInSAR) technique, it is possible to measure the changes of the surface elevation. These measurements are shown by a series of colored bands, the so-called fringes or interferogram as shown in Fig. 1, as the Angarpathra, Godhar, and Bastacola mines of Jharia coal field (JCF). 1 The DInSAR technique also has limitations of accuracy and reliability of specific area deformation. The Persistent Scatterer Interferometric Synthetic Aperture Radar (PSInSAR) technique shows relatively better results for deformation of the land surface and with higher accuracy.
The Long short-term memory (LSTMs) are the variants of recurrent neural networks (RNNs) mainly used for time series forecasting. The LSTMs are categorized into Vanilla, Bidirectional, and Stacked LSTMs. The LSTMs have been used to predict the Land subsidence in (Li et al. 2020;Qiao et al. 2020;Mubashar et al. 2021;Pu et al. 2018) and Bidirectional LSTMs are used in (Shen et al. 2021;Qu et al. 2019). The LSTMs are efficiently predicting the Land subsidence from the given datasets with better accuracy. It is better to use the Bidirectional LSTMs because it processes the data in both forward and backward time order. The Vanilla and Stacked LSTMs are the variants of LSTMs, where Vanilla LSTMs reduces the memory usage four times better than the LSTMs. So that we can minimize the additional memories during evaluation, and it also speeds up the training and testing. The stacked LSTMs are performing the operations simultaneously using a hierarchy of hidden layers, and it evaluates the complex data very easily with improved accuracy (Gangopadhyay et al. 2018;Liang et al. 2018).
Our research is based on the processing and analysis of free available SAR datasets as SENTINEL-1A SAR data. These SAR Datasets are available on a specific web portal in archive format for the researchers' but lagging by a few months to be available some times. The modified PSInSAR technique is giving large Excel sheet of movement of location termed as PS point and pictorial view with coded colour to identify subsidence on map. The Vertical movement of the PSs near GNSS control stations ( Fig. 2) have been validated with processed GNSS data output. Here an excel sheet of 26,000 rows and 100 columns had been generated but we ran a sample of 14 rows data sets to study the LSTM module for 1-year prediction. The prediction is for the future so the output of the processed SAR data in Excel format is transferred to the LSTM module for further prediction at 12 days' intervals for a year. This prediction is giving alarming information to deal with settlement and inhabitants residing in the vicinity and saving the life and economy of the nation.
The contributions of this article are summarized as follows: • We develop a modified PSInSAR processing Technique by changing the range, subset area, pixel value, boxcut filtering on SAR images. This processing technique helps in increasing the accuracy of the deformation velocity. The outcome of the approach is in the form of colour coded subsidence locations, referred as Permanent Scatterers (PSs). In our experiment, we considered SAR images of Sentinel-1A of JCF and converted them into PSs • The movement of PSs are validated with the processed GNSS Data. This strategy helps in assessment of accuracy of deformation obtained by modified PSIn-SAR. After validation, the PSs of higher coherence are selected for model testing and converted into numerical format and stored in excel using SAR data processing tool called SARPROZ. • The excel sheet consists of 26,000 rows and 100 columns generated using SARPROZ tool and we use a sample of 14 rows datasets to predict the land subsidence using LSTM. • Train and test the data for accurate predictions of the land surface subsidence using Vanilla and Stacked LSTMs. • We present the comparisons of observed and predicted values of land subsidence in JCF of mining as well as GNSS-based locations. Finally we also present the predicted values of land surface subsidence for one year. • All the predicted locations are finally validated with physical site visits accordingly. The prediction showed the subsidence rate in Nai-dunia basti in JCF is alarming as 93.8 mm /year whereas Digwadih and Godhar showed the critical rate as 82 mm/year.
The remaining sections of this article is arranged as follows. In Sect. 2, we review some of the existing but related approaches for land subsidence prediction in Mining areas. In Sect. 3, we provide the preliminaries used in the problem formulation. In Sect. 4, we describe the proposed model in detail. In Sect. 5, we present the experimental results of the proposed algorithm. This paper is concluded in Sect. 6.

Related work
In recent years Jharia Coal field (JCF) has witnessed large number of land subsidence and coal fires. Underground activities and coal fires are main cause of such incidents. Nearly 150 such cases have been reported in the past decade. In (BCCL 2008), Master plan prepared by Bharat Coking Coal Limited (BCCL) in association with central mine planning and design institute (CMPDI) also signifies huge losses in BCCL due to mine fire and induced subsidence in rural as well as urban area of JCF. The researchers have used conventional DInSAR techniques for monitoring long-term land subsidence phenomenon and achieved the deformation by analyzing the fringes obtained by SAR images Strozzi et al. 2001;Lanari et al. 2004;Raucoules et al. 2003;Amelung et al. 2000). However, the conventional DInSAR techniques have limitations in terms of (1) very small spatial baseline (\200 m) (2) baseline dependent accuracy of external DEM (3) no reduction of atmospheric phase. To overcome the limitations associated with conventional DInSAR techniques, a first-generation time series InSAR (Advanced DInSAR) technique was introduced by (Ferretti et al. 2001). Many researchers have worked on advanced DInSAR techniques and developed different approaches to deformation analysis (Hooper et al. 2004;Hooper and Zebker 2007;Hooper 2008;Crosetto et al. 2005;Mora et al. 2003;Schmidt and Bürgmann 2003). Many studies have been taken out internationally for land subsidence monitoring using Advance DInSAR techniques (Dong et al. 2013;Qin and Perissin 2015;Chatterjee et al. 2015Chatterjee et al. , 2016Ishwar and Kumar 2017;Engelbrecht et al. 2011;Gupta et al. 2014;Przyłucka et al. 2015).
In order to overcome the limitations of second generation advanced DInSAR techniques, the development of the Persistent Scatterer Interferometric Synthetic Aperture Radar (PSInSAR) technique took place to detect land deformation at the millimeter level. The PSInSAR technique is the geodetic SAR processing technique that uses two or more SAR images to generate maps of topography or deformation of the Earth's surface (Prati et al. 2010;Hooper et al. 2012;Mura et al. 2016).
The PSInSAR technique has been applied in the coal filed or nearby areas to detect land deformation using C-Band and L-Band SAR Data (Abdikan et al. 2014;Thapa et al. 2016). The short wavelength C-band is more suitable to detect slow velocity subsidence if the optimal baseline is maintained, whereas long wavelength L-band can effectively detect rapid velocity subsidence (Chatterjee  Abdikan et al. 2011;Yue et al. 2011). Its limitation is as less number of PSs achieved due to exclusion of partially correlated scatterers for the analysis resulting in some information losses in the vicinity. The authors in (Perissin and Wang 2011) have developed a modified PSInSAR approach which is also called Persistent Scatterer Interferometry (PSI). In PSI the correlated scatterers are also included along with permanent scatterers to increase the point target density in the highly susceptible area for decorrelation (Sefercik and Soergel 2014). PSI is more reasonable than the Advance DInSAR time series for recognizing most stable scattered pixels where pixel properties don't fluctuate with time and radar look point (Crosetto et al. 2016;Mura et al. 2016).
A study conducted by Central Ground Water Board (CGWB), in Lucknow city, indicated land subsidence is likely to occur due to over exploitation of groundwater in the next 15-20 years if immediate step to increase recharge is not taken at some of the localities of Lucknow in Uttar Pradesh: such as Narhi, Charbagh, Rajajipuram and Gomtinagar regions may see land subsidence by 2026 (Trivedi 2020).
The authors in (Zhou et al. 2017) have studied the spatial-temporal analysis of land subsidence caused by groundwater pumping from 2010 to 2015 in the Beijing plain using the SBAS InSAR technique. 69 interferograms generated using 47 TerraSAR images were utilized to investigate the land subsidence where long haul groundwater over exploitation and the use of shallow metropolitan space has prompted land subsidence. The highest yearly land subsidence rate was 146 mm/year from 2011 to 2015. The study between the SBAS InSAR results and the ground leveling measurements demonstrated that the SBAS InSAR results accomplished an accuracy of 2 mm. This research work is aimed at the study of the feasibility of the modified PSInSAR technique with C-band SAR data for finding the slow surface deformation caused by coal mine fire and underground mining activities in JCF. Also, a multi-temporal analysis of SAR images of ENVISAT ASAR has been carried out for monitoring and mapping of temporal land subsidence of the area under study. The modified PSI technique has proven its ability to detect land subsidence over the vegetated and rural areas. It also resolves low spatial density of permanent scatterers by considering partially correlated scatterers as permanent scatterers (PSs) and extracting information from these PSs. The study has been focused on detecting continuous slow rate subsidence of fifteen major sites of JCF. The imaging techniques also reduce the safety risk and decrease the expenses that are inherent in conventional methods due to extensive fieldwork. In this study, the SAR images acquired by SENTI-NEL-1A of the European Space agency have been processed by SARPROZ Software for deformation analysis of JCF, Dhanbad, India. Further Prediction analysis has been carried out for the determination of vertical shifting of the objects which will be helpful to the safe planning of the projects. In (Iwanec et al. 2016; sar 2019), authors have performed theoretical studies how to predict the subsidence above the single and multi-seam longwall mines. However, the work has certain limitations and the prediction is totally based on the physical conditions.

Preliminaries
In this section, we provide the preliminaries used in the proposed work. Long Short-Term Memory (LSTM) is used to time series forecasting (TSF). LSTM model that is used for univariate TSF problem is Univarite LSTM. While predicting the future values using the past observations through single series of observations and a model is a complicated issue.
As per earlier discussions, the LSTMs operate on sequential data and increasing the number of layers increases the levels of abstraction overtime on the input data. The network depth is more important than the number of memory cells considered for a layer. The LSTM with multiple layers is treated as Stacked LSTM. These upstream layers always provide the sequential output instead of a single output value. Each layer in Stacked LSTMs consider a 3D input for its memory cell and produce a single value as 2D array in an output. Each layer of the stacked LSTM is a chain-like processing framework shown in Fig. 3. The cell state in the Fig. 3 is running on top with minor interactions to just flow unchanged data. This layer can remove or add information to the cell state using the two gates such as and È. The r layer provide the values between 0 or 1 and it is described in Eq. (1).
The process framework performs mainly four operations to produce the final outcome. First, the forget gate starts removing the useless informations form the input gate as shown in Eq. 3. It provides input data to the sigmoid function (rð. . .Þ) such as previous layer output i.e. h tÀ1 and hidden layers feature information x t . The f t maps the previous layer cell state output i.e. C tÀ1 to the current cell state C t .
where W f is the weight matrix and b f is the bias vector of the forget gate. The r is computed as shown in Eq. (4).
In the next state, the LSTM performs two operations rð. . .Þ and tanh ð. . .Þ simultaneously to update the current cell state C t as shown in Eqs. (5, 6), respectively. The rð. . .Þ and tanh ð. . .Þ produces the i t and C t , respectively using previous layer output i.e. h tÀ1 and hidden layers feature information x t .
where W i and W c are the weight matrix of input gate and state update vector, respectively, and b i and b c are the bias vectors of input gate and state update vector, respectively. The hyperbolic tangent function is represented using tanh and it is computed as shown in Eq. (7).
tanhðaÞ ¼ e a À e Àa e a þ e Àa Finally, the current state of LSTM has been updated using previous two operations as shown in Eq. (8).
where is point-wise scalar multiplication of the vectors. In the last phase, the LSTM produces the output through output gate by using the current cell state information. The outputs such as O t and h t are computed as shown in Eqs. (9) and (10).
In general, the LSTM is comprised of a hidden layer followed by an output layer. Depends on the types of TSF problem, several LSTM techniques are existed in the literature among Vanilla and Stacked LSTMs are more popular.
In this work, we tested two univariate LSTM model such as (1) Vanilla LSTM and (2) Stacked LSTM.

Vanilla LSTM and stacked LSTM
While making the prediction through Vanilla LSTM, it uses an output layer and a single hidden LSTM layer. It supports the sequence data as an input. The LSTM model reads single time step in each time, unlike the Convolution neural networks (CNN) and the data represented in state format while learning the model. The LSTM can extend with multiple hidden layer which makes the deeper model called stacked LSTM and which reflects the deep learning model. Increasing the number of layers in the LSTM will increase the prediction accuracy. But, increasing the too many layers also increase the complexity as well the computational time. The additional layers recombine the representations to provide the new combinations of the predictions with increased abstraction. Instead of giving importance to the memory cells of a layer, considering the depth of network is more important. An LSTM produces the two dimensional output by taking the three dimensional input to the system while performing the learning process. This issue can be addressed by taking the output of each LSTM layer at each time stamp and set the input data with a return sequences ¼ True argument. This allow us to have three dimensional output from hidden LSTM layer as input to the next. Here we build both Vanilla LSTM and Stacked LSTM for predicting the next one Year land subsidence and compared both model with their loss and accuracy. These four steps discussed earlier will repeated over multiple LSTM layers as shown in Fig. 4.

Proposed work
The proposed model is illustrated and depict using flow chart shown in Fig. 5, which is primarily partitioned into four parts including data collection (along with masking), data augmentation, Training and Evaluation models. Using Vanilla and Stacked LSTMs, predicting land surface subsidence values for one year.

Data collection
A network of fifty eight Global navigation satellite system (GNSS) survey points established in JCF to find the subsidence magnitude during 11 phases of the survey at quarterly intervals. The GNSS data which we have collected has been processed with respect to BASE Station data. The errors during post-processing have been minimized by Precise Point Positioning. Archives data of the same time period has been made available by International GNSS Service (IGS) to reduce the ephemeris. Some important parameters are shown in Table 1 and locations are pointed on a map shown in Fig. 6. The measurement of land surface subsidence for some of the locations in JCF using GNSS are shown in Fig. 7. SAR images of SENTINEL-1A of the same period acquired from the European Satellite Agency. The 67 SAR Images of 12 days Temporal Resolution is processed in SARPROZ  Software to obtain the persistent Scatterer in the JCF. During processing of SAR data by PSInSAR technique a huge excel sheet with a lot of information has been produced in terms of deformation magnitude of temporal baseline of SAR images coherence, Displacement velocity cumulative displacement etc. Ultimately the result of deformation will be up to the last date of the SAR image used. An excel sheet of 26,000 PSs obtained with more than 0.5 temporal coherence value after the processing. A sample of 22 locations has been collected covering a few of the GNSS survey points' locality of JCF. The GNSS processed result was used to validate the occurrence of subsidence shown by the PSInSAR technique. For further prediction the excel sheet of 14 locations is converted into Dataframe with normalized data (using pandas and sklearn python packages) which inputs the LSTM. The JCF is a large coal mines area in India with 19.4 billion tonnes of available coal. Since 1916, this area suffered a coal bed fire and consumed nearly 37 millions tons of coal. It results in water and air pollution and also land subsidence in the cities of JCF (Pai and Carr-Wilson 2018). We consider various parameters from the JCF including the land type (Agricultural or Barren), transportation facility (Road, Rail, or others), Displacement velocity, depth of the seam and the SAR image availability from 03-10-2016 to 28-12-2018 on every twelve alternative days using the remote sensing. The land subsidence near to the river area is high as compared to the other areas during rainy seasons. In general, the land subsidence patterns have more diversity among the various locations around the Jharia coal field, and great challenges to predict the short-term land subsidence. The missing values during the data collection are masking based on the method used in (Che et al. 2018).

Data augmentation and preparation
For any machine learning methods including stacked LSTMs, training dataset is one of the important considerable factors. The best training dataset will be generated through data augmentation and preparation, and the process is summarized using Fig. 8. The primary goal of the data augmentation is to amplify training data by dividing  original data with overlaps to minimizes the complexity of the computations. In this work, we consider the timewarping data augmentation approach to prepare the data for training process. We operated with various warping ratios to show the variability of the synthetic training data (Rashid and Louis 2019). With this augmentation process, the training data increases 4-fold with various temporal length.
The data collected and augmented from the remote sensing is one dimensional, where the LSTM requires 3D input data. In data preparation phase, the data initially split in to multiple short sub-sequences and then reshape subsequences. We illustrate the data preparation through an example for better understand. The learning process of LSTMs follows a sequence of inputs to an output. An Example of Univariate time series: {6, 8, 10, 12, 13, 16, 18, 20, 22, 24, 26}. Split these series in to multiple input/ output samples. Here, the input can contains more that one step (in this example we consider three and represented using X) and an output as one step (treated as Y) as shown below: X Y ½6 8 10 12 ½8 10 12 13 ½10 12 13 16 ::: ::: ½20 22 24 26 In this way, the input data split and reshare the data for providing the input to training or testing modules. The input data categorized into training and test data. The training data used in Stacked LSTM learning models where as the test data used for predicting the results.

Training and evaluation model
The training data augmented and prepared for input to the Stacked LSTM. The over-fitting is avoided by adopting the dropout strategy after each LSTM for better generalization. The dropout rate is approximately 10% (Wei 2020

Experimental results
In this section we evaluate the prediction accuracy of the land subsidence through Vanilla LSTM and Stacked LSTM. Initially, we discuss the model construction followed by the results analysis.

Model construction
In this work, we collect the data from 14 locations around Jharia coal field from Jharkhand state in India between 3rd Oct 2016 to 28th Dec 2018, i.e. 817 days. We consider 14 environmental conditions and each condition 817 samples are collected. All these conditions are segmented with the length of 128 and 60% overlap. In this work, we use 64, 32, and 32 hidden units in LSTM layer 1, layer 2 and layer 3, respectively. The input layer has 128 units which is equal to input sample dimension. We can fit the training dataset once the complete model is defined. For the whole dataset, 70% is used for training and remaining used for testing. The Root Mean Square Error (RMSE) is an indicator to evaluate the performance of training model.

Results analysis
In this section, we compare the simulation results of Stacked LSTM and Vanilla LSTM using various metrics such as accuracy, land subsidence, and etc. The proposed approach used regression model, so Root mean squared error is a good measure of accuracy.

Root mean square errors
It is measured as the difference between values predicted by a model and the values observed and it is denoted using R m . It is calculated using Eq. (11) where X i is the actual data andX i is the predicted data of data set i, and N indicates the number of samples.

Root mean absolute percentage errors
It is measured as the difference between values predicted by a model and the values observed and it is denoted as R p . It is calculated using Eq. (11) where X i is the actual data andX i is the predicted data of data set i, and N indicates the number of samples. The comparison results of running example dataset in terms of R m and R p are presented in Table 2 The comparison between the observed and predicted train and test land subsistence (reduced level) values of Location L1-L14 with Vanilla LSTM and Stacked LSTM model are presented in Figs. 9, 10 and 11. We observe that the prediction accuracy of both the models is between 80%

Hyperparameters' influence
In this work, two hyperparameters are highly influenced on the accuracy and Time consumption which are (a) #of LSTM layers and (b) Size of the input. We conduct experiments on these hyperparameters and plot the results in Figs. 12 and 13.
In Fig. 12, we evaluate the accuracy by considering the two hyperparameters. In Fig. 12a, we assess the accuracy of fourteen data labels by increasing the # of LSTM layers in the stacked LSTM. We notice that increasing the number of layers also increases the accuracy. Similarly, we evaluate the accuracy by increasing the input data's size, and the results are plot in Fig. 12b. Here, we observe that increasing the input size also increases the accuracy. In Fig. 13, we evaluate the computation time per epoch by considering the two hyperparameters. In Fig. 13a, we assess the time computation of fourteen data labels by increasing the # of LSTM layers in the stacked LSTM. We notice that increasing the number of layers also increases the computation time slightly. We conclude that increasing the layers in the stacked LSTM also increases the computational time. Similarly, we evaluate the time computation by increasing the input data's size, and the results are plot in Fig. 13b. Here, we observe that increasing the input size also increases the computational time. In the proposed work, we used three layers to balance the prediction accuracy and computational time.

Conclusion
In this work, we have presented a new scientific approach for monitoring, validation and prediction of mining induced land subsidence in Jharia Coalfield using modified PSInSAR, GNSS and Recurrent Neural Networks. We have used the modified PSInSAR technique to collect the land subsidence value in an Excel sheet as well as colour coded pictorial views of the subsided locations. Some of PSs in the vicinity of GNSS survey points have been selected for prediction modelling. The cumulative displacement and deformation velocity during datasets period and during prediction period have been analysed and the suitable precaution mobel may be developed. We have used two variants of RNN are Vanilla LSTM and Stacked LSTM. We have collected 67 datasets pertaining to land subsidence at 14 various locations in JCF at an interval of 12 days. To perform prediction of land subsidence, we have trained, tested, and validated the predictive models by splitting datasets into 7:2:1. Finally, we have predicted the land subsidence for one year i.e, next 30 predictions in the interval of 12 days and demonstrated the prediction deformation values of all the 14 locations. All the predicted locations are finally validated with physical site visits accordingly. The prediction showed the subsidence rate in Nai-dunia basti in JCF is alarming as 93.8 mm /year whereas Digwadih and Godhar showed the critical rate as 82 mm/year. status of subsiding land vulnerable to roof collapse in the Jharia coalfield, India, as obtained from shorter temporal baseline c-band dinsar by smaller spatial subset unwrapped phase profiling. Int J Remote Sens 37 (1)