Using the Bentley MicroStation environment to program calculations of predicted ground subsidence caused by underground mining exploitation

The paper presents a new concept of creating a program for calculating ground subsidence caused by underground mining extraction, which is substantially different than previously used solutions. Instead of compiling the calculation algorithm to an executable file, the whole application algorithm has been written in the Bentley MicroStation graphic environment’s interpreted language. This graphic environment is used in mining for keeping mining maps and the format is used for storing finished and designed mining exploitation data sets. The paper describes S. Knothe’s theory of the influence of mining exploitation on the terrain surface in the context of the calculation methodology used for terrain subsidence computation. Based on the theory of influence, the structure of the created MVBA (MicroStation Visual Basic for Applications) algorithm for interpolating subsidence contours is described. The new application is called uDEFO. Calculation results from the program are compared with calculation results from software currently used in coal mines.


Introduction
Underground exploitation of mineral resources causes many negative effects on the terrain surface. It directly causes the partial or full degradation of civil and technical infrastructure, and causes changes in ground water conditions, which lead to changes in the natural environment [1,2]. One of the basic effects of underground exploitation is ground deformation caused by the extraction of mineral deposits that can be seen on the surface as a subsidence basin [3]. There are many theories that describe the influence of extraction and make it possible to calculate the shape of the subsidence basin and distribution of its deformation [3,4,5] on the basis of the shape of the extracted deposit and its parameters. The calculation of deformation process factors has been significantly simplified by the use of computers since as early as the mid-1970s [6]. Since then, several other programs for forecasting post-mining terrain deformation have been created [7,8]. All of these applications were complicated to use and included complex functions which often duplicated the functionality of other computer systems used in coal mines. They required the user to spend substantial time to learn their functionality. However due to their high cost, these software packages have been only installed in small quantities, one or two per mine. This is sufficient for preparing mining plans but limits access to only selected employees of survey and geology departments. In 20 th Century was increasingly common for CAD (Computer Aided Design) software to be used for keeping survey and geological documentation, which is slowly replacing analog records. For this purpose paper records are converted into electronic records using CAD software plugins [9]. Specialist systems for processing digital maps have also been developed [10,11]; however, this type of software causes a problem with the integration of parcel data in CAD systems with geometrical data in ground deformation forecasting software. This is often solved at the lowest level of integration by exporting data and importing it into another system. Despite applying simplifications that facilitate automation of the process to the second level of integration, which means the systems connect to each other and directly exchange data [12] (without an exchange file), there is still the complication of keeping the data in both systems simultaneously. The results of forecasts from computation programs are most often imported into CAD Software to prepare cartographic documentation and are also stored in two different systems.
To solve this problem the authors propose a new solution which involves the use of a third level of data integration. This involves creating an application that processes data within the environment of the platform which manages spatial information. In underground mining in Poland, many mines use Bentley MicroStation software to store so-called hybrid maps of mining excavations (mining-geology maps). This software is used, among others, in the coal mines of the former KHW S.A. or KGHM PM S.A. Therefore, the MicroStation platform was chosen to implement calculations of deformation using the MVBA language. The designed application will eliminate the process of data exchange between systems and limit the possibility of human error.

Forecasting the impact of exploitation on the ground surface
The problem of predicting continuous deformations has been the subject of research since the beginning of the twentieth century. Many theories have been developed that made it possible to calculate depressions and deformations of the terrain. In the past, methods based on empirical formulas and axioms were used; currently, methods based on models of stochastic environments are commonly used, and research is increasingly being conducted on the use of models based on continuous environments and neural networks [13]. In Poland, the currently used theory of forecasting is Stanisław Knothe's theory, which is based on the stochastic model of J. Litwiniszyn [14]. In this model, changes in the surface area caused by exploitation are of random nature and depend on many factors. Not all of these factors can be clearly defined, and the forecast is only an approximation of reality. In his theory, S. Knothe adopted the following assumptions: the deposit is horizontal or almost horizontal (α ≤10÷15 o ), the rock mass, in which the process takes place is homogeneous and incompressible, the description refers to the final process of displacement and deformation, horizontal displacement of points and inclination of the basin profile are proportional [14]. These assumptions were preceded by a series of analyses of subsidence basins that appeared over large fields of exploitation. These observations led to the following conclusions: significant influences extend only a certain distance from the edge of the field; the maximum subsidence occurs directly above the exploitation at a distance larger than the radius of impact; the point of inflection occurs above the edge of the exploitation or above the exploitation; subsidence equal to 50% of the maximum subsidence occurs directly above the edge of exploitation [14].

Parameters of the calculation process and methodology of calculating subsidence
Four parameters are of fundamental importance to the land subsidence calculations in S. Knothe's theory [4]. The first is the depth of exploitation, H, which, depending on the elevation of the deposit in relation to the land elevation, exploitation speed and the type of rock mass above, determines the possibility of discontinuous deformations, simultaneous occurrence of continuous and discontinuous deformations or the occurrence of continuous deformations only. Depth H is crucial for the calculation of the extent of deformation, in which influences may appear. The second parameter is the thickness of the exploited deposit, g (the height of the post-exploitation gap), on which the maximum depth of the settlement basin depends. The third parameter, a, describes how the post-exploitation void is filled, which determines the percentage conversion of the thickness of the selected deposit in relation to the maximum settlement of the area. It can take values from about 0.3 in the case of using backfill up to approx. 0.8 in the case of collapse. The fourth parameter is the extent angle of the main influences, , which depends on the construction of the rock mass. These four factors are used in the calculation of the parameters that describe the course of the subsidence basin. The relationships between these parameters are subsequently defined. The maximum value of ground subsidence is defined by the formula [4]: The r parameter is a parameter of the dispersion of influences and determines the profile of the depression basin, which in theory extends to infinity [4]: It is also called the radius of influence range, which does not mean that there can be no subsidence at a distance farther than r [6]. According to S. Knothe, at a distance greater than r, the model assumes the occurrence of decreases lower than 0.61% Wmax [4]. These parameters and indicators are valid for the so-called flat state of deformation, i.e. for a crosssection of a subsidence basin formed above a large field. A large field means a field whose dimensions are greater than 2r [6]. The shape of the basin in the perpendicular profile of the exploitation basin is shown below in Fig. 1. In the case of exploitation in the form of a polygon, the terrain subsidence at a specific coordinate is given by the formula [6]: The course of terrain subsidence in the basin profile perpendicular to the edge of a large field of exploitation [6].
where: qk (x, y) is the angle included between the axis x of the assumed coordinate system and the straight line passing through the calculation point and the k-th vertex of the exploitation parcel ql + 1 = q1 xk,yk -coordinates of the k-th vertex of the exploitation parcel,

Fig. 2.
Sketch of symbols for a field in the form of any polygon [6].

Determination of surface deformation in the extraction of inclined deposits
The above formulas that relate to determining settlements for horizontal deposits can also be used for calculations in the case of inclined deposits with a relatively small angle of inclination (<10 o ). The subsidence basin above inclined deposits is also characterized by asymmetry, which is bigger the higher the inclination angle is. For slight slopes, this asymmetry of the distribution of subsidence relative to the center of the basin can be disregarded. The basin, however, is shifted towards the direction of the decline (downwards slope). If the angle and azimuth of the dip are known, the calculation activities are as follows.
Determine the value of deviation vector v and vertical thickness of the chosen layer of the deposit g' from the formulas [6]: where: α -slope angle of the deposit ξ -angle of deviation for surfaces, which according to S. Szpetkowski takes the value of 0.7α [6].
Next, the calculation point is shifted by the vector v in the direction of the incline and for the new values x and y the calculations of subsidence are carried out according to formulas (3) and (4), taking into account the new value g when calculating the value Wmax.
For mining cases with greater inclination, the extraction field should be divided into narrower horizontal strips and calculations should be made for each of them. The values obtained should be added up for each calculation point [6].

Programming the algorithm of total subsidence calculation
The theoretical description of the rock mass depression presented above is the basis for its implementation in the form of a computational algorithm. The authors decided to program in Bentley MicroStation, which has the built-in MicroStation Visual Basic for Applications (MVBA) programming language [15]. This is a Microsoft scripting language dedicated to implementation in applications for the Windows operating system. The version for MicroStation has been expanded by a set of objects that represent the geometry of the vector object of the *.dgn file and objects that represent tools that process geometry [16]. This enables reading of data about exploitation parcels directly from the map as well as saving the results of calculations in new layers of the same map. The terrain subsidence calculation application was named uDEFO.

Sources and data types used in uDEFO
The application receives the parcel's geometry directly from the mining map file, whereas the descriptive data must be entered individually for each parcel on a separate drawing layer. The description should be located inside the parcel and should contain seven values entered in seven subsequent lines. In order, these are the name of the deposit, the name of the parcel, the average depth of the parcel, the coefficient of exploitation a, the average thickness of the parcel, the azimuth of the parcel, and the angle of the parcel's decline. Due to the fact that maps of mining excavations are still created in the form of two-dimensional projections of space on a plane, it is necessary to enter text values of parameters of the deposit's geometry such as the angle and azimuth of the decline or the depth of deposition.
The basic complex type of data in the application is the object class "Parc", which contains data about a single parcel. Particularly important is the field "elev_xxx", which is used to record information on whether the depth of the parcel is calculated from the vertices or will be given by the user. In this data type the Wri() table is included, which is the set of coordinates of the plot points x, y, z. The next two types of composite data are the "Side" and "Square" classes, which are used to interpolate isolines of subsidence from the calculation grid of subsidence values. The value of the subsidence is calculated in points and requires transformation into an isoline plot of the spatial distribution of the terrain subsidence with a description of the isolines.

Processing of data in uDEFO
The uDEFO program is run in the MicroStation platform environment, either from the command line (text) or graphically using the VBA Project Manager menu. The first stage of the program's operation is loading of data. In the project file, the user selects the data to be calculated with the MicroStation tools (parcels with their descriptions). The launched application has a "Load Parcels" key, which can be used to load the selected data into the application. The list of loaded parcels is displayed so the user can check the number of parcels and their parameters. If necessary, the user can also edit the parcel text attributes.
The third step is performed by the program without the user's participation. The program calculates the range of influence of the exploitation based on the boundaries and parameters of the loaded parcels. Then, based on the range of impact limits, the program builds a rectangle for the future grid of calculation points that covers the entire deformation range. The user can select the density of the grid of calculation points.
In the fourth stage, the program calculates the plane coordinates of all points in the calculation grid. When loading data, attention must be paid to determining the azimuth and decline of the deposit and its mean depth.

Determining the azimuth, the decline, and mean depth of the parcel
In the case of two-dimensional geometry files or when all vertices had a z = 0 coordinate, it was assumed that the average depth was specified in the description of the parcel, which is read and loaded by the application. The second case is full three-dimensional mining maps of deposits, the vertices of which have individual z coordinates. In this case, the program uses a procedure that calculates the azimuth and the fall angle of the parcel. Due to the fact that the number of vertices is usually greater than three, the plane that is the best approximation of the area determined by all vertices is calculated. The calculations are performed in the slope procedure according to the least squares method. The following matrices A and L are created (Fig 4):

Fig. 4. Matrices A & L.
Using the matrices, the value X is calculated with the following formula: where: Δz -difference between the elevation of points (0,0) i (a1,a2), which is calculated with the formula: (10) v -length of the stretch (0,0) (a1, a2) The calculated values are used to calculate the vector v of the offset of the calculation points according to the formula (5). This value is then divided into components parallel to the axes of the coordinate system and saved to the corresponding fields of the variable type Parc. For the discussed case, in which the vertices are described by three coordinates, the procedure computes the arithmetic mean of the z coordinates and assigns it to the Hparcel field with the average depth of the deposit. For the correct operation of the program, it is necessary that the z values of the parcels' elevations be absolute altitudes (relative to sea level).

Editing data and the computation factor
Loading data, in addition to displaying them in the menu, also activates the function of editing the text parameters of the parcel. Selecting a parcel from the list loads the data from the fields contained in the selected line into the text boxes below. These values can be edited in these fields and they are updated when "Enter" is clicked. In addition, two options are available to facilitate the editing of data. The "Set for All" button allows the same value for the tgB parameter to be set for all parcels at once. In the case of a three-dimensional parcel for which the average elevation value is not given in the text description but results from the calculation of the average height based on the coordinates of the vertices, the "Set" button in the "Terrain Elev" frame is activated. Vertex heights are values above sea level. Entering the land elevation results in the calculation of the depth of the parcel. The last element of the data editing frame is the "Save parcels" button, which completes the editing and prepares the data for the calculation. If the data is correct, the last part of the interface is activated, in which there is a "Grid density" field for setting the distances between the grid calculation points. The second "Contour Step" field determines the vertical distance between the calculated contour lines. Setting these values completes the entering and editing of data. The program interface is presented in Fig. 5.   Fig. 5. uDEFO main application window

Generating the resulting spatial distribution of subsidence
The second phase of data processing begins when the user presses the "Draw Subsidence" button to calculate the value of subsidence at each point based on formula (3). The results of calculations in binary form are saved temporarily on the hard drive. Clicking the "Draw Subsidence" button also calls the "ww_parc" procedure, which is responsible for plotting the contours of terrain subsidence. Drawing the contours takes place in the following stages: -Determination of mesh parameters. For all parcels selected for calculation, the extreme values of x and y coordinates are searched for by means of the P_min and P_max procedures. The resulting range is extended by three times the range of influence; -Calculation of the value of subsidence in grid points. The Grid procedure creates a file of points which have 'x' and 'y' coordinates that define their position in the plane and a 'z' coordinate that represents the amount of subsidence. With a few loops, the coordinates of each subsequent point are calculated based on the bottom left corner, grid density and the number of rows and columns. In the same loops, subsidence is calculated after determining the position. The calculation of values for a single parcel is done by running the procedure "w_pkt". This process is preceded by the calculation of the values of Wmax and r; in the case of an inclined parcel, the vector V offset of the calculation point is also determined. The procedure "w_pkt" contains a loop, in which the integral of the formula (3) is calculated for each edge of the current parcel. The integral is computed numerically by the procedure. The whole process is repeated for all parcels, and the final value is the sum of all calculated subsidence values. The last task of the procedure is to search for the largest and smallest values of depressions, which are used to determine subsequent values of contours; -Creating a file of squares of the interpolation grid. The next step after creating the point file is the "rys_ww" procedure. Once it has been started, the "p_kwd" procedure is called, which creates a file of squares. Due to the complexity of the contour interpolation process, it was necessary to create a file describing the neighborhood between points. For this purpose, an additional "Square" type was created, the structure of which is illustrated in Figure 6. The position in the file and the square number describes the number "n". The fields preceded by the letter "k" contain numbers of "neighbors" in a grid of squares. The fields N, S, E and W -are "Side" variables which contain the coordinates of the vertices and the variable, in which information about the intersection with the contour is recorded. Points inserted as vertices are taken from the points file based on the calculated position from the number of rows and columns. In the case of squares on the edge of the grid, the numbers of adjacent squares on the edge side take the value 0; -Creating contours in a *.dgn file. The task of creating graphic elements representing contours is performed by a layer procedure that is called in the "rys_ww" procedure several times. On the basis of the minimum and maximum values of subsidence, all the values of the contours H are calculated at intervals specified by the user-entered step. In the first stage of operation, the layer procedure retrieves subsequent squares from the file and the condition for intersecting the side of the grid with the contour is checked. This intersection occurs when "h" is between the subsidence values at the ends of the side. This condition is checked by the "check_side" procedure. If no square satisfying the condition is found, the procedure ends. Otherwise, the found square is marked as the initial square (by giving the aforementioned field the value of 2), and the first two points of intersection are entered into the temporary file with contour lines. The coordinates of the intersection points are calculated from the ratio based on the value of z of the vertices and H. This process is performed by the "interp" procedure. The sides through which the contour line passes are marked by giving the value 1 to the "ww" field on the given side. One of the "neighbors" at the sides for which the intersection was found is also remembered. In the second stage, a loop is executed which loads the memorized neighbor and checks subsequent 4 sides for intersections. When checking, the side through which the contour has already passed is omitted, i.e. the side joint with the previous square. When the intersection condition for a given side is met, the coordinate values of the intersection point are calculated and the point itself is saved to the contours point file (Fig 7). The side is marked as being intersected by the contour. A "neighbor" corresponding to this side is remembered. The actions in the loop are repeated until the "neighbor" is the square that is marked as the starting one. The program does not introduce solutions for cases, in which the contour is not closed because it goes beyond the boundaries of the grid. The extents of the grid created within the program exceed the range of occurrence of subsidence, therefore an uncompleted contour should not occur. The next exception is the situation, in which the contour line passes through the vertex of the square. This creates a necessity to search for three "neighbors" and creates an infinitely large number of exceptions in the case of subsequent transitions through the corners. This problem was solved by moving points that could generate exceptions. Before searching for intersection points, all points are checked in all squares. If any point is found whose z value is equal to the level of any contour, the z value is increased by 0.0001 m. This does not cause significant differences in the geometry of the line but eliminates the possibility of error. When all points of the contour are stored in the file in the proper order, creating a graphic object in the file is done by assigning a variable of type "CurveElement", an array of points taken from a temporary file. "CurveElement" is a curve passing through given points, thanks to which the contour is not a broken line. The last task of the procedure is to reset the condition variables in all squares in order to prepare for the creation of the next contour.

Tests verifying the calculations of subsidence
In order to verify the results of the uDEFO program, the EDN program by Jan Białek was used [7]. The data for analysis was compiled on the basis of mining maps from the "Murcki-Staszic" black coal mine. The calculations were performed with the assumption of complete terrain subsidence. In addition to isolines, the parameters of 10 randomly selected points were summarized in the values of parameters. For the calculation, 4 excavation walls with a thickness of 1.8 meters were selected. The exploitation coefficient a was set as 0.8, and angle of draw tgB was assumed to be 2.0 for this exploitation; the exploitation depth was set to 610 m. The results are shown in table 2 and Fig. 8.  0.01 6.08% The average differences were about 6.1% between the forecasts carried out by the EDN program and the uDEFO program, both of which were conducted as part of a larger series of tests with 4 objects. This value shows that the calculations with the developed program are reliable. The difference between the calculated values of subsidence in these programs is close to zero. However, these differences may be a result of different ways of storing values with inadequate precision. In the uDEFO program, the decimal part of variables is not rounded. The exported values from uDEFO have been rounded after importing them to an Excel spreadsheet. EDN software export data files with rounded values of their decimal part.

Summary and conclusions
As a result of the work carried out, the assumed solution was obtained in the form of a program that calculates terrain subsidence directly in the graphic environment of the MicroStation platform. The most important achievement, however, was the development of an original method of contour interpolation. The developed solution maximally uses the possibilities offered by the MicroStation CAD environment, and the authors maximally used the features of the described phenomenon, i.e. the deformation of the terrain in the form of a subsidence basin created as a result of underground exploitation. The ready-made curvilinear geometry was obtained from the CAD program, which made it very easy to draw subsidence isolines. This knowledge about the process of calculating subsidence values has significantly simplified the process of interpolating contours. Isolines of subsidence cannot go beyond the range of exploitation influences, which means that they do not intersect with the boundary of the computational grid. This significantly simplifies the interpolation of the contours of subsidence. For a small number of parcels, the program works surprisingly quickly despite the fact that it is interpreted (not compiled) code. The program can be extremely useful in the geological/surveying sections of mines, where it can be used for calculations related to the forecasting of underground and soil water states [17], height changes of watercourses, sewers and waterworks. In fact, it can be used wherever information about changes in the height of the mining area is necessary.
In the future, the program could be extended with a module supporting inclined parcels exploited in complicated geological conditions [18,19]. However, this would require the use of an algorithm dividing the parcel into smaller zones with new parameters g and H. In addition, it would be necessary to program tasks such as calculating the remaining deformation indicators (deformations, inclinations and horizontal displacements) and to introduce a time function into the calculations that would make it possible to calculate intermediate states of subsidence.
Hereby, we would like to thank Bentley Systems Inc. -Poland for financing the publication of this article. The uDEFO application was created as part of Paweł Owsianka's unpublished master's thesis project [20].