- Chn Energy Menxi Coal-Chemical Co Ltd., Ordos, China

The probability integral method is one of the most widely used methods for predicting surface subsidence induced by underground mining in China. In its parameter calculation, the least square algorithm is commonly employed for fitting parameters. However, in the process of fitting parameters, the results are easily affected by ill-conditioned normal matrices and the interference of outliers, resulting in divergent problems. To solve these problems, the principle of robust ridge estimation was introduced in this paper, and a parameter calculation model for the probability integral method based on this principle (hereafter referred to as the established model) was established. Besides, a parameter calculation experiment with manual intervention was conducted in combination with engineering examples. The results demonstrate that the parameter calculation method based on robust ridge estimation can suppress the interference of outliers, overcome the problem of ill-conditioned matrix, and ensure the effectiveness and reliability of parameter estimation results. Compared with the conventional least squares method, the robust ridge estimation method demonstrates greater accuracy in predicting surface subsidence parameters, which validates its rationality and accuracy in underground mining engineering. The research findings provide technical support for obtaining similar parameters for surface subsidence in mining areas and hold significant engineering application value.

## 1 Introduction

China is the world’s largest coal producer. Extensive underground coal mining is prone to causing large-scale surface subsidence, leading to problems such as vegetation damage, soil degradation, and lowered underground water level in the coal mine subsidence areas, which has brought considerable damage to the living environment in these areas. To address these issues, it is necessary to predict the extent and scope of surface damage caused by underground mining in advance, and then adopt targeted mining technology with low ecological damage. The work of predicting surface damage in advance is known as mining-induced subsidence prediction (MSP).

At present, scholars at home and abroad have conducted extensive research on MSP and have put forward numerous MSP models. In 1954, Litwiniszyn (Litwiniszy, 1956) introduced the stochastic medium theory which considered rock mass displacement as the stochastic movement of countless small unit particles. The theory laid the foundation for modern MSP methods. In 1993, Shu and Bhattacharyya (Shu and Bhattacharyya, 1993) established an empirical prediction model by investigating the relationship between underground strata and surface subsidence movements. The model is presented in the form of graphs and tables, which can be used to estimate the maximum subsidence, tilting, and horizontal deformation of the surface. In 1995, by analyzing the surface subsidence induced by the mining of shallow thick coal seams in Raniganj coalfields in India, Singh and Yadav (Prediction of subsidence due to, 1996) proposed a viscoelastic model for surface MSP and validated the feasibility of this model through the mining conditions at two other local coal mines. In 2003, Ambrožič et al. (Ambrožič and Turk, 2003) employed artificial neural networks for MSP. This method is based on a substantial number of observation data and uses variables to represent the data of surface subsidence in the process of prediction. In 2014, Ren et al. (Ren et al., 2014) first proposed the generalized influence function method that utilizes computer simulation. The method can not only be applied to more complex engineering geological conditions but also effectively express its stress-strain relationship during mining. In 1959, Liu and Liu (China University of Mining and Technology, 1981) first translated and introduced the stochastic medium theory into China, bringing advanced technology and mature experience to MSP in the country. In 1963, Zhou and Yu (Zhang, 2010) established the negative exponential profile function of the subsidence basin based on the analysis of numerous measured data. In 1965, the probability integral method was first proposed, and later Chinese scholars, including Baoshen Liu and Guohua Liao, made indelible contributions by introducing the method to the domestic context, thus laying a solid theoretical foundation for the quantitative calculation of MSP in China.

MSP models can be broadly categorized into three types, i.e., curve prediction methods, influence function methods, and profile function methods (He et al., 1991; Guo, 2019). Among them, the probability integral method, which is most widely applied in China, is a typical influence function method. The probability integral method, based on the stochastic medium theory, gets its name from the inclusion of the probability integral in the prediction formula of movement and deformation. It is currently an important method used in China for predicting surface movement and deformation of mining-induced subsidence and is one of the methods for MSP as specified in the reference (State Administration of Work Safety, 2017). The accuracy of MSP methods primarily depends on the selection of prediction parameters. The least square fitting algorithm that boasts simple calculation is usually used for calculating the prediction parameters, and the accuracy of the calculated parameters can basically meet the engineering requirements. The least squares estimation possesses desirable properties in parameter estimation, and when the error follows the normal distribution, it is unbiased, consistent, and effective among all unbiased estimation classes. However, it has two problems: first, in the presence of many independent variables, including approximately linear related variables, the parameter value it estimated deviates notably from the true value; second, when the observed value is contrary to the normal distribution assumption and outliers are present in the data, the least squares estimation can be interfered, and the deviation of a single observed value may have a significant impact on the parameters (Sui, 1994; Wu, 2009; Wang et al., 2012).

Regarding the first problem, Shu and Bhattacharyya (Wang et al., 2012) proposed ridge estimation, which is a kind of biased estimation of compressibility designed to reduce the mean square error. It can improve the ill-condition of the normal matrix and stabilize the parameter solution. As for the second problem, Ambrožič et al. and Ren et al. (Wu, 2009; Dong-Sheng et al., 2023) introduced robust estimation, which, by selecting the appropriate equivalent right, helps overcome the difficulty of parameter calculation posed by model bias and the presence of outliers. However, these two methods can only address one problem, not both simultaneously. By combining the advantages of ridge estimation and robust estimation, the principle of robust ridge estimation proposed in this paper can resist both the influence of ill-conditioned normal matrices on the parameter calculation results, and the impact of outliers or gross errors on them, ensuring the validity and reliability of the results. It can provide technical support for obtaining estimated parameters of surface subsidence in similar mining areas.

## 2 Robust ridge estimation

### 2.1 Principle of robust ridge estimation

The introduction of robust ridge estimation is in response to the fact that, when outliers are present in the observed values, and the coefficient matrix *A* of the equation shows an ill-conditioned tendency, neither *LS* estimation nor *LS* ridge estimation can deal with this ill-conditioned result (Sui, 1994; Zhou et al., 2020). Thus, it is necessary to mitigate the ill-condition of the normal matrix through ridge estimation. The fundamental idea of ridge estimation is to exchange an appropriate increase in bias for a significant reduction in variance, thereby reducing the error and improving the calculation accuracy of the sample. The basic calculation method of robust ridge estimation is as follows:

Let the observation equation be:

where *L* is the *n*-dimensional observation vector; *A* is the *X* is the *t*-dimensional parameter vector; and *n*-dimensional error vector,

The corresponding error equation is:

The robust ridge estimation solution of the parameter, obtained from the principle of robust ridge estimation (Zhou et al., 2020; Lian et al., 2021), is defined as:

where *K* is the ridge parameter; and

From Eq. 3, it is evident that the determination of the ridge parameter and the equivalent weight matrix are necessary conditions for obtaining the solution of robust ridge estimation. Thus, the key to robust ridge estimation lies in selecting appropriate ridge parameters and equivalent weight matrices.

### 2.2 Determination of ridge parameters

Comparative study on the determination methods of ridge parameters is one of the current hot topics. Several methods are widely used for the determination of ridge parameters, including the ridge trace method, the generalized cross-validation (GCV) method, the double-h formular method, the L-curve method, among others (Wang et al., 2012). Despite its remarkable applicability, the ridge trace method is too arbitrary and has no strict theoretical basis. The GCV method is theoretically capable of selecting the optimal ridge parameter, but sometimes the change of the GCV function is too smooth or even divergent, so it is difficult to determine its minimum value. Although the double-h formular method is simple in calculation and flexible in application when determining the ridge parameters, its effect is not obvious when the coefficient matrix of the normal equation is severely ill-conditioned. In contrast, the L-curve method, an extremely rigorous method to study ridge parameters in theory, is characterized by accurate determination and good applicability (Aerospace Research, 2018; Cwiakala et al., 2020; Jiang et al., 2020; Lian et al., 2021). Considering its superiority, the L-curve method is selected to determine the ridge parameters in this paper.

According to the regularization theory, the robust ridge estimation criterion of Eq. 1 is as follows:

where *k*;

With *K* will get many various (

Let

The maximum curvature can be obtained by calculating the maximum value of the above equation, and the corresponding maximum point is the desired one. Then, the robust ridge parameter *k* corresponding to the point can be calculated.

### 2.3 Determination of equivalent weight function

The robust effect of robust ridge estimation mainly depends on the equivalent weight function. Different equivalent weight functions lead to different robust estimation models and consequently different robust effects (Lawrence and Marsh, 1984; Wu, 2009; Li S. et al., 2017; Dong-Sheng et al., 2023). The IGG robust scheme proposed by Professor Zhou is to modify the weight matrix of the least squares parameter estimation solution, that is, to replace the prior weight with the equivalent weight, so that we can use the equivalent weight to reconstruct the robust estimation solution and the unit weight mean error (Zhou et al., 2020). In this paper, the IGG robust scheme was adopted. The form of the equivalent weight is defined as:

where

However, the probability that the absolute value of an error is greater than the mean square error, double of the mean square error, and three times of the mean square error is 31.7%, 4.5%, and 0.3% respectively, and most measurement specifications in China stipulate that double of the mean square error is the limit error. Considering the above fact, we take

## 3 Methods

### 3.1 Traditional parameter calculation models

The probability integral method mainly involves five predication parameters, namely, surface subsidence coefficient q, horizontal movement coefficient b, tangent of major influence angle tanβ, major influence propagation angle θ, and deviation of inflection point (including the offset of strike inflection point S_{1}, the offset of dip inflection point S_{2}, the offset of uphill inflection point S_{3}, and the offset of downhill inflection point S_{4}). When conducting robust ridge estimation for parameter determination, q, tanβ, θ, S_{1}, S_{2}, S_{3}, and S_{4} should be obtained firstly by using the measured observed values of surface subsidence. Subsequently, the horizontal movement parameter b can be obtained by iterative fitting of the seven obtained parameters, along with the observed values of horizontal movement of the surface points.

According to the reference (Shuaiying et al., 2021), it can be found that the movement and deformation prediction model of the strike main section is given by:

where

The prediction model for the movement and deformation of the inclined main section is as follows:

where

The essence of the traditional curve fitting method is an estimation method, which takes the least squares principle as the criterion to calculate parameters according to the measured subsidence values and horizontal movement values on the main section. The basic principle is that the target variable

With *n* pairs of observed values (

In other words, a set of parameters B is selected to minimize the sum of squares of deviations between the fitting curve and the measured results.

### 3.2 Parameter calculation model based on robust ridge estimation

The prediction parameters of the probability integral method can be obtained from the measured data of surface movement and deformation at any points along any direction. When the arranged observation stations are conventional, the least squares curve fitting method is used; when they are non-conventional, that is, when they include a series of scattered points, the least squares surface fitting method is employed. However, when the normal matrix of the least squares is ill-conditioned or outliers are present in the observed data, distorted parameter results are likely to occur. The least squares surface fitting method based on robust ridge estimation can not only effectively resist the interference of outliers and ill-conditioned normal matrix, but also apply to the parameter calculation of arbitrarily-shaped or incomplete observation station data (Wang et al., 2012).

According to the probability integral method’s model for surface movement and deformation at any point, the subsidence value *W* at any point can be expressed as a function of the measuring point coordinates and the estimated parameters, namely,

Select initial parameter values “q_{0}, tanβ_{0}, θ_{0}, S_{10}, S_{20}, S_{30}, and S_{40},” and linearize them as follows:

Then, the general form of the error equation is given as:

where

The matrix form of Eq. 10 is:

By determining the ridge parameter *k* from Eq. 5 and the equivalent right from Eq. 6, the robust ridge estimation solution of the corrections of the seven parameters can be obtained, as shown in Eq. 15:

After the results are obtained, the results of the previous calculation are used as the initial values for the next iteration, and the iterative calculation is carried out until the difference of the parameter estimates of the last two calculations meets the iterative convergence accuracy. At this point, the robust ridge estimation solution for parameter estimates can be obtained.

## 4 Result and discussion

Engineering examples in this paper are based on the data of a surface movement observation station in references (Chen, 2010; Li Shuaixin et al., 2017). The observation station is located above the auxiliary 271 working face. Since the observation station records data during both single-seam mining and multi-seam mining, and more detailed data in both strike and dip directions are provided during the former, this paper uses the data of the observation station during single-seam mining. The mining details of the working face corresponding to the observation station are listed in Table 1.

To test the effectiveness of the established robust ridge estimation model, firstly, the least squares method and the robust ridge estimation algorithm were used to fit the parameters respectively. Subsequently, a random error with a standard deviation of 10 mm was generated and added to the observed values using the “random” command in MATLAB software. Afterwards, the least squares method and the robust ridge estimation algorithm were used again to fit the parameters respectively. Finally, a comparative study on the two results of parameter calculation was conducted in both vertical and horizontal directions.

1) The results of fitting the parameters using the original measured data are shown in Table 2.

2) The results of parameter calculation from data fitting after manual intervention are given in Table 3

A comparative analysis on Tables 2–4 reveals the following:

1) When no manual intervention is involved in the observed values, the results of parameter calculation from the robust ridge estimation theory and the traditional least squares estimation theory are essentially the same.

2) When some outliers and random errors are artificially added to the observed values, the parameters obtained using the traditional least squares estimation theory show significant deviations, while the parameters obtained using the robust ridge estimation theory differ insignificantly from those obtained without manual intervention (Table 4).

**Table 4**. Absolute value of the difference between data fitting and parameter calculation results before and after manual intervention.

## 5 Conclusion

1) Traditional observation data are prone to problems such as outliers and ill-conditioned matrices, which tend to occur when using the least squares method to fit the parameters. In other words, when the conventional probability integral method is utilized to obtain the surface deformation parameters, the parameters are prone to considerable errors, resulting in distorted results.

2) Applying the robust ridge estimation theory to the process of obtaining parameters using the probability integral method can automatically eliminate the interference of ill-conditioned normal matrices and outliers.

3) The parameter calculation experiment with manual intervention demonstrates that the established parameter calculation model for the probability integration method based on the robust ridge estimation theory, in comparison with the conventional least squares method, has smaller errors in predicting surface subsidence parameters before and after manual intervention. By adopting the theory for parameter calculation, the predication parameters obtained through the probability integral method boast better effectiveness and reliability.

## Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

## Author contributions

JL: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Writing–original draft.

## Funding

The authors declare that no financial support was received for the research, authorship, and/or publication of this article.

## Conflict of interest

Authors JL was employed by Chn Energy Menxi Coal-Chemical Co Ltd.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

Aerospace Research (2018). *Studies Conducted by S.X. Li et al on Aerospace Research Recently Reported (A three-dimensional robust ridge estimation positioning method for UWB in a complex environment)*. Defense and Aerospace Week.

Ambrožič, T., and Turk, G. (2003). Prediction of subsidence due to underground mining by artificial neural networks. *Comput. Geosciences* 29 (5), 627–637. doi:10.1016/s0098-3004(03)00044-x

Chen, Y. (2010). *Study on surface movement law of deep mining in Kailuan mining area*. Henan, China: Henan Polytechnic University.

China University of Mining and Technology (1981). *Coal mine strata and surface movement*. Beijing, China: Coal Industry Press.

Cwiakala, P., Gruszczynski, W., Stoch, T., Puniach, E., Mrocheń, D., Matwij, W., et al. (2020). UAV applications for determination of land deformations caused by underground mining. *Remote Sens.* 12 (11), 1733. doi:10.3390/rs12111733

Dong-Sheng, H., Cheng, X. K., Ya-Fei, Z., and Li, T. (2023). LIAN Xu-gang Parameter inversion of mining subsidence probability integration prediction method based on space-air-ground integrated monitoring. *COAL Eng.* 55 (1), 81–86.

He, G., Yang, L., and Jia, F. (1991). *etc. Mining subsidence*. Beijing, China: China University of Mining and Technology Press.

Jiang, Y., Misa, R., Tajduś, K., Sroka, A., and Jiang, Y. (2020). A new prediction model of surface subsidence with Cauchy distribution in the coal mine of thick topsoil condition. *Arch. Min. Sci.* 65 (1), 147–158. doi:10.24425/ams.2020.132712

Lawrence, D. K., and Marsh, C. L. (1984). Robust ridge estimation methods for predicting u. s. coal mining fatalities. *Commun. Statistics - Theory Methods* 13 (2), 139–149. doi:10.1080/03610928408828669

Li, S., Li, G., Wang, L., Zhou, Y., Peng, Y., and Fu, J. (2017a). A three-dimensional robust ridge estimation positioning method for UWB in a complex environment. *Adv. Space Res.* 60 (12), 2763–2775. doi:10.1016/j.asr.2017.10.040

Li, S., Li, G., Wang, Li, Zhou, Y., Peng, Y., and Fu, J. (2017b). A three-dimensional robust ridge estimation positioning method for UWB in a complex environment. *Adv. Space Res.* 60 (12), 2763–2775. doi:10.1016/j.asr.2017.10.040

Lian, X. G., Wu, Y. R., Ge, L. L., Du, Z., and Liu, X. (2021). DInSAR monitoring of surface subsidence by fusing sentinel-1A and -1B data to improve time resolution in a mining area. *Can. J. Remote Sens.* 47 (4), 596–606. doi:10.1080/07038992.2021.1952554

Litwiniszy, j (1956). Application of the equation of stochastic processes to mechanic of loose bodies. *Arch. Mech. Stosow.*, T8.

Prediction of subsidence due to coal mining in Raniganj coalfield, West Bengal, Indi. *Int. J. Rock Mech. Min. Sci. Geomechanics Abstr.*,1996,33(2).

Ren, G., Li, G., and Kulessa, M. (2014). Application of a generalised influence function method for subsidence prediction in multi-seam longwall extraction. *Geotechnical Geol. Eng.* 32 (4), 1123–1131. doi:10.1007/s10706-014-9787-y

Shu, D. M., and Bhattacharyya, A. K. (1993). Prediction of sub-surface subsidence movements due to underground coal mining. *Geotechnical Geol. Eng.* 11 (4), 221–234. doi:10.1007/bf00466365

Shuaiying, P., Shengwu, Q., Ming, W., and Li, G. (2021). Mining subsidence prediction for multi-seam and non-rectangular goafs based on probability integral model: a case study from China. *Arabian J. Geosciences* 14 (13), 1231. doi:10.1007/s12517-021-07620-3

State Administration of Work Safety, (2017). *Regulations on mining under buildings, water bodies and railways and coal safety pillars in main roadway*. Beijing, China: China Coal Industry Publishing House. (in Chinese)).

Wang, M., Guo, G., Wang, L., et al. (2012). Prediction parameters calculation of probability integration method based on ridge estimate. *Coal Min. Technol.* 17 (2), 17–19.

Wu, R. (2009). Study of robust estimation parameters model for probability integral method. *Coal Technol.* 28 (9), 167–168.

Zhang, Y., Discussion on the determination method of a and b coefficients of negative exponential function method.*Tongmei Technol.*, 2010,4 : 32-38.

Keywords: robust ridge estimation, probability integral method, parameter determination, illconditioned normal matrix, outliers

Citation: Lou J (2024) A parameter inversion method for the probability integral method based on robust ridge estimation. *Front. Earth Sci.* 11:1330163. doi: 10.3389/feart.2023.1330163

Received: 30 October 2023; Accepted: 27 November 2023;

Published: 07 June 2024.

Edited by:

Yongqiang Zhou, Chinese Academy of Sciences (CAS), ChinaReviewed by:

Xiaoding Xu, China University of Mining and Technology, ChinaShuyang Yu, Hohai University, China

Hongru Li, Tongji University, China

Copyright © 2024 Lou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Jie Lou, loujie2023@163.com