# An Improved Knothe Time Function Model in the Prediction of Ground Mining Subsidence by Using the Kalman Filter Method

Authors: Shifang Zhang; Jin Zhang
Journal: Geofluids
 (2022)

---

## Abstract

The Knothe time function is a classical method in predicting the ground mining subsidence. Nevertheless, it does not take the observation data into account in the prediction process. The Kalman filter method can solve this issue at large extent. Taking a coal mining work face of Xishan Coalfield as an example, this research compares the performance of the traditional Knothe time function and that of the improved Knothe time function by using the Kalman filter method. The comparison results show that through an improvement by using the Kalman filter method, theRMSE is improved from 133.4 mm to 78.3 mm; ME, from 91.9 mm to 3.1 mm; and the relative error, from 8.1% to 5.7%. Meanwhile, the improved model has good astringency. These verify that the improved model has higher accuracy and reliability. Hence, this research presents an effective method in predicting ground subsidence of mining area by improving the Knothe time function using the Kalman filter method.


---

## Body

## 1. Introduction

Ground subsidence due to coal mining is a complicated four-dimensional time-space process [1, 2]. This process leads to various forms of damage to subsurface facilities, such as buildings, roads, and water facilities [3]. Moreover, it relates a series of problems including environmental deterioration, land desertification, and economic loss [4]. Consequently, it has both theory and theoretical values for accurately predicting the ground subsidence process due to coal mining, which can eliminate or prevent these damages and problems largely [5, 6].Time function method presented by Knothe in 1953 is a widely used model to predict the subsidence process due to coal mining [7, 8]. As a classical model, it is developed and improved by many researchers in ground subsidence prediction. Cui et al. [1] proposed the time function model parameters in predicting progressive surface subsidence above coal mining. Hu et al. [2, 9] determined the parameters of the time function model through the probability integral method. Taherynia et al. [10] used the Knothe time model and Geertsma influence function to predict the subsidence over oil and gas fields. Zhang et al. [3] presented an improved two-parameter time function model to predict the ground mining subsidence in Heze City of China. Cheng et al. [4] established an improved time function model by analogical reasoning from a phenomenological perspective. Besides the Knothe time function, many other prediction models were also developed by the researches [11–15]. Most of these models can predict the whole subsidence process completely, while they rarely consider the observation results at different phases [16, 17]. Still, these models can provide the state equation for the Kalman filter method.The Kalman filter method is one of the most widely used technologies in the data processing field [18]. Even if the accurate simulation model cannot be acquired fluently, the Kalman filter can still predict the status of the process in the past, at present, and in the future [19]. Due to its powerful ability, the Kalman filter is widely used in many military and civilian fields, such as navigation, positioning, surveying, electron, control, tracking, surface deformation, and sensor network [20–25]. Furthermore, the Kalman filter method can predict the status of the process with incorporating the observation data simultaneously. However, it has difficulty in determining the state equation.In this research, the aim is to present an improved Knothe time function model for the prediction of ground mining subsidence by using the Kalman filter method, which is beneficial for the environment and ecology protection of the mining area.


## 2. Study Area and Data Sources

### 2.1. Study Area

This research selects the ground surface region over a coal mining workface of Tunlan Minefield, Xishan Coalfield in Shanxi Province, as the study area. It is located in the central Shanxi Province and north-central China (Figure1).Figure 1
Location and elevation distribution of the study area.With a geographical location between 112°05′E and 112°07′E and 37°52′N and 37°53′N, the study area has the elevation between 1105.7 m and 1266.4 m. There are a series of observation points along the strike section and incline section to monitor the subsidence process during coal mining activities. The direction of the strike section is from southwest to northeast, as to that of the incline section, from northwest to southeast. Figure 1 shows both the strike section and incline section distributed in undulating topography.Tunlan Minefield is located in the northwest section of Taiyuan Xishan Coalfield and the east wing of Malan Syncline. The workface is in the right flank of lower group of South Fourth Panel of Tunlan Minefield. The mining method of the work face is fully mechanized mining, and the roof management method is full span falling method.


### 2.2. Data Sources

The main data sources used in this research are the coal mining status and observation data.The coal mining status includes the parameters used in ground subsidence prediction, such as influence radius, mining depth, mining thickness, mining progress, subsidence coefficient, and inclination angle.The observation data mainly refers to the observation points along the strike section and incline section. This research mainly predicts the observation points along the strike section. The elevation distribution and the subsidence status of these points are presented in Figure2.Figure 2
Elevation distribution and subsidence status of the observation points along the strike section.Figure2 shows that there are 43 points along the strike section. With a horizontal distance of about 1100 m, these points distribute in the undulating topography. For example, A8 has the lowest elevation value of about 1100 m, whereas A36 has the highest elevation value of about 1225 m.As to the ground subsidence status for these points, some observation data are lost during the monitoring process, such as A14 and A30. The subsidence profile is in accord with the ground subsidence rule due to coal mining generally, although there are some exceptions. For example, the subsidence values have evident undulations, such as for A20 and A43. These exceptions may result from complicated factors including ground surface conditions and coal mining status.


## 2.1. Study Area

This research selects the ground surface region over a coal mining workface of Tunlan Minefield, Xishan Coalfield in Shanxi Province, as the study area. It is located in the central Shanxi Province and north-central China (Figure1).Figure 1
Location and elevation distribution of the study area.With a geographical location between 112°05′E and 112°07′E and 37°52′N and 37°53′N, the study area has the elevation between 1105.7 m and 1266.4 m. There are a series of observation points along the strike section and incline section to monitor the subsidence process during coal mining activities. The direction of the strike section is from southwest to northeast, as to that of the incline section, from northwest to southeast. Figure 1 shows both the strike section and incline section distributed in undulating topography.Tunlan Minefield is located in the northwest section of Taiyuan Xishan Coalfield and the east wing of Malan Syncline. The workface is in the right flank of lower group of South Fourth Panel of Tunlan Minefield. The mining method of the work face is fully mechanized mining, and the roof management method is full span falling method.


## 2.2. Data Sources

The main data sources used in this research are the coal mining status and observation data.The coal mining status includes the parameters used in ground subsidence prediction, such as influence radius, mining depth, mining thickness, mining progress, subsidence coefficient, and inclination angle.The observation data mainly refers to the observation points along the strike section and incline section. This research mainly predicts the observation points along the strike section. The elevation distribution and the subsidence status of these points are presented in Figure2.Figure 2
Elevation distribution and subsidence status of the observation points along the strike section.Figure2 shows that there are 43 points along the strike section. With a horizontal distance of about 1100 m, these points distribute in the undulating topography. For example, A8 has the lowest elevation value of about 1100 m, whereas A36 has the highest elevation value of about 1225 m.As to the ground subsidence status for these points, some observation data are lost during the monitoring process, such as A14 and A30. The subsidence profile is in accord with the ground subsidence rule due to coal mining generally, although there are some exceptions. For example, the subsidence values have evident undulations, such as for A20 and A43. These exceptions may result from complicated factors including ground surface conditions and coal mining status.


## 3. Methodology

The probability integral method can predict the subsidence results to determine the typical points; further, the Knothe time function model is used to provide the state equation, which can be improved by using the Kalman filter method. The final method is the key innovation in this research.


### 3.1. Probability Integral Method

Developed from stochastic medium theory, the probability integral method can predict any point’s subsidence due to the mining activities of the working face [26, 27]. Referencing the subsidence computation equations, the inclination, horizontal movement, horizontal deformation, and curvature can be predicted accordingly. Among the observation points along the strike section, five typical points are selected from these prediction results, as shown in Table 1:Table 1
Prediction results for the five typical points.
PointsSubsidence (mm)Inclination (mm)Horizontal movement (mm/m)Horizontal deformation (1/km)Curvature (mm/m)A20-6286402-0.85-0.013A24-12602.87200-3.17-0.048A28-1419-0.57-26-0.99-0.015A29-1390-0.71-35-0.2-0.003A36-1505-0.162-1.53-0.023In Table1, A20 has the highest inclination and horizontal movement values; A24 has the highest horizontal deformation and curvature values; A28 has the x-coordinate of 219.384 m, about one time of influence radius (216.667 m); A29 has the highest subsidence value through observation (-1519.865 mm); A36 has the highest subsidence value through prediction (-1505 m).The prediction of the subsidence process is conducted for the five typical points in this research.


### 3.2. Knothe Time Function Model

According to the Knothe time function model, the relationship between the ground subsidence valueWt and time t can be defined as follows [1, 3]:
(1)Wt=W0φt,where
(2)φt=1−e−c∗t.In equations (1) and (2), W0 is the maximum subsidence; c is the time influence parameter, which is used to describe the subsidence process in time influenced by geological conditions and mining progress [28–30]. Parameter c can be computed as the following [17]:
(3)c=−vL1ln0.02.In equation (7), v is the advancing rate of the working face; L1 is the critical mining dimension.Through investigation,L1 is 100 m in this research. The mining progress is not at a stable velocity. The propulsion meters are various for different months, which is shown in Table 2.Table 2
Propulsion status of the working face in each month of the mining process.
MonthDay (d)Propulsion meters (m)Advancing rate (m/d)MonthDay (d)Propulsion meters (m)Advancing rate (m/d)2018.073169.62.252019.073173.62.372018.083188.82.862019.083199.23.202018.093096.83.232019.093077.62.592018.103155.21.782019.103181.62.632018.113076.02.532019.113088.02.932018.123164.82.092019.123175.22.432019.013186.42.792020.013164.82.092019.022884.83.032020.022991.23.142019.033179.22.552020.033182.42.662019.0430106.43.552020.043061.62.052019.053188.02.842020.053141.61.342019.063086.42.88In Table2, the advancing rate is computed by dividing days by the propulsion meter. The advancing rate varies remarkably in different months. For example, it is 3.23 m/d in Sep. 2018, while it is 1.78 m/d in Oct. 2018. To simplify the computation process, the advancing rate is regarded as a constant value for the whole mining process. As the total propulsion meters are 1819.2 m, and the total observation period is 763 days. Thus, the whole advancing rate can be computed by dividing the total observation days by the total propulsion meter, which is about 2.4 m/d. Then, the time influence parameter c can be computed, whose value is 0.0215 d-1.In this research, the monitoring of the observation points is conducted at 19 phases from 2018/6/13 to 2020/7/15, about two years. The mining statuses during the 19 phases are computed in Table3 as follows.Table 3
Mining statuses during 19 phases for the observation points.
PhaseDay (d)Propulsion meters (m)Advancing rate (m/d)PhaseDay (d)Propulsion meters (m)Advancing rate (m/d)2018/6/130002019/8/53587.232.492018/7/62313.470.592019/9/53196.133.102018/8/32864.722.312019/10/1540104.152.602018/8/242160.152.862019/11/153186.122.782018/10/2461159.592.622019/12/153080.392.682019/1/1583195.072.352020/1/153170.172.262019/3/1559167.722.842020/4/1591237.852.612019/4/153194.083.032020/5/254064.351.612019/5/153095.783.192020/7/15518.050.162019/7/147134.192.86In Table3, the column “day” is the duration between this phase and the last phase. The column “propulsion meter” is the propulsion progress during these days, which is computed according to the advancing rate in Table 2. The “day” is divided into several pieces to accord with the advancing rate in Table 2, and the propulsion meters are computed in every piece. The final propulsion meters are acquired by adding the pieces of the propulsion meters together. The column “advancing rate” is accrued by dividing “day” by the “propulsion meter” column.From Table3, we can see that the advancing rate has some undulation in different periods. It is low at the early and last phase when the mining activities start and end.Using 19 phases, the working face can be divided into 18 elements. For theith mining element, its mining period is ti and the advancing rate is vi, so the mining length is vi∗ti. For the moment t, the ground subsidence of the mining element from 1 to n can be computed as follows:
(4)W1t=W0φt,W2t−t1=W0φt−t1,W3t−t1−t2=W0φt−t1−t2,⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯..Wnt−t1−t2−⋯−tn−1=W0φt−t1−t2⋯−tn−1.According to the superposition principle, the subsidence value of any point at timet can be expressed as
(5)Wx,t=φtWx−Wx−v1t1+φt−t1Wx−v1t1−Wx−v1t1−v2t2+φt−t1−t2Wx−v1t1−v2t2−Wx−v1t1−v2t2−v3t3+⋯+φt−t1−t2−⋯−tn−1Wx−v1t1−v2t2−⋯−vn−1tn−1−Wx−v1t1−v2t2−⋯−vntn,where
(6)Wx=W02erfπrx+1.Using equations (5) and (6), the ground subsidence process can be predicted accordingly.


### 3.3. Improved Knothe Time Function Model Using the Kalman Filter Method

Based on equation (2) and equation (5), the state equation of the Kalman filter method can be transformed as [31]
(7)Wx,tk=Wx,tk−1+1−e−ctk−tk−1e−ctk−1Wx−Wx−v1t1+e−ctk−1−t1Wx−v1t1−Wx−v1t1−v2t2+e−ctk−1−t1−t2Wx−v1t1−v2t2−Wx−v1t1−v2t2−v3t3+⋯e−ctk−1−t1−t2⋯−tn−1Wx−v1t1−v2t2−⋯−vn−1tn−1−Wx−v1t1−v2t2−⋯−vntn.Using equation (7), the ground subsidence can be predicted at any time for the typical observation points. Meanwhile, the prediction results can be rectified by using the improved method (namely, Kalman filter method after here).The equations for the prediction process are(8)Wk−=A∗Wk−1+B∗Uk,(9)Pk−=A∗Pk−1∗A′+Q.Equation (8) is a simplified format of equation (7); Wk− is the predicted subsidence value at time k; A is the state transition matrix; Wk−1 is the rectified subsidence value at time k−1; Uk and B are the control vector and its coefficient matrix, respectively. In equation (9), Pk− is the priori error matrix at time k; Pk−1 is the posterior error matrix at time k−1; Q is the covariance of the state equation.The equations for the rectified process are(10)Kk=Pk−∗Ak′∗Ak∗Pk−1−∗Ak′+R−1,(11)Wk=Wk−+Kk∗yk−Ak∗Wk−,(12)Pk=I−Kk∗Ak∗Pk−.In equations (10), (11), and (12), Kk is the gain matrix; Ak is the coefficient matrix of the observation equation; R is the covariance of the observation equation; yk is the observation vector; the meanings of other symbols can be acquired from the above equations.According to the observation status,R is 16 mm2; through a series of experiments, Q can be determined as 25 mm2. Based on the Kalman filter method, the predicted and rectified values of the ground subsidence at different phases can be computed using the equations from (8) to (12) in an iterative way.


## 3.1. Probability Integral Method

Developed from stochastic medium theory, the probability integral method can predict any point’s subsidence due to the mining activities of the working face [26, 27]. Referencing the subsidence computation equations, the inclination, horizontal movement, horizontal deformation, and curvature can be predicted accordingly. Among the observation points along the strike section, five typical points are selected from these prediction results, as shown in Table 1:Table 1
Prediction results for the five typical points.
PointsSubsidence (mm)Inclination (mm)Horizontal movement (mm/m)Horizontal deformation (1/km)Curvature (mm/m)A20-6286402-0.85-0.013A24-12602.87200-3.17-0.048A28-1419-0.57-26-0.99-0.015A29-1390-0.71-35-0.2-0.003A36-1505-0.162-1.53-0.023In Table1, A20 has the highest inclination and horizontal movement values; A24 has the highest horizontal deformation and curvature values; A28 has the x-coordinate of 219.384 m, about one time of influence radius (216.667 m); A29 has the highest subsidence value through observation (-1519.865 mm); A36 has the highest subsidence value through prediction (-1505 m).The prediction of the subsidence process is conducted for the five typical points in this research.


## 3.2. Knothe Time Function Model

According to the Knothe time function model, the relationship between the ground subsidence valueWt and time t can be defined as follows [1, 3]:
(1)Wt=W0φt,where
(2)φt=1−e−c∗t.In equations (1) and (2), W0 is the maximum subsidence; c is the time influence parameter, which is used to describe the subsidence process in time influenced by geological conditions and mining progress [28–30]. Parameter c can be computed as the following [17]:
(3)c=−vL1ln0.02.In equation (7), v is the advancing rate of the working face; L1 is the critical mining dimension.Through investigation,L1 is 100 m in this research. The mining progress is not at a stable velocity. The propulsion meters are various for different months, which is shown in Table 2.Table 2
Propulsion status of the working face in each month of the mining process.
MonthDay (d)Propulsion meters (m)Advancing rate (m/d)MonthDay (d)Propulsion meters (m)Advancing rate (m/d)2018.073169.62.252019.073173.62.372018.083188.82.862019.083199.23.202018.093096.83.232019.093077.62.592018.103155.21.782019.103181.62.632018.113076.02.532019.113088.02.932018.123164.82.092019.123175.22.432019.013186.42.792020.013164.82.092019.022884.83.032020.022991.23.142019.033179.22.552020.033182.42.662019.0430106.43.552020.043061.62.052019.053188.02.842020.053141.61.342019.063086.42.88In Table2, the advancing rate is computed by dividing days by the propulsion meter. The advancing rate varies remarkably in different months. For example, it is 3.23 m/d in Sep. 2018, while it is 1.78 m/d in Oct. 2018. To simplify the computation process, the advancing rate is regarded as a constant value for the whole mining process. As the total propulsion meters are 1819.2 m, and the total observation period is 763 days. Thus, the whole advancing rate can be computed by dividing the total observation days by the total propulsion meter, which is about 2.4 m/d. Then, the time influence parameter c can be computed, whose value is 0.0215 d-1.In this research, the monitoring of the observation points is conducted at 19 phases from 2018/6/13 to 2020/7/15, about two years. The mining statuses during the 19 phases are computed in Table3 as follows.Table 3
Mining statuses during 19 phases for the observation points.
PhaseDay (d)Propulsion meters (m)Advancing rate (m/d)PhaseDay (d)Propulsion meters (m)Advancing rate (m/d)2018/6/130002019/8/53587.232.492018/7/62313.470.592019/9/53196.133.102018/8/32864.722.312019/10/1540104.152.602018/8/242160.152.862019/11/153186.122.782018/10/2461159.592.622019/12/153080.392.682019/1/1583195.072.352020/1/153170.172.262019/3/1559167.722.842020/4/1591237.852.612019/4/153194.083.032020/5/254064.351.612019/5/153095.783.192020/7/15518.050.162019/7/147134.192.86In Table3, the column “day” is the duration between this phase and the last phase. The column “propulsion meter” is the propulsion progress during these days, which is computed according to the advancing rate in Table 2. The “day” is divided into several pieces to accord with the advancing rate in Table 2, and the propulsion meters are computed in every piece. The final propulsion meters are acquired by adding the pieces of the propulsion meters together. The column “advancing rate” is accrued by dividing “day” by the “propulsion meter” column.From Table3, we can see that the advancing rate has some undulation in different periods. It is low at the early and last phase when the mining activities start and end.Using 19 phases, the working face can be divided into 18 elements. For theith mining element, its mining period is ti and the advancing rate is vi, so the mining length is vi∗ti. For the moment t, the ground subsidence of the mining element from 1 to n can be computed as follows:
(4)W1t=W0φt,W2t−t1=W0φt−t1,W3t−t1−t2=W0φt−t1−t2,⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯..Wnt−t1−t2−⋯−tn−1=W0φt−t1−t2⋯−tn−1.According to the superposition principle, the subsidence value of any point at timet can be expressed as
(5)Wx,t=φtWx−Wx−v1t1+φt−t1Wx−v1t1−Wx−v1t1−v2t2+φt−t1−t2Wx−v1t1−v2t2−Wx−v1t1−v2t2−v3t3+⋯+φt−t1−t2−⋯−tn−1Wx−v1t1−v2t2−⋯−vn−1tn−1−Wx−v1t1−v2t2−⋯−vntn,where
(6)Wx=W02erfπrx+1.Using equations (5) and (6), the ground subsidence process can be predicted accordingly.


## 3.3. Improved Knothe Time Function Model Using the Kalman Filter Method

Based on equation (2) and equation (5), the state equation of the Kalman filter method can be transformed as [31]
(7)Wx,tk=Wx,tk−1+1−e−ctk−tk−1e−ctk−1Wx−Wx−v1t1+e−ctk−1−t1Wx−v1t1−Wx−v1t1−v2t2+e−ctk−1−t1−t2Wx−v1t1−v2t2−Wx−v1t1−v2t2−v3t3+⋯e−ctk−1−t1−t2⋯−tn−1Wx−v1t1−v2t2−⋯−vn−1tn−1−Wx−v1t1−v2t2−⋯−vntn.Using equation (7), the ground subsidence can be predicted at any time for the typical observation points. Meanwhile, the prediction results can be rectified by using the improved method (namely, Kalman filter method after here).The equations for the prediction process are(8)Wk−=A∗Wk−1+B∗Uk,(9)Pk−=A∗Pk−1∗A′+Q.Equation (8) is a simplified format of equation (7); Wk− is the predicted subsidence value at time k; A is the state transition matrix; Wk−1 is the rectified subsidence value at time k−1; Uk and B are the control vector and its coefficient matrix, respectively. In equation (9), Pk− is the priori error matrix at time k; Pk−1 is the posterior error matrix at time k−1; Q is the covariance of the state equation.The equations for the rectified process are(10)Kk=Pk−∗Ak′∗Ak∗Pk−1−∗Ak′+R−1,(11)Wk=Wk−+Kk∗yk−Ak∗Wk−,(12)Pk=I−Kk∗Ak∗Pk−.In equations (10), (11), and (12), Kk is the gain matrix; Ak is the coefficient matrix of the observation equation; R is the covariance of the observation equation; yk is the observation vector; the meanings of other symbols can be acquired from the above equations.According to the observation status,R is 16 mm2; through a series of experiments, Q can be determined as 25 mm2. Based on the Kalman filter method, the predicted and rectified values of the ground subsidence at different phases can be computed using the equations from (8) to (12) in an iterative way.


## 4. Results

The prediction curves acquired by three different ways are analysed at different phases and for five typical points. Thus, the error curves are analysed accordingly. Finally, the values of the error indexes are computed and compared for the five typical points.


### 4.1. Prediction Results

Based on the time function method and Kalman filter method, the prediction results of the two methods and the rectified results of the Kalman filter method can be computed for the typical observation points at different phases. The distribution curves of these results are shown in Figure3.Figure 3
Prediction curves at different time phase for the typical observation points.
(a)
A20(b)
A24(c)
A28(d)
A29(e)
A36In Figure3, the rectified values of the Kalman filter method are taken as the reference values. Thus, the prediction values have a sharp difference for the five typical observation points. For example, the prediction curves of some observation points are approximate to the rectified curves, such as A28 and A29 points, while other points have evident sharps compared with the rectified curves, such as A20, A24, and A36.In addition, the prediction results acquired from the Kalman filter method have much better astringency than that from the time function methods, especially for the observation points A20, A24, A29, and A36.Generally speaking, the prediction curves from the Kalman filter method are much closer to the rectified curves than that from the time function curves, although occasionally the curves of the time function have better performance, such as that for A28, A29, and A36 at some phases.


### 4.2. Error Distribution

The prediction curves in Figure3 present the prediction quality of the Kalman filter method and time function method by comparing the rectified curves. To quantify the quality of the prediction results of the two methods, the error is computed by taking the rectified values as the real value.Through computing the difference between the predicted values and the rectified values, the error distributions of the two methods for the five observation points at different phases are presented in curve form as shown in Figure4.Figure 4
Error distribution curves at different time phase for the typical observation points.
(a)
A20(b)
A24(c)
A28(d)
A29(e)
A36In Figure4, the “Knothe” refers to the Knothe time function method. Figure 4 shows that the error values of the prediction results by the Kalman filter usually have lower error values than that by time function method, although it has higher absolute values occasionally at a few phases, such as A28 and A29 at 2019/1/13 and A36 at 2019/4/13.In addition, the prediction values by the Kalman filter are gradually approximated to 0 mm at last, whereas those by time function do not have the astringency characteristics. Moreover, the prediction results by time function method have system errors for some points, such as A20, A24, and A29.


### 4.3. Error Index Computation

Figure4 presents the error distribution curves of the Kalman filter method and time function method. To compute the exact error values, the mean error (ME) and root mean square error (RMSE) indexes are adopted in this research. Other error indexes, such as mean absolute error and standard deviation error, have similar characters to ME or RMSE. The ME index represents the difference of the mean values between the predicted and rectified results, which is an important index to evaluate the prediction quality. The RMSE index represents the square root of the ration of the square of the deviation between the predicted value and the rectified value to the number of observations. The RMSE index is sensitive to extralarge or extrasmall errors, so it is good at reflecting the precision of the prediction results.The ME and RMSE values of the two methods are computed for the five typical observation points, which is shown as in Figure5.Figure 5
Error index computation for different observation points.Figure5 shows that all ME values acquired by the Kalman filter are lower than those by the time function method. As to the RMSE, except for A28, the values for the Kalman filter method are much lower than those for the time function method, especially for A20 and A24 points. Except for A29, other ME values are positive. The ME and RMSE values vary largely for different points: A20 has the highest value; then for A24 and A36, their values are at the middle level; A28 and A29 have the lowest values. Therefore, the prediction results have the highest quality for A28 and A29 points, the middle quality for A24 and A36 points, and the worst quality for A20 points.Through computing the average values of error indexes for the five points, theME values for the time function and Kalman filter method are 91.9 mm and 3.1 mm, respectively, and the RMSE values are 133.4 mm and 78.3 mm, respectively. Hence, the prediction results by the Kalman filter have higher accuracy than those by time function method, especially for ME index.


## 4.1. Prediction Results

Based on the time function method and Kalman filter method, the prediction results of the two methods and the rectified results of the Kalman filter method can be computed for the typical observation points at different phases. The distribution curves of these results are shown in Figure3.Figure 3
Prediction curves at different time phase for the typical observation points.
(a)
A20(b)
A24(c)
A28(d)
A29(e)
A36In Figure3, the rectified values of the Kalman filter method are taken as the reference values. Thus, the prediction values have a sharp difference for the five typical observation points. For example, the prediction curves of some observation points are approximate to the rectified curves, such as A28 and A29 points, while other points have evident sharps compared with the rectified curves, such as A20, A24, and A36.In addition, the prediction results acquired from the Kalman filter method have much better astringency than that from the time function methods, especially for the observation points A20, A24, A29, and A36.Generally speaking, the prediction curves from the Kalman filter method are much closer to the rectified curves than that from the time function curves, although occasionally the curves of the time function have better performance, such as that for A28, A29, and A36 at some phases.


## 4.2. Error Distribution

The prediction curves in Figure3 present the prediction quality of the Kalman filter method and time function method by comparing the rectified curves. To quantify the quality of the prediction results of the two methods, the error is computed by taking the rectified values as the real value.Through computing the difference between the predicted values and the rectified values, the error distributions of the two methods for the five observation points at different phases are presented in curve form as shown in Figure4.Figure 4
Error distribution curves at different time phase for the typical observation points.
(a)
A20(b)
A24(c)
A28(d)
A29(e)
A36In Figure4, the “Knothe” refers to the Knothe time function method. Figure 4 shows that the error values of the prediction results by the Kalman filter usually have lower error values than that by time function method, although it has higher absolute values occasionally at a few phases, such as A28 and A29 at 2019/1/13 and A36 at 2019/4/13.In addition, the prediction values by the Kalman filter are gradually approximated to 0 mm at last, whereas those by time function do not have the astringency characteristics. Moreover, the prediction results by time function method have system errors for some points, such as A20, A24, and A29.


## 4.3. Error Index Computation

Figure4 presents the error distribution curves of the Kalman filter method and time function method. To compute the exact error values, the mean error (ME) and root mean square error (RMSE) indexes are adopted in this research. Other error indexes, such as mean absolute error and standard deviation error, have similar characters to ME or RMSE. The ME index represents the difference of the mean values between the predicted and rectified results, which is an important index to evaluate the prediction quality. The RMSE index represents the square root of the ration of the square of the deviation between the predicted value and the rectified value to the number of observations. The RMSE index is sensitive to extralarge or extrasmall errors, so it is good at reflecting the precision of the prediction results.The ME and RMSE values of the two methods are computed for the five typical observation points, which is shown as in Figure5.Figure 5
Error index computation for different observation points.Figure5 shows that all ME values acquired by the Kalman filter are lower than those by the time function method. As to the RMSE, except for A28, the values for the Kalman filter method are much lower than those for the time function method, especially for A20 and A24 points. Except for A29, other ME values are positive. The ME and RMSE values vary largely for different points: A20 has the highest value; then for A24 and A36, their values are at the middle level; A28 and A29 have the lowest values. Therefore, the prediction results have the highest quality for A28 and A29 points, the middle quality for A24 and A36 points, and the worst quality for A20 points.Through computing the average values of error indexes for the five points, theME values for the time function and Kalman filter method are 91.9 mm and 3.1 mm, respectively, and the RMSE values are 133.4 mm and 78.3 mm, respectively. Hence, the prediction results by the Kalman filter have higher accuracy than those by time function method, especially for ME index.


## 5. Discussion

This research presents a new model to predict subsidence due to coal mining by incorporating the Kalman filter into the Knothe time function method. Through comparison, the new model has higher accuracy than the time function method. TheME is improved from 91.9 mm to 3.1 mm, and the RMSE improved from 133.4 mm to 78.3 mm. Moreover, the error is approximate to 0 mm at last, which has good astringency characteristics. As to the Knothe time function model, it may have system errors for some points, such as A20, A24, and A29. Hence, the new model presented in this research improves the performance of the time function largely, which has significance in predicting ground subsidence due to coal mining.Compared to previous similar researches, Cui et al. [1] deemed that the average relative error is 8% by using time function method; Zhang et al. [3] improved the Knothe time function by considering the velocity, acceleration, and the proper value of c, and the average relative standard error was improved from 23% to 4.8%; Cheng et al. [4] improved the Knothe time function model by analogical reasoning, and the improved model has average relative standard deviation of 4.9%, far lower than 23.1% of Knothe time function. In this research, the time function is improved by incorporating the Kalman filter, and the relative error is improved from 12.2% to 7.2%, the ME from 91.9 mm to 3.1 mm, and the RMSE from 133.4 mm to 78.3 mm. The typical point A20 has ground rise at first, which violates the time function rule. It may be due to undulating topography. By only computing the other four points, the relative error will be improved from 8.1% to 5.7%. The performance of the time function is similar to that in Cui et al.’s research [1]. These aspects prove that the improved model has higher accuracy and reliability.This research conducts subsidence prediction by the two methods by regarding the surface as flat topography. In fact, it is undulating. The undulating topography can be incorporated in future researches. Besides the ground subsidence, there are other deformation factors, such as inclination, horizontal movement, horizontal deformation, and curvature. These factors can also be used in predicting the mining progress. Moreover, there are other Kalman filter methods, such as extended Kalman filter [32, 33], unscented Kalman filter [34–36], and total Kalman filter [37, 38]. These methods can be widely used in prediction models.


## 6. Conclusions

Taking a coal mining work face of Xishan Coalfield in Shanxi Province of China as an example area, this research presents an effective method in predicting ground subsidence of mining area by improving the Knothe time function using the Kalman filter method. Meanwhile, the performance of the traditional Knothe time function and its improved version is compared in different aspects. The research results can be concluded as follows:(1)
Incorporating probability integral method, time function method, and Kalman filter method, an improved method is presented in predicting ground subsidence due to coal mining. Probability integral method is adopted to predict the final subsidence status and select typical observation points; the time function method provides the state equations for the Kalman filter method; based on the probability integral method and time function method, the Kalman filter method predicts the ground subsidence by assimilating the observation values. This research presents an effective method in predicting ground subsidence of mining area by improving the Knothe time function using the Kalman filter method(2)
The performance of the two prediction methods is compared through prediction curves, error curves, and error indexes of the five typical observation points. Through comparison, the prediction results by the Kalman filter method have remarkably higher accuracy than those by the time function method. In addition, the prediction results by the Kalman filter are gradually approximated to 0 mm at last, whereas those by time function do not have this astringency characteristic. Moreover, the prediction results by the time function method have systematic errors for some points. Thus, the prediction results by the Kalman filter method have higher reliability



---

*Source: 3567808-2022-06-20.xml*