While promoting the development of industry and human civilization around the world, coal inevitably leads to many environmental problems, one of which is ground subsidence caused by mining. For a long time, the research on ground movement caused by mining mainly focuses on the initial period, active period and weakening period of ground point movement, because the main movement and deformation occur in this period, and it is generally believed that the ground has been basically stable after the weakening period (Vervoort 2016; Cai et al. 2023; Qin and Li 2014). However, more and more studies have shown that the process of ground subsidence after mining can last for months to decades or even hundreds of years (Yi 2016; Duo et al. 2022; Zhang 2017). Many countries have security and environmental problems caused by residual subsidence. For example, more than 20 years after the end of mining in a mine in Shanxi, China, a large range of cracks appeared in the village houses near the old mining area (Li 2014). In the United States, the residual subsidence of mining forms a large area of subsidence depression, resulting in many former ranches flooded (Tajduś et al. 2021). Many coal mines are mined in mountainous areas. For example, about 1/3 coal mines in China are located in mountainous areas; 90% of the main producing areas of Russian coal are located in eastern and western Siberia. The terrain in mountainous areas is complex and changeable, and the influence of residual subsidence caused by mining is far greater than that in plain areas (Sun 2018; Shen 2017; He 2003). Therefore, for the safety of human production and life and the protection of the mining environment, it is necessary to research the residual subsidence of underground mining in mountainous areas. The subsidence development of ground points after mining is shown in Fig. 1.

Ground point subsidence development and subsidence velocity curve

In order to monitor the residual subsidence of the mining area, the traditional method is level survey and traverse measurement. Although the traditional monitoring method has high precision, it has high working intensity and poor real-time data. For the mountainous terrain, it is not only more difficult to measure, but also the measurement accuracy will be greatly decreased. With the development of new technologies, especially the development of remote sensing technology, many scholars have used InSAR technology to monitor the long-term residual subsidence after mining (Modeste et al. 2021; Shi et al. 2020; Liu et al. 2014; Wang 2015; Xu et al. 2017). InSAR technology can monitor the weak subsidence of the ground, and the monitoring time span is large and the precision is high, which makes InSAR have unique advantages in monitoring the micro-deformation of the ground. However, InSAR technology not only has a long satellite observation period, but also has serious time and geometric decorrelation, which reduces the reliability of InSAR monitoring. (Yang et al. 2017). With the development and improvement of GNSS system, GNSS online monitoring technology has been more and more used in ground monitoring of mining areas. GNSS online monitoring technology has high monitoring accuracy and rapid and flexible network deployment. It can monitor the ground subsidence movement of mining areas in real time and can accurately predict potential geological disasters. However, when the GNSS satellite signal is interfered or occluded, the monitoring accuracy and reliability will be greatly decreased (Xue 2012). For the monitoring of residual subsidence in underground mining in mountainous areas, which is mostly located in non-urban areas covered by cultivated land or vegetation, with the characteristics of large terrain fluctuation, complex terrain, high occlusion and strong reflection. The traditional monitoring method is difficult to measure and consumes human and financial resources, which is not suitable for the monitoring of residual subsidence in underground mining in mountainous areas. InSAR and GNSS technologies are also affected by the above characteristics and thus affect the accuracy of monitoring, but GNSS technology can select the location with better observation conditions to arrange monitoring points. Therefore, this research uses GNSS technology to monitor the residual subsidence of mining in mountainous areas online.

The movement of overlying strata and ground caused by mining is a complex space–time process. The residual subsidence of mining area as the dense subsidence of rock mass continues to appear on the ground, and its movement law is consistent with the mining subsidence law during the ground movement period (Wang and Zhang 2005). For the simulation and prediction of residual subsidence in mining area, there are empirical equation method (Wang and Deng 2015; Bai et al. 2009; Guo et al. 2021), prediction theory based on time function (Jia 2018; Li et al. 2022), viscoelastic prediction model based on rock mechanics (Zhu et al. 2014), similar material simulation (Xu et al. 2015), numerical simulation experiment (Yi 2016), simulation prediction based on data processing (Chen 2019; Wang 2020; Zhang and Zhang 2021). Through comparative analysis, it is found that the methods based on rock mechanics and numerical simulation involve many physical and mechanical parameters and are not convenient for application, but they can explain the causes of rock strata and ground movement and fragmentation from the mechanical mechanism. The parameter determination of empirical equation method is easy and practical. The dynamic process prediction method based on time function and data processing has few model parameters and is easy to determine, but it is seriously separated from the actual mining process (Cui et al. 2022). In addition, most of the existing mathematical models of residual subsidence in mining areas are based on plain mining areas, and the maximumed mathematical models hinder the research of residual subsidence in underground mining in mountainous areas.

2.1 GNSS observations and processing

GNSS online monitoring of mining subsidence in the mining area is a real-time continuous monitoring of the subsidence deformation of the ground of the mining area in the time and space domain by using professional GNSS monitoring equipment. Its working principle is carrier phase difference technology. The base station transmits its carrier observation and station coordinate information to the monitoring station in real time through the communication link. The monitoring station receives the carrier phase of the GNSS satellite and the carrier phase from the reference station, and forms the phase difference observation value for real-time processing, which can give centimeter-level or even millimeter-level positioning results in real time, and the user can obtain the mining subsidence information in real time through the client. The principle of GNSS online monitoring is shown in Fig. 2.

Principle of GNSS online monitoring

GNSS observations are inevitably affected by the complex environment of mountainous areas, resulting in a large number of gross errors. The residual subsidence process of mining is a small deformation, so it is necessary to effectively process the original GNSS observations with gross errors. The gross errors of GNSS observations in mountainous areas mainly come from external periodic pulse interference, sensor noise and observation process noise. In order to eliminate the above errors, this paper successively uses median average filtering and Kalman filtering to denoise the original GNSS observations, and then uses the built-in Smooth Spline function in MATLAB to fit the filtered discrete data.

2.1.1 Median average filtering

The median average filtering is to remove a maximum value and a minimum value in the collected subsidence data, and then average the remaining data. It can effectively overcome the fluctuation interference caused by accidental factors and has a good inhibitory effect on periodic interference.

2.1.2 Kalman filter

Kalman filtering is a recursive algorithm for dynamic data processing. It combines the dynamic information of the target with the observation results to suppress the influence of noise and obtain an estimated value. The filtering body is divided into two processes: system state prediction and system state update. The GNSS observations of residual subsidence will be updated in real time, and Kalman filter conforms to this characteristic (Chen 2019). Therefore, Kalman filter can minimize the interference of noise to obtain more accurate values.

1. 1.

   System state prediction equation

The system state prediction equation is shown in Eq. (1).

where \(X(t)\) is the system state at time \(t\), \(A\) is the state transition matrix, \(B\) is the system model parameters, \(U(t)\) is the control quantity of the system at \(t\) time, \(W(t)\) is the process noise.

1. 2.

   System state update equation

The system state update equation is shown in Eq. (2).

where \(P(t)\) is the measured value at time \(t\), \(H\) is the parameter of the measurement system, \(V(t)\) is the observation noise.

2.1.3 Smooth spline function

Smoothing Spline is used for function fitting of discrete data after filtering. The fitting principle is to find a fitting function \(f(t)\) for all functions with second-order continuous derivatives, which can minimize the value of Eq. (3) and \(\text{RSS}\) can be understood as a penalty coefficient.

where \(\lambda\) is a constant called the smoothing coefficient.

2.2 Residual subsidence time series model in underground mining in mountainous areas

In order to build a time series model of residual subsidence in underground mining in mountainous areas to describe or predict the residual subsidence of single ground point, at least the relationship between the residual subsidence of ground point and time in mountain area is needed, which is the core idea behind the method. The probability integral method is one of the empirical equation methods for residual subsidence. The Kelvin model of rock mechanics can well reflect the creep process of rock during residual subsidence (Zhu et al. 2014). For this reason, this section constructs the relationship between residual subsidence of single ground point in mountainous area based on probability integral method, and then uses Kelvin model to construct the relationship between residual subsidence and time.

2.2.1 Calculation of residual subsidence at arbitrary point of ground in underground mining in mountainous areas

The probability integral method believes that the ground movement caused by mining is a random event, and the profile equation of the subsidence depression caused by mining can be expressed as an integral equation of the probability density function (Hu 2019). The same coordinate system as the probability integral method is used for the ground movement and subsidence in the mountain area. The essence of the calculation of ground movement and subsidence in mountainous or hilly areas is to add the influence of mining slip caused by topographic changes on the basis of the flat probability integral method. In the case of known topography, geology and mining conditions, the subsidence of arbitrary point on the ground of mountainous or hilly areas within the mining influence range can be calculated (He 2003).

The line passing through the inflection point of the downhill boundary of the inclined main section of the mined-out area and parallel to the strike is the x-axis, and the line passing through the inflection point of the left boundary of the main section of the mined-out area and parallel to the inclined direction is the y-axis. The subsidence calculation equation of arbitrary point in the ground movement period in the mountain area \(W^{\prime}(x,y)\) is as follows (He 2003):

where \(W(x,y)\) is the subsidence of arbitrary point \((x,y)\) on the plain land under the same geological and mining technical conditions. \(\alpha^{\prime}_{{_{x,y} }}\) is the dip angle of point \((x,y)\) on the terrain trend ground. \(\varphi\) is the ground tilt direction angle of point \((x,y)\).\(D_{x,y}\) is the ground characteristic coefficient of point \((x,y)\). \(P[x]\)、\(P[y]\) are the landslide influence function of point \((x,y)\) on the main section of tendency and inclined. \(A\)、\(P\)、\(T\) are the parameters of the landslide influence function。\(W_{\text{m}}\) is the maximum subsidence under the geological mining conditions. \(l\) is the calculated length of maximumed mining. \(L\) is the calculated length of inclined working face。

Then the calculation equation of arbitrary point \((x,y)\) subsidence \(W(x,y)\) is:

where \(W_{{a_{y} }}\) is the subsidence of \((x,y)\) point at the projection point \(a_{y}\) inclined to the main section. \(C_{\text{xm}}\) is the full mining coefficient of the main section.

The ground movement caused by underground mining tends to be stable during the residual subsidence period. However, a large number of studies have shown that bulking is a characteristic of the increase in the fracture volume of intact rock mass (Palchik 2015). There are a large number of voids in the mined-out area and the overlying strata due to the bulking of the caving rock, the fracture of the rock strata and the separation of the strata. Under the action of long-term load, it is theoretically possible to be compacted and to the original rock state, that is, there is a ground maximum residual subsidence coefficient \(q_{\text{c}} = 1 - q\).However, the possibility of caving broken rock being compressed back to the original state is very small. Therefore, the ground maximum residual subsidence coefficient can be further expressed as (Cui et al. 2020, 2022):

where \(q\) is the subsidence coefficient; \(n\) is a coefficient related to the compaction degree of caving broken rock, referred to as rock coefficient.\(0 < n < 1\), when \(n = 1\), the caving broken rock is compressed to the original rock state.

Thus, the maximum residual subsidence \(W_{{\text{c}}} (x,y)\) at arbitrary point \((x,y)\) on the ground is as follows:

Mining residual subsidence is the continuation of traditional mining subsidence. Substituting Eq. (11) to Eq. (4), The maximum residual subsidence \(W^{\prime}_{\text{c}} (x,y)\) of arbitrary point in underground mining in mountainous areas is as follows,

2.2.2 Residual subsidence time series process of arbitrary point in underground mining in mountainous areas

The Kelvin model in rock mechanics can well reflect the creep process of rock in the residual subsidence stage of mining (Zhu 2014). In the residual subsidence stage, with the increase of time, the broken rock mass and cavity in the mined-out area are compacted to stability, and the residual subsidence tends to the maximum residual subsidence \(W^{\prime}_{\text{c}} (x,y)\). Therefore, the residual subsidence time series function of arbitrary ground point \((x,y)\) is:

where \(\lambda\) is the influence parameter of creep equation.

2.3 Inversion of model parameter

Suppose that GNSS observations \(W_{{\text{c}_{i} }} (i = 1,2, \cdots ,m)\) of residual subsidence at arbitrary point m days are obtained within time \((t_{0} ,t_{1} , \cdots ,t_{n} )\). Under the condition that the geological conditions of the mining area are known, \(m\) equations can be obtained according to the functional relationship between the residual subsidence of the ground point and the time. Let rock coefficient and creep equation influence parameter \(P = [n,\lambda ]\) as model parameters, according to the criterion of minimizing the sum of squared errors, the process of inverting the optimal parameter \(P\) can be expressed as an unconstrained optimization problem.

where \(\text{SSE}\) is the sum of squared errors, \(\hat{W}_{\text{c}}\) is the residual subsidence model value of ground point \(i\).

Establish fitness function,

The PSO (particle swarm optimization) algorithm is easy to implement, has high precision and fast convergence speed, can find the optimal value of the parameters, and can reduce the parameter error of the predicted model. Therefore, the PSO algorithm is proposed to invert the optimal parameters. Considering that the algorithm is sensitive to the initial data and parameters, it is easy to fall into the local optimum. Therefore, this paper increases the randomness method so that the particles have more opportunities to jump out of the local optimum, that is, on the basis of the lithology classification of overlying strata above the working face and the provisions of the literature, multiple different initial positions are used to run the algorithm to find the global optimal solution of the parameters.

2.4 Prediction of residual subsidence in underground mining in mountainous areas

After the optimal model parameter \(P\) of the residual subsidence time series model of underground mining in mountainous areas at the ground point is obtained by the least squares inversion optimized by the PSO algorithm, the residual subsidence of the period \((t_{n - m} , \cdots ,t_{n} , \cdots )\) can be predicted according to the Eq. (16).

Similarly, the residual subsidence of the remaining ground points can also be calculated by the above method. Figure 3 is the flow chart of the residual subsidence model.

Establishment process of residual subsidence model of arbitrary point on the ground of underground mining in mountainous areas based on GNSS online observation

3.1 Study area and GNSS observations acquisition

The study area is located in the Yangquan mining area in eastern Shanxi Province, China (Fig. 4). The topography of the mining area is complex and crisscrossed by gullies. The ground is mainly covered by farmland and deciduous broad-leaved forest. The terrain in the mining area is high in the northwest, low in the east and south, and the maximum relative height difference can reach 179 m. In this case, the high occlusion and strong reflection characteristics of the mining environment will cause serious interference to GNSS observations.

Location and landform of the study area

There are two working faces (P01 and P02) in the mining area, which are excavated from south to north on March 1st, 2019 and from north to south on June 27th, 2019 respectively, and the mining ends on April 5th, 2020 and March 1st, 2020 respectively. The mining depth is deep (about 340–550 m) and the mining height is large (about 3–8 m). Considering the location of coal mining face and the influence of topography and vegetation, six GNSS monitoring stations and one base station were installed at six different positions on the ground of the region on August 22nd, 2019, and the observations began at the end of August 2019. The positional relationship between these two working faces and seven GNSS stations is shown in Fig. 5. Figure 6 is data acquisition.

GNSS monitoring stations location

Data acquisition base station and monitoring station. a GNSS base station b GNSS monitoring station

The data collected by six GNSS monitoring stations in this experiment are shown in Table 1.

3.2 Data processing

Based on all the data collected by the six GNSS monitoring stations, the relative subsidence data of each station are obtained based on the first observed subsidence data, as shown in Fig. 7. The results show that there are a lot of gross errors in the original GNSS observations, with an average deviation of 47 mm and a maximum deviation of 118 mm.

The original observation subsidence data sequence of GNSS monitoring stations during the experiment

3.2.1 Median average filtering and Kalman filtering

The median average filtering eliminates the occasional pulse interference (such as occasional human interference) in the GNSS original observation subsidence data, and suppresses the periodic interference (such as electric current, temperature). On this basis, Kalman filter further eliminates the interference and influence of process noise and sensor noise. After two kinds of filtering, the absolute error of the data is smaller, the data is more stable, and the average relative deviation is only 5 mm. Taking the observations of the monitoring station P-3 as an example, the subsidence values before and after filtering are compared, and the results are shown in Fig. 8.

Filter denoising processing of GNSS original observation subsidence data. a Comparison before and after filtering b The details of comparison before and after filtering

3.2.2 Smoothing spline function fitting and changing time benchmark

After filtering and denoising, the GNSS observations fluctuate around the true value on average. In order to study the residual subsidence, the Smoothing Spline function of MATLAB is used to fit the subsidence curve closer to the real value according to the discrete data. By averaging the data of one day, the time benchmark is unified to one day, as shown in Fig. 9.

Smooth processing of GNSS observation subsidence data. a Comparison before and after smoothing b The details of the comparison before and after smoothing

3.3 Dynamic residual subsidence time series model of underground mining in mountainous areas based on GNSS online observation

3.3.1 Determination of time series model parameters of residual subsidence in underground mining in mountainous areas

The high-precision DOM (Digital Orthophoto Map) and DEM (Digital Elevation Model) obtained by the UAV (Unmanned Aerial Vehicle) in the early stage of the experiment were used to extract the topography and geomorphology of the six monitoring stations in the study area. According to the processed GNSS observations, the residual subsidence of all monitoring points were extracted on the basis of the definition of residual subsidence, as shown in Eq. (17). Table 2 shows the slope aspect and residual subsidence time information of six monitoring points.

where \(i\) is the subsidence time; \(W_{i}\) is the subsidence at time \(i\); \(k\) is the beginning time of entering the residual subsidence.

From Table 2, it can be seen that the four ground GNSS monitoring stations in the study area has entered the residual subsidence stage, but it is found that it is reactivated by external environmental factors during the monitoring period, and the ground subsidence becomes larger and no longer meets the definition of residual subsidence. Therefore, based on the residual subsidence before activation, this paper studies the residual subsidence stage of ground monitoring points in the study area.

According to the mining conditions of the mining area and the geological data such as the rock movement observation report of the working face, all the parameters in Eq. (13) are determined. Some parameters in the model are shown in Table 3.

More specifically, this experiment uses the observations of four GNSS monitoring stations (P-1, P-3, P-4, P-5) in the mining area during the residual subsidence period to invert the model parameters and verify the accuracy. In order to avoid the influence of systematic errors in the experiment, the same dynamic prediction duration (30 days) is selected for the measured residual subsidence with different durations, that is, the observation value from the beginning of residual subsidence to 30 days before the end of residual subsidence is used to invert the optimal parameter \(P_{i} = [n_{i} ,\lambda_{i} ],i = 1,3,4,5\) of the model in this paper, and the dynamic residual subsidence prediction is carried out from 30 days before the end of residual subsidence to the end of residual subsidence. We set the particle size, the maximum number of iterations, and the maximum particle velocity of the PSO algorithm to 50,100, and 10, respectively. Finally, the optimal parameter results of the four GNSS monitoring stations are shown in Table 4.

3.3.2 Prediction of residual subsidence using time series model of residual subsidence in underground mining in mountainous areas

After determining the model parameters of the time series model of residual subsidence in underground mining in mountainous areas at the above four monitoring points, the Eq. (16) can be used to predict the residual subsidence value of the ground of the monitoring point from the beginning to the end of the residual subsidence period based on the model in this paper. The result is shown in Fig. 10. From the diagram, it can be seen that the residual subsidence prediction curve based on the residual subsidence time series model of underground mining in mountainous areas basically coincides with the measured value curve in the inversion parameter phase and the dynamic prediction phase, which is in line with the mining residual subsidence law of the mining area.

The comparison between the predicted value of the GNSS monitoring station based on the model in this paper and the measured value in the inversion parameter phase and the dynamic prediction phase. The black dotted line indicates the beginning period of residual subsidence, and the orange dotted line indicates the prediction period of residual subsidence. a P-1 monitoring station b P-3 monitoring station c P-4 monitoring station d P-5 monitoring station

3.4 Accuracy assessment of residual subsidence prediction

The average error and median error of the predicted values of the four stations in the study area based on the residual subsidence time series model of underground mining in mountainous areas in the inversion parameter phase, the dynamic prediction phase and the measured values are shown in Table 5.

From the Table 5, it can be seen that the average error of the residual subsidence prediction results of the four monitoring station observation data based on the residual subsidence time series model of underground mining in mountainous areas established in this paper is small in the inversion parameter phase and the dynamic prediction phase. The maximum average error in the inversion parameter phase is 5.3 mm, and the minimum average error is 2.0 mm; in the dynamic prediction phase, the maximum average error is 4.3 mm and the minimum average error is 0.1 mm. In addition, the maximum error of the predicted value in the inversion parameter phase is 6.5 mm, and the minimum value is 2.3 mm; in the prediction dynamic phase, the maximum error is 4.6 mm, and the minimum is 0.2 mm.

In addition, it can be seen from Tables 2 and 5 that the larger the ratio of the inversion parameter period to the dynamic prediction stage of the residual subsidence prediction process, the higher the prediction accuracy of the residual subsidence. On the contrary, the smaller the ratio is, the lower the prediction accuracy of residual subsidence is.

It shows that the time series model of residual subsidence in underground mining in mountainous areas established in this experiment has good adaptability, high correlation of nonlinear fitting, and can realize high precision prediction of residual subsidence time series at any point on the ground of mountain area.

In order to demonstrate that the residual subsidence time series model of underground mining in mountainous areas proposed in this paper has strong applicability, the data of the study area in Reference (Cui et al. 2020) (Lingxin mining area in Ningxia, China) are used for experimental analysis. According to Reference (Cui et al. 2020), the residual subsidence observed by the GNSS LX 2–4 in the study area was 24 mm in the first year (December 29, 2015-December 29, 2016), 37 mm in the second year (December 29, 2016-December 29, 2017), and 22 mm in the third year (December 29, 2017-July 28, 2018). According to this, the residual subsidence of the ground point is shown in Table 6.

Combined with the literature and the observation data of working face rock movement in Lingxin mining area of Ningxia, all the known parameters in Eq. (13) are determined. Some parameters in the model are shown in Table 7.

According to the functional relationship between the residual subsidence of the ground point and the time, three equations can be obtained. At this time, only two parameters \(n\) and \(\lambda\) (rock coefficient, creep equation influence parameters) need to be determined. Therefore, it is feasible to optimize the model parameters by PSO based on three sets of residual subsidence. Similarly, we set the particle size, the maximum number of iterations, and the maximum particle velocity of the PSO algorithm to 50,100, and 10, respectively. Finally, the optimal model parameter \(P = [n,\lambda ] = [0.5,0.0008]\) of the monitoring point LX 2–4, according to Eq. (16), the results of the residual subsidence prediction value in the inversion parameter phase are shown in Table 8.

The residual values of the model values and the measured values of the LX 2–4 are shown in Table 9.

The results show that the residual subsidence prediction value of the monitoring point LX 2–4 is not much different from the actual value in the inversion parameter phase, and the fitting effect is preferable. The average absolute value of the residual is 7.3 mm, and the maximum residual is 12 mm in the first year. However, it is impossible to predict the dynamic residual subsidence of the mining area. The main reason is that there are too few data samples and lack of original data of literature experiments. In addition, the terrain of the mining area is flat, and the model is still applicable. Only the ground characteristic coefficient (\(D\)) in Eq. (12) is assigned to 0, which proves that the model is applicable to both mountain and plain mining areas.

1. (1)

   Due to the complex and changeable terrain and dense vegetation in mountainous areas, it is more and more difficult to monitor the residual subsidence in underground mining in mountainous areas. The leveling measurement is time-consuming and labor-intensive. The InSAR is affected by the characteristics of high occlusion and strong reflection due to the undulating terrain and dense vegetation cover in mountainous areas, and its monitoring accuracy is greatly reduced. The GNSS online monitoring technology has the characteristics of high monitoring accuracy, rapid and flexible network layout and strong applicability. This paper comprehensively compares the advantages and disadvantages of the above monitoring technologies, and selects GNSS online monitoring technology as the optimal means to monitor the residual subsidence of underground mining in mountainous areas.

2. (2)

   Most studies have not explained the processing of GNSS original observations in detail. The gross error elimination of GNSS original observations is an indispensable step. The gross error of GNSS monitoring in underground mining in mountainous areass is more serious. This paper proposes a method of using Median average filtering, Kalman filtering and Smoothing Spline function, and describes in detail. In addition, experiments show that this method can effectively eliminate errors in GNSS original observations.

3. (3)

   Cui et al. (2020) proposed the coefficient of residual subsidence in mining area related to the compaction degree of caving broken rock. Zhu et al. (2014) obtained the function of residual subsidence and time based on the Kelvin model of rock mechanics. The residual subsidence process is the continuation of the traditional mining subsidence. In view of the blank of the time series model of residual subsidence in underground mining in mountainous areas, based on the subsidence model of underground mining in mountainous areas proposed by He. (2003), this paper introduces the correlation coefficient of compaction degree of caving broken rock and the Kelvin model function of residual subsidence and time, and constructs the residual subsidence time series model of arbitrary point in underground mining in mountainous areas. The GNSS observations of Yangquan and Ningxia mining areas in China prove that the model has high predicted accuracy and strong applicability.

4. (4)

   The traditional method of least squares inversion model parameters is often difficult to achieve and easy to fall into local circulation. In this paper, the PSO algorithm is used to optimize the least squares to invert the optimal parameters. The algorithm is simple and has the ability of global search, and the model parameters obtained are more accurate.

1. (1)

   Due to its undulating terrain and high occlusion and strong reflection environment in mountainous areas, GNSS monitoring equipment will be affected by uncontrollable environmental and human factors during the observation process. Therefore, compared with the flat mining area, the GNSS original observations in the mountainous area has more gross errors. Through median average filtering, Kalman filtering and Smoothing Spline function, the denoising processing of GNSS original observations in mountainous areas can be effectively realized.

2. (2)

   Based on the probability integral method and the coefficient of residual subsidence stage related to the degree of compaction of caving broken rock, the residual subsidence relationship of single ground point in mountainous area is constructed. Based on the Kelvin model of rock mechanics, the relationship between residual subsidence and time is studied, and the time series model of residual subsidence in mountainous area is proposed.

3. (3)

   Through the GNSS observations of the residual subsidence period of the four GNSS monitoring stations in the study area, the residual subsidence time series model of underground mining in mountainous areas is applied, and the PSO algorithm is used to optimize the least squares fitting to find the optimal model parameters. The residual subsidence of the four monitoring points was predicted, including a 30-day dynamic residual subsidence prediction, and the accuracy was evaluated. The results show that the time series model of residual subsidence in underground mining in mountainous areas can accurately describe the residual subsidence of mining area under mountain terrain, and can reasonably and accurately predict the residual subsidence.

4. (4)

   It is found that the model is also applicable to the case of flat mining area, which proves that the model has high predicted accuracy and strong applicability. However, due to the lack of original data, the number of samples in the verification experiment is small and the dynamic residual subsidence prediction cannot be carried out. This is a regret of this research.