Research on time series InSAR monitoring method for multiple types of surface deformation in mining area

Interferometric Synthetic Aperture Radar (InSAR) has become the primary remote sensing detection method for monitoring significant surface subsidence and deformation caused by underground coal mining. Time Series InSAR has been intensively researched in mining-caused deformation monitoring for its ability to provide surface deformation time series, with Small Baseline Subset InSAR (SBAS-InSAR) and Permanent Scatter InSAR (PS-InSAR) being two classical methods. To obtain the incline and curvature in mining areas and assess the building damage, this paper exploits the SBAS-InSAR and PS-InSAR techniques to estimate surface incline and curvature, following the principle of directional derivatives. The surface beneath the villages and the industrial square near the 7221 Grout-filled working face in the Huaibei mining area of Anhui Province, China, was analyzed. The estimated surface incline and curvature with SBAS-InSAR were within the threshold values specified for Class I damage, while PS-InSAR showed <1% of measurements that exceeded this threshold. This paper provides a method for monitoring the surface deformation beneath the buildings in the mine area. Meanwhile, the dynamic analysis of the damaged buildings provides a basis for determining the influence range with grout-filled working faces.


Introduction
China is a large producer and consumer of mineral resources, and its large-scale development and exploitation have made significant contributions to the national economy and social development. However, it also seriously damages the geological environment of mines and directly affects the effective use of national resources and the everyday life safety of residents in mining areas (Bell et al. 2000). For the disaster monitoring method of surface deformation caused by mining, geodetic measurement (Sedlak and Sgem 2015;Wesolowski 2016;Gao et al. 2014;Li et al. 2018) and other means need to spend a lot of human resources and material resources. The Interferometric Synthetic Aperture Radar (InSAR) is an active microwave remote sensing technique for deformation monitoring (Bhattacharya and Mukherjee 2017;Moreira et al. 2013;Mastro et al. 2020). InSAR uses two SAR images acquired along similar orbits to detect changes occurring in ground topography. The key advantages of InSAR over traditional ground-based techniques are costeffective and wide coverage. With development in recent decades, today's InSAR mainly includes Digital Elevation Model (DEM) production (Zhang, Wang, and Chen 2012), differential InSAR (D-InSAR) (Wang et al. 2009), and time-series InSAR (TS-InSAR) deformation monitoring (Fan et al. 2018). InSAR deformation monitoring has been widely used to monitor surface deformation, such as mining (Yang et al. 2020;Sui et al. 2020), landslides (Castaneda et al. 2009), volcanic eruptions (Olivier et al. 2019;Comerci et al. 2015), glacial movements (Zhou et al. 2011;Li et al. 2009), spatial changes in sea level (Tang et al. 2021), and urban ground subsidence Ding et al. 2021;Li et al. 2016).
Small Baseline Subset InSAR (SBAS-InSAR) and Permanent Scatter InSAR (PS-InSAR) have been widely used as two classical TS-InSARs. Accurate SBAS-InSAR data provide the theoretical basis for environmental restoration, safe construction and can compensate the need for manpower and the limitation of GPS for monitoring deformation of a large, sinking open-pit mine in cold and high altitude Tianshan (Du et al. 2021). New timeseries InSAR technique combined with ascending and descending orbit data to analyze the surface deformation after mining at the French and German border. The high accuracy of subsidence and uplift data identifies more deformation areas and compensates for the insufficient number of level survey datum points (Samsonov et al. 2013). The surface deformation values calculated by SBAS-InSAR were validated using GPS data at a mine site in Seoul. The support vector regression algorithm used the validated surface deformation values to obtain high accuracy surface deformation predictions (Shi et al. 2020). Surface subsidence diagram and monitoring villages near mining areas in India using the modified PS-InSAR have good practical results for monitoring the slow settlement of villages with high accuracy (Kumar et al. 2020). However, the above results are applied to the calculation of subsidence monitoring of the mine surface, while the estimation and analysis of the surface incline, curvature, and other deformation values are lacking.
Results got by InSAR are the projection of line-of-sight deformation to vertical deformation. The surface subsidence value of the mine area is obtained by this method, which cannot be directly applied to the evaluation of the level of structural damage of the surface buildings. Building damage assessment around the Fengfeng mine area was mainly calculated using the results of Sentinel-1A images, and it computed the conversion formula between surface subsidence values and building deformation values for the evaluation of building damage (Diao et al. 2018). The assessment of building damage around the Panji coal mine was approximated using satellite image calculations combined with the field deployment of conductors and using the surface subsidence velocity and slope as evaluation factors to classify the surface building damage level (Li et al. 2019). Since most mining areas are far from cities, PS points are scarce, and the time-space decoherence of mining features is severe. The application of PS-InSAR technology is difficult (Maghsoudi et al. 2018). Xing increased the number of PS points by widely distributing artificial angle reflectors in the mine area and applied PS-InSAR technology to mine deformation monitoring (Xing et al. 2013). Therefore, remote sensing technology should reasonably extract "surface area" deformation data in mining areas and establish a building damage evaluation system based on surface deformation variables. The paper mentions the Cross-Reference Table of Surface Structural Damage Levels and Ground Deformation, used as a building damage level and evaluation index (Diao et al. 2018). The reference table is based on the actual calculated inclination value, curvature value, and horizontal movement value of the ground surface-the threshold value as a basis for judging damage to the building. This evaluation system is essential to analyze the damage to the mine buildings during the mining of the working face.
In this paper, with the help of Sentinel-1A, 21 images data from October 11, 2017, to June 22, 2018, SBAS-InSAR (Berardino et al. 2002Fan et al. 2018) was used to select the pixels with coherence greater than 0.18 near the 7221 Grout-filled working face of a mine in Huaibei Mining Group, Anhui Province. The line-of-sight deformation was projected onto the vertical deformation to obtain the time series deformation diagram of the 7221 working face in this period. Meanwhile, using PS-InSAR (Hooper et al. 2004(Hooper et al. , 2007Ferretti et al. 2000), PS points within the area of Gaochangying village, Houlou Gaojia village, and Industrial Square maintenance belt near the 7221 Grout-filled working face of the mine were selected to estimate the secular subsidence values of each PS point during the mining period. With the help of directional derivatives and geodesy, each pixel or PS point is assumed to be a ground deformation observation conductor point. The calculated diagram of the inclination and curvature values corresponding to each period of the time-series sedimentation diagram is obtained. Referring to the Cross-Reference Table of Surface Structural Damage Levels and Ground Deformation (Diao et al. 2018), building damage levels can be evaluated for villages and Industrial Square buildings in the study area. The active period of surface deformation and influence occurring during the mining time of the 7221 working face mined by grout filling technology is analyzed.

SBAS-InSAR methods
The SBAS-InSAR combines small baseline interference pairs by setting a specific time baseline threshold, and inverse performs the optimal deformation phase sequence according to the minimum parametric criterion to limit temporal and spatial decoherence. Suppose that N + 1 SAR images gained in the same region in the time series (t 0 , t 1 , …, t n ) have been acquired. After setting a specific time baseline threshold, M Interference pairs can be formed, and M satisfies: SBAS-InSAR requires aligning all images to the same super master image and generating M differential interferograms images with multi-look processing, topographic differencing, and phase unwrapping before performing a parametric inversion. Now, consider two SAR images acquired at moments t A and t B (assuming t A < t B ) with differential j-interferogram. Without considering errors such as atmospheric delay, the interferometric phase composition of the pixels with coordinates (x, r) in the azimuthal and distance directions is: where (t B , x, r) and (t A , x, r) are the phases of the corresponding SAR images at time t B and t A , d(t B , x, r) and d(t A , x, r) are the line-of-sight deformations at moments t B and t A to moment t 0 . To obtain a smooth and continuous solution, it is necessary to convert the phase in Eq.
(2) to the form of the average phase velocity multiplied by the time baseline, i.e., where v j (x, r) denotes the average phase velocity of the j-interferogram. Thus, the phase of the j-interferogram can be expressed as the integral of the phase velocity across the periods: So far, we can establish the observation equation connecting the phase velocity and M differential interference phases for each period of (t 0 , t 1 , …, t n ): where B is the M × N coefficient matrix containing information on how the interferometric pairs are composed, v is the N × 1 vector containing the phase velocity of each period, and is the M × 1 vector containing the differential interferometric phase. If B is rank-deficient, Singular Value Decomposition (SVD) can obtain the generalized inverse matrix of B and then, obtain the minimum parametric solution of the deformation velocity.

PS-InSAR methods
The basic idea of PS-InSAR is to screen out points that are only slightly affected by spatial and temporal decoherence in long-time series. The scattering characteristics of these points should remain stable. Ferretti first screened a series of PS candidate points by amplitude departure index DA (Hooper et al. 2004(Hooper et al. , 2007Ferretti et al. 2000). This index is constructed considering that when the phase stability is high, the dispersion of the amplitude can better reflect the distribution of the phase. After aligning the N + 1 view SAR images covering the same area to the same master image, the D A of a pixel is calculated as follows: where σ A and m A denote the standard deviation and average value of the time series amplitude of the image A pixel. The smaller the D A value represents, the better the amplitude stability of the pixel. The pixel with D A < 0.25 is usually selected as the PS candidate point with better phase stability.
After the PS candidate points are identified, the Delaunay triangle network is constructed using these points. Assume that the existing DEM has roughly removed the topographic phase, the relative differential i-interferogram between the PS point pairs (also called "arc segments") formed by the adjacent PS points x and y in the differential interferogram of view Δ i x,y,int can be expressed as: ,Δ i x,y,Δh and Δ i x,y,res denote the relative linear deformation phase, residual elevation phase, and other residual phases between the two points. The relative atmospheric delay, orbital error, nonlinear deformation, and other noise are mainly included in Δ i x,y,res . B i ⊥ ,t i denote the spatial vertical baseline and the temporal baseline of the differential i-interferogram. ,R i , denote the wavelength, sensor-to-target distance, and radar incidence angle.Δv x,y denotes the relative deformation velocity between two points.Δh x,y denotes the relative residual elevation between the two points.
Atmospheric phase and orbital errors are strongly correlated in space, and residual terrain phase is also spatially correlated, so modeling for PS point pairs can effectively weaken their effects. To solve for Δv x,y and Δh x,y in Eq. (7), the time domain coherence | | | x,y | | | needs to be maximized: where N denotes the total number of PS point pairs. The periodogram method is generally used to solve Eq. (8) by first giving approximate upper and lower bounds for Δv x,y and Δh x,y then iterative sampling. Since the approximate upper and lower bounds of Δv x,y and Δh x,y are required, it is essential to know the study area. After estimating Δv x,y and Δh x,y for each arc segment, a particular path integral can obtain the absolute deformation.
Using the low and high-frequency characteristics of atmospheric phase and orbital error in space and time, respectively, the atmospheric phase, orbital error, and nonlinear deformation in the residual phase can be separated by spatial and temporal filtering. Then, the atmospheric phase and orbital error on each pixel of each interferogram can be estimated by spatial interpolation. After evaluating and removing the atmospheric phase and orbit errors, the underestimation in analyzing the phase time domain stability has been mitigated. The final PS point can be selected based on the time domain coherence and amplitude dispersion simultaneously, and the inversion is performed again to complete the measurement of the deformation field finally.

Surface incline and curvature deformation calculation
The magnitude of surface incline and curvature deformation can be calculated based on the vertical deformation values of pixels and PS points. When considering the effect of horizontal movement, the effect of the obtained line of sight deformation on settlement can be converted to settlement perpendicular to the surface direction by the following equation: where DLOS denotes the line-of-sight deformation value, DV denotes the vertical deformation value, DE denotes the east-west deformation value, DN denotes the north-south deformation value, θ is the radar beam incidence angle, and α is the satellite heading angle.

inversion for surface incline with SBAS-InSAR
In the conventional ground deformation calculation of the measured data, the incline between two observation points is the ratio of the difference between the sinking of two points and the horizontal distance between two points after deformation. The incline is between two observation points, not on a particular observation point. By treating each pixel in the SBAS-InSAR deformation diagram as an approximate surface deformation observation line point, the vertical deformation value of the pixel is equivalent to the subsidence value measured by the ground observation point in the corresponding period. The distance between the centers of adjacent pixels in the SBAS-InSAR deformation diagram is equivalent to the distance between adjacent observation points. Therefore, like the calculation method of surface deformation of conventional measured data, the surface incline value between two pixels can be obtained by dividing the differencing vertical deformation of adjacent pixels by the distance between the center points of adjacent pixels. The incline is the first-order derivative of the subsidence to the horizontal distance from the surface. In this paper, the incline on a particular pixel should be calculated by an interval centered on that pixel. It is assumed that the x and y axes of the plane are oriented parallel to the north and east of the sedimentation map. The least-squares method is used to calculate the tilt value along the coordinate axis near any point (m,n) in the subsidence diagram. The direction along the x positive semi-axis is as an example: where i x (m,n) denotes the incline value of the point (m,n) along the positive direction of the x-axis, f x (m,n) is the first-order partial derivative in the direction of the x-axis at the point (m,n) of the subsidence diagram. R x is the distance between two adjacent pixels.
For the surface incline in any direction, the surface moving basin can be considered as a two-dimensional surface. The incline of the surface along the coordinate axis direction is the same size and opposite direction of the first-order partial derivatives of the sinking in the x and y-axis directions. From the definition of the first-order directional derivative, it can be similarly deduced that: where i(β,m,n) denotes the incline of the (m,n) point along the angle of β direction.
Maximum incline for a point on the surface and its direction: where i′(m,n) denotes the absolute value of the maximum incline of the point (m,n), which is called the maximum incline value, β′(m, n) is the maximum inclined azimuth of the point (m,n), "|*|" means modulo the vector " *".

inversion for surface curvature with SBAS-InSAR
Curvature is generated by the uneven surface incline and reflects curvature of the surface convexity. In mining subsidence observation, the surface curvature is the ratio of the difference in the incline of two adjacent line segments to the sum of the lengths of the two-line segments before deformation. The curvature is the second-order derivative of the subsidence to the horizontal distance from the surface. The curvature of a point on the ground in any direction can be calculated according to the following formula.
where K(β, m, n) denotes the curvature of the angle along β direction at the point (m,n), with positive values for convexity on the surface and negative values for depression on the surface. f xx (m, n), f yy (m, n), and f xy (m, n) are the second-order partial derivatives of the sink at the point (m,n) along the x-axis, the second-order partial derivatives of the sink at the point (m,n) along the y-axis and the second-order mixed partial derivatives of the sink at the point (m,n) along the x-and y-axes. The second-order partial derivative is similar to Eq. (10). For the maximum curvature and its direction at a point j on the surface, the solution is as follows: where β 1, β 2 , … β j are all values of β that allow Eq. (14) to hold, K′(m, n) is the maximum curvature at point j, β j is the orientation angle of the maximum curvature value at point j.

inversion for surface incline and curvature with PS-InSAR
The PS points' subsidence values are used to obtain the incline and curvature values of the ground surface in this paper. The PS points are equated to the geodetic surface deformation observation conductor points. The author uses the permutation principle to obtain the surface tilt between two PS points by traversing any pairs of PS points whose differential distances meet the requirements. The incline is the first-order derivative of the subsidence to the horizontal distance from the surface. In this paper, it is considered that the incline value between a particular pair of two PS points should be calculated by using the two PS points as an interval endpoint. As in Fig. 1, the two PS points A, C.
It is assumed that the four points A, B, C, D are the surface PS points extracted by PS-InSAR technology. Due to the influence of mining of the nearby coal seam workings, subsidence occurs at the four points A, B, C, D in space to the points A 1 , B 1 , C 1 , D 1 , and the specific subsidence value is W A , W B , W C , W D . The PS points A, B, C, D are assumed to be conductor points A, B, C, D , and the model diagram for calculating the incline's surface and curvature is shown in Fig. 1.
Calculate the surface incline value at two points A, C.
where i AC is the incline value of two points A, C . W A and W C are the subsidence value of two points A, C . L A 1 C 1 is the spatial distance of two points A 1 , C 1 . Translate A 1 C 1 upward AA 1 so that the point A 1 coincides with the point A and C 1 reaches the point C 2 . Since L CC 2 ≪ L AC 2 and L AC 2 =L A 1 C 1 , this can be approximately considered as L AC ≈ L C 1 A 1 .
The calculated incline value in the interval and the endpoints of the interval (the two PS points of the interval) are two observation conductor points and an observation interval in the deployed ground deformation observation conductor. This paper selects two observation intervals using the permutation principle, containing four observation lead points (four PS points). The two observation intervals are required to have an identical endpoint (i.e., two PS points overlap), thus forming a curvature value calculation interval, as shown in Fig. 1 for three points A, B, D . In mining subsidence observation, the surface curvature is the ratio of the difference in incline value of two adjacent line segments to the sum of the lengths of the two-line segments before deformation.
Calculate the surface curvature value at three points A, B, D.
where K A∼B∼D is the three-point A, B, D curvature value. i AB is the two-point A, B incline value, and i BD is the two-point B, D incline value, which can be calculated by Eq. (15). L AB , L BD are the values of horizontal distance between A, B , B, D.

Ground deformation for building damage level evaluation
The Cross-Reference 1 3 deformation index for evaluating the damage of buildings within the interval. The curvature value of the interval can be compared with the surface curvature threshold in Table 1 as an indicator to evaluate building damage within the interval. A diagram can also be made to observe the spatial distribution of PS points in the interval exceeding the specified threshold to obtain a diagram of the spatial distribution of building vulnerability areas.

Study area
The study area is in Guoyang

Analysis of SBAS-InSAR measurements
A total of 21 views of Sentinel-1A images were selected from October 11, 2017, to June 20, 2018, and the time threshold was set to 24 days to obtain the time-space baseline diagram shown in Fig. 3. The subsidence diagram calculated by SBAS-InSAR, as shown in Fig. 4a-g, shows that one month after mining the 7221 grouted working face, subsidence appeared on the surface south of the working face. Because of the minor effect of coal seam mining on surface subsidence, the subsidence area is only near the working face of 7221. From Fig. 4h-o, although the mining of 7211 working face has reached the active period, the speed of surface subsidence has accelerated, and the range of surface subsidence basin is expanding. Since the grouting activity of the grouting holes in the working face, the surface subsidence of the nearby villages is well controlled. The influence of coal seam mining on surface subsidence has still not reached the nearest Houlou Gaojia village. The overall Houlou Gaojia village has no apparent subsidence. From the last five subsidence diagrams in Fig. 4p-t, it can be found that the coal seam mining at the 7221 working face and the grouting activity of the grouting hole make the surface activity intense. A sizeable area of the surface subsidence phenomenon, the subsidence range, has entered the interior of Houlou Gaojia village. The village is affected by the mining at the 7221 working face. However, the analysis of the severity of damage to the village buildings cannot be directly obtained from the analysis in Fig. 4 Conclusion.

Analysis of PS-InSAR measurements
Again, the data were selected from October 11, 2017, to June 20, 2018, with 21 views of Sentinel-1A images to obtain the time-space baseline diagram shown in Fig. 5.  As it is sparsely distributed on the calculated PS points near the working face, the value available for reference is low. The focus is on selecting the PS points within the maintenance belt of Gaochangying village, Houlou Gaojia village, and the Industrial Square for analysis. The calculation shows that 257 PS points were obtained from Houlou Gaojia Village, 3,186 PS points were obtained from Gaochangying Village, and 1,696 PS points from Industrial Square. Figure 6 shows that the west side of Houlou Gaojia Village was subject to minor subsidence by mining the 7221 working face. The village center area showed slow surface changes, while the east side of the village showed strong subsidence and a trend of change from the village's east side to the center. Minor subsidence also occurred along the southwest side of Gaochangying village along the 7221 working face, with minor surface deformation as the grouting work proceeded. The surface on the northeast side of Gaochangying Village has been in the process of subsidence, which is sharp and has a tendency toward the village center area.

Multiple types of surface deformation based on SBAS-InSAR
With the help of directional derivative and geodesy, the subsidence diagram obtained by SBAS-InSAR was further calculated and analyzed to obtain the maximum incline value and maximum curvature value calculation diagram during the mining cycle of 7221 grouting working face, as shown in Figs. 7 and 8.
The area showing blank in Fig. 7 is the area of lower coherence, such as water bodies. Evaluate the threshold value of 3 mm/m for surface incline division between Class I and Class II buildings in the mine area. The buildings in Gaochangying, Houlou Gaojia, and the Industrial Square maintenance belt near the 7221 working face are within the prescribed threshold of Class I damage. To analyze in more detail the effect of  The images period selected for this paper contains the mining cycle of work face 7221, which was mined on December 6, 2017, until June 3, 2018, when the work face stopped mining. From the analysis of Fig. 7a-g, it is concluded that the working face 7221 was mined from south to north, and the grouting holes did not start grouting activities. The southwestern part of Houlou Gaojia Village, which is the closest to the working face during this period, was also not affected by mining, which is consistent with the observation of the subsidence diagram. In (h)-(m) in Fig. 7, the working face was mined for nearly two months, and the two grouting holes in the south of the working face started grouting for nearly one month. Some pixels with surface incline values exceeding 3 mm/m appeared around the surface. For the maintenance belt of Houlou Gaojia Village, closer to the working face, the incline value exceeding the threshold value specified for Class I damage still did not appear. It can be concluded from Fig. 7n-q that the 7221 working face continues to be mined to the north, and the grouting holes in the north are grouted. The incline value becomes more extensive in the western part of Houlou Gaojia village. The northeastern village maintenance belt of Gaochangying and Houlou Gaojia gradually reaches the incline value of secondary damage because of the influence of unknown workings mining, and the buildings should be monitored and protected with emphasis. From the analysis of Fig. 7r-t, it is concluded that the working face mining reached the last month, the three grouting holes in the north exceeded one month, and the surface movement changed dramatically. The southwestern part of Houlou Gaojia Village is in the dynamic range of surface movement. As it is affected by surface movement and deformation, the southwestern Houlou Gaojia Village buildings should be monitored and protected.
The blanks in Fig. 8 show low coherence and curvature values of -0.05 to 0.05 (× 10 -3 /m) pixels within the building protection line, which are not considered. Since the curvature takes positive or negative values to represent surface raise or depression, the images of curvature values calculated using SBAS-InSAR subsidence diagrams were divided into positive and negative pixel curvature values. Gaochangying, Houlou Gaojia, and Industrial Square are all within the threshold of direct damage. To analyze the surface deformation and building damage of overlying solo rock grout filling mining in more detail, the author will add three segments to the color banding of building curvature values within the first level of damage.
From Fig. 8a-l, we can see that only a part of the area north of Gaochangying village appears to have a high curvature value, which is far away from the mining of 7221 working face and not caused by the mining of 7221 working face. In Fig. 8m-t, the surface movement is more active in the west and north of Houlou Gaojia Village. The buildings in the maintenance belt of the village have not reached the first level of destruction, although the surface subsidence is caused by the mining of the 7221 working face. Due to the exploitation of working face 7221 and the grouting action of three grouting holes in the north, the surface subsidence and uplift occurred. Surface deformation easily damaged the buildings during this period. For the critical buildings in the southwestern part of the village, wires should be deployed to shorten the observation period and conduct high-precision and comprehensive observation, to strengthen the critical protection of the buildings in this area.

Multiple types of surface deformation based on PS-InSAR
Due to the large number and dense distribution of PS points within the village, it was impossible to conduct cumulative subsidence, and subsidence velocity analysis of PS points to use the data fully. This paper proposes a method based on the principle of permutation and combination, assuming PS points as surface deformation observation conductor points. Make comprehensive observation and calculation of all PS points' geographic location and time series subsidence values in the village, selecting pairs of PS points with a spatial distance of 3 m ~ 30 m to form the incline value calculation interval. Due to the small number of PS point pairs with calculated incline values exceeding the 3 mm/m threshold during the initial work face mining, the value available for analysis is low. Therefore, this paper selects 2017-10-11-2018-02-08, 2017-10-11-2018-03-04, 2017-10-11-2018-04-09, 2017-10-11-2018-05-03, 2017-10-11-2018-05-27, 2017-10-11-2018-06-20 for 6 periods of incline calculations of the time-space distribution of overrun PS point pairs for the diagram. Set the diagram out period of 24 days for intense surface deformation. The technology can still reduce the period of comprehensive observation and analysis to 12 days. The time-space distribution of the calculated incline value exceedance point pairs for PS within the village maintenance belt is shown in Fig. 9.
Analysis of Fig. 9 shows  Table 2 can be obtained. They correspond to the spatial distribution of PS points within the village maintenance belt to calculate the incline value exceedance point pairs Fig. 9a to f for six periods. The incline calculation interval composed of 854 pairs of PS points meeting the distance requirement is calculated in Houlou Gaojia Village. There are fewer inclination calculation intervals exceeding the limit, and the proportion of the inclination calculation interval exceeding the limit in all inclination calculation intervals does not exceed 0.23%. A total of 25,933 incline calculation intervals composed of PS point pairs satisfying the distance requirement were calculated in Gaochangying village. Less than 70 inclination calculation intervals exceeded the limit, and the proportion of the inclination calculation intervals that exceeded the limit in all inclination calculation intervals did not exceed 0.24%. Seventeen thousand four hundred one incline calculation intevrvals composed of PS point pairs satisfying the distance requirements were calculated for the Industrial Square. Less than 50 incline calculation intervals exceeded the limit, and the proportion of the incline calculation intervals that exceeded the limit in all incline calculation intervals did not exceed 0.28% (Table 2).

3
The results of the incline calculation interval that meet the requirements are used as known data to calculate the curvature values within the curvature calculation interval for three critical buildings in Houlou Gaojia Village, Gaochangying Village, and Industrial Square. The dates as the incline calculation interval were chosen to calculate the maximum curvature values for buildings in these three village maintenance belts. The PS point pairs at the endpoints of the curvature calculation interval exceeding 0.2 (× 10 -3 /m) were plotted to analyze the spatial distribution in time of the PS point pairs exceeding the calculated curvature values at the surface. The spatial distribution of the curvature calculation interval endpoints PS point pairs for the calculated curvature values exceeded within the maintenance belt of Houlou Gaojia Village is shown in Fig. 10. The spatial distribution of curvature value exceeded the PS point pairs diagram of the exact dates selected in Fig. 10. From the analysis of (a) to (f) in Fig. 10, the calculated curvature values in the curvature calculation interval of Houlou Gaojia Village exceeded the threshold in a small number and were concentrated near the village isolation zone in the village's northeast and near the village square. In the village's southwest, the nearest area to the 7221 working face does not have the calculated curvature value exceeding the threshold, which is less affected by the mining of the 7221 working face and has higher building stability.
The number of PS point pairs of curvature calculation interval endpoints in the maintenance belt of Houlou Gaojia Village is small. Analyze the proportion of the curvature calculation interval that exceeds the limit among all curvature calculation intervals in Houlou Gaojia Village. The total number of curvature calculation intervals in Houlou Gaojia Village is 7,255. The number of curvature calculation intervals exceeding the limit does not exceed 30, and the proportion does not exceed 0.39% (Table 3).
The spatial distribution of the curvature calculation interval endpoints PS point pairs for the calculated curvature value exceeded within the maintenance belt of Gaochangying village is shown in Fig. 11. The number of exceedances in the curvature calculation interval is large in Gaochangying village. However, the distribution is relatively concentrated in three areas: the square on the west side of the village, around the buildings on the northeast side of the village, and near the maintenance belt on the village's northwest side. The buildings are easily damaged in these three areas, so monitoring and protection should be strengthened. The spatial distribution of the curvature calculation interval of the exceedance in other areas within the village is relatively scattered. The reasons for the exceedance should be analyzed by field inspection against the spatial location distribution in the diagram.
To specifically analyze the proportion of the exceeded curvature calculation intervals in the total number of curvature calculation intervals within the maintenance belt of Gaochangying village, Table 4 was obtained. The total number of curvature calculation intervals in Gaochangying Village is 752088. The number of curvature calculation intervals exceeding the limit does not exceed 7000, and the proportion does not exceed 0.91%. Since a curvature calculation interval needs to be expressed with the help of three PS   To specifically analyze the proportion of curvature calculation intervals that exceed the limit in the total number of curvature calculation intervals in Industrial Square, Table 5 was made. The total number of curvature calculation intervals in industrial plazas is 417,551. The number of curvature calculation intervals exceeding the limit does not exceed 2,500, and the proportion does not exceed 0.59%.

Analysis of SBAS-InSAR surface subsidence monitoring
The SBAS-InSAR selects the spatial-temporal baseline threshold to remove the interferogram with poorer interference-effect compared with the D-InSAR. SBAS-InSAR can reject the pixels with coherence lower than the threshold value and ensure the accuracy of subsidence image metadata to analyze surface subsidence. Meanwhile, the SBAS-InSAR also improves the application of surface subsidence in mining areas for the problems of rapid surface subsidence and significant changes in the backscattering properties of the features. In this paper, SBAS-InSAR is chosen to analyze the dynamics of the surface subsidence problem caused by the mining of the 7221 grouted working face. The results obtained by SBAS-InSAR (Fig. 4) can accurately reflect the temporal and spatial extent of surface subsidence initiation and rapid subsidence in the mining area, which provides a reasonable data basis for analyzing surface dynamics. Then, with the help of the theory of directional derivatives, each pixel is assumed to be a "surface deformation observation line point" so that the maximum incline and curvature value of the mining area are calculated in pixel units. As shown in the figures (Figs. 7 and 8), the maximum incline value of the mine surface, the dynamic change of the curvature value, and the period and spatial distribution of the surface subsidence and uplift phenomenon occurred drastically during the mining process of 7221 grouting face. The method provides suggestions for the location and encryption of the observation period of the surface movement observation wires during the mining process at the working face, and the construction of the isolation zone for the key protection buildings in the mining area.
Low resolution caused by multi-look average filter leads to the edge length of 12.75 m for each pixel in the subsidence diagram calculated by SBAS-InSAR. Due to the distance and direction of the satellite, the azimuthal spatial resolution is low, which  point. Therefore, the subsidence value obtained is smaller than the calculated value observed by the actual deployment of the surface movement observation wire, as can be obtained from the analysis in Fig. 13. For the surface beneath the village's deformation monitoring analysis with high coherence and many PS points, the inversion analysis method based on the PS point pair incline value calculation interval and the curvature value calculation interval proposed in this paper can be selected.

Analysis of PS-InSAR mining surface monitoring
This paper proposes that the PS points in the mine area are assumed to be surface deformation observation conductor points so that the theoretical model of incline value and curvature value calculation and analysis can be established. The following conclusions were obtained from the monitoring and analysis of the buildings' surface within the village's maintenance belt near the 7221 grouting work face: (1) The PS points in the maintenance belt of the mining village are assumed to be surface deformation observation conductor points. The distance requirement of adjacent surface deformation observation conductor points is set so that many inclines and curvature calculation intervals can be obtained and meet the requirements. Considerable data damage evaluation of buildings in mining areas and spatial distribution analysis of building vulnerable areas can effectively provide critical building damage evaluation and protection suggestions.
(2) The incline calculation interval and curvature calculation interval obtained by this method are not directional, and the resultant data are redundant, so how to use the data set more fully and rationally is a current problem that needs to be solved. (3) Since the principle of PS-InSAR is to select features with high coherence for subsidence analysis, the solution for applying PS-InSAR to the surface near the mine working face is to deploy artificial angle reflectors, which requires a lot of human and material resources. Otherwise, the number of PS points near the mine working face cannot meet the calculation requirements. Therefore, in this paper, only the PS points within the village maintenance belt are selected for analysis, and the surface PS points near the working face cannot be analyzed.

Conclusion
In this paper, 21 scenes of Sentinel-1A images were selected as the data source, and the starting and ending times included the mining time of 7221 grout working face in a mine of Huaibei Mining Group Anhui Province. SBAS-InSAR and PS-InSAR were used to analyze and calculate the data, focusing on deformation monitoring and building damage analysis of the ground surface within the village maintenance belt of the mining area. The specific conclusions are: From the time series subsidence diagram (Fig. 4), it can be found that the surface subsidence affected by the Grout-filled mining technology is reduced because of the grouting activity during the working face mining process. The active period when the surface activity is intense appears later but lasts longer. The mining of 7221 grouting face has less impact on the nearest Houlou Gaojia village. In contrast, the mining of 3233 working face of a mine west of Gaochangying village affected the severe surface deformation of Gaochangying village and Industrial Square.
The subsidence results calculated by SBAS-InSAR are used as the data source. The pixels are assumed to be surface deformation observation conductor points for calculating the maximum incline and maximum curvature value of the mine surface. It can be dynamically observed from the diagram that, due to the mining of 7221 grouting face, the surface deformation extends outward from the working face area to the inside of the village maintenance belt of Houlou Gaojia Village. The calculation results of this theoretical model can also reflect the spatial distribution of active surface deformation areas in the village more accurately and provide suggestions for monitoring and protecting the building deformation in this area.
The subsidence results of PS points calculated by PS-InSAR are used as the data source. The PS points are hypothetically used as surface deformation observation conductor points to conduct the calculation interval of surface incline and curvature calculation interval analysis. It was concluded that the proportion of PS point pairs exceeding Class I damage within the maintenance belt of Houlou Gaojia Village, Gaochangying Village, and Industrial Square around the mine area were all less than 1% of the total calculated results. The few areas that exceed the threshold requirements of the I level of damage have a regular spatial distribution. The buildings are very little affected by the mining movement of the working face during the monitoring cycle. The results can lay conductors in intense surface movement and buildings easily damaged and carry out encrypted observation. Guidance can be provided for monitoring surface time series deformation of buildings in the maintenance belts of villages in mining areas and the critical protection of buildings.