Determination of parameters for the prognosis of rock mass deformation with the use of the cellular automaton method

The article describes the basics of construction and operation of a simple computational model built on the basis of the cellular automata theory designed to simulate rock mass subsidence caused by underground mining exploitation. The construction of the computational model as a deterministic finite cellular automaton has never been described in the literature in the context of the choice of method parameters. In the context of the construction of the model, an exemplary methodology for the selection of model parameters was presented on the example of the results of the geodetic observation of the consequences of underground mining exploitation.


Introduction
As a result of underground mining exploitation in the rock mass, a post-mined void is created. Huge pressure and the mass of rock layers causes the emptiness to be naturally filled. The process of filling the post-exploitation void leads to the creation of rock mass deformation. When certain conditions are met, continuous deformations occur in the form of a subsidence trough on the surface of the mining area. This type of deformation in contrast to discontinuous deformations is predictable with much greater accuracy.
The problem of description of phenomena occurring in the rock mass as a result of mining operations was and still is the subject of research by numerous authors [1,2]. In the period preceding the development of computer techniques, the study of rock mass movement was usually performed on models of equivalent or loose materials.
In model studies, two main principles of similarity are adopted. The principle of geometric similarity is based on rescaling the dimensions of the real fragment of rock mass. In turn, the principle of mechanical similarity is based on [3]:  geometrical compatibility between model displacements and displacements determined from geodetic measurements,  compatibility between forces affecting individual model parts and their actual counterparts  if the time factor is also taken into account, maintaining a constant ratio of model deformation velocity to the speed of actual deformations. Numerous studies have led to the creation of a group of methods for forecasting the deformation of the mining area, which include [4]:  methods based on empirical formulas, where deformation indices were determined on the basis of various graphs, patterns and nomograms [5], IOP Conf. Series: Earth and Environmental Science 261 (2019) 012045 IOP Publishing doi: 10.1088/1755-1315/261/1/012045 2  methods which uses geometric-integral theory to describe the distribution of deformation [4,6]  methods basing on models of continuous centre where the state of stresses and displacements defines a system of differential equations and the equation of state dependent on the adopted model [3,7],  methods based on models of stochastic centers used for the first time by J. Litwiniszyn [8],  methods that uses the artificial neural networks [9]. The development of computer techniques and attempts to describe nonlinear phenomena, not recognized by the widely used Budryk-Knothe theory, has revived in recent years the emergence of new possibilities in the field of forecasting the deformation of rock mass through the use of the theory of cellular automata [10]. The usefulness of the method has already been demonstrated on both theoretical and real examples [11,12]. The method adopts the principle of geometrical similarity, however, its characteristics correspond to the characteristics of stochastic centers models [8].
The article presents the basics of construction and characteristics of the computational model as a spatial deterministic finite automaton. With reference to the presented information, the method of proceeding in choosing the parameters of the method was discussed. Sample methodology of selecting the parameters of the method was presented on the example of the analysis of the effects of real mining exploitation.

The computational model as a deterministic finite cellular automaton
A cellular automaton is a numerical model of a dynamic system. Dynamic system means one in which component elements move one another. In contrast to the Monte Carlo method, which in general is a mathematical method of simulating a stochastic center, in a cellular automaton, the simulated evaluation process of a given phenomenon may be of a random or deterministic nature. The history of both methods dates back to the 1940s. The mathematician of the Hungarian origin of J. Von Neuman is widely regarded as the author of the method. The real development of the method took place about 20 years later thanks to E. Codd. Due to his work, mobile automata have become practically applicable [13]. Countless observations of the effects of various cellular automata have made their basic classification. The fundamental one was made in 1983 by Stephen Wolfram [14]. It shows that there are cellular automata distinguishing by certain features.
Cellular automata are widely used in many fields of science. In the area of mining area protection, the possibility of practical application was shown for the first time in the world due to the works of P. Sikora. [10]. However the originator of the possibility of using this method in this area was T. Niemiec [15].
Each automaton defines several basic parameters. These include, first of all: cellular grid, cellular neighborhood, boundary conditions and transition function. The elements listed below are briefly discussed with reference to a computational model that maps a fragment of a rock mass. The rock mass model as a cellular automaton arises through its discretization in the form of a regular grid of cells that are of the same size and shape and adhere closely to each other (figure 1).  Figure 1. Spatial grid of the cellular automaton with determined number of rows (W) and columns (D x S). For individual cells can be ascribed real dimensions: length (Dk), width (Sk) and height (Wk) [16].
The cell grid can be two or three-dimensional. In any case, it has a limited size. Cells are assigned real-size dimensions: length Dk, height Wk and, in the case of spatial variant, length Sk. In this way, it is possible to reproduce by scaling the mining exploitation in a specific fragment of the rock mass (figure 2).
Cells store numeric data. Cells representing the selected seam are assigned a value equivalent to the volume of the resulting post-mining void. The way of filling the void is taken into account. Ultimately, the void size assigned to the cell is equal to the product of the thickness of the seam g and the exploitation rate a. In the sense of the method, this is understood as the initiating subsidence. The mapping of the deck or multiple decks in the described way determines the initial conditions for the operation of the simulation.
The operation of the automaton consists in the simultaneous evaluation of all the cells of the machine. Each cell in the machine changes its state depending on the condition of the cells in its socalled cellular neighborhood. The initiator of the changes will be a different states of cells mapping the post-exploitation void. It is assumed that the process of rocks collapsing to the resulting post-mining void is simulated. In the simplest terms, it is assumed that only the force of gravity is responsible for this phenomenon. In the sense of the presented method, it allowed to define the so-called cell neighborhood. According to the definition of a cellular automaton, it specifies cells that are in the close vicinity of a given cell considered in a discrete time of evaluation of the model within the scope of which data exchange will take place. The following figure illustrates the arrangement of cells being in the vicinity of the base cell located in the row below (figure 3). These cells have been marked in green. The cellular neighborhood determines the direction of the model's evaluation. In this case, it is directed towards the row reflecting the surface of the mining area. Thus, the end of the simulation will take place when all the cells will be evaluated and the sum of subsidences on the surface of the model will be equal to the sum of the subsidences assigned at the beginning of the simulation. This means that the deformation distribution is lossless.
The final characterization of the model depends on the assumed transition function -the algorithm being the heart of the method. For the assumed cellular neighborhood, the basic characteristic of the subsidence propagation from the base cell is assumed to be in accordance with the figure below (figure 4).  Figure 4. Characterization of the transition algorithm from the base cell to cells from the determined cellular neighbourhood [16].
The above picture (figure 4) presents the initial deformation distribution algorithm from the base cell for the simulation of horizontal displacements. The basic parameter of the algorithm is the socalled main transition P. It determines the ratio of rock masses moving to the cells from the neighborhood in the vertical direction to the volume moving to the remaining cells. It is a parameter that is particularly important in the simulation of horizontal displacements -not discussed in a minor article. In other words, the parameter Pϵ(0,1) specifies the part of the subsidence assigned to the base cell that will be passed to the cell lying directly above it.
In this way, the distribution simulation is performed for each cell of the automaton. The same cellular neighborhood and the same transition algorithm are assumed (defined before starting the simulation). The deformation distribution is registered in the entire grid of the automaton. Simulation of subsidence using the described method leads to results consistent with the characteristics of the loose center model of J. Litwiniszyn [8].
Numerous observations of the effects of the method allowed to develop the initial characteristics allowing its practical application (1) [16].
where: TX;TYmaximal inclination in case of a full subsidence trough along X and Y axis. agmaximal subsidence in case of a full subsidence trough, Hdepth of mining exploitation, A -adjustment parameter due to the value of the main transition parameter P. As in the case of the geometrical-integral theory of Knothe's theory, a number of geological properties of the rock mass defining its ability to create subsidence trough with a given slope are  [6]. In the cellular automata method the parameter of the maximum slope of aT is defined analogously to the tgβ parameter of Knothe's theory (equation 2).
In contrast to the tgβ parameter, the aT parameter is variable depending on the depth of mining exploitation H. This confirms the results of Z. Rogusz's research [17] on the variability of the tgβ parameter depending on the depth of mining operations.
The maximum slope parameter aT can be determined from the criterion of matching the theoretical subsidence trough to the trough resulting from geodetic measurements.
Next, the methodology of selecting the basic parameters of the method to describe the deformation of the rock mass is presented.

Selection of basic parameters of the method for the simulation of rock mass vertical displacements
In order to perform simulations of rock mass subsidences, for circular mining and geological conditions, the basic parameters of the method should be determined. Assuming the previously described cellular neighborhood and the transition function, it is necessary to define at least:  cell dimensions mapped in reality Wk, Sk and Dk,  the value of the maximum slope parameter aT,  the size of the main passage P.
Choosing the size of cells is a complex task. The adopted cell size will affect:  the size of the mapped rock mass for the number of columns, rows and levels adopted,  the model resolution determining the accuracy of the post-mined void mapping and obtained simulation results,  the efficiency of numerical calculations. When determining the size of a cell grid, must remember to ensure an unfettered deformation distribution. In other words, the cell grid size should be larger than the expected extent of deformation on the surface of the model.
The smaller the size of the cells we accept, the higher the resolution of the obtained results will be. However, this leads to increased demand for the amount of computer's RAM (random access memory). It is considered reasonable to accept the size of the cells adapted to the distance between the points of the measuring line and to the dimensions of the building objects to be protected.
Next, an example of how to proceed in the selection of method parameters for a part of a rock mass in the area of done mining exploitation is presented. The example applies only to the simulation of rock mass subsidence as an indicator characterized by the smallest random dispersion [18,19].
In the years 1974-75, the "Dębieńsko" mine exploited the longwall 4a in the 326/5 deck with the longitudinal system with the roof collapse to the average height of 1.4 m. The deck was located at an average depth of about 155 m. The average inclination of the layers is about 10 °. The ratio of the dimensions of the longwall and the depth of mining operations allows to accept the conditions of a flat state of deformation.
In the area of the southern edge of the longwall, a measuring line was established on the surface of the mining area. The average distance between the measuring points was approx. 15 m. The contours of the longwall against the background of the measuring line are shown in figure 5.
Maximal subsidence measured on the line was equal 1050 mm. A full subsidence trough was created on the land surface.
To determine the basic parameters of the method, i.e. the maximum slope ratio aT and the ceiling control coefficient a, width Sk and cell length Dk equal to 15 m were assumed.  The graph below (figure 6) presents the values of the R 2 determination coefficient determined by comparing the profile measured with calculated for different sizes of aT parameter. The results obtained indicate a good fit in each case. However, for the parameter size aT equal to 1.8 and 2.1, the R 2 coefficient of determination was the largest.
The graph below (figure 7) presents the differences between the standard deviation obtained from calculations (for the points of the measurement line) and the size of the standard deviation determined from the results of geodetic measurements.  Big differences between the profiles of subsidence troughs result from the failure to take into account the so-called operating edge. Then it was taken into account geometrically by appropriate reduction of the operating field. The parameters were determined again -in the sense of the least squares method -finally assuming the parameter values: aT = 2.1, a = 0.75 and operating edge   [m]. The contours of a suitably reduced panels' size of the operating edge are shown in figure 5 with a dashed line. Figure 9. Profiles of the subsidence trough slope determined from measurements and numerical calculations for the best values of parameters a, aT and operational edge obr in the sense of the least squares method.
Taking into account the operating edge significantly improved the quality of adjustment of the subsidence troughs ( figure 9). In the end, the size of the parameters should be considered as final in this case.
It should be remembered that when using the criterion of parameter selection, it is always necessary to simulate the distribution of vertical displacements. A single simulation is numerically effective compared to other methods, eg finite element methods, however, multiple simulations are a significant overload for the computer. In the case of using CA3D software (copyrights by P. Sikora) appropriate optimization techniques are used and the processor is multithreaded.
However, some differences can still be noticed. They result from the failure to take into account the deviation of influences observed on the surface due to mining exploitation of inclined seams. Unfortunately, this is a complex issue and not yet even theoretically solved for the case of the spatial model with characteristics different from the flat model [16].

Summary
The article describes the basics of construction and operation of a simple computational model built on the basis of the cellular automaton designed to simulate rock mass depressions caused by underground mining exploitation. The construction of the computational model as a deterministic finite cellular automaton has already been described in the literature but never in the context of the choice of method parameters.
The computational model described here presents a fragment of a rock mass through its discretization in the form of a spatial grid of cells of equal dimensions and shapes that closely adhere to each other. By assigning dimensions to the cell mapped in reality, the method adopts the principle of geometric similarity.
For the practical application of the method, mathematical relationships, developed by P. Sikora [10,16], binding the most important parameters of the method are used.
It is worth noting that the characteristics of the model coincide with the characteristics of the loose medium of J. Litwiniszyn [8]. On the other hand, the principle of geometric similarity makes it possible to determine parameters from the results of geodetic measurements of underground rock mass exploitation. In analogy to the method of S. Knothe 11 mass to create subsidence trough with a certain slope are included in one parameter -aT. However, unlike the tgβ parameter, the value of the maximum slope parameter aT varies with the depth of mining exploitation (equation 2).
The maximum slope parameter aT as well as the exploitation parameter a can be determined from the criterion of matching the theoretical subsidence trough to the trough resulting from geodetic measurements. The methodology of such conduct is presented in the article on the example of the real mining exploitation of the longwall 4a in the 326/5 seam. Taking into account the possibility of the exploitation edges, included in the model in geometric form, made it possible to obtain a match with a coefficient of determination R 2 =0.993 and an mean error of 34 mm. The parameters determined in this way could be the basis for forecasting the impact of another mining operation in this part of the rock mass.
Despite the simplicity of building the computational model as well as the transparency of parameters, there are still things to be solved. An example is the development of practical characteristics of the model for the simulation of vertical and horizontal displacement in the case of inclined seams. Including this factor would certainly increase the accuracy of the forecast.
Finally, it should be recognized that the method has significant potential for further development. An additional possibility of including nonlinear properties of summing up the influences and discontinuities in the structure of the rock mass to the deformation distribution (in a direct way) makes the method an excellent complement to S. Knothe's geometrical-integral method commonly used in Poland.