## 2 Research Methodology  
### Principle of the Probability Integral Method  
The probability integral method is a mining subsidence mobile deformation prediction method based on stochastic medium theory, which considers that the surface subsidence belongs to random events [20]. It proposes to divide the whole mining area into many tiny units, and the surface subsidence is the sum of the effects of each unit on the surface superimposed on each other and can be solved directly by the integration method [21; 22; 23].  
Take the origin o of the unit coordinate system o-xyz as the center of the mining unit, and the unit subsidence basin-wide expression on the neighboring domain dGA (area microproduct) of point A(x,y,z) is  
\[W_{e}=f\Big{(}x^{2}\Big{)}f\Big{(}y^{2}\Big{)}=\frac{1}{r^{2}}e^{-\pi\frac{x^{2 }+y^{2}}{\pi^{2}}} \tag{1}\]  
In the formula presented, \(f\big{(}x^{2}\Big{)}\) and \(f\big{(}x^{2}\big{)}\) denote the probability density function that follows a normal distribution; \(r\) denotes main influence radius; and the calculation formula is  
\[r=\frac{H}{\tan\beta} \tag{2}\]  
In the formula presented, \(H\) denotes the mining depth of the working face along the strike; \(\tan\beta\) denotes the tangent of major influence angle.  
\(W_{e}\) divide the entire mining area to obtain the expression for the subsidence full-basin integral:  
\[W(x,y)=W_{0}\int\int_{s}W_{e}dS \tag{3}\]  
In the formula presented, \(dS\) is the area of coal seam micro-original mining unit, and \(W_{0}\) is the maximum subsidence value.  
Therefore, the formula model for calculating subsidence at any point on the surface is  
\[W(x,y)=W_{0}\times C_{xm}\times C_{ym} \tag{4}\]  
In the formula presented, \(C_{xm}\) and \(C_{ym}\) denote \(x\) (strike) and \(y\) (inclination) directions, respectively; \(W_{0}\) denotes the maximum subsidence value, which is calculated by the following formula:  
\[W_{0}=m\times q\times\cos\alpha \tag{5}\]  
In the formula presented, \(m\) denotes mining thickness in the vertical tangential direction; \(q\) denotes subsidence coefficient; and \(\alpha\) is the dip angle of coal seam.  
\[C_{xm}=\frac{1}{\sqrt{\pi}}\times\left(\int_{-\sqrt{\pi}\times\frac{y}{\pi}}^{ \infty}e^{-\lambda^{2}}d\lambda-\int_{-\sqrt{\pi}\times\frac{y-1}{\pi}}^{ \infty}e^{-\lambda^{2}}d\lambda\right) \tag{6}\]  
\[C_{ym}=\frac{1}{\sqrt{\pi}}\times\left(\int_{-\sqrt{\pi}\times\frac{y}{\pi}}^ {\infty}e^{-\lambda^{2}}d\lambda-\int_{-\sqrt{\pi}\times\frac{y-1}{\pi}}^{ \infty}e^{-\lambda^{2}}d\lambda\right) \tag{7}\]  
In the formula presented, \(d\lambda\) denotes the unit thickness of the coal seam; \(l\) is the calculated length of the working face along the strike; \(L\) is the calculated length of the tendency working face; and their calculation formulas are  
\[l=D_{3}-s_{3}-s_{4} \tag{8}\]  
\[L=(D_{1}-s_{1}-s_{2})\times\frac{\sin(\theta_{0}+\alpha)}{\sin\theta_{0}} \tag{9}\]  
In the formula presented, \(D_{1}\) is the working face tendency length; \(D_{3}\) is the working face strike length; and \(\theta_{0}\) is the propagation angle of extraction; and \(s_{1}\), \(s_{2}\), \(s_{3}\), and \(s_{4}\) are deviation of the inflection points of the downhill boundary, uphill boundary, left boundary, and right boundary, respectively.  
Similarly, the horizontal movement of any point on the surface is given by  
\[U(x,y)_{\varphi}=U_{ax}\times C_{ym}\times\cos\varphi+U_{ay}\times C_{xm} \times\sin\varphi \tag{10}\]In the formula presented, \(\varphi\) is the angle value from the positive counterclockwise of the \(x\)-axis to the specified direction; and \(U_{ax}\) and \(U_{ay}\) are the horizontal movement values of the projective point positions on the principal section of the strike and tendency, respectively.