## 2 Materials and Methods
### Dynamic Deformation Prediction and GIS Overlapping Analysis Models  
The two main functions were developed related to MIDA, all of which are similar to the actual requirements of the mining engineer/enterprise. One is the PIM, and the other is the GIS spatial overlapping analysis model. The former uses the dynamic prediction method (DPM) of probability integral, which considers the correlation between the process of surface subsidence and the spatial position and time. The DPM is the final state of surface subsidence caused by underground coal mining, namely PIM, multiplied by the time effect function. The PIM is based on the random medium theory and has been extensively used in China. The PIM can be used to illustrate the prediction model of an arbitrary point on the surface, and the formula of unit mining surface subsidence basin is determined according to Equation (5):  
\[W_{e}(x)=\frac{1}{r}e^{\pi\frac{x^{2}}{r^{2}}} \tag{5}\]  
Considering the 3D characteristics in Figure 11, the surface subsidence of any point A(\(x,y\)) due to the unit mining B(\(s,t\)) can be expressed as Equation (6):  
\[W_{e}(x,y)=\frac{1}{r^{2}}e^{\pi\frac{(x-s)^{2}+(y-t)^{2}}{r^{2}}} \tag{6}\]  
Figure 10: Procedure of collecting and storing the entity data. (**a**) The interactive editing interface and (**b**) the established GIS database of house entity and query.  
When the mining is O\({}_{1}\)CDE, the subsidence prediction formula of point A caused by the whole mining is as follows:  
\[W(x,y)=W_{0}\int_{0}^{l}\int_{0}^{L}\frac{1}{r^{2}}e^{\pi\frac{(x-x)^{2}+(y-t)^{ 2}}{r^{2}}}\,\mathrm{d}t\mathrm{d}s \tag{7}\]  
where \(W_{0}=mq\cos\alpha\) is the maximum subsidence at supercritical mining, \(r=H_{0}/\tan\beta\) is the main influence radius, \(l=D_{3}\)-\(S_{3}\)-\(S_{4}\) is the strike length, and \(L=(D_{1}\)-\(S_{1}\)-\(S_{2})\cdot\sin(\theta_{0}+\alpha)/\sin\theta_{0}\) is the dip length.  
To use this model, it is essential to determine the prediction parameters \(q\), \(\tan\beta\), \(S_{1}\), \(S_{2}\), \(S_{3}\), \(S_{4}\), and \(\theta_{0}\). The sample data for our case study are listed in Table 1. Based on the subsidence formula, aside from the surface subsidence (referred to as w from now on), other deformation values of any point A(\(x,y\)) comprising slanting, curvature, horizontal movement, and displacement along the strike direction (90\({}^{\circ}\)) and dip direction (0\({}^{\circ}\)) for the surface (\(X\), \(Y\)) can be deduced, hereinafter referred to as (i90, i0), (k90, k0), (u90, u0) and (s90, s0), respectively. After that, the deformation isolines for each surface deformation can be obtained using a certain interpolation calculation method based on its corresponding scattered data points. Interested readers can refer to [1] for more detailed PIM information. When we accomplish the dynamic deformation prediction calculation represented by deformation isolines, a subsequent GIS overlapping spatial analysis based on clipping operation can be conducted, which requires the intersection between the damage areas depicted by closed deformation contours and the GIS entities, including point, line, or polygon geometry. Based on the intersected results, the new GIS entity data layers can be produced hereafter. It should be noted that the two main functions used Visual Studio 2017 C# and were encapsulated in the server-client as a dynamic link library (DLL) format, which was called when the Cesium web-based framework triggered the MIDA engine.