Distributed Fuzzy Clustering Analysis of Time-Lapse Electrical Resistivity Tomography for Water Inrush Monitoring in Coal Mines

Abstract

The majority of water inrush accidents in coal mines are caused by mining engineering activities. To avoid water inrush accidents, the Time-lapse Electrical Resistivity Tomography (TLERT) is applied to monitor water migration in fractured zone. A great challenge for the application of TLERT monitoring is the huge and numerous time series data sets generated by monitoring systems, which are difficult to process manually. This research proposed a distributed fuzzy clustering algorithm based on kernel function estimation to analyze TLERT images automatically. The resistivity date can be classified with different cluster centroids. The fuzzy c-means algorithm was chosen to display resistivity change. The algorithm was validated using a floor water inrush model. The results indicate that the water migration in the fractured zone can be monitored automatically and the edge of the resistivity changing area can be shown clearly.

1. Introduction

Among the major accidents in coal mines in China, the number of occurrences of mine water inrush accidents and associated deaths is second only to that of gas accidents, and the economic losses caused by it have always been the first in various coal mine disasters. Most water inrush accidents are caused by mining engineering activities. The traditional geophysical exploration in coal mines uses advance detection, but the results can only reflect the distribution of primary water flowing structures such as hidden water sources and faults in the roof, floor or side of the working face. During the mining process, because of the damage to and movement of the mining strata [1,2,3], the water sources hidden in the roof, floor and sidewall of the working face will migrate. The cracks and breaks in the roof, floor and the coal seam will form real-time changing water flowing channels, which will bring serious difficulties to the prediction and prevention of water hazards.

Time-lapse Electrical Resistivity Tomography (TLERT) is a geophysical dynamic monitoring technology that performs 2D/3D data acquisition at a certain time interval with the installation of several electrodes, which can be used to monitor the electrical changes of the underground medium in real time. Authors such as Liu et al. (2009), Zhang et al. (2009) and Gao and Meng (2011) have used Electrical Resistivity Tomography (ERT) to dynamically detect the overburden deformation and fracture [4,5,6]. Liu et al. (2009) applied ERT monitoring to the mine water inrush model test [7]. Nevertheless, the large amount of monitoring data is difficult to process manually. Most of the existing inversion algorithms are based on the idea of iterative fitting, which is time-consuming and easily influenced by artificial factors, limiting the application of TLERT in mine water inrush monitoring.

Aiming at the inversion algorithm of Time-lapse Electrical Resistivity Tomography, Daily and Ramirez (1992) proposed the ratio inversion, which uses the data ratio of adjacent times to highlight the abnormal part of the image [8]. However, time-lapse electrical resistivity images are likely to cause false anomalies due to inversion errors, and these false anomalies stand a good chance of magnifying image differences. Loke (2001) proposed to use the initial data inversion results as a prior model to invert the data of each monitoring period, which improved the calculation speed [9]. Subsequently, Oldenborger et al. (2007) and Miller et al. (2008) conducted an in-depth study of this method [10,11]. LaBrecque and Yang (2001) used a difference algorithm to invert the difference of resistivity at different times in the same place [12]. The idea is to calculate the background model by the conventional inversion method, and afterwards to invert multiple sets of monitoring data, respectively, using the differential inversion method. The model is used as a prior model for all monitoring data and is used to calculate differential images, and has achieved good results. In recent years, new inversion methods such as joint delay inversion, synchronous delay resistivity inversion, and 4D “active time” constrained resistivity inversion have emerged [13,14,15,16]. In China, Liu et al. (2013) and Liu et al. (2012) have applied inversion methods such as the smooth constraint, the inequality constraint and the space structure constraint to the data processing of mine floor water inrush monitoring and grouting effect monitoring, which have achieved good results [17,18].

Machine learning is such a powerful tool that we have seen tremendous success in areas such as Computer Vision, Speech Recognition, and Natural Language Processing [19,20,21,22,23]. The same is true in the field of anomaly detection. Wu et al. (2018) used automatic identification and classification of Palomar Transient Factory astrophysical objects in GLADE [24]. Xiao et al. (2018) applied the SMK-means which is achieved by Mini Batch K-means based on simulated annealing algorithm for anomalous detection of massive household electricity data, which can give the number of clusters and reduce the number of iterations and improve the accuracy of clustering [25]. Fuzzy clustering is an un-supervised machine learning algorithm that can realize image analysis automatically and its application effect in ERT as well. Ward et al. (2014) classified the resistivity data of three-dimensional bedrock inter-face by using the fuzzy recognition algorithm [26]. Chambers et al. (2015) used the fuzzy clustering method to analyze images of dynamic monitoring of groundwater level changes in quarries, and delineated the water level change interface clearly [27]. In this research, a distributed fuzzy clustering algorithm based on kernel function estimation is applied in TLERT image processing to monitor underground water migration in coal mines.

2. Methodology

2.1. Time-Lapse Electrical Resistivity Tomography

Time-lapse Electrical Resistivity Tomography (TLERT) can continuously observe the potential or apparent resistivity response caused by the change of electrical properties. The electrical distribution of dielectrics can be continuously imaged by the inversion calculation. As is shown in Figure 1, two electrodes (A, B) are used to supply power to the medium, and another pair of electrodes (M, N) are used to measure the potential difference. The high resistivity body has the ability to repel current. The higher the resistivity, the stronger the repelling ability. In contrast, a low resistance body has the ability to attract current; the lower the resistivity, the stronger the attraction. Each current line leads to the B electrode according to the path with the least resistance encountered, that is, it obeys the so-called minimum energy principle. Therefore, the distribution of an underground resistivity anomaly can be calculated by measuring the potential difference.

2.2. Fast Inversion Algorithm of Time-Lapse Electrical Resistivity Tomography

Time-lapse resistivity inversion imaging aims to reconstruct the underground resistivity distribution with inversion of the data at each monitoring time. The traditional resistivity inversion method is time-consuming and cannot meet the needs of the real-time processing and interpretation of mine water inrush monitoring data. Based on forward modeling, we adopted the smooth-constrained least squares algorithm developed by Akca to inverse the TLERT monitoring data [28]. Using the inversion results of the previous time as a prior model for the inversion of monitoring data at the next time, the calculation speed has been significantly improved.

The specific calculation steps can be shown as follows:

(1)

Establishing an initial model. If the time is t₀, then the initial model is set as a uniform model, and the resistivity value is the arithmetic mean of the measured apparent resistivity data; otherwise, the final inversion result of the apparent resistivity at the previous moment will be used as the initial model for inversion calculation.

(2)

Calculating the apparent resistivity theoretical value and the objective function value of a given model, and then calculating its Jacobian matrix. According to the inversion principle, the Jacobian matrix needs to be recalculated every time the model parameters are corrected, which is time consuming. Therefore, the Jacobian matrix is recalculated only when the model correction error is greater than 5%.

(3)

Singular decomposition. This method is used to solve the normal equations to obtain the correct value of model parameters, and then the modified model is calculated.

(4)

Performing a forward calculation on the corrected model, and comparing the apparent resistivity value with the measured apparent resistivity value. When the following two conditions are not met at the same time, step (2) is continued: (i) when the root mean square error is less than 3%; (ii) the maximum number of iterations reaches 30 times.

In the inverse calculation of the Time-lapse Electrical Resistivity Tomography, the matrix to be solved is often nonlinear and ill-conditioned, and the final solution is not unique. The program uses a smooth constrained least squares inversion algorithm, and the correct value of model parameters is [28,29]:

$$\mathrm{\Delta }{m}_i={({\mathit{J}}^T{W}_d^T{W}_d\mathit{J}+\lambda L)}^{-1}({\mathit{J}}^T{W}_d^T{W}_d\mathrm{\Delta }d-\lambda L{m}_{i-1})$$

(1)

where $\mathrm{\Delta }{m}_i$ is the correct value of model parameters, $\mathit{J}$ is the Jacobian matrix, also known as the sensitivity matrix; $W_d$ is the weighting matrix of the original data; $\lambda$ is the regularization factor (the initial value of regularization in this program is the logarithmic standard deviation of the measured data, and the minimum value should not be less than 1% of the initial value per iteration); $L$ is a five-point finite difference Laplace operator, which is calculated by the function “delsq” in MATLAB. $\mathrm{\Delta }d$ has the following relationship with $W_d$:

$$\mathrm{\Delta }d=\ln {\rho }_s^{true}-\ln {\rho }_s^f$$

(2)

$$W_d=diag(\frac{1}{\sqrt[4]{{\rho }_s^{true}}})$$

(3)

where ${\rho }_s^{true}$ is the measured apparent resistivity value and ${\rho }_s^f$ is the apparent resistivity value of the forward calculation.

According to the correct values of the model parameters, the updating formula can be obtained.

$$m_{i+1}=m_i+\mathrm{\Delta }{m}_i=m_i+{(J^T{W}_d^T{W}_{d}J+\lambda L)}^{-1}(J^T{W}_d^T{W}_d\mathrm{\Delta }d-\lambda L{m}_{i-1})$$

(4)

The data fit is measured by:

$$R_{MS}=\sqrt{{(W_d\mathrm{\Delta }d)}^T(W_d\mathrm{\Delta }d)/N}$$

(5)

2.3. Automated Image Analysis

The fast inversion of monitoring data can be realized with the smooth-constrained least squares algorithm, whereas the abnormal boundary cannot be recognized explicitly and automatically as before. As an unsupervised classification method, distributed fuzzy clustering has been applied to resistivity image analysis successfully. In this paper, the distributed fuzzy clustering algorithm is applied to analyze the inversion results and resistivity difference data. The position and number of centroids are determined according to the peak value of the probability density curve. The anomaly and the positive and negative change areas of resistivity difference are delineated automatically with the fuzzy clustering C-means.

2.3.1. Fuzzy C-Means Clustering

Fuzzy clustering analysis is a clustering method based on the similarity between data objects. According to the characteristics of resistivity, this paper uses the fuzzy C-means clustering algorithm for monitoring data analysis. The core idea of the algorithm is to divide the data into several clusters. The data differs greatly between clusters, and the data similarity between the same clusters is higher. Finally, the classification correlation is represented by the fuzzy membership value, the non-similarity measure is defined by the Euclidean distance, and the objective function is constructed to iteratively relocate the data. The data is moved repeatedly from one cluster to another for classification and partition to search for optimal membership matrix and clustering center. When the objective function reaches the minimum, the iteration terminates.

The set $\rho$ of N resistivity data is: $\rho =\{{\rho }_1,{\rho }_2,\cdots ,{\rho }_N\}$ and ${\rho }_j$ represents the number j data in the data set. The matrix $U=[u_{ij}]$ is obtained by solving the membership matrix of sample data set $\rho$. U is a matrix of $C\times U$. C denotes that the sample data set $\rho$ is classified into c categories. $U_{ij}$ denotes the degree of the membership of the data ${\rho }_j$ in the data set $\rho$ under the i-th classification. The degree of the membership $u_{ij}$ should satisfy the following conditions:

(1)

The sum of all membership degrees of the given arbitrary data ${\rho }_j$ is 1, as follows.

$${\displaystyle \mathbf{\sum }_{i=1}^c{u}_{ij}=1,\forall j=1,2,\cdots ,N}$$

(6)

(2)

Each cluster c contains at least one point with a non-zero membership, but cannot include all points with membership degree 1.

$$0<u_{ij}<1,1\le i\le c;1\le j\le N$$

(7)

The fuzzy matrix $U^{(0)}$ is randomly generated based on the above criteria and is used for iterative initialization. The formula for solving the membership degree u is:

$$u_{ij}={\left[{\displaystyle \mathbf{\sum }_{j=1}^c{\left(\frac{d_{ik}}{d_{jk}}\right)}^{\frac{2}{m-1}}}\right]}^{-1},\forall i\in \{1,\cdots ,c\};\forall k\in \{1,\cdots ,N\}$$

(8)

where $m\in [1,\infty )$ is a fuzzy index or a weighted index; $d_{ik}=‖x_k-v_i‖$ represents the Euclidean distance between the data ${\rho }_k$ and the cluster centroid $v_i$. Cluster centroid collection is represented by the formula $V=\{v_1,v_2,\cdots ,v_c\}(v_i\in V,i=1,2,\cdots ,c)$, where $v_i$ represents the i-th cluster centroid. The $v_i$ calculation formula is as follows:

$$v_i=\frac{{\displaystyle {\sum }_{k=1}^N{u}_{ik}^m\cdot x_k}}{{\displaystyle {\sum }_{k=1}^N{u}_{ik}^m}}$$

(9)

The process of obtaining the cluster centroid and membership matrix by fuzzy C-means clustering algorithm is an iterative process. The specific iterative process is as follows:

(1)

Given the cluster type number c, initialize each cluster center according to the iterative convergence conditions (6) and (7).

(2)

Repeat the following steps until the data membership is stable.

(a)

The membership function is calculated according to Formula (8) using the current cluster centroid.

(b)

The centroid is recalculated according to Formula (9) using the current membership function.

When the iteration converges, various types of centroids and membership values are obtained, and the division of fuzzy clustering is completed.

Since the initial value of the Fuzzy C-Means (FCM) algorithm needs to be randomly set, this randomness cannot guarantee that the global optimization can be achieved every time. That is to say, only the initial cluster centroid setting has global characteristics, and the clustering result can reach the global optimal. Due to the large amount of monitoring data and the time-consuming iterative process, a distributed fuzzy clustering method is proposed. The method firstly uses the kernel density estimation to obtain the global cluster centroid value, and then uses the fuzzy clustering C-means to cluster, thus overcoming the defect that the fuzzy cluster C-means is sensitive to the initial value and can easily fall into the local optimal solution. The calculation result can achieve global optimization, and it skips the iterative optimization process, which can reduce the calculation time.

2.3.2. Kernel Density Estimation

The kernel density estimation is used to create a histogram to represent the distribution of data. According to the difference of probability density, a symmetric weighting function summation (kernel function) is applied to each data, instead of simply assigning each point to each interval. In the kernel density estimation, the set of N resistivity data is given as $\rho =\{{\rho }_1,{\rho }_2,\cdots ,{\rho }_N\}$. The general form of the quasi density function can be estimated as:

$${\widehat{f}}_{\sigma }(\rho )=\frac{1}{N}{\displaystyle \sum _{i=1}^N{K}_{\sigma }(\rho -{\rho }_i)}$$

(10)

where $K_{\sigma }(\rho -{\rho }_i)$ is the kernel function, $\sigma$ is the bandwidth, $\rho$ is the independent variable to be estimated, $\widehat{f}$ is the probability density function, N is the number of resistivity values, and ${\rho }_i$ is the i-th resistivity value. The selection of kernel function needs to satisfy the conditions of attenuation, normalization and symmetry. Gauss function, trigonometric function, dual-bandwidth function and Epanechnikov function are commonly used, but the selection of kernel function has little influence on the estimation effect, while the influence of bandwidth is obvious [30]. Small bandwidth will produce many small peak distributions, while large bandwidth will produce fewer peaks and smoother curves, as shown in Figure 1. In theory, the optimal choice of bandwidth is to make the kernel density estimation function closest to the kernel density function of the actual data (Figure 2). The bandwidth can be selected by establishing an optimization mode. However, since the kernel density function of the actual data is unknown, it is very difficult to establish the optimization model directly. Therefore, based on the normal reference criterion, the normal distribution is used as the reference distribution to select the optimal bandwidth, which is called the embedding method [31].

$$\widehat{\sigma }={\left(\frac{4{\sigma }^5}{3N}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$5$}\right.}$$

(11)

where $\sigma$ represents the standard deviation of all resistivity values.

In this paper, Gauss function was chosen to be the kernel function:

$$K_{\widehat{\sigma }}(\rho -{\rho }_i)=\frac{1}{\widehat{\sigma }\sqrt{2\pi }}{e}^{-\frac{{(\rho -{\rho }_i)}^2}{2{\widehat{\sigma }}^2}}$$

(12)

The final calculated probability density function is:

$${\widehat{f}}_{\widehat{\sigma }}(\rho )=\frac{1}{N}{\displaystyle \sum _{i=1}^N\frac{1}{\widehat{\sigma }\sqrt{2\pi }}{e}^{-\frac{{(\rho -{\rho }_i)}^2}{2{\widehat{\sigma }}^2}}}$$

(13)

2.3.3. Fuzzy Clustering C-Means Based on Kernel Density Estimation

For a set of random variables (the data set studied in this paper is the resistivity value or the resistivity difference), the cluster centroids are pre-estimated by the initial approximate probability density distribution of the data set. The kernel density estimation method is used to find and determine the cluster centroids, and the probability density function ${\widehat{f}}_{\widehat{\sigma }}(\rho )$ is used to analyze the data set. Each peak value in the probability density function represents the probability estimation of the data in total data. In the calculation of ${\widehat{f}}_{\widehat{\sigma }}(\rho )$, the resistivity corresponding to the peak value of the probability density curve is approximately equal to the mean value of the sample data. Therefore, the number of peaks represents the appropriate number of clusters to packet data, which is the number of categories.

Assuming that the bandwidth in the kernel function K is suitable for the data set, the number of peaks represents the number of clusters, and the ${\rho }_i$ of each peak that corresponding to the cluster is its cluster centroid.

The specific steps of program design are as follows:

(1)

Importing the resistivity data that needs to be inverted, then performing preprocessing, and then mapping the 2D resistivity matrix into a 1D vector;

(2)

Calculating the variance of resistivity, and then the standard deviation of resistivity is obtained;

(3)

Calculating the optimal bandwidth using Equation (11);

(4)

Using the Formula (13) to calculate the probability density function and discretizing the probability density. Determining the peak value and the number, so as to determine the cluster centroid and the number of categories of resistivity data;

(5)

Finding the membership degree corresponding to each item of resistivity data, and then forming an N × K membership matrix, where N represents the number of resistivity data, and K represents the number of categories;

(6)

Finding the maximum value of each column of the membership matrix, which is the maximum membership degree corresponding to each resistivity value;

(7)

Assigning the category attribute of the resistivity value according to the maximum membership degree of each resistivity value;

(8)

Outputting class matrix and inversing 1D matrix to 2D matrix;

(9)

Outputting the results and drawing the clustered contour map.

It is worth noting that in the actual calculation process, there are often many small peaks (the peak value is too small, indicating that the clustering range is small, there is no research significance), which results in more classification categories. To solve this problem, the centroids with small probability density are filtered out by setting the threshold value. The disturbance caused by the fluctuation of a very small range of data is reduced, and the threshold value is set as 10% of the mean value of the probability density of all data in the same data sample. The centroid whose probability density is less than the threshold is abandoned.

3. Cases Study

3.1. Identification of Bedrock Interface

The smoothness-constrained inversion technique is applied in the fast inversion algorithm of the monitoring data, which results in the indistinct appearance of anomalous boundaries. Nevertheless, the resistivity change is subtle over time usually and delineating anomaly boundary accurately is important to monitor water migration. To address this problem, the clustering was applied to analyze the ERT resistivity data and it was proved effective in the case of identification of the bedrock interface.

The test site is located in the campus of the China University of Mining and Technology, Xuzhou, China. The site is a landform unit of piedmont alluvial plain. A 28-electrode dipole-dipole array with dipole separations of 4 m was used for detecting the bedrock interface (Figure 3). In order to validate the availability of the method, an exhaustive geotechnical engineering investigation (Figure 4) was made before this test, which indicates that the depth of bedrock interface ranges between 11 m and 13 m.

When the data collection has been finished, the smooth-constrained least squares algorithm is firstly used to inverse resistivity distribution (Figure 5). The x coordinate represents the horizontal distance, and the h coordinate represents the depth. Based on the resistivity difference, the bedrock interface can be distinguished at the depth of 10~14 m. We can see that this bedrock interface identification process is inaccurate. The results of data processing vary with the processors and there are unavoidable errors compared with the real depth of the bedrock interface.

The fuzzy clustering algorithm was applied to analyze inversion resistivity. The resistivity values are classified into two categories with kernel density estimation. The resistivity classification diagram (Figure 6) can be drawn using a distributed fuzzy clustering algorithm. The x coordinate represents the horizontal distance, and the h coordinate represents the depth. It can be seen that the depth of bedrock interface is between 11 m and 13 m, which is consistent with the actual situation. The fluctuation of bedrock interface can also be seen from the map, which is almost consistent with the bedrock depth and its extension fluctuation in the engineering geological profile drawn by borehole data. Therefore, the application of distributed fuzzy clustering algorithm has a significant effect on bedrock interface delimitation.

3.2. Water Inrush Dynamic Monitoring Model of the Coal Seam Floor

The confined water will migrate through the water flowing fractured zone in the course of mining. The resistivity of the fractured zone is obviously different from that of surrounding rocks before and after the water filling, which makes water migration monitoring possible using Time-lapse Electrical Resistivity Tomography. The relevant geo-electric model of the coal seam floor confined water migration can thus be established. Before working face mining, the model can be considered as a level layered media model contained surrounding rocks, aquiclude and aquifer when time is t₀. The time series of t₁, t₂, t₃ and t₄ represent the changing process of water flowing fractured zone which will break through the aquiclude at time t₄. Most coal seams are high-resistance, according to the research; the underground whole-space model can be reduced to a surface half-space model [32] (Figure 7).

The software “ELRIS2D” written by Akca is used for the forward and inverse calculation of the water inrush dynamic monitoring model [28]. The results show that the inverse location and range of the aquiclude and the water flowing fractured zone are in accordance with those in the real model, while the inverse range of the aquifer is larger. The water migration in the fractured zone can be observed approximately (Figure 8).

In order to make the edge of the water flowing fractured zone clear, we analyze the inverse resistivity images with distributed fuzzy clustering. For this study, the clustering is utilized to the log resistivity ratio rather than inversion resistivity to highlight subtle change over time. The kernel density of resistivity change logarithm is estimated first. After the optimal bandwidth is obtained with the embedding method, we can find three centroids in the probability density curve. Then the logarithm of resistivity change can be divided into three classes. The blue zone represents the resistivity decrease, which matches the water flowing fractured zone. The yellow and white zone shows less and no change of resistivity, respectively. Compared with inverse results, the distributed fuzzy clustering graph of log resistivity ratio plot can reveal resistivity change more clearly. The water migration underground can be monitored automatically with the blue zone changing over time (Figure 9).

4. Conclusions

A distributed fuzzy clustering algorithm based on kernel function estimation is applied to analysis of TLERT images. The distributed fuzzy clustering method is used to analyze the inversion results and resistivity difference data in view of the problem that the boundary of the resistivity image is fuzzy and the anomaly delineation is easily affected by human factors. According to the peak value of the probability density curve, we can determine the position of the center of mass and the number of classes, and delineate the range of abnormal bodies. This method is an unsupervised machine learning algorithm. The difficulty of analyzing a huge amount of monitoring data automatically can be solved. Meanwhile, the edge of an abnormal resistivity area can be shown more clearly compared with a conventional apparent resistivity contour map. The distributed fuzzy cluster analysis is carried out for the resistivity difference at different times of the floor confined water uplift. Compared with the results without cluster analysis, the location of the uplift zone delineated by this method is more consistent with the actual situation. The prospect of distributed fuzzy clustering analysis of Time-lapse Electrical Resistivity Tomography for water inrush monitoring in coal mines is an important contribution.