An improved method of three-dimensional displacement field generation in mining areas with a single InSAR pair

ABSTRACT An improved method combining a probability integral method and interferometric synthetic aperture radar (InSAR) for retrieving three-dimensional (3D) displacement field caused by coal mining is proposed. In this paper, fully considering dip angles of coal seams, we improved the method of obtaining 3D deformation using the linear relationship between horizontal displacement and surface tilt. Compared with the former method ignoring dip angles of coal seams, the improved method is more accurate in retrieving 3D displacement in mining areas. The improved method is first validated with a simulated coal seam with a dip angle of about 25°. The root-mean-square errors (RMSEs) in the Up-Down, East-West and South-North directions are 23.4, 33.0 and 8.0 mm, respectively. Finally, real SAR data are used to retrieve the 3D displacement fields of the Number 15235 working face of the Jiulong coal mine with a dip angle of about 13°. Due to the lack of field observations to evaluate horizontal displacement accuracy, Up-Down displacement is compared with leveling data. The RMSEs of the improved and former methods are 6.4 and 7.1 mm, respectively. The real data also indicate that the improved method is more accurate. In summary, the improved method is effective for obtaining a 3D displacement field, especially in the inclined coal seams mining area.


Introduction
After underground coal mining, the original stress balance of the rock mass around the mining area is destroyed, and the stress is redistributed until a new balance is reached (Deng et al., 2014). During this process, surface movement and deformation are common. To protect ground buildings, railways and other infrastructures, ensure the human safety and guide underground mining work safely, it is important and necessary to monitor ground surface deformation. At present, commonly used methods for monitoring surface displacement are traditional geodetic approaches, such as leveling and GPS surveying. Although the precision can reach the millimeter level, there are some disadvantages: (1) heavy labor requirements, high cost and uncertain human factors; (2) good weather conditions are necessary; and (3) low-density points are obtained.
Differential interferometric synthetic aperture radar (D-InSAR) has proven to be a feasible method for monitoring land deformation, such as the deformation caused by earthquakes (Atzori et al., 2009), landslide (Cascini et al., 2013), water exploitation (Motagh et al., 2008), volcano movement (Jung, Lu, Won, Poland, & Miklius, 2011) and coal mining (Liu et al., 2014;Przylucka, HErrera, Graniczny, Colombo, & Bejar-Pizarro, 2015), etc. Compared with common monitoring methods, D-InSAR technology has several advantages, including wide spatial coverage, high sampling rate, all-time and all-weather capability, low-labor intensity, low cost and high accuracy (HEarrera et al., 2007;Wang, Yu, et al., 2018). However, D-InSAR technology can only measure one-dimensional displacement along the radar line of sight (LOS), and surface subsidence obtained using D-InSAR technology neglects the influence of horizontal movement. Surface displacement in the mining area has characteristics of small range and large horizontal displacement. If horizontal displacement is ignored, D-InSAR accuracy in monitoring surface subsidence will inevitably be affected.
At present, five InSAR-based methods have been developed to measure surface 3D displacement. (1) The ascending and descending D-InSAR and Multiple Aperture InSAR (MAI) method (Jung, Lu, Won, Poland, & Miklius, 2011) uses ascending and descending SAR images to obtain two different LOS displacements and then uses the MAI to obtain two different azimuth deformations. Finally, the 3D deformation is estimated using the least squares method.
(3) The multi-sensor and multi-track SAR method (Gray, 2010;Ng, Ge, Zhang, & Li, 2012;Wright, Parsons, & Lu, 2004) uses SAR images from different platforms or different tracks obtained from the same area and time to calculate 3D displacements. (4) The GPS and InSAR method is described in references (Catalao, Nico, Hanssen, & Catita, 2011;Gudmundsson, Sigmundsson, & Carstensen, 2002;Hu, Li, Sun, Zhu, & Ding, 2012). (5) Finally, methods based on a prior model are described (Gernhardt, Cong, Eineder, Hinz, & Bamler, 2012;Kumar, Venkataramana, & Hogda, 2011;Li et al., 2015;Zheng, Deng, Fan, & Huang, 2018). These methods use the mining subsidence theory, glacier surface movement model and earthquake models as prior knowledge establish the relationship between horizontal and vertical displacements to retrieve the 3D displacement. In the above five methods, (1) to (3) need at least SAR images of two orbits, and the accuracy of the south-north displacement calculated from the MAI and image matching is low. The high placement cost of GPS units in method (4) restricts the density of the GPS net and reduces spatial resolution. In order to obtain the temporal evolution characteristics, some efforts have been made for retrieving 3D deformation time series. For example, Casu, Manconi, Pepa, & Lanari (2011) retrieved the range and azimuth displacement time series by adopting Pixel-offset (PO-) SBAS method, and the accuracy is on the order of 1/30th of a pixel for both range and azimuth directions. However, the approach only obtained 2D displacement time series. (Pepe, Solaro, Calo, & Dema, 2016) obtained 3D displacement time series using multiplatform SAR data by adopting minimum acceleration algorithm. To a certain extent, method of combining multiplatform/multi-sensor SAR data to obtained 3D deformation time series still exist some limitations, such as difficulties in SAR images acquisition and low South-North accuracy (Ozawa & Ueda, 2011). The mining subsidence model has been used in method (5) to assist InSAR in retrieving 3D displacement from a single InSAR pair. The prior model (Li et al., 2015) assumed that the coal seam was horizontal, while, in reality, the coal seams usually have dip angles. Therefore, compared with the horizontal coal seams, the surface deformation caused by inclined seams exploiting will show different characteristics: (1) the relationship between surface horizontal movement and tilt is different; (2) the horizontal displacements along the dip direction are not symmetrical. Therefore, the 3D deformation calculated by ignoring the dip angles in mining area will undoubtedly influence reliability and reduce the calculation accuracy.
To address these described deficiencies, an improved method for retrieving 3D displacement of mining area by integrating probability integral method and D-InSAR was developed. First, the displacement along the LOS direction in the mining area is calculated by D-InSAR, and then according to the probability integral model, a 3D decomposition equation is established based on the relationship between the vertical subsidence and horizontal movement in the inclined coal seam. Finally, the vertical subsidence and horizontal displacement are calculated. Based on experiments with simulated and real SAR data, the method can effectively improve the accuracy of 3D displacement and provide support for the early warning of construction damage in the mining area.

Relationship between vertical subsidence and horizontal displacement
Generally, the mining direction follows the strike direction, which is perpendicular to the inclined direction, to extract coal easily. Figure 1 illustrates the profile of a mining subsidence basin after the coal mine has been exploited in the strike direction. At this point, the coal seam can be regarded as horizontal.
In Figure 1, D 1 is the strike length of the working face, W m is the maximum subsidence value and the arrow on the right side of the coal seam indicates the mining direction.
Underground mining can use parameters, for example, mining depth, mining thickness, dip angle of coal seam and the size of the working face, to predict surface subsidence based on a model (Fan, Figure 1. Profile of the mining subsidence basin in the strike direction. Lu, & Yao, 2018;Wang, Li, et al., 2018). The probability integral model, which was constructed in the 1960s using random media theory, is widely used in China (He, Yang, Ling, Jia, & Hong, 1991). In this model, the horizontal displacement along profile in the strike direction is proportional to the surface tilt i.
where b and r are the horizontal displacement coefficient and major influence radius, β and H are the major influence angle and average mining depth, respectively. Considering surface inclination is the first-order derivative of subsidence, we obtain (2) Figure 2 illustrates the vertical subsidence profile (Fan, Gao, Yang, Deng, & Yu, 2015) after the coal mine has been exploited in the inclination direction. The surface subsidence from inclined coal seam mining is equivalent to the superposition of horizontal coal seams AC and BD, called the equal influence principle (He et al., 1991). In Figure 2, H is the average mining depth, D 2 is the inclination width of working face, θ 0 is the mining influence propagation angle, and α is the dip angle of the coal seam. (Fan, Gu, Qin, Xue, & Chen, 2014) provides the final relationship between horizontal displacement and vertical subsidence of an inclined coal seam.
where θ 0 is the mining influence propagation angle and a parameter in the probability integral model. Generally, θ 0 = 90°-0.5α. Here, the W(y)·cot(θ 0 ) is different with the former research.
Calculating 3D displacement fields integrating probability integral method and InSAR Assuming that the strike direction and inclination direction of the coal seams are approximately South-North and East-West (if not, we can project strike and inclination directions to South-North and East-West directions), the following conditions are true: the dimension of the geocoded LOS deformation map is n × m and H (i, j) and W (i, j) (i= 1, 2 . . . n; j= 1, 2 . . . m) are the average mining depth and subsidence of the pixel (i, j), respectively. According to Equations (2) and (3), horizontal displacements (U E , U N ) along the South-North and East-West direction are: where ΔN, ΔE are the spatial resolution of the LOS deformation map in South-North and East-West directions, respectively. To calculate U E and U N simply, Equation (4) can be expressed as: From Equation (5), U E and U N can be expressed by the subsidence. If we know the subsidence of every pixel, horizontal movement U E and U N can be obtained. According to radar imaging principles, the LOS displacement is the synthesis of the Up-Down, East-West, and South-North displacement. The resulting expression (Pepe & Calo, 2017) is: where LOS, α and θ are the LOS displacement, azimuth angle of the SAR orbit and incident angle, respectively. Then, by substituting Equation (5)   C3ði; jÞ ¼ In Equation (7), the number of equations is (n-1) × (m-1), whereas the number of unknowns is n × m and the equation is under-determined. Considering the small deformation at the edge of the mining subsidence basin, we assume the horizontal displacements of these pixels are negligible. Thus, the LOS displacement can be expressed as: LOSði; mÞ ¼ cosðθÞ Á Wði; mÞði ¼ 1; 2; :::; n À 1Þ LOSðn; jÞ ¼ cosðθÞ Á Wðn; jÞðj ¼ 1; 2; :::; mÞ Combining Equations (7) and (8), the number of unknowns is n × m and the number of equations is also n × m, so the equation has a unique solution. At the same time, the coefficient matrix of the equation is the upper triangular matrix and the equation can be solved directly by back substitution. The flow chart of 3D displacement retrieving method proposed in this paper is shown in Figure 3.

3D displacement field simulations
To verify the effectiveness of the proposed method, we simulate 3D displacement fields using Mining Subsidence Prediction System (MSPS) (Wu, Wang, & Wang, 2012;Wu & Zhou, 1999) software developed by the China University of Mining and Technology. The size of the simulated working face is 650 × 250 m, mining thickness is 5 m, subsidence constant is 0.6, and dip angle of the coal seam is 25°. The strike and inclination directions are in the South-North and East-West directions, respectively. In the inclination direction, the coal seam is low in the east. According to the geological and mining conditions for this mine, the simulation parameters are as follow. Average mining depth H is 558 m, major influence angle tangent tanβ is 2, horizontal movement coefficient b is 0.3, mining influence propagation angle θ 0 is 77.5°. The simulated 3D displacement fields are shown in Figure 4(a-c). As shown in Figure 4, the maximum Up-Down displacement is about 2.1 m, South-North horizontal displacement is almost symmetrical and maximum value is about 0.65 m. The horizontal displacement in the East direction is larger in terms of deformation dimension and value, the maximum value is about 1.2 m and maximum horizontal displacement in the West direction is about 0.6 m. Clearly observed in the simulated results, when a coal seam has a dip angle, the horizontal displacement along the inclination direction is not symmetrical, but instead toward the lower coal seam. From this observation, the theoretical and practical basis for constructing the improved method in this study is provided.
To project the simulated 3D displacement to the LOS direction, we first use two-pass D-InSAR to process two RADARSAT-2 images obtained on 17 January and 5 March 2016 to generate the LOS deformation map. Based on the azimuth angle of the SAR orbit (349.1377°) and incidence angle (35.5085°), we project the simulated 3D displacement into the LOS direction (Figure 4(d)).

Results and accuracy analysis
After obtaining the LOS displacement, we use Equations (5)-(8) to calculate the 3D displacement fields, which are shown in Figure 5. The retrieved 3D displacement fields are in good agreement with the simulated ones.
We compared the retrieved and simulated 3D displacements, and RMSEs in the Up-Down, East-West and South-North directions are 23.4 mm, 33.0 mm and 8.0 mm, respectively. The main reason of errors is that we made some simplification in this study. In Figure 2, the horizontal displacement coefficients for two hypothetical horizontal coal seams AC and BD are b down and b up , respectively, with influence radius of r down and r up , respectively (He et al., 1991). Here we simplified the horizontal displacement coefficient and mining influence radius as b and r in Equation (3).
The SAR system is insensitive to the South-North displacement because of the near-polar orbits (Wright et al., 2004). However, accuracy for the South-North direction is higher in this study. The reasons are as follows: (1) In this simulated experiment, we did not add random noise to the simulated LOS displacement; (2) For the improved method, component W(y)•cot(θ 0 ) is added to the East-West displacement; therefore, the accuracy of Up-Down displacement directly affects the East-West displacement. In order to explain the reason (2) more clearly, comparisons are made between simulated and retrieved East-West and Up-Down deformation, which are showed as Figure 5(d) and (e). We can see that Figure 5(d) have great relevance with (e), which can be the evidence for reason (2). Due to the characteristics of horizontal movement error along the inclination direction, the accuracy of horizontal displacement in the middle area of the mining subsidence basin will be lower than that in other regions using the proposed method.
We also retrieve 3D displacement fields using the former method (Li et al., 2015) and the results are showed as Figure 6. Comparing Figure 6 and the simulated 3D displacement in Figure 4, retrieved 3D displacement using the former method also keep consistent with the simulated ones, while, we find that the Up-Down displacement is smaller than that of Figure 5(a) because C1 in Equation (7) increases from 6.38 to 6.50 and C1 is a denominator. East-West displacement changes considerably, the East-West horizontal displacement is almost symmetrical, in large contrast to the simulated results. We also evaluate the accuracy, the RMSEs along the Up-Down, East-West and South-North directions are 89.8 mm, 127.9 mm and 26.2 mm, respectively. Compared with the improved method, the accuracy is lower, which indicates that considering the dip angle of coal seam is very necessary and can improve accuracy.

Study area and SAR data processing
The study area is located on the 15,235 working face of the Jiulong mine in Handan City, Hebei Province, which is shown with the black box on the left of Figure 7. The red arrow on the working face is the mining direction. The coal seam is 4.5 m thick, the subsidence constant is 0.4, coal seam dip angle is about 13°(the mining depth in the East is deeper than that in the West), and the size of the working face is 935 m × 142 m (strike length by inclination width). According to the geological and mining data summarized from other working faces for the mine, the parameters required for retrieving the 3D displacement fields are as follows. Average mining depth H is 740 m, major influence angle tangent tanβ is 1.79, horizontal movement coefficient b is 0.25, and mining influence propagation angle θ 0 is 85°.
In this experiment, five repeat-pass RADARSAT-2 images are used with the HH (horizontal transmit and horizontal receive) polarization mode and pixel resolution in the range and azimuth directions of 2.66 m and 2.91 m, respectively. Because image acquisition was from 6 November 2015 to 5 March 2016, coherence loss caused by vegetation coverage is greatly reduced, and the accuracy of surface information extraction ensured. During D-InSAR processing, two-pass D-InSAR are applied. Every two images in the acquisition time constitute interferometric pairs (Table 1). Finally, we add fourtime period LOS displacements to acquire the total LOS displacement ( Figure  8 (a)) from 6 November 2015 to 5 March 2016. The pixel resolution in the South-North and East-West directions are 3 m and 3.75 m after Geocoding. In the D-InSAR processing, Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) data (90 m × 90 m) are used to remove the topographic phases, the adaptive filtering method is applied to weaken the effect of noise on the interference phase and the minimum cost flow method is used to unwrap the differential interferogram.

Results
After obtaining the LOS displacement, the 3D displacement fields were obtained by Equations (5)-(8), which are shown in Figure 8(b-d). From Figure 8(b), we can see that there were two mining subsidence areas. Due to a lack of mining information and surface measurement, we only analyze the land deformation caused by 15,235 working face exploitation. For this mining area, the maximum vertical subsidence was about 0.2 m, and the maximum horizontal movements along the East- West and South-North were about 0.1 m, respectively. In the center of mining subsidence basin, the horizontal displacement along East-West and South-North directions is to zero, which is accord with mining subsidence law (He et al., 1991). It should be noted that the land deformation generated along the East-West direction was better than that of South-North direction, and the main reason is that the SAR system is insensitive to the South-North displacement in the current near-polar orbit of SAR satellites.

Discussion
To monitor surface deformation, two observation lines were set up along the strike direction (line AA′, 50 points) and inclination direction (line BB′, 38 points), and the leveling measurement was conducted when the RADARSAT-2 images were captured every 24 days. The point distribution is shown as Figure 8(b).  The leveling is carried out in accordance with the fourgrade leveling standard, so it is feasible to use leveling data to verify the calculation accuracy of the improved method. Figure 9 presents the comparison of the retrieved Up-Down displacement using the improved method (red line) and leveling (blue line) results, which are in good agreement. From Figure 9 we can see that the maximum land subsidence values along the strike and inclination directions were about 0.18 m, the retrieved Up-Down displacements were smaller than the leveling results. In the strike and inclination directions, the maximum absolute errors were 25.7 mm and 23.5 mm, respectively, and the RMSE of all the points was 6.4 mm.
In the former study (Li et al., 2015), the 3D displacement calculation method was derived in horizontal coal seam exploitation. In order to make a comparison, we ignored the dip angle of the coal seam to extract the Up-Down displacements which are shown in Figure 9 (green line). The results of two observation lines are also in good agreement with leveling. However, compared with the results generated by the improved method, the deviation from leveling was greater. In the strike and inclination directions, compared with leveling observations, the maximum absolute errors were 31.5 mm and 27.3 mm, respectively, and the RMSE of all points was 7.1 mm. Due to small Up-Down displacements in this area, the accuracy is not improved obviously. However, from the real SAR and the simulated experiment, we can speculate that when the deformation or dip angle is large, the retrieved results calculated by the improved method will be better than that of the former method obviously.
In reality, coal seams always have dip angles. Therefore, compared with the former method, the improved method is more practical. In addition, the proposed method only needs two SAR image of the same obit to calculate 3D deformation, so it is simple and can reduce the monitoring cost. It should be noted that in the improved method, we simplify horizontal displacement coefficient, mining influence radius and surface tilt in Equation (4), which could reduce the accuracy of retrieving 3D displacement. However, from the simulated and real SAR experiments, simplified improved method is feasible to obtain 3D displacement caused by inclined seam exploitation.

Conclusion
Considering differences in the surface movement between inclined and horizontal or near-horizontal coal seams, this study derived equations for retrieving 3D displacement fields using the linear relationship between horizontal displacement and surface tilt.
Using these equations a method for retrieving 3D displacement fields using a single InSAR pair was optimized. The resulting accuracy is higher compared with the former method, particularly for coal seams with larger inclination angles.
The method was verified using simulated data and real SAR images. In the simulated experiment, the maximum absolute errors for the Up-Down, East-West and South-North directions were 100 mm, 140 mm and 40 mm with RMSEs of 23.4 mm, 33.0 mm and 8.0 mm, respectively. Using the former method, the maximum absolute errors for the Up-Down, East-West and South-North directions were 400 mm, 560 mm and 140 mm with RMSEs of 89.8 mm, 127.9 mm and 26.2 mm, respectively. Using real SAR data, the horizontal displacement cannot be evaluated because GPS data were unavailable; therefore, we analyzed vertical subsidence alone. Using the proposed method, the maximum absolute error for the retrieved vertical subsidence and leveling was 25.7 mm with an RMSE of 6.4 mm. The former method had a maximum absolute error of 31.5 mm and RMSE of 7.1 mm. Both the simulated and real data prove that the method in this study is effective for inverting 3D displacement fields of mining areas, with an accuracy higher than the former method.

Disclosure statement
No potential conflict of interest was reported by the authors.