An Improved GNSS and InSAR Fusion Method for Monitoring the 3D Deformation of a Mining Area

Accurate measurement of the surface three-dimensional deformation fields in a mined-out area is essential for understanding the law of mining subsidence and guiding the safe production and construction of mining areas. In this study, we proposed an improved fusion method to monitor this three-dimensional deformation. We used the Helmert variance component estimation (HVCE) method and the back-propagation neural network (BPNN) to fuse the global navigation satellite system (GNSS) and interferometric synthetic aperture radar (InSAR) measurements. Based on this improved fusion method, we measured the three-dimensional deformation in the West Second Mining Area of Jinchuan, Jinchang City, Gansu Province, China, from March 22, 2019, to June 8, 2020. The root mean square errors (RMSEs) based on our method for the three directions of east-west (E-W), north-south (N-S), and up-down (U-D) were 21.4 mm, 7.6 mm, and 34.2 mm, respectively. These RMSE values are lower than those obtained from previous methods. The results show that the spatial distribution of the three-dimensional deformation follows the law of mining subsidence.


I. INTRODUCTION
Although mineral resources have made an essential contribution to the rapid social and economic development of the world, they have also brought a series of geological and security problems in mining areas [1][2][3][4][5][6]. The surface deformation of the mined-out area is one of the most severe manufactured geological disasters in a mining area [7][8][9][10]. Therefore, three-dimensional deformation monitoring of mining areas helps to understand the processes and mechanisms of surface deformation caused by mining [11][12][13][14][15].
In recent decades, with the continuous development of surveying and mapping technologies, significant progress has been made in monitoring methods and means of surface deformation caused by mining. Traditional global navigation satellite system (GNSS) measurements can accurately measure three-dimensional surface deformation, which are point-based measurements. However, GNSS points are relatively sparse and have a low spatial resolution, which is insufficient to show detailed surface deformation [16][17][18][19][20].
As an increasing number of spaceborne synthetic aperture radar (SAR) sensors have been launched into space, interferometric SAR (InSAR) measurements have achieved rapid development in the deformation monitoring of mining areas in the past 25 years [21][22][23][24][25][26][27]. InSAR measurements can measure large-scale surface deformation with broad spatial coverage, day and night observations, and short revisit intervals, which have been widely used [28][29][30]. However, InSAR measurements are the projection of the threedimensional spatial deformation in the satellite line-of-sight (LOS) direction, which is a deformation in a single direction. In terms of two-dimensional and three-dimensional deformation monitoring, multitracks and multiplatform SAR image monitoring results have often been used for joint calculations [31][32][33]. However, the satellite flies in polar orbit, resulting in InSAR measurements that are insensitive to deformation in the north-south direction [34][35][36][37].
There has been interest in combining GNSS and InSAR measurements to monitor three-dimensional deformation owing to the complementary nature of the measurements [38][39][40]. A combination method is to apply GNSS and In-SAR separately to comprehensively analyze the surface deformation, such as the study of spatial and temporal deformation evolution in Yellowstone [41,42]. In the other method, the sparse GNSS points are interpolated into the same spatial grid of InSAR, and then three-dimensional deformation fields are obtained by building a fusion model. In [43], the authors applied Markov random field-based regularization and a simulated annealing algorithm to estimate the three-dimensional deformation fields of the Reykjanes Peninsula, southwest Iceland. An improved method named analytical optimization in [44] was suggested to combine GNSS and InSAR observations and then used to estimate the three-dimensional deformation fields of Southern California. All these methods need to accurately determine a priori variances of the GNSS and InSAR observations in the fused model. In general, it is of great significance to determine the exact weights for each kind of measurement. To achieve this purpose, [45][46][47][48][49] introduced the Helmert variance component estimation (HVCE) method into the fusion of GNSS and InSAR measurements to estimate their a posteriori variances. It has been demonstrated that the HVCE method can determine accurate weights for each type of measurement without any prior information. However, since the performance of the HVCE method depends on the redundancy of observations, the problem of negative variances may occur, resulting in failure to fuse properly [50,51]. The back-propagation neural network (BPNN) was proposed by Rumelhart in 1986 [52] and is one of the most widely used artificial neural networks in machine learning models [53][54][55]. Unlike the traditional linear regression model, the BP neural network has no strict requirements for data distribution. Therefore, the results of using this neural network to deal with negative variances are usually robust.
In this study, an improved method that uses the HVCE-BPNN method to obtain three-dimensional deformation is proposed. It is well known that the BPNN has a strong learning ability to compensate for the disadvantages of the HVCE method. The proposed method was used to obtain the three-dimensional surface deformation caused by mining of the 1580 working face in the Jinchuan West Second Mine, Jinchang City, Gansu Province, China. According to the validation of different fusion schemes along with GNSS and InSAR measurements, this method is useful for monitoring three-dimensional ground surface deformation in a mining area.

II. STUDY AREA AND DATA SOURCES
The study area is the West Second Mining Area of Jinchuan, located in Jinchang City, Gansu Province, China (Fig. 1). The West Second Mining Area is located in the middle part of the Hexi Corridor and the northern foot of the Qilian Mountains at high elevations ranging from 1700 to 1331 m [56]. The geographic coordinates of the study area are between 102°7′46″-102°8′25″ E and 38°29′16″-38°29′45″ N.
In this study, the exploitation range is 420 m long from north to south and 230 m wide from east to west, with a total area of 43 512 m 2 . The West Second Mining Area has been affected by underground mining and F8 fault activity for a long time, resulting in a large area of surface collapse, and its scope continues to expand [57].
Both GNSS and InSAR measurements are space geodetic methods capable of all-weather observations. The GNSS observation data used in this study were obtained from GNSS receivers in the mining area. As shown in Fig. 1, 139 GNSS observation points were uniformly laid on the surface of the exploitation range, monitored from March 22, 2019, to June 8, 2020, once every 2 months. The image data used in this study are a total of 67 frames of single look complex (SLC) SAR images, which are generated in interferometric wide swath (IW) mode from the Sentinel-1A satellite operated by the European Space Agency (ESA). The SAR images were cropped into a region of interest (ROI) to improve the processing efficiency. A zone of approximately 1.36 km × 1.20 km focusing on the mining subsidence area is then selected in this study. The coverage of Sentinel-1A is shown in Fig. 2, and the related parameters are listed in detail in Table 1.

A. INSAR GEOMETRIC TRANSFORMATION
Satellite-borne InSAR technology is a common method for monitoring surface deformation. The disadvantage of this measurement is that it can only monitor deformation in the LOS direction, and the satellite LOS deformation represents the projection of the three-dimensional surface deformation [58]. The LOS deformation of a satellite is composed of east-west (E-W), north-south (N-S), and up-down (U-D) deformations. Fig. 3 shows the deformation transformation relationship between the InSAR LOS and the three components of East-North-Up (ENU), and the spatial geometric relationship is given by (1).
represents the distortion of the satellite's LOS, and = [ ] represents the unknown deformation vector in E-W, N-S, and U-D. = [ ] represents the unit projection vector of the three components, where = − , = , and = . Furthermore, represents the incidence angle of the satellite, and represents the azimuth angle.

B. HVCE-BPNN METHOD FOR GNSS AND INSAR FU-SION
The acquisition of three-dimensional deformation fields requires the GNSS and InSAR observations to be appropriately weighted. The HVCE-BPNN method can balance the contributions of different grouped data and provide a posteriori variance for each pixel. In the following, we provide the deduction of HVCE-BPNN for the multisource data fusion method.
Assume that there are groups of observations that are not related. The indirect adjustment model can be expressed as: where is the vector of observation residuals; is the design matrix; ̂ is the vector of the deformation, and represents the vector of observations. If the variances of the observations are known, a weighted least square (WLS) adjustment can be performed to determine the vector of the deformation ̂= ( (1/ 1 2 ⋯ 1/ 2 ) is the weight matrix. Moreover, the correlation between the observations was not considered in this case.
According to (3), we can see that the weight matrix is vital for data fusion. However, as mentioned above, it is difficult to accurately determine the variances of the GNSS and InSAR measurements. HVCE was first proposed by Helmert to iteratively determine the variances or weight matrices of the observa-tions [59]. In applying this method, the observations were first classified into several groups according to their statistical properties. In this case, as the accuracies of the U-D and horizontal directions were significantly different, the GNSS measurements were divided into two groups. Furthermore, the InSAR measurements, which were regarded as homogeneous, were categorized into one group.
In the current study, we assumed that 1 , 2 , and 3 are the three groups of observations, and 1 , 2 , and 3 are the corresponding weight matrices. Then, we obtain the first solution . In the first adjustment, 1 , 2 , and 3 can be assigned arbitrarily positive values, such as identity matrices. The relationship between the variance components and observation residuals can be written as: The variance components are then used to determine the new weight matrices where is an arbitrary constant. The new weight matrices ̂ thus determined are then used in (4) to update the WLS solution. The process is repeated until the following condition is satisfied: (10) However, the weight matrices ̂ from the final iteration by HVCE are not always satisfied, which may result in negative variances [60]. Based on this, an improved method using a back-propagation neural network is proposed to avoid the occurrence of negative variance in HVCE.
The BPNN is a multilayer feed-forward neural network and an error inverse propagation learning algorithm. This algorithm is the most widely used, intuitive, and easy to understand in artificial neural networks [61]. In this section, a BPNN with a three-layer structure is designed. As shown in Fig. 4, the three-layer BPNN consists of one input layer, one hidden layer, and one output layer. Let , , and denote the unit numbers of the input layer, hidden layer, and output layer, respectively. In the given task, five sets of observations (GNSS deformation in three directions, ascending and descending InSAR deformation in the LOS direction) were defined as inputs, and the outputs were the three-dimensional deformation results. The hidden layer unit number, , was set up according to the following empirical equation: = √ + + (11) where represents a constant belonging in the range [1,10].
In designing a BPNN, the selection of the activation function should not be ignored. This study used the Purelin function expressed as follows: ( ) = (12) where represents the weight input of a neural unit.
Meanwhile, the choice of the BPNN learning algorithm is related to the efficiency of network learning [62]. This study uses the Levenberg-Marquardt (LM) algorithm, and the weight increment ∆ is expressed as follows: where represents the identity matrix; represents the user-defined learning rate, and ( ) represents the Jacobian matrix.

C. DATA PROCESSING
The observation data were obtained by the GNSS static network observation mode under the WGS-84 coordinate system, and the geodetic coordinates and geodetic heights were generated using Hi-Target Geomatics Office (HGO) software. The SP3 postaccuracy ephemeris data released by the International GNSS Service (IGS) are used to provide satellite coordinate information to maximize the accuracy of the solution in the baseline solution step. The small baseline subsets (SBAS) algorithm is a time-series InSAR method based on D-InSAR that overcomes the traditional InSAR problems of spatiotemporal decorrelation and atmospheric interference [63]. In this study, the ascending and descending surface deformations in the LOS direction were obtained using the SBAS-InSAR method. To minimize the orbital error in the interferogram, the orbital parameters of the SLC image data are provided with the postaccurate orbit data (POD) provided by the ESA. Moreover, the SAR images were resampled from 5 m × 20 m (range resolution × azimuth resolution) to a ground surface resolution of 20 m × 20 m using multilook processing. The 30 m SRTM data were then used as the external digital ele-vation model (DEM) for the processing of SBAS-InSAR [64].
It can be seen from the above that the single-point GNSS data and the InSAR surface monitoring data are different types of data. The single-point GNSS data must be converted into the surface, which can be fused with the surface monitoring data of InSAR. Thus, a kriging interpolation of the semivariogram is used to interpolate GNSS points into a surface with the same pixel resolution as the InSAR measurements for data fusion, and the east, north, and up are positive directions, as shown in Fig. 5(a, b, c). The LOS deformation fields of the surface of the study area obtained by the ascending and descending SAR are presented in Fig.  6(a, b).

FIGURE 6. Ascending and descending deformation of InSAR measurement in LOS directions: (a) ascending cumulative deformation of InSAR in LOS direction; (b) descending cumulative deformation of InSAR in LOS direction.
To extract the three-dimensional deformation of the ground surface, the results shown in Figs. 5 and 6 were used as the basis for data fusion. Based on the InSAR geometric transformation and Table 1, the initial error Equation (2) can be expressed as: According to the equations deduced in Section B, iterative weighting is performed to output convergent and divergent pixels using the classical HVCE method. The BPNN established by the convergent pixels was then used to calculate the divergent pixels. In this study, there were 3 080 pixels in each deformation map, of which 2 472 were convergent pixels and 608 were divergent pixels. This paper uses the MATLAB neural network toolbox to construct the BPNN model. 70% of the converged pixels are used as the training set, 15% as the validation set, and 15% as the test set. The training progress is shown in Fig. 7; the best validation performance is shown in Fig. 8, and the fitting results for each kind of dataset are shown in Fig. 9. Figs. 8 and 9 demonstrate the great performance of the BPNN model when fusing GNSS and InSAR measurements. Finally, a detailed data processing flowchart using the HVCE-BPNN method is shown in Fig. 10, and the complete three-dimensional deformation fields of the mining area can be obtained.

IV. RESULTS AND DISCUSSION
To evaluate the accuracy of the HVCE-BPNN method in calculating three-dimensional surface deformation results with GNSS and InSAR, the following four schemes are proposed to study the effects of different methods on the accuracy of three-dimensional deformation. Among them, schemes 1 to 3 are the traditional methods, and scheme 4 is the method proposed in this paper. The specific schemes are as follows.
Scheme 1: LOS deformation of ascending and descending orbits of InSAR are substituted into (1), and the E-W and U-D deformation results are calculated by the method that ignores the N-S deformation [65,66]. Scheme 2: We proposed inferring three-dimensional deformation by combining the N-S deformation result from GNSS interpolation with the InSAR deformation result [67]. The N-S deformation of GNSS and ascending and descending orbit LOS deformation of InSAR are substituted into (1). Moreover, (1) can be expressed as: represents the deformation vector in the N-S direction.
Scheme 3: The three-dimensional deformation of GNSS and ascending and descending orbit LOS deformation of InSAR are used in this scheme. Based on the principle of the LS method, is assigned as the identity matrix, and then the three-dimensional deformation results are calculated. Scheme 4: The same data as scheme 3 are used in this scheme. Different from scheme 3, the HVCE-BPNN method belongs to the posterior weight determination method, and its purpose is to carry out reasonable weighting for different types of data.
To assess the accuracies of the three-dimensional deformation calculated by different schemes, the observed data of 139 GNSS points were considered real values, and the threedimensional deformation results were compared during the 2019.03. 22-2020.06.08 period. The three-dimensional deformation errors are shown in Fig. 11, and RMSEs of different schemes between the three-dimensional deformation results and GNSS points are shown in Table 2. The results of scheme 1 are very similar to those of scheme 2 in the E-W direction, as shown in Fig. 11(a). However, the U-D direction shows larger RMSEs in scheme 2 than in scheme 1, which are listed in Table 2. This occurs because the N-S deformation of the GNSS directly replaces the fusion N-S deformation, which ignores the N-S deformation component of the InSAR measurement. In addition, both schemes 3 and 4 have much smaller RMSEs than the others in the E-W and U-D directions because the threedimensional deformation of GNSS and multiorbit deformation of InSAR are considered. However, the improvement in the N-S accuracy of the two schemes is not significant, as InSAR measurements are insensitive to deformation in the N-S direction. The projection vectors of the ascending and descending orbits of the Sentinel-1A satellite in this direction are only -0.154 8 and -0.159 9, respectively.

FIGURE 11. Three-dimensional deformation errors of different schemes: (a) E-W deformation errors of observation points; (b) N-S deformation errors of observation points; (c) U-D deformation errors of observation points.
As expected, scheme 4 has the smallest RMSEs in all three directions, and the HVCE-BPNN method obtained threedimensional deformation results that are more reasonable than the others, although they are quite consistent in general. Compared to the results of scheme 3, improvements of approximately 2.2 mm, 0.5 mm, and 3.5 mm can be achieved by the results of scheme 4 in E-W, N-S, and U-D, respectively. This demonstrates that the HVCE-BPNN method is helpful for reasonable weight determination and improvement of calculation accuracy. Moreover, our method is limited mainly by the density of the GNSS. Because the observation points away from the exploitation range are sparse, the results are mainly dominated by InSAR, and if we can obtain denser GNSS observation points, the accuracy is expected to be further improved.
The HVCE-BPNN method was used to calculate the threedimensional surface deformation, and the E-W direction result obtained by three-dimensional decomposition is presented in Fig. 12(a). The deformation in the N-S direction obtained by the three-dimensional decomposition is presented in Fig. 12(b), and the deformation in the U-D direction obtained by the three-dimensional decomposition is presented in Fig. 12(c). The directions to the east, north, and up were positive. Our experiment demonstrates the advantages of the three-dimensional deformation results, showing that the large deformation is concentrated in the southwest area of the exploitation range. The maximum deformations in the E-W, N-S, and U-D directions were 228 mm, 297 mm, and 190 mm, respectively, during 2019.03. 22-2020.06.08.
In addition, the GNSS and InSAR measurements are fused in the time domain using the HVCE-BPNN method. The period of InSAR was chosen to be 60 days to harmonize in the time domain. Four points were selected in the vicinity of the mined-out area (i.e., 4-4 in the east, 4-7 in the south, 6-8 in the west, and 8-5 in the north). The three-dimensional time series deformation results are shown in Fig. 13. Extracting the three-dimensional results of the time series, we found that substantial deformation occurred with a continuous trend of deformation after August 2019. In Fig. 13(a), it can be seen that 4-4 and 6-8 show a greater trend of variation, while 4-7 and 8-5 show a more stable performance. Conversely, as shown in Fig. 13(b), 4-7 and 8-5 show a greater trend of variation, while 4-4 and 6-8 show a more stable performance. As shown in Fig. 13(c), all four points show a sinking trend in the time series. Among them, the points far from the deformation area have small deformation variables, and those near the deformation area have large deformation variables. When studying the law of surface subsidence caused by mining, the three-dimensional space problem of surface deformation is usually to analyze the deformation along with the strike and dip line points. The distribution of topography is also an important factor affecting the spatial distribution of three-dimensional deformation on the ground surface of the mining area. Therefore, the high-precision topography of the mining area obtained by UAV photogrammetry technology and the deformation points in the strike and dip directions were extracted, as shown in Fig. 14.
The three-dimensional decomposed deformation results were compared with the GNSS data from points in the strike and dip directions during the observation period, as shown in Fig. 15. The E-W deformation is shown in Fig. 15(a) and 15(d). The E-W horizontal movement of most points mainly occurred from west to east. The deformation in the N-S direction is shown in Fig. 15(b) and 15(e). In the N-S direction, the deformation mainly moved from south to north. Namely, the horizontal displacement of the mining area mainly moved from southwest to northeast. The horizontal displacements of the E-W and N-S directions are caused not only by underground mining but also by two other factors. The first factor is that there is an F8 fault near the exploitation range, which destroys the ore integrity and reduces the mechanical properties of the rock mass. The second reason is the influence of terrain. In the horizontal direction, there is a significant difference in elevation on the surface, and the terrain trend is low in the northeast direction but high in the southwest direction, as shown in Fig. 14. Therefore, there is a clear trend of moving from high altitude to low altitude, which represents moving from southwest to northeast. The U-D deformation is displayed in Fig. 15(c) and 15(f). The maximum deformation values in the strike and dip directions are S11 and D5, respectively, instead of the points near the center of the mined-out area. Furthermore, the subsidence value of the southwestern region is larger than that of the northeastern region, which is consistent with the law of mining subsidence.
Overall, with regard to the E-W and N-S deformation results, the western region of the surface subsidence area moves toward the east, and the eastern region of the surface subsidence area moves toward the west, and the southern region of the surface subsidence area moves toward the north, and the northern region of the surface subsidence area moves toward the south. The E-W and N-S deformation values in the center of the surface subsidence area were close to zero. The magnitude of the horizontal deformation from the center to both sides first increased and then decreased, and the center of the surface subsidence area was located above the mined-out area but not directly above the area that is being mined. This is because of the activity of the F8 fault and the hysteresis of surface subsidence caused by mining underground working faces. The three-dimensional deformation results obtained by our method are consistent with the minedout area and are in agreement with the law of mining subsidence.

V. CONCLUSIONS
Accurate information on the three-dimensional deformation of the mined-out area surface is of great significance for the production and construction of a mining area. In general, it is beneficial to fuse GNSS and InSAR observations to obtain the surface deformation results in three directions. This study presents a method that is proposed to optimally determine the weights of GNSS and InSAR observations in threedimensional deformation results. The significance of the proposed method is that the BPNN has a better learning ability, which can improve the HVCE method. The analysis of the 4 schemes shows that our new method can provide more reliable results than the previous fusion methods in threedimensional directions, and the RMSEs of the HVCE-BPNN method in three directions are 21.4 mm, 7.6 mm, and 34.2 mm. The proposed method was validated by analyzing the three-dimensional deformation results of the observation points in the strike and dip lines. The acquisition of threedimensional deformation fields by this method can be utilized both to provide a basis for safe production in local mining areas and to serve as a guide to the population to take appropriate measures in a timely manner to protect personal safety and property.
In addition, the introduction of a learning algorithm in the HVCE-BPNN method also increases the complexity and runtime of the proposed method. Therefore, our method can better serve the study of deformation in mining areas. On the other hand, some methods to optimize the learning algorithms more easily require further investigation.