Definition of the form of coal spontaneous combustion source as the inverse problem of geoelectrics

The paper reviews the method of determining the shape and size of the coal self-heating source on coal pit benches and in coal piles during mining of coal by the open method. The method is based on the regularity found in the 1970s of the previous century and related to the distribution of potential of the natural electrical field arising from the temperature in the vicinity of the center of self-heating. The problem is reduced to the solution of inverse ill-posed problem of mathematical physics. The study presents the developed algorithm of its solution and the results of numerical simulation.


Introduction
Spontaneous combustion of coal is one of the main causes of fires in open and underground development of coal seams.Geological and mining factors play a significant role in the process of spontaneous combustion of coal.In case of open mining the tectonic disturbances of coal benches and enclosing strata can be attributed to the first factors when heavily fragmented and cracked coal can lead to increase in the porosity and surface area of coal exposed to oxidation, as a result even at normal temperatures of the ambient air the heat is generated.The order of mining operations in relation to the strike of seams, method of conducting drilling and blasting operations, parameters of the system of mining, the height and time of bench standing can be attributed to the second group of factors.The total influence of all factors shows itself in formation of the loose accumulation of coal, in some cases exceeding the critical mass, which ensures accumulation of heat due to oxidation of coal than its consumption by convection, conduction, thermal radiation and other causes.
For control and diagnostics of the processes of self-heating and spontaneous combustion of coal it is required to promptly assess the power of the heat source, its size and shape.These parameters allow the proper positioning of the wells for injection of antipyrogens in the rock mass and their necessary volume in order to retard coal oxidation process.The theoretical basis of the process of coal self-heating control and definition of the form of the heat emitting source is detected in the 1970s steady relationship between the emergence of pockets of heat and the formation on the surface of coal accumulations of changes in the potential appearing in the process of self-heating of quasi-stationary electric field.There-fore, the problem of determining the shape of the source by measuring the field potential on the surface of the accumulation is very important.
From the point of view of applied mathematics, the problem of determining the shape and dimension of the sources of fields having different nature (electrical, gravitational, thermal, magnetic or other) is the inverse and ill-posed problem of mathematical physics [1,2].Most of the existing theoretical methods and numerical algorithms are devoted mainly to the solution of linear inverse problems, where the unknown function is a co-factor of some well-known expressions characterizing the environment [3][4][5][6].Much less papers are devoted to development of algorithms for solving non-linear inverse problems [7][8][9][10][11].

Formulation and solution of the direct problem
Let us assume that the natural electric field (NEF) is generated by a horizontal cylindrical source of current with arbitrary but constant cross section.We assume as well that the length of the source and depth in the lower layer of two-layer homogeneous isotropic space with plane-parallel boundaries are known (Fig. 1).In this case, the potential of the point source will be determined by the known formula where C -coefficient characterizing the force of the electric current source; M -the point of measurement of the field potential on the earth's surface.
To obtain the computing formula for the value of the cylinder potential, we integrate the expression (1) by the volume of cylinder where Р -the point located in the cylinder volume.
We pass in the integral (2) to non-dimensional variables by the general formula , whereby w we understand coordinates of the point of the field М and the integration variables of cylinder Р.Then, (2) will take the form Then, to simplify the form of the formulas we will not draw the horizontal line at dimensionless values.Let's integrate (3) where ) To calculate the double integral, we define the boundary of the region as a function ) (ϕ ρ in the polar coordinate system.The transition to the latter is performed in the standard way: Then (4) will assume the form . Thus, the problem of determining the form of the source of the field from measurements of the potential on the day surface has been reduced to the Fredholm-Uryson integral equation of the first kind (hereinafter referred to as EI) with respect to the unknown function where -nonlinear kernel; * U -values of potential measured on the day surface (NEF).The formula (5) allows us to compute the potential of NEF of the cylinder in any point on the surface of earth.Further on, to solve the inverse problem of restoration of the cylinder form, we will use only the central axis of symmetry 0 = M y .

Solution of the inverse problem
For the definition of function we consider the problem of finding the minimum of regularized functional of A.N. Tikhonov [1,2,14].
-the right part of the IE, which in practice is set experimentally, and for test problems is determined by solving a direct problem using formulas (4) with the addition of a random correction simulating inaccuracies in full-scale measurements; α is the regularization parameter.
We shall seek the minimum of the functional ( 6) by the method of conjugate gradients [13], whose general iterative scheme has the form where , q -iteration index; k -minimization step; I -linear combination of the derivatives (gradients) in the previous steps; coefficient is determined by one of the formulas [13], / ρ Φ -the Frechet derivative of the functional (6).Let's consider some aspects of numerical realization of formulas (5 -7).The integrand of the kernel R is a smooth monotonically increasing function on the interval , where -fixed consistent values.To calculate the integral in (5), we use the 5-point Gauss formula.
To calculate the integrals and the derivative in ( 6), we define a grid with respect to the variables, M x , φ.The step of the grid with respect to the variable φ will be equal to ) and in case of the non-linear IE, the resulting errors have a significant effect on the computational process.Therefore, we approximate the derivative ϕ ρ′ by a central difference relation of the second order of accuracy . Due to the periodicity of the function ) (ϕ ρ this approximation can also be used for boundary points.In fact, let the function be defined on an interval ] ; 0 [ Λ , where Λ is the period of the function ; then this equality will be true , where . Thus, the derivative on the borders of the grid will be equal to , where the number of the array of values is indicated in square brackets ] [ j ρ ; K-is the number of elements in the grid with respect to the variable φ.Note that if the region has a symmetry about the Ox axis, then this derivative will be zero. Taking all the above into account, we obtain a discrete analogue of the expression (6) of the form -the difference scheme indicated above; µ σ, -are the coefficients of the Simpson formula.
To implement the scheme (7), we calculate the derivative of ( 8) with respect to the variable , and expressions TT[t] are implemented in the following way: ; for the rest four values the following expressions are used: To find the step of minimizing k, we will use the "golden section" method for the objective function ( ) . Numerical experiments have shown that this function is unimodal on the interval 0 ≥ k .Since it is impossible in advance to determine the interval at which the minimum point is contained, we adhere to the following algorithm.

Algorithm
Step 1.We set the starting interval for the change of k: and the error of the "golden section" method ε << ε 1 .
Step 2. Using the "golden section" method, we define the minimum point of the objective function -k min .We calculate the auxiliary value and go to step 2. If "no", then the value of k min is the required minimum of the objective function.
The regularization parameter α will be determined, following the studies [13][14][15], by the following iterative way: ( ) , while in the first step of the conjugate gradient method we set 0 ) 0 ( = α , and on the second step - . In the works [9,10], a saw-tooth effect of "blurring" an approximate solution is observed, depending on the choice of the regularization parameter α.In case of a linear IE, as these studies show, it is possible to select such a constant value of the parameter const q = α ) ( that the indicated 'saw-tooth' formation will not be noticeable.In case of a nonlinear IE, such a value of the α parameter could not be chosen, so to neutralize the "blurring" after each step of the conjugate gradient method, the values were smoothed using the arithmetic mean of two neighboring values ] The 1 st International Innovative Mining Symposium Considering the above solution method, it can be noted that it depends on the choice of the starting values and the patterns of variation of the following parameters ) (q γ , ) (q α , ε, ε 1 , ρ (0) .Numerical experiments have shown that the choice of concrete realization of each of them has a significant effect on the efficiency of solving the inverse problem.
Here is an example of the implementation of a numerical algorithm.
Let us consider a test region in the form of a cylinder 100.0 m long and a cross-section in the form of an ellipse, the dimensionless equation of which in the polar coordinate system has the form where we assume z M = 150.0;a = 130.0; b = 70.0 m.To simulate errors in full-scale measurements, we will add a random correction to the potential values calculated from formula (2) as follows: ( ) , where rand are uniformly distributed numbers on the interval [0; 1], U ex -the exact calculated values of the potential according to the formula (4).
The choice of the selected parameters will be the following:  ( ) We will show in Figures 2 and 3 the graphic result of restoration of the cylinder crosssection form after 200 iterations.The numerical method of interpreting geophysical data of measurements of the potential of a quasi-stationary electric field is developed in the article in order to determine the form of the heat release source.Let's give practical recommendations on the use of the developed methodology.
A coordinate grid is constructed on the surface with a step of about 25-50 meters.In each node of this grid, the field potential is measured.The contour map is constructed, along which the elongation line of the proposed source is determined.More accurate measurements in the amount of at least 100 are taken along this line.These data will be used as initial data.If detailed measurements are impossible, then an approximate interpolation polynomial is constructed, along which the required values of the potential are determined.Further, with the help of the developed program, the form of the section of the cylindrical heat source is defined.

Fig. 1 .
Fig. 1.The scheme of the containing space and cylindrical field source

1 , 2 .
the step by the variable M x -All integrals in (6) will be calculated by the Simpson formula.In the work[16]  in case of the linear IE, the derivative ϕ ρ′ is approx- imated by the right difference formula.Numerical experiments have shown that with a relatively small (about 100) number of nodes ] [ j ϕ , the accuracy of such an approximation is low (

Step 3 .
We verify the truth of the inequality 1 , 0 < ∆ .If "yes", then increase the right boundary of the minimum search interval 1 + = b b

Fig. 2 .Fig. 3 .
Fig. 2. Restoration of the area and value of the potential after 200 iterations and smoothing only after the first iteration of the conjugate gradient method )