Underground Goaf Parameters Estimation by Cross ‐ Iteration with InSAR Measurements

: Determining the geographic location and spatial distribution of underground goaf is of great significance for the prevention of mining subsidence hazards and the detection of illegal min ‐ ing. However, traditional goaf detection techniques mainly focus on geophysical methods that are labor intensive, have low efficiency, and are expensive. Due to the large range and off ‐ site monitor ‐ ing capability of interferometric synthetic aperture radar (InSAR) techniques, research on goaf loca ‐ tion detection based on InSAR measurements has been increasing. This paper proposes a new method for locating underground goaf based on cross ‐ iteration and InSAR measurements. Firstly, the functional relationship between the geometric parameters of the goaf and the line of sight (LOS) deformation retrieved by InSAR techniques is constructed. Then, the three initial model parameters of the probability integration method (PIM) are determined by mining geological conditions. Fi ‐ nally, the cross ‐ iteration method is used to determine the parameters to characterize the spatial lo ‐ cation of underground goaf. The experimental results show that the average relative errors of the simulated experiment and the real experiment are 1.5% and 5.1%, respectively, and the inverted goaf parameters are in good agreement with the real values. Moreover, the proposed method only requires the main lithology of the overlying rock in the goaf and does not depend on the accuracy of PIM model parameters. Therefore, this method has engineering application value for the detec ‐ tion of goaf lacking actual measurement data or that caused by illegal mining.


Introduction
Underground coal mining is the main way to obtain coal resources. The underground cavities formed by the mining of coal seams are called goafs. With the large-scale mining of coal, a large number of unknown goafs mainly caused by illegal mining activities have been generated in the world. In addition, the destruction or loss of ancient coal seam mining records make it difficult to know the actual situation of many abandoned goafs. These unknown goafs may not only endanger people's lives and property and destroy surface infrastructure, but also cause landslides, surface subsidence, and other geological disasters. Therefore, determining the geographic location and spatial distribution of underground goafs is of great significance for the early warning and management of geological disasters. Traditional underground goaf detection technologies focus on geophysical techniques, such as transient electromagnetic method [1,2], microgravity measurement [3], and remote thermal IR surveying [4]. These methods can roughly determine the spatial distribution and geographic location of the goaf. However, they are time-consuming, labor-intensive, costly, and have a small detection range.
Interferometric synthetic aperture radar (InSAR) is a new technology that can use the interferometric phase information of SAR image pairs to obtain the surface deformation in the line of sight (LOS) direction [5]. InSAR technology can provide wide area coverage, high spatial resolution, and high monitoring precision of surface deformation fields [6,7]. It is gradually being used in many fields, such as the monitoring of volcanic activities [8,9], permafrost [10][11][12][13], landslides [14,15], and earthquakes [16][17][18][19]. In addition, InSAR techniques have many successful applications in critical infrastructure monitoring fields [20][21][22]. Moreover, with the continuous development of InSAR technology, a wealth of SAR sensors and a large number of historical archived images already exist, and free SAR image acquisition methods are constantly enriched. Compared with traditional goaf detection technology, InSAR technology has the unique advantages that provide the possibility for precise detection and positioning of goaf in a large area [23].
At present, there are few studies using InSAR monitoring results to locate goaf. In 2013, Hu et al. [24] proposed a differential InSAR (DInSAR)-based illegal-mining detection system (DIMDS), which can determine the center geographic coordinates of the underground goaf by combining the acquired LOS deformation and the mining subsidence law. However, the rough center coordinates of the goaf are not sufficient to accurately assess potential geological disasters. In order to obtain more effective goaf geometric information, Yang et al. [25] have constructed the functional relationship between the LOS deformation and the eight parameters of underground goaf based on the probability integration method (PIM). Then, the SAGA hybrid algorithm [26] was used to invert the unknown goaf parameters. However, the complexity of the function model will affect the efficiency of parameter inversion. Du et al. [27] extracted feature points on the sinking curves to participate in the calculation of the coal seam boundary and mining depth. Subsequently, Xia et al. [28] proposed a goaf positioning method that can calculate more parameters and is suitable for inclined coal seams. Similar to the method proposed by [27], this method ignores the influence of horizontal deformation and obtains incomplete goaf parameters. More importantly, these methods [25,27,28] all require prior model parameters of the PIM that are usually obtained from measured data. If the actual measurement data are insufficient, large PIM model parameter errors may be introduced, and the positioning accuracy of underground goafs may be seriously affected. Therefore, reducing the number and accuracy requirements of PIM model parameters plays an important role in the precise positioning of the unknown goaf.
In summary, this paper presents a cross-iteration method for accurate inversion of full underground goaf parameters. This method reduces the input of model parameters of the PIM and is not significantly affected by the accuracy of the required model parameters. We first modify the PIM model to reduce the input model parameters. Then, we constructed the functional relationship between the LOS deformation and the geometric parameters of the goaf and determined the three initial model parameters of the PIM related to the geological mining conditions. Finally, the geometric parameters of the goaf were inverted by a cross-iteration method. The reliability of the method was verified through simulation experiments and real experiments.
The remaining structure of the paper is as follows. The principle of the cross-iteration method is described in Section 2. Experiments and results are shown in Sections 3 and 4, respectively. We state the discussion in Section 5, and the conclusion is presented in Section 6.

Principle of the PIM
The PIM model is a widely used mining subsidence prediction model developed on the basis of the random medium theory [29,30]. This model divides the goaf into many small mining units and regards the impact of mining on the ground surface as the sum of the mining impacts of all small units [31]. Supposing that there is an underground goaf caused by long-wall mining, in order to determine the geographic location and spatial distribution of this underground goaf, a set of parameters is usually required. This set of parameters usually obtained through actual observations is regarded as the geometric parameters of the goaf. As is shown in Figure 1, the geometric parameters of the goaf are given by a true boundary length , true boundary width , coal seam mining height , dip angle , geodetic coordinate of the origin , , and azimuth angle of mining direction , i.e., , , , , , , , . As long as the mining direction and related parameter values are determined, the vertical deformation , and horizontal deformation , , in direction at a point , can be predicted based on the PIM model [32].
, , cos sin / In Equations (1) and (2), cos is the maximum subsidence value of the ground; is the subsidence factor; and are the vertical subsidence on the main section of the strike and dip directions, respectively (see Figure 2); and are the horizontal deformations on the main section of the strike and dip directions, respectively. Furthermore, and sin /sin are the calculated boundary length and width of the goaf, where (see Figure 2 b) is the influence propagation angle of mining, , , , and (see Figure 2) are the offsets of inflection points in strike-back, strike, down-dip, and up-dip directions.
/ tan , / tan , and / tan are the main influence radius in down-dip, up-dip and strike directions, respectively, where tan , tan , and tan are the tangents of the major influence angles in the three directions.
(see Figure 2 a) is the average mining depth, and (see Figure 2b) are the mining depths in down-dip and up-dip directions, and is the horizontal displacement constant.
In addition to the geometric parameters of the goaf , , , , , , , , there are a set of model parameters , , tan , tan , tan , , , , , that depend on the geological mining conditions [31]. Model parameters of the PIM are usually determined based on actual measured data and will vary with different geological mining conditions. The tangent of the major influence angle depends on the main properties of the overlying rock [25,31]. It can be approximately assumed that the tangents of the major influence angle in three directions are equal, i.e., tan tan tan . The offsets of the inflection points 1,2,3,4 are only used to obtain the calculated boundary length l and width of the goaf, and and are negatively correlated with the actual boundary length and width of the goaf, respectively. In order to simplify the PIM model, we can assume that and . Based on the above assumptions, the model parameters are modified to , , tan , , , , which not only simplifies the PIM model, but also reduces the difficulty for model parameter acquisition.

The Function Model linking InSAR LOS Measurement to Goaf Parameters
Using InSAR technology to monitor surface deformation can obtain millimeter-level or even higher precision LOS deformation fields. Based on the imaging geometry of SAR, the InSAR-derived deformation can be regarded as the comprehensive projection of the north-south, east-west, and vertical deformation components , , in the LOS direction.
cos ⋅ sin sin ⋅ sin cos ⋅ where is the incidence angle of the SAR sensor; is the heading angle. The PIM model can predict the 3D displacement of the surface settlement point and fully includes the required spatial distribution parameters of the goaf. Therefore, substituting Equations (1) and (2) into Equation (7), we can construct the functional relationship between the LOS deformation and the PIM model.
According to Equation (8), when the model parameters are known, the non-linear optimization algorithm can be selected to reverse the geometric parameters of the goaf . Moreover, using the large-area surface deformation data monitored by InSAR to directly invert fewer geometric parameters of the goaf inevitably reduces the calculation efficiency. Therefore, selecting an appropriate subsampling method not only effectively avoids data redundancy, but also saves sufficient data information. In this paper, an adaptive quadtree subsampling method [33,34] is selected to process deformation points. Compared with the traditional quadtree method, this method can divide irregular polygonal areas and reduce the influence of null values in InSAR data.

Model Resolution via Cross-Iteration
The PIM model is complex and nonlinear, which makes it difficult to directly determine the value of the model parameters . For mining areas lacking actual measurement data, model parameters of mining areas with similar geological structure conditions or model parameters of different working faces in the same mining area are usually used as approximate replacements. Thus, large model parameter errors may be introduced, which will result in lower inversion accuracy. In addition, a set of empirical parameters can also be determined according to the geological mining conditions, which may also seriously affect the inversion accuracy. Therefore, a cross-iteration method for the inversion of goaf parameters is proposed. This method only needs the empirical values of three initial model parameters , , tan to obtain the goaf location results with high accuracy. Cross-iteration is an iterative method of cross-updating between grouped parameters, through the optimization algorithm to invert the pending parameters of each iteration. Since the magnitude of the subsidence factor , the horizontal displacement constant , and the tangent of the major influence angle tan is small, the value interval is easily determined. First, we can obtain the empirical values of these three model parameters , , tan according to the main properties of the overlying rocks in the goaf. Then, we must select an improved genetic algorithm (GA) [32] to invert the goaf parameters , , , , , , , and the remaining model parameters , , . The fitness function of the improved GA can be expressed as min , is the predicted LOS deformation based on Equation (8). The processing steps of the improved GA are as follows: By repeatedly running GA for times, solutions 1,2, ⋯ , and corresponding fitness values 1,2, ⋯ , can be obtained. Then, the RMSES of each parameter were calculated based on the GA-derived resolutions. Subsequently, we eliminated solutions with an error of more than twice the corresponding RMSES. Finally, the reciprocal of the fitness values 1,2, ⋯ , were used to weight the remaining solutions 1,2, ⋯ , to obtain the unknown parameters , i.e., The influence propagation angle of mining , the offsets of inflection points , , the dip angle , the azimuth angle of mining direction , and the geodetic coordinate of the origin , are less affected by , , tan . Therefore, these seven parameters are used as known to invert the remaining parameters using the improved GA. In this way, two iterations were carried out, and a crossover process was completed. When the ratio of the difference between the mining depths obtained in two adjacent iterations to the previous mining depth is less than the threshold, the geometric parameters of the goaf are relatively stable. Thus, the cross-iteration is terminated. The specific process of crossiteration (see Figure 3) is as follows.
(1) First Iteration : Determine the empirical values of the three initial model parameters  , , tan  and  invert  the  remaining  unknown  parameters  , , , Second Iteration: Treat , , , , , , as known and update the remaining parameters to , , tan , , , , .
(3) Iteration Termination Judgment: When the condition / 0.05 is satisfied, the iteration is terminated, and the target parameter value is , , , , , , , . If the termination condition is not satisfied, the parameter values , , tan obtained in the second iteration are used to repeat the two steps (1) and (2) until it is satisfied / 0.05.

Overall Processing
The proposed method in this paper is essentially based on the constructed function model to perform a cross-iterative inversion to obtain the required parameters. Therefore, we first had to build a function model. Prior to this, InSAR data processing was required to obtain the cumulative LOS deformation within the target time range. Then, we linked the modified PIM model shown in Section 2.1 with the cumulative LOS deformation to construct the inversion model of the goaf parameters. Subsequently, the empirical values of the initial model parameters were determined according to the mining geological conditions of the study area. At this point, we had completed the preparatory work before the cross-iteration processing. Finally, we performed the inversion of the goaf parameters according to the cross-iteration steps shown in Section 2.2.2 and output the results. In this way, we completed the goaf positioning process as shown in Figure 3.

Simulated Experiment
A simulated experiment was carried out. At first, a ground displacement field was simulated. We assumed that there was an underground goaf with a coal seam mining height of 3 m, a boundary length of 600 m, a boundary width of 150 m, a dip angle of 12°, a mining direction angle of 100°, and an average mining depth of 500 m. The overlying rock of the assumed goaf was mainly medium sandstone and limestone. The large mining depth and the type of overlying rock determined the continuity of surface deformation. Subsequently, we simulated the 3D displacements of the surface caused by the mining of the goaf. Assuming that the incidence angle and heading angle of the SAR sensor were 33.66° and −10.39°, the 3D displacements were converted to the simulated LOS deformation observations according to Equation (7). The spatial resolution of the simulated LOS deformation is 5 m × 5 m. As can be seen in Figure 4, the LOS deformation basin was approximately symmetrical, the largest deformation point was located in the center of the basin, and the largest deformation value was 0.54 m. Moreover, the relative positional relationship between the simulated goaf (the black polygon in Figure 4) and the LOS deformation basin in the same plane conformed to the characteristics of mining subsidence.

Study Area
Fengfeng coalfield is located in the southwest of Handan City, Hebei Province, China. The fold that controls the overall morphology and structure of the mining area is the Gushan-Zishan anticline; the Zishan anticline in the northern section has an axis of NNE and the Gushan anticline in the southern section has an axis close to SN [35]. The Gushan-Zishan anticline divides the mine into two parts, the eastern monoclinic and the western syncline [36]. The coal-bearing strata of the Fengfeng mining area belong to the Permo-Carboniferous, and they include six to seven mineable coal seams [37]. The Shanxi Formation of the Permian is mainly composed of siltstone, fine gray sandstone, and coal seams, with one mineable coal seam [35]. The Taiyuan Formation of the Carboniferous is mainly composed of gray to dark gray siltstone, sandstone, limestone, and coal seams, with five to six mineable coal seams [35]. Long-term and large-scale coal mining has resulted in a large number of abandoned underground goafs that have caused many related disasters [38]. Therefore, the accurate determination of these unknown old goafs can protect people's lives and property and allow sustainable development of the mining area. In this study, we selected the surface subsidence area caused by mining the 132610 working face (the white polygon in Figure 5) as the study area. From the literature [28], the coal seam mining height is 4.5 m, the boundary length is 319 m, the boundary width is 165 m, the dip angle is 31°, the mining direction angle is 169°, and the average mining depth is 774 m.

SAR Data Sets and InSAR Processing
Twenty C-band Sentinel-1A images covering the study area were selected to monitor the surface deformation. The basic parameters of the Sentinel-1A data set were shown in Table 1 The acquisition time interval was largely consistent; thus, we considered the actual goaf parameters to still be applicable. The surface deformation caused by mining develops rapidly and has a large magnitude. Therefore, we chose the two-connection network construction method to reduce the impact of time decoherence. As shown in Figure 6, a total of 37 interferometric pairs were selected.  In order to obtain a reliable cumulative deformation caused by the mining of the 132,610 working face, we performed SBAS-InSAR processing [39]. GAMMA software was used to process the selected interference pairs. First, we used a network-based method [40] for image registration. Then, the DEM data of SRTM 1 arc second (30 m) were applied to remove the topographic phases of the interferograms, and we performed a 5 × 1 multilook operation on the range and azimuth directions to obtain rectangular pixels. Third, the adaptive filter based on the local fringe spectrum [41] was used to further suppress phase noise in the interferograms and generate filtered coherence maps. Afterwards, we masked the points with coherence less than 0.2 and used the minimum cost flow unwrapping method [42] to calculate the ambiguity of the differential interferometric phases. Finally, a polynomial model [43,44] was selected to eliminate the residual phase ramps caused by the atmospheric or baseline modeling error in the differential interferograms. Subsequently, we performed spatio-temporal filtering to remove atmospheric delay components. The cumulative deformation sequence from 17 June 2015 to 7 March 2016 is shown in Figure 7. Because the coal roof is composed of limestone and siltstone, and the mining depth is large (an average of 774 m), continuous deformation of the surface was caused. This was consistent with the time series deformation shown in Figure 7, which proved the accuracy of the InSAR monitoring results from the side.

Inversion of Goaf Parameters
Based on the simulation of ground displacement, the inversion of goaf parameters was conducted. We first constructed the functional relationship between the simulated LOS deformation and the geometric parameters of the goaf. Subsequently, the proposed cross-iteration method was used to estimate the geometric parameters of the goaf. The improved GA was used to retrieve unknown parameters in the cross-iteration processing. More specifically, we repeated the GA for 100 times, with a population size of 100, cross fraction of 0.95, and maximum generations of 500. After eliminating gross errors, the final results were obtained by weighting the remaining solutions with the reciprocal of the fitness. It was noteworthy that the empirical values of the subsidence factor , the horizontal displacement constant , and the tangent of major influence angle tan needed to be determined in advance. Since the main lithology of the overlying rock in the simulated goaf was medium-hard, the initial model parameters were obtained based on [45], i.e., 0.72, 0.23, tan 2.1. In the cross-iteration procedure, the value ranges of the parameters were fixed. The value range of each parameter was determined by [45], the common value range and the geometric characteristics of the LOS deformation basin.

Accuracy Evaluation
To visually show the effect of the proposed method on the spatial positioning of the goaf, the simulated underground goaf and the positioning result are depicted in a threedimensional space, as shown by the black and purple rectangles in Figure 8. It is easy to see that the shapes of the two goafs are similar, but the average mining depth from the inversion is greater than the simulated average mining depth. We projected the two goafs onto a two-dimensional plane. It can be seen from Figure 4 that the gap between the two is not obvious, and they roughly coincide with each other.
To quantitatively evaluate the inversion accuracy of the goaf, we calculated the absolute error and relative error between the inverted and simulated geometric parameters of the goaf (see Table 2). It can be seen that the maximum absolute error was -44.6 m, and the minimum absolute error was 0.03°. There was a significant difference between the absolute error of each parameter. The main reason is that the magnitudes and dimensions of the parameters were not exactly the same, making it difficult for the absolute error to meet the overall requirements for accuracy evaluation. Since the relative error ignores the unit, the accuracy of the parameters can be evaluated uniformly. Table 2 shows that the relative error of goaf parameters was concentrated between 0.0% and 8.9%, and the average relative error was 1.5%. This indicates that the inversion accuracy of the goaf parameters was high, which further proves the reliability of the proposed method.

Inversion of Goaf Parameters
According to Equation (8), the construction of the function model needed to determine the local coordinates of the InSAR monitoring deformation. The ground resolution of the geocoded deformation data was 14.324 m × 15.286 m. Therefore, we chose the latitude and longitude of the center of the pixel as the geographic coordinates of the deformation data and converted them to geodetic coordinates. Then, the geodetic coordinates were converted to local coordinates by determining the origin geodetic coordinates and azimuth angle of mining direction. In the cross-iteration process, the local coordinates of the LOS deformation point changed with the iteration of the origin geodetic coordinates and azimuth angle of mining direction.
In this experiment, we used the adaptive quadtree subsampling method to process deformation points. The calculated LOS deformation after adaptive quadtree subsampling is shown in Figure 9. Before the cross-iteration, the empirical values of the three initial model parameters were determined by [45], i.e., 0.6, 0.3, tan 1.95. According to the subsampling LOS deformation data and the initial model parameters, the proposed method was used to invert the geometric parameters of the underground goaf. We used the improved GA to retrieve unknown parameters in the cross-iteration processing. In each iteration, we repeated the GA for 100 times, with a population size of 100, cross fraction of 0.95, and maximum generations of 500. The results are shown in Table 3.

Results and Accuracy Evaluation
We comprehensively analyze the reliability of the cross-iteration method from two perspectives. First, we drew the underground goafs in Figure 10 based on the inverted and actual geometric parameters of the goaf to qualitatively evaluate the positioning results. As is shown in Figure 10, the inverted goaf is not located in the center of the surface subsidence basin but is biased towards the side of the largest sinking point. This is consistent with the general rule of mining subsidence.
Then, we respectively calculated the relative error and absolute error between the inversion value and the measured value of the goaf parameters to quantitatively test the effect of the cross-iteration method. Due to the lack of complete measured data, we did not perform precision analysis on the geodetic coordinate of the origin 4,054,325m, 514,206m obtained from the inversion. In Table 3, the maximum relative error is 10.7% of the coal seam mining height . From the perspective of absolute error, the coal seam mining height was only overestimated by about 0.48m. In a comprehensive consideration, the accuracy of the coal seam mining height was acceptable. The relative error of the average mining depth was 9.1%, similar to the result of the simulation experiment. Furthermore, the difference between the inverted values and measured values of the remaining goaf parameters was small. The average relative error of the six parameters of the goaf, shown in Table 3, was 5.1%, which meets the actual requirements for the location of the goaf. Then, we compared the goaf parameters obtained by the proposed method with those obtained by Xia et al. [28] to further verify the reliability of the cross-iteration method. It can be seen in Table 3 that the geometric parameters of the goaf, except for the average mining depth, have a higher accuracy. The sources of the errors in the inverted geometric parameters of the goaf are discussed in Section 4.3. Figure 10. Comparison between the inversion of the goaf geometry and the actual goaf geometry. The gray polygon is the actual goaf, and the light blue polygon is the inverted goaf.

Impact of the Cross-Iteration Method
The cross-iteration method uses iteration to improve the accuracy of parameter inversion. Therefore, we first analyzed the impact of iteration on the results. We compared the accuracy of iterative and non-iterative inversion results in real experiments. As shown in Tables 3 and 4, the absolute and relative errors of the inversion parameters without iteration have increased significantly. The average relative error of the six parameters of the goaf increased from 5.1% to 15.5%. The inversion results show that cross-iteration can improve the accuracy of the goaf parameters. Simulation experiments can independently verify the influence of each impact condition on the results. Hence, simulation experiments can be designed to better test the performance of the cross-iteration method. We choose the simulated goaf described in Section 3.1, and assumed the PIM model parameters, i.e., 0.6, 0.3, tan 2.0, 82.2, 30, 30. Subsequently, we predicted the surface LOS deformation by Equation (8). We constructed the functional relationship between the predicted LOS deformation and the PIM model and used the initial model parameters to invert the geometric parameters of the underground goaf. The results obtained by the cross-iteration method are listed in Table 5. In order to verify the impact of the cross-iteration method, we calculated the relative error and absolute error of each geometric parameter of the goaf. The accuracy of this experiment can be evaluated according to the magnitude and error of the goaf parameters. As shown in Table 5, when there is no error in the LOS deformation, the absolute error and the relative error are both concentrated on the average mining depth . The larger error proves that the proposed cross-iteration method had the greatest impact on the average mining depth . Although the accuracy of the average mining depth was lower than the other goaf parameters, the relative error of within 10% can meet the requirements of goaf positioning. The remaining geometric parameters of the goaf were roughly in agreement with the simulated values. In addition, the average relative error of the geometric parameters of the goaf was 3%. Therefore, the cross-iteration method obtained the goaf parameters with high overall accuracy, which meets actual engineering needs. Table 5. Comparison between the inversion results without errors and with random errors in deformation.

Deformation Parameters
Error-free

Influence of Deformation Error
InSAR technology solves the displacement vector of the ground object in the LOS direction by separating the interferometric phase [46]. Although we performed the operations as described in Section 3.2.2, the unwrapped deformation phase still contained terrain residuals, orbital errors, unwrapping errors, atmospheric effects, and noise. We tested the influence of LOS deformation error on the positioning results of the goaf through a simulation experiment. We used the predicted LOS deformation described in Section 4.3.1. In order to approximate the actual solution, we added −5cm to 5cm random errors to the predicted LOS deformation. Furthermore, we used the precise values of the three initial model parameters to participate in the cross-iteration so that the impact of LOS deformation errors on the positioning results could be effectively evaluated. Afterwards, a cross-iteration experiment was carried out on the LOS deformation with random errors. The inversion results and relative errors of goaf parameters are shown in Table 5. Table 5 shows the accuracy of the coal seam mining height and dip angle obtained by the predicted LOS deformation with a decrease in random errors. The absolute error of the two parameters changed from 0.13 to 0.35 and from 0.45 to 1.05, respectively. As the magnitude of the two parameters is small, the absolute error did not clearly change. The relative error of the two parameters changed significantly, increasing by 7.4% and 5%, respectively. Moreover, the accuracy of the remaining goaf parameters was not distinctly affected by the LOS deformation error. In particular, the average mining depth , which was greatly affected by the cross-iteration method, was not obviously affected by the LOS deformation error. The relative error range of the inverted goaf parameters based on LOS deformation with random errors was 0% to 11.7%, and the average relative error was 4.4%, while the relative error range of the inverted goaf parameters based on error-free LOS deformation was 0% to 8.6%, and the average relative error was 3.0%. The mean relative error and relative error range of the goaf parameters obtained by error-free LOS deformation were both smaller than the former. This shows that the reliability of the goaf parameters was negatively affected by the accuracy of the InSAR-derived deformation. Therefore, using accurate surface LOS deformation can improve the inversion accuracy of the cross-iteration method.

Effect of the Initial Model Parameter Error of the PIM
For goafs that are unknown or lack actual measurement data, the proposed method determines the empirical values of the initial model parameters of the PIM through geological mining conditions, which inevitably contains certain errors. In this section, simulation experiments were conducted to test the influence of model parameter error distribution on the inversion accuracy of goaf parameters. First, we chose the predicted errorfree LOS deformation as in Section 4.3.1. Afterwards, we sequentially added -50% to 50% error with an interval of 10% to the three initial model parameters. The detailed method is to select one of the initial model parameters to add −10% to −50% and 10% to 50% error, and the other two parameters do not have an added error. By selecting different initial model parameters that required the addition of errors, a total of 30 sets of parameters can be obtained. Finally, we used the predicted LOS deformation and each set of initial model parameters to perform cross-iteration to invert the geometric parameters of the goaf. The experimental results are shown in Figure 11.
It can be seen from Figure 11 that the accuracy of the initial model parameters of the PIM had no significant effect on the results of the cross-iteration processing. When the error contained in the subsidence factor and tangent of major influence angle tan increased, the boundary length , boundary width , azimuth angle of mining direction , and geodetic coordinate of the origin , were largely consistent with the simulated values (red line in Figure 11). When the error contained in the horizontal displacement constant increased, the boundary length , boundary width , and origin ab-scissa of the local plane coordinate system had no obvious change. However, the accuracy of dip angle , the origin ordinate of the local plane coordinate system , and the azimuth angle of mining direction had a slight downward trend, but there was no significant difference between these two parameters and the simulated values. The calculated average mining depths were similar to the results in Section 4.3.1. Therefore, we believe that the accuracy of the average mining depth was not greatly affected by the initial model parameters. Afterwards, we calculated the average relative error of the goaf parameters, except the geodetic coordinate of the origin , , when different initial model parameters included errors. As shown in Figure 12, the maximum average relative errors of the coal seam mining height and dip angle were 8.4% and 8.9%, respectively. We found that the horizontal displacement constant had a greater impact on the coal seam mining height than other initial model parameters, and the subsidence factor and tangent of major influence angle tan had a greater impact on the dip angle than the horizontal displacement constant . Furthermore, we considered the absolute error to more fully analyze the impact of the initial model parameters. When the average relative error was equal to 10%, the coal seam mining height and dip angle changed by 0.3 m and 1.2°, respectively. This was sufficient to prove that the relative error of less than 10% did not have much influence on the inversion values of the coal seam mining height and dip angle . In addition, the initial model parameters, including errors, had little effect on the geometric parameters of the goaf except for the coal seam mining height and dip angle (see Figures 11 and 12).

Impact of the Tangent of the Main Influence Angle on the Average Mining Depth
We analyzed the reasons for the large average mining depth error from the principle of the PIM model. Equations (3) and (5) show that the average mining depth was only used to calculate the main influence radius in strike direction , i.e., / tan . Furthermore, we found that the tangent of the major influence angle tan was overestimated in the cross-iteration processing. If the main influence radius in strike direction is regarded as a constant, the tangent of the major influence angle tan has a positive correlation with the average mining depth . Based on the simulated experiment described in Section 4.3.1, we calculated the ratio of to tan in each iteration. The ratios were 251.7, 252.8, and 250.2, respectively. This shows that the ratio of the average mining depth H to the tangent of major influence angle tan was close to the true value 250, so the assumption that is a constant was established. Moreover, the magnitude of the two parameters was quite different; a slight fluctuation of the tangent of the major influence angle tan will cause a large change in the average mining depth . Therefore, the main reason that the average mining depth was overestimated is that the tangent of the major influence angle tan was overestimated. Solving this problem will be the subject of our future research.

Influence of Subcritical Mining
The ratio of goaf size to mining depth determines the mining degree. When / , / 1.2 ∼ 1.4, the mining degree is subcritical, and when / , / 1.2 ∼ 1.4, the mining degree is critical or supercritical, where and are the true boundary length and width of the goaf and is the average mining depth. The ratio of the true boundary length and width of the goaf to the average mining depth in the real experiment is 0.4 and 0.2, respectively. This shows that the mining degree of the 132610 working face is subcritical in the strike and dip directions. The PIM model generally overestimates the 3D surface deformation induced by subcritical mining, because the subsidence factor for subcritical mining is smaller than the subsidence factor under the critical or supercritical mining. As analyzed in Section 4.3.3, the increase of the subsidence factor has no obvious effect on the geometric parameters of the goaf. Therefore, subcritical mining is probably not the main factor affecting the accuracy of goaf parameters. The results of the real experiment roughly conform to this point of view.

Conclusions
This paper proposes a cross-iteration method to determine the spatial location of the underground goaf. The proposed method reduces the required PIM model parameters to three and does not significantly depend on their accuracy. We verified the reliability of the proposed method through experiments. The average relative errors between the geometric parameters of the goaf inverted by the simulation experiment and the real experiment and the real value were 1.5% and 5.1%, respectively.
Analyzing the influence of various factors on the inverted goaf parameters, it was found that the cross-iteration method was the main reason for the overestimation of the average mining depth, but the overestimated value generally did not exceed 10% of the true value. The deformation error monitored by InSAR caused a negative impact on the goaf positioning results, thus ensuring that the accuracy of the InSAR monitoring can improve the accuracy of the goaf positioning. Furthermore, the main reason for the overestimation of mining thickness and coal seam dip in real experiments and simulation experiments may be that the LOS deformation and initial model parameters contained errors, but the overestimated values were not significant. However, there was no obvious linear correlation between the magnitude of the error included in the initial PIM model parameters and the accuracy of the geometric parameters of the goaf. In addition, we should also consider the inapplicability of the PIM model in subcritical mining, which may also be one of the reasons that affect the positioning accuracy of the goaf. Therefore, the improvement of the cross-iteration method and the consideration of the impact of subcritical mining will be the main directions of our future research.