Retrieving Three-Dimensional Large Surface Displacements in Coal Mining Areas by Combining SAR Pixel Offset Measurements with an Improved Mining Subsidence Model

: Interferometric synthetic aperture radar (InSAR) technology can obtain one-dimensional surface displacements in the radar line of sight (LOS). In the ﬁeld of mining subsidence, large 3D movements often occur at the same time, and hence InSAR derived one-dimensional LOS displacements can hardly reﬂect the actual surface motion in mining areas. To realize the monitoring of three-dimensional large surface displacements in mining areas, a method for monitoring three-dimensional large surface displacements in mining areas that combines SAR pixel offset tracking (OT) and an improved mining subsidence model is proposed in this article. First, a new functional relationship between surface subsidence and horizontal movement combined with the characteristics of the overburden rock stress and the deformation characteristics of the fractured rock mass in coal mining areas is established. Then, a three-dimensional surface deformation model is established based on the proposed relationship between surface subsidence and horizontal movement and the radar projection equation, and ﬁnally, the optimal parameters of the deformation model are inverted iteratively using LOS deformation results obtained by OT method to retrieve the three-dimensional large displacements of the surface. The signiﬁcant advantage of the method proposed in this article is that it can accurately acquire the 3D large surface displacements using only two SAR amplitude images with the same imaging geometry. To verify the accuracy and reliability of the proposed algorithm, two scenes of high-resolution spotlight TerraSAR-X images are used in this paper to conduct a three-dimensional surface displacement monitoring experiment on a working panel in the Daliuta mining area in Shaanxi Province, China, based on the proposed method. Experimental monitoring results show that the maximum surface subsidence is approximately 4.5 m, and the maximum horizontal movements in the strike and dip directions are approximately 1.4 m and 1.2 m, respectively. Using GPS measurements to verify the monitoring results, the root mean square error (RMSE) of vertical subsidence is 6.8 cm, and the RMSE of horizontal movement is 7.1 cm. Compared with those in the original mining subsidence model, the accuracies of vertical subsidence and horizontal movement in the proposed model are increased by 28.2% and 37.5%, respectively, which proves the reliability and accuracy of the proposed method.


Introduction
Coal has been a basic energy source and an important industrial raw material for all countries in the world, and it has strongly supported the development of the world economy. However, large-scale and high-intensity coal mining often causes severe and large-gradient deformation on the Earth's surface during a short period, which in turn not only causes the ground surface to sink but also could result in damage to the buildings, infrastructure and roads above the goaf. Monitoring of coal mining induced large-gradient three-dimensional deformation on the Earth's surface is of great significance to accurately assess the impacts of coal mining activities, adjust mining plans in mining areas, and predict potential geological disasters during mining.
Traditional monitoring techniques for mining subsidence include triangulation, leveling and GPS, but they have disadvantages such as high cost, heavy workload, and sparse monitoring point density [1,2]. Interferometric synthetic aperture radar (InSAR) has been proven to be a powerful tool for monitoring surface displacements [3,4] with the advantages of all-time, all-weather, wide coverage (tens of kilometers to hundreds of kilometers), high precision (centimeters to millimeters) and high spatial resolution (tens of meters to submeters). It overcomes the shortcomings of the above mentioned traditional monitoring techniques, and has been successfully used to monitor mining subsidence [5][6][7][8][9]. Coal mining often causes large-gradient surface displacements during a short period, and the surface displacement gradient can exceed the maximum detectable deformation gradient of InSAR [10,11], which impedes InSAR to determine the correct surface displacements. Yang et al. [12] proposed using the interpolation multiview method to increase the detectable deformation gradient, and Zhao et al. [13] used the full resolution method of a single interferogram to increase the detectable deformation gradient. However, the extent for these two methods to increase the detectable deformation gradient is limited and the problem how to use InSAR to obtain large-gradient surface displacements remains unsolved. To obtain large-gradient surface displacements in mining areas, Zhao et al. [14] used SAR pixel offset tracking (OT) method to obtain surface subsidence signals in a mining area with an accuracy of approximately 1/25 pixels. Huang et al. [15] managed to recover the surface displacements in the Daliuta mining area in Shaanxi Province, China; the maximum LOS displacement was found to be approximately 3.6 m, and the RMSEs in the strike and dip directions were 0.143 m and 0.108 m, respectively. Note that this method obtained vertical displacements only with an assumption of zero horizontal movement.
Conventional InSAR can only provide one-dimensional surface displacements in the direction of the radar's line of sight (LOS) in the absence of multi-tracks of SAR images, but the surface displacements caused by underground mining often occur in both vertical and horizontal directions, and hence InSAR LOS results cannot fully reflect the actual surface motion in mining areas, which may even cause misinterpretation in some cases. To derive three-dimensional (3D) surface displacements from InSAR observations, multi-tracks of SAR images and/or a dense GNSS network are often required [16][17][18][19][20][21][22]. On the other hand, the existing mining subsidence theory [23] reveals that there is a functional relationship between vertical and horizontal movements caused by coal mining. Therefore, it can be used to determine 3D surface displacements from 1D LOS displacements. Fan et al. [24] and Yang et al. [25] used OT method to obtain large-gradient surface displacements in the LOS direction and then combined them with a mining subsidence model to obtain largegradient 3D surface displacements. It should be noted that the existing mining subsidence model has two shortcomings: (i) it is assumed that the horizontal movement is proportional to the tilt in the existing mining subsidence model, indicating that horizontal displacements are equal to zero in the critical and supercritical areas, but that assumption does not always hold true [26][27][28]; (ii) the predicted curves of vertical and horizontal displacements from the existing mining subsidence model converge too fast at the boundaries; that is, the prediction errors of the existing mining subsidence model are greater near the goaf boundaries than those away from the boundaries [29]. These two deficiencies lead to low accuracy of the 3D surface displacements through the fusion of the existing mining subsidence model and the OT method. Chen et al. [30] presented a method that integrates OT and clastic medium theory to resolve 3D large gradient displacements caused by underground coal mining. This method solves the above mentioned problem (i), but the problem (ii) has not been solved yet.
To solve the abovementioned problems and obtain a higher accuracy, in this paper, an improved mining subsidence model is proposed based on previous study [30]. It is combined with the OT method to retrieve precise three-dimensional surface displacements in mining areas, with a working panel in the Daliuta mining area in Shaanxi Province, China, being taken as an example.

Mining Subsidence in the Original Mining Subsidence Model
The mining subsidence model based on the random medium theory is one of the most widely used mining subsidence prediction methods. In this theory, the main rock mass of the subsidence area is considered as a discontinuous medium, and its medium model is similar to the granular medium model. It is assumed that, although the rock is deformed under mining, the total volume remains unchanged [23]. As shown in Figure 1a, we assume that xOy is the ground coordinate system and vO 1 u is the coal seam coordinate system. According to the principle of mining subsidence in the existing mining subsidence model, in horizontal or near-horizontal coal seam mining unit B(u, v), the vertical and horizontal displacements of any point A(x, y) on the ground surface can be written as follows [23]: where W e (x, y) is the surface subsidence caused by mining unit B(u, v), r is the main influence radius, r = H/tanβ. b, H and tanβ represent the coefficient of the horizontal movement, the mining depth of the extracted unit element, and the tangent of the major influential angle, respectively. U e (x, y) and I e (x, y) are the unit horizontal movement and tilt deformation caused by exploiting unit B(u, v). B is the proportionality coefficient between U e (x, y) and I e (x, y).
Equation (1) is established according to the theory of granular medium mechanics in discontinuous medium [23]. In this theory, the granular medium is composed of medium particles similar to sand particles or relatively small rock masses, and the particles have completely lost contact with each other and can move around. The movement of a granular medium is characterized by the random movement of particles, and the movement of a large number of granular media is regarded as a random process [23]. As shown in Figure 1b, in the flat bottom of the mobile basin, the stress of the mining overburden along the bedding direction is relatively small, which is basically in line with the theoretical assumption of the existing mining subsidence model. Above the coal pillars in the goaf, the rock mass is subject to strong tensile stress directed towards the center of the goaf. Under the action of tensile stress, the mining overburden above the coal pillar produces horizontal displacements towards the center of the mined-out area and sink in the vertical direction, which leads to the deviations between the predicted surface subsidence and the measured values. In addition, if there is a certain thickness of loose layers (such as alluvial layers and residual layers) overlying the coal seam, cracks are prone to appear under the action of tensile stress due to the low compressive and tensile strength of the loose layers. A more significant plastic flow towards the center of the mined-out area can be generated, which increases the amount of subsidence above the coal pillar and make it more difficult to predict at the boundaries. In addition, in mining areas with high phreatic water levels, the surface subsidence in the mined-out area causes the groundwater level to rise, and the collapsibility of loess can also result in prediction biases at the boundaries. Therefore, the subsidence curve predicted by the existing mining subsidence model (as Equation (1) converges too fast at the edges (as shown in Figure 1b) and has a large deviation from the actual subsidence curve. subsidence curve predicted by the existing mining subsidence model (as Equation (1)) converges too fast at the edges (as shown in Figure 1b) and has a large deviation from the actual subsidence curve. Equation (2) shows the proportionality between horizontal movement and tilt according to the stochastic medium theory model (SMTM), which holds true only in a limited number of cases. This is because the SMTM assumes that the total volume of the rock mass remains unchanged with coal mining but the assumption of the SMTM is inconsistent with the actual movement of the rock mass. The actual rock mass is anisotropic and discontinuous (with cracks), and thus its volume varies during mining due to the bulk expansion and rheology of the rock mass; the rock mass should be considered as a clastic medium [30].

Modeling 3D Surface Displacements in Coal Mining Areas with Improved Mining Subsidence Model
In this paper, phenomenological theory [31] is invoked to address the problem of converging too fast at the edges in Equation (1), which refers to the laws obtained by summarizing experimental results instead of using its internal mechanism. Phenomenological theory is the generalization and refinement of experimental phenomena and it has the function of describing and predicting physical phenomena, but no explanatory function. W(x) is the subsidence curve under semi-infinite mining condition, i.e., when time t → ∞ and in x > 0, all the coal seam is exploited and in x < 0, the unmined coal seam remains. Surface deformation is mainly concentrated in the range of 2r above the mining boundary, this is the origin of parameter r; β is the angle between the horizontal and the line connecting the boundary point of the major influence range (i.e., −r or r) and the mining boundary point. (e) Plane-coordinate system for coal mining. Note that (i) B(u, v) represents one unit element, A(x, y) represents one surface point influenced by exploiting unit element B(u, v), D 1 and D 3 are the extraction length and width of the working panel, respectively. (ii) parameters r, β, and H represent the major influence radius, major influence angle and mining depth. (iii) W i e is the vertical displacement caused by exploiting unit element, parameters θ and α represent the propagation angle of the coal mine and coal seam inclination, respectively. s 1 and s 3 are the offsets of the inflection points at the lower and upper boundaries of the goaf in the main section. Equation (2) shows the proportionality between horizontal movement and tilt according to the stochastic medium theory model (SMTM), which holds true only in a limited number of cases. This is because the SMTM assumes that the total volume of the rock mass remains unchanged with coal mining but the assumption of the SMTM is inconsistent with the actual movement of the rock mass. The actual rock mass is anisotropic and discontinuous (with cracks), and thus its volume varies during mining due to the bulk expansion and rheology of the rock mass; the rock mass should be considered as a clastic medium [30].

Modeling 3D Surface Displacements in Coal Mining Areas with Improved Mining Subsidence Model
In this paper, phenomenological theory [31] is invoked to address the problem of converging too fast at the edges in Equation (1), which refers to the laws obtained by summarizing experimental results instead of using its internal mechanism. Phenomenological theory is the generalization and refinement of experimental phenomena and it has the function of describing and predicting physical phenomena, but no explanatory function. This method can avoid the exploration of the complex mechanical properties and mechanical processes of the internal movements and deformation of the rock layer, and an appropriate fitting can be obtained as long as the appropriate parameters are determined. According to the phenomenological theory [31], the expression form of the surface unit sinking basin can be changed to a combination of two unit sinking basins with different main influence radii according to a certain proportion, and then the surface subsidence prediction formula is established according to the superposition method, which can be expressed as [23]: where W i e (x, y) is the subsidence caused by mining unit B(u, v) that is obtained by the improved mining subsidence model, and W 1 e (x, y) and W 2 e (x, y) represent the subsidence of the unit corresponding to the changed main influence radius r 1 and r 2 , respectively. ρ is the scale factor, r 1 and r 2 represent the main influence radius of the improved mining subsidence model, where r 1 = H/tanβ 1 and r 2 = H/tanβ 2 .
According to the superposition principle, the surface point A(x, y) caused by the mining of the entire working panel can be expressed as: where m is the coal seam thickness, q a is the sink factor of improved model, α is the coal seam dip, l = D 3 −s 4 , L = (D 1 −s 2 )cosα, and D 3 and D 1 represent the lengths of the strike and dip directions, respectively. s 1 , s 2 , s 3 and s 4 are the offsets of the inflection points at the lower and upper boundaries of the goaf in the main section and the left and right boundaries of the main section, respectively. According to Chen et al. [30], in the x-z plane, the horizontal movement of the surface caused by unit mining can be expressed as: where θ is the propagation angle of coal mining. Combining Equation (6) with Equation (2), the horizontal movements of point A on the surface can be written as: where U EW (x, y) and U NS (x, y) represent horizontal movements in the easting and northing directions, respectively; and I EW (x, y) and I NS (x, y) represent the east-west and north-south tilts, respectively. Since tilt is the first derivative of the vertical displacements [23], the corresponding tilts (I EW (x, y) and I NS (x, y)) of point A(x, y) in the easting and northing directions can be expressed as: where ∆x and ∆y are the distances between adjacent points in the easting and northing directions, respectively. In Equations (5)-(8), parameters H, D 1 , D 3 , ∆x and ∆y are constants of the coal mining area. Therefore, 3D surface displacements can be derived according to Equations (5) and (7) when parameters q a , tanβ, θ, b a , s and ρ are known. Note that the tangent of the major influential angle in the improved mining subsidence model is changed to two, so tanβ needs to be determined by Equation (9): Remote Sens. 2021, 13, 2541 6 of 19

Estimation of Model Parameters using OT Results
From Section 2.2, it is clear that three-dimensional movements can be retrieved only if the parameters q a , tanβ, θ, b, s and ρ are determined. Here, an iteration method is proposed to determine the parameters q a , tanβ, θ, b, s and ρ using OT results as follows: (1) Obtain the initial vertical subsidence. Assuming that there is no horizontal movement in the initial settlement field, the initial vertical subsidence can be determined. The vertical subsidence of any pixel (i, j) can be expressed as: where ϕ represents the incident angle of the radar, W (i, j) is the vertical subsidence of pixel (i, j), and LOS(i, j) is the LOS-direction deformation of pixel (i, j) obtained by OT method.
(2) Determine the initial parameters of the model based on the genetic algorithm (GA) and fruit fly optimization algorithm (FOA). The increased number of the parameters in the improved mining subsidence model may lead to a decrease in the inversion accuracy of the existing intelligent optimization algorithm. Therefore, the GA+FOA integrated method is used to invert the parameters q a , tanβ, θ, b a , s and ρ. The GA is a random search algorithm that can simulate natural evolution processes. The GA is derived from the ideas of "genetic selection" and "natural selection" [32]. The main function of the GA is to compare the optimization problem to a population, encode the parameters to be sought as chromosomes, perform the steps of selection, crossover, and mutation in an iterative manner, and finally generate chromosomes that meet the target requirements. The GA has high versatility and strong robustness and is suitable for solving multi-parameter models; its disadvantage is that the local optimization ability is insufficient. The fruit fly has a keen sense of smell and the algorithm tries to simulate this [33]. The FOA has a very keen sense of smell and vision and can find the source of a scent tens of kilometers away [34]. When it is close to the source of a scent, its keen vision can find the location of food and companions and then it quickly flies towards the target [35]. The advantage of the FOA is that the optimization accuracy is high and the algorithm is easy to be implemented by computer programming. Its disadvantage is that when the step length is set too large, the local optimization ability decreases, and if the setting is too small, the global optimization ability is poor. The GA+FOA integrated method combines the excellent global optimization ability of the GA and the good local optimization ability of the FOA. It can also achieve higher inversion accuracy for the improved mining subsidence model with an increased number of parameters. Based on the principle of the minimum RMSE between the predicted and true values (here, the true values are GPS displacements), the output parameters are the initial optimal parameters.
(3) Calculate the initial three-dimensional surface displacements. Based on the initial optimal parameters output in Step (2), the three-dimensional initial deformation of the ground surface is estimated according to Equations (5) and (7).
(4) Obtain the new LOS deformation field. Based on the three-dimensional surface displacements obtained in Step (3), the following equation can be used to obtain the simulated LOS deformation LOS'(i, j) [36]: where δ represents the heading angle of the radar and U EW (i, j) and U NS (i, j) represent the horizontal movement of pixel (i, j) in the easting and northing directions, respectively. Now, we can calculate the absolute residual of the LOS-direction deformation, which can be expressed as Threshold judgment and loop iteration. We set a threshold ε for the absolute residual ∆LOS(i, j), here set to 10 mm. (i) If ∆LOS(i, j) > ε, return to Steps (1)-(4) to start a new iteration; (ii) if ∆LOS(i, j) ≤ ε, then output the current parameters as the optimal parameters of the model. Set the maximum number of iterations to 100.
(6) Output the optimal parameters of the model and obtain the final three-dimensional surface displacement field. Repeat Steps (1) (use the LOS' as new input for the next iteration) through (5) until the closure error ∆LOS(i, j) < ε, thereby retrieving the parameters q a , tanβ, θ, b a , s and ρ and deriving the time-varying three-dimensional deformation fields from Equations (5) and (8).
Since this algorithm is established by combining OT and improved random medium theory (as introduced in Section 2.2), this algorithm is referred to as OT-IRMT. The flow chart of the OT-IRMT model is shown in Figure 2. The IRMT model (yellow box in Figure 2) is incorporated and utilized to analytically separate the vertical and horizontal components (Equation (7) in Section 2.2). The OT results (LOS dataset in Figure 2) are used to estimate the parameters of IRMT model. The significant advantage of the OT-IRMT is that it can accurately acquire the 3D large surface displacements using only two SAR amplitude images with the same imaging geometry. new iteration; (ii) if ∆ ( , ) ≤ ε, then output the current parameters as the optimal parameters of the model. Set the maximum number of iterations to 100. (6) Output the optimal parameters of the model and obtain the final three-dimensional surface displacement field. Repeat Steps (1) (use the LOS' as new input for the next iteration) through (5) until the closure error ΔLOS(i, j) < ε, thereby retrieving the parameters qa, tanβ, θ, ba, s and ρ and deriving the time-varying three-dimensional deformation fields from equations (5) and (8).
Since this algorithm is established by combining OT and improved random medium theory (as introduced in Section 2.2), this algorithm is referred to as OT-IRMT. The flow chart of the OT-IRMT model is shown in Figure 2. The IRMT model (yellow box in Figure  2) is incorporated and utilized to analytically separate the vertical and horizontal components (Equation (7) in Section 2.2). The OT results (LOS dataset in Figure 2) are used to estimate the parameters of IRMT model. The significant advantage of the OT-IRMT is that it can accurately acquire the 3D large surface displacements using only two SAR amplitude images with the same imaging geometry.

Study Area and Data Used in this Study
The study area is located in the Daliuta coal mining area, Shaanxi Province, China. The study area has thick loose layers, shallow coal seams, thin bedrock, and a small mining depth to mining thickness ratio, and the thickness of the coal seam is approximately 6.94 m. The average mining depth in the mining area is 235 m, and the inclination of the coal seam is 1°~3°, which is a near-horizontal coal seam. The strike length of the working panel is approximately 4547 m, the dip length is approximately 300 m and the strike is super-critical mining. The small mining depth to thickness ratio and the rapid advancement of the working panel make the ground movement in this area large. Ground movement due to coal mining has unique features, including the fast propagation of coal mining influences, discontinuous surface displacements, and a large probability of catastrophes such as stepped collapses, hillside slip and ground fissures (shown in figures 3b and

Study Area and Data Used in this Study
The study area is located in the Daliuta coal mining area, Shaanxi Province, China. The study area has thick loose layers, shallow coal seams, thin bedrock, and a small mining depth to mining thickness ratio, and the thickness of the coal seam is approximately 6.94 m. The average mining depth in the mining area is 235 m, and the inclination of the coal seam is 1 •~3• , which is a near-horizontal coal seam. The strike length of the working panel is approximately 4547 m, the dip length is approximately 300 m and the strike is supercritical mining. The small mining depth to thickness ratio and the rapid advancement of the working panel make the ground movement in this area large. Ground movement due to coal mining has unique features, including the fast propagation of coal mining influences, discontinuous surface displacements, and a large probability of catastrophes such as stepped collapses, hillside slip and ground fissures (shown in Figure 3b,c). A single reference real-time kinematic (RTK) GPS was used to collect the coordinates of 71 ground stations along the strike (direction of the line formed by the intersection of the coal seam with a horizontal plane) and dip (direction of the line formed by the intersection of the coal seam with a vertical plane perpendicular to the strike of the feature) directions in the coal mine working panel (shown in Figure 3a).   Table 1, and the schematic diagram of the mining panel is shown in Figure 3a.  Table 1, and the schematic diagram of the mining panel is shown in Figure 3a.  [37]. The search window in the OT computation process was set to 64 × 64 pixels, and the oversampling factor to 2. The OT-derived surface displacements in the LOS direction were then employed to model 3D displacements with the method presented in Section 2.2. According to the algorithm flow shown in Section 2.3, 71 points (these points are colocated with the GPS measurements) along the strike line A-A' and dip line B-B' were extracted for parameter inversion for the improved mining subsidence model.
The GA+FOA integrated parameter method was selected to invert the model parameters: (i) set the initial range of parameters (as shown in Table 2); (ii) use the GA to perform parameter inversions 20 times independently; (iii) take the maximum and minimum values of each parameter inversion as the new optimization range; (iv) use the FOA to perform parameter inversions 20 times independently; and (v) the average values are output as the optimal parameters of the model. Table 2 shows the initial parameter ranges and the inversion results, whilst Figure 4 shows the distribution of parameter inversion values. Based on the optimal inversion parameters of proposed model and Equations (5) to (8), the vertical subsidence and the horizontal movement are obtained (as shown in Figure 5).  Figure 5e,f shows that the maximum horizontal movement in the strike direction was greater than that in the dip direction and that the horizontal movements in the dip direction were nearly symmetric. In addition, Figure 5e,f show that the horizontal movements were slight near the subsidence center, the maximum horizontal movement occurred at the edges of the mining area, and the horizontal movements decreased outside of the mining area.

Performance Assessment
To assess the performance of the OT-IRMT model, we compare the 3D surface displacements with GPS derived displacements over the 71 GPS points on the working panel ( Figure 3). Figure 6a,c show that the vertical displacements predicted by the OT-IRMT model were closer to GPS derived vertical displacements than those predicted by the original mining subsidence model, especially for the points at edges, e.g., Points 6-11 and 43-46 in Figure 6a and Points 5-8 and 19-24 in Figure 6c.
It is clear in Figure 6b that the horizontal displacements along Profile A-A' were equal to zero in the area of critical or super-critical extraction for the original mining subsidence model, whilst the horizontal displacements predicted by the OT-IRMT model agreed well with GPS ones. For Profile B-B', the OT-IRMT model was also in a better agreement with GPS measurements at the edges than the original model (Points 1-8 and 20-25 in Figure  6d). In terms of the RMSE of the differences between the predicted displacements and GPS-derived ones, the accuracies of vertical and horizontal displacements obtained by the OT-IRMT model increased by 28.2% and 37.5%, respectively, compared with those obtained by the original model. Figure 6e shows that the OT-IRMT derived vertical displacements were in good agreement with GPS measurements, with an RMSE of 6.8 cm (equivalent to 1.7% of the maximum subsidence). Figure 6f shows that the RMSE of the differences between the OT-IRMT derived horizontal movement and GPS measurements was 7.1 cm, which was 6.7% of the maximum horizontal displacement.
Chen et al. [30] and Yang et al. [38] combined the OT method and the original mining subsidence model to monitor the three-dimensional deformation in a coal mining area. The method proposed by Yang et al. [38] is based on the assumptions that (i) the horizontal movement at the edge of the mining area is zero and (ii) the horizontal movement is proportional to the tilting deformation, and the three-dimensional deformation model is solved by back-substitution method. However, as described in Section 2.1, the hypothetical model may not be consistent with the actual ground movements, and if there were uncertainties in the initial LOS-direction deformation, the back-substitution method would accumulate errors and thus lead to poor performance. The method proposed by Chen et al. [30] was mainly to address the problem of zero horizontal movements in the area of critical or super-critical extraction, but the issue of the fast convergence at the edges remained unsolved. Table 3 shows the comparisons among the offset tracking-clastic medium theory (OT-CMT) model proposed by Chen et al. [30], the OT-IRMT model proposed in this paper, and the alternative offset tracking-single amplitude pair (AOT-SAP) model proposed by Yang et al. [38]. Compared with the OT-CMT method, the OT-IRMT model improved 6.8% and 18.4% for the vertical and horizontal displacements, respectively. Compared with the AOT-SAP model, the OT-IRMT model improved 60.4% and 51.7% in the vertical and horizontal directions, respectively. The comparisons suggest the superiority of the OT-IRMT model. Note that the improvement of the OT-IRMT model compared with the OT-CMT method occurred mainly at the edges of the mining area, and hence the improvement was marginal. In contrast, the improvement of the OT-IRMT model compared with the AOT-SAP model was much greater, this was because the assumption of the AOT-SAP model that the horizontal movements and the tilt deformation are proportional during the modeling process may not hold true in the coal mining.
To assess the performance of the OT-IRMT model, we compare the 3D surface displacements with GPS derived displacements over the 71 GPS points on the working panel ( Figure 3). Figure 6a,c show that the vertical displacements predicted by the OT-IRMT model were closer to GPS derived vertical displacements than those predicted by the original mining subsidence model, especially for the points at edges, e.g., Points 6-11 and 43-46 in Figure 6a and   Figure 6c.

Discussion
The inversion accuracy of the OT-IRMT model parameters directly determines the accuracy of the OT-IRMT derived 3D surface displacements, and hence an orthogonal experiment was performed to examine the sensitivity of model parameters. The orthogonal experiment method is an experimental method to study the degree of the influence of multiple factors on a model [39]. This factor-equilibrium method can use a limited number of experiments to analyze the sensitivity of the model to the parameters. In this paper, the OT-IRMT model parameters, i.e., q a , tanβ 1 , tanβ 2 , θ, s and ρ were selected as the factors that influence the model accuracy. Table 4 shows the settings of the model parameters in the orthogonal experiment. Taking the three-dimensional directions of the ground motion, namely, the vertical, easting and northing directions, the differences between the measured and estimated values were used as indicators in the orthogonal experiment to analyze the sensitivity of the model to the parameters. Since the horizontal movement coefficient b only affects the horizontal motion, it was not included in the experiment.
Taking the geological mining conditions of the mining area as an example, the thickness of the coal seam at the simulated working panel was 5 m, the mining depth was 200 m and the coal seam inclination angle was 0 • . The lengths of the strike direction D 3 and the dip direction D 1 were 400 m and 200 m, respectively. The overlying rock lithology was medium hard. Seventy observation points were arranged at 10 m intervals along the strike and dip directions (as shown in Figure 7).

Discussion
The inversion accuracy of the OT-IRMT model parameters directly determines the accuracy of the OT-IRMT derived 3D surface displacements, and hence an orthogonal experiment was performed to examine the sensitivity of model parameters. The orthogonal experiment method is an experimental method to study the degree of the influence of multiple factors on a model [39]. This factor-equilibrium method can use a limited number of experiments to analyze the sensitivity of the model to the parameters. In this paper, the OT-IRMT model parameters, i.e., qa, tanβ1, tanβ2, θ, s and ρ were selected as the factors that influence the model accuracy. Table 4 shows the settings of the model parameters in the orthogonal experiment. Taking the three-dimensional directions of the ground motion, namely, the vertical, easting and northing directions, the differences between the measured and estimated values were used as indicators in the orthogonal experiment to analyze the sensitivity of the model to the parameters. Since the horizontal movement coefficient b only affects the horizontal motion, it was not included in the experiment.
Taking the geological mining conditions of the mining area as an example, the thickness of the coal seam at the simulated working panel was 5 m, the mining depth was 200 m and the coal seam inclination angle was 0°. The lengths of the strike direction D3 and the dip direction D1 were 400 m and 200 m, respectively. The overlying rock lithology was medium hard. Seventy observation points were arranged at 10 m intervals along the strike and dip directions (as shown in Figure 7).  The orthogonal table L 27 3 6 was calculated with 6 factors and 3 levels, as shown in Table 5. After the 27 tests shown in Table 5, we obtained the analysis results of the orthogonal test (shown in Table 6). Through multifactor analysis of variance, the coefficient of determination R 2 of the vertical monitoring accuracy in this experiment was 0.813. This means that the independent variables, namely, the parameters q a , tanβ 1 , tanβ 2 , θ, s and ρ, can explain 81.3% of the change in the dependent variables. In other words, if the independent variables remained unchanged, the degree of the change in the dependent variables would be reduced by 81.3%. Parameters q a , tanβ 1 , tanβ 2 , θ and s had significant impacts on the vertical displacements (p < 0.05), and ρ did not have a significant impact on the vertical displacements (p > 0.05). For the monitoring accuracy of horizontal motion in the easting and northing directions, the coefficients of determination R 2 were 0.892 and 0.863, indicating that the independent variables q a , tanβ 1 , tanβ 2 , θ, s and ρ were able to explain 89.2% and 86.3% of the change in the monitoring accuracy in the easting and northing directions, respectively. Parameters tanβ 1 and ρ had significant impacts on the easting and northing monitoring accuracy (p < 0.05), and q a , tanβ 2 , θ and s did not have significant impacts on the easting and northing monitoring accuracy (p > 0.05).
To determine the sensitivity of the 3D surface displacements to each parameter, the single-factor range comparison method was used to analyze the results of this experiment. Figure 8 shows the parameter trends with the order of sensitivity of each parameter as the horizontal axis and the error in the fit between the measured and predicted values as the vertical axis. It can be seen in Figure 8 that for the vertical direction, when parameter q a was at different levels, the range of the fitting error was the largest, indicating that it had the most significant impact on the vertical monitoring accuracy. The range of the fitting error of parameter ρ was the smallest, suggesting that it had the smallest impact on the vertical monitoring accuracy. For horizontal movements in both the easting and northing directions, parameter ρ had the most significant impact, and parameter s had the smallest impact. To determine the sensitivity of the 3D surface displacements to each parameter, the single-factor range comparison method was used to analyze the results of this experiment. Figure 8 shows the parameter trends with the order of sensitivity of each parameter as the horizontal axis and the error in the fit between the measured and predicted values as the vertical axis. It can be seen in Figure 8 that for the vertical direction, when parameter qa was at different levels, the range of the fitting error was the largest, indicating that it had the most significant impact on the vertical monitoring accuracy. The range of the fitting error of parameter ρ was the smallest, suggesting that it had the smallest impact on the vertical monitoring accuracy. For horizontal movements in both the easting and northing directions, parameter ρ had the most significant impact, and parameter s had the smallest impact. To quantitatively analyze the impact of parameter changes on the accuracy of the 3D surface displacements in the mining area, the range of empirical values corresponding to the nature of the overburden was used for each parameter [40] combined with the abovementioned orthogonal test parameter sensitivity ranking results. The more sensitive parameters such as qa, tanβ1, tanβ2, θ and ρ were changed with a step of 0.01 or 0.1, the least To quantitatively analyze the impact of parameter changes on the accuracy of the 3D surface displacements in the mining area, the range of empirical values corresponding to the nature of the overburden was used for each parameter [40] combined with the abovementioned orthogonal test parameter sensitivity ranking results. The more sensitive parameters such as q a , tanβ 1 , tanβ 2 , θ and ρ were changed with a step of 0.01 or 0.1, the least sensitive parameter s was changed with a step of 1, and the estimated values in the vertical, easting and northing directions were calculated. Figure 9 shows the fitting errors.
It is clear in Figure 9 that the changes in q a , θ and s had great impacts on the monitoring accuracy in the vertical direction but had marginal impacts on the monitoring accuracy in the easting and northing directions. As shown in Figure 9a, when the error of q a reached 13% of the given true value, that is, the maximum error of the given range, the error in the vertical fitting accounted for 7.3% of the maximum subsidence while the easting and northing monitoring accuracy was almost unaffected. sensitive parameter s was changed with a step of 1, and the estimated values in the vertical, easting and northing directions were calculated. Figure 9 shows the fitting errors.
It is clear in Figure 9 that the changes in qa, θ and s had great impacts on the monitoring accuracy in the vertical direction but had marginal impacts on the monitoring accuracy in the easting and northing directions. As shown in Figure 9a, when the error of qa reached 13% of the given true value, that is, the maximum error of the given range, the error in the vertical fitting accounted for 7.3% of the maximum subsidence while the easting and northing monitoring accuracy was almost unaffected. In addition, as shown in Figure 9f, ρ had a great impact on the monitoring accuracy in the easting and northing directions but had limited effects on the vertical accuracy. When the error of ρ reached 125% of the given true value, which was the maximum error of the given range, the errors in the easting and northing accounted for 28.9% and 40.4% of the maximum horizontal displacements, respectively. This conclusion is consistent with the results of the orthogonal experiment. Tanβ1 and tanβ2 had impacts on the vertical, easting and northing monitoring accuracy. Note that tanβ1 had a greater impact on the vertical monitoring accuracy than on the easting and northing monitoring accuracy, which is consistent with the orthogonal experimental results. When the error of tanβ1 reached 30% of the given true value, which was the maximum error of the given range, the errors in the vertical, easting and northing accounted for 6.2%, 1.5% and 2.1% of the maximum vertical and horizontal displacements, respectively. When the error of tanβ2 reached 46.7% of the given true value, which was the maximum error of the given range, the errors in the vertical, easting and northing fittings accounted for 8.4%, 2.9% and 3.8% of the maximum vertical and horizontal displacements, respectively. Note that although in the analysis of variance tanβ2 showed limited impacts on the monitoring accuracy in the easting and northing directions, its parameter sensitivity ranked the third in the range analysis, indicating that it still had to some extent impact on the monitoring accuracy in these two directions. This suggests that it is not rigorous to use multifactor analysis of variance to draw conclusions on the influence of parameters in the model, and it is necessary to combine range analysis to assess the sensitivity of the parameters. Based on the above analysis, when a parameter (except for ρ) fluctuated within the given range, the fitting error did not exceed 10% of the maximum deformation value. In this experiment, the parameters qa, tanβ1, tanβ2, θ and s were selected according to the lithology of the overlying rock to meet the requirements of mining subsidence monitoring accuracy, that is, no more than 10% of the maximum deformation value. On the other hand, it is necessary to In addition, as shown in Figure 9f, ρ had a great impact on the monitoring accuracy in the easting and northing directions but had limited effects on the vertical accuracy. When the error of ρ reached 125% of the given true value, which was the maximum error of the given range, the errors in the easting and northing accounted for 28.9% and 40.4% of the maximum horizontal displacements, respectively. This conclusion is consistent with the results of the orthogonal experiment. tanβ 1 and tanβ 2 had impacts on the vertical, easting and northing monitoring accuracy. Note that tanβ 1 had a greater impact on the vertical monitoring accuracy than on the easting and northing monitoring accuracy, which is consistent with the orthogonal experimental results. When the error of tanβ 1 reached 30% of the given true value, which was the maximum error of the given range, the errors in the vertical, easting and northing accounted for 6.2%, 1.5% and 2.1% of the maximum vertical and horizontal displacements, respectively. When the error of tanβ 2 reached 46.7% of the given true value, which was the maximum error of the given range, the errors in the vertical, easting and northing fittings accounted for 8.4%, 2.9% and 3.8% of the maximum vertical and horizontal displacements, respectively. Note that although in the analysis of variance tanβ 2 showed limited impacts on the monitoring accuracy in the easting and northing directions, its parameter sensitivity ranked the third in the range analysis, indicating that it still had to some extent impact on the monitoring accuracy in these two directions. This suggests that it is not rigorous to use multifactor analysis of variance to draw conclusions on the influence of parameters in the model, and it is necessary to combine range analysis to assess the sensitivity of the parameters. Based on the above analysis, when a parameter (except for ρ) fluctuated within the given range, the fitting error did not exceed 10% of the maximum deformation value. In this experiment, the parameters q a , tanβ 1 , tanβ 2 , θ and s were selected according to the lithology of the overlying rock to meet the requirements of mining subsidence monitoring accuracy, that is, no more than 10% of the maximum deformation value. On the other hand, it is necessary to set the range of ρ to [0.19, 0.69] (indicated by the dashed line in Figure 9f). From the analysis above, when we use the model in this paper to obtain the 3D surface displacements in mining areas, if the vertical subsidence error is large, we can focus on adjusting the parameters q a , tanβ 1 , tanβ 2 , θ, and s, which have greater impacts on the vertical monitoring accuracy, to achieve high monitoring accuracy. When the errors of horizontal movements in the easting and northing