A three-dimensional electromagnetic perspective detection method for the internal structure of underground coal mine working faces

Abstract

The present invention discloses a three-dimensional electromagnetic perspective exploration method for the internal structure of an underground coal mine working face. The method arranges boreholes in a broken line shape, with the upper and lower rows serving as transmitting and receiving boreholes, respectively. Electromagnetic wave perspective detection is performed using a one-shot, multiple-receiver method to obtain perspective data. During three-dimensional observation data processing, the detected space is divided into grid cells. Electromagnetic wave rays pass through several grids from the transmitting point to the receiving point, and the length and absorption attenuation coefficient of the rays in each grid are calculated. Based on the absorption attenuation coefficient, a K-means clustering algorithm is used to perform anomaly segmentation, enabling three-dimensional exploration of anomalies within the working face. The present invention achieves the positioning of anomalies in the height direction of the coal seam, solving the problem of spatial three-dimensional positioning of anomalies. It can present the three-dimensional spatial distribution of uneven anomalies within the working face, helping coal mining to move further towards unmanned and intelligent operations.

Description

Technical Field

[0001] The invention belongs to the technical field of geophysical exploration and relates to a three-dimensional electromagnetic perspective exploration method for the internal structure of an underground working face in a coal mine. Background Art

[0002] Current coal mine face internal structural exploration techniques, due to limited construction conditions, typically employ a two-dimensional (2D) approach. This involves working within two existing tunnels within the working face. The resulting results primarily interpret the planar position of target geological bodies within the working face, but there is no robust method for assessing the spatial distribution and scale of anomalous geological bodies. With the advent of the information age, the development of transparent and intelligent mines places greater demands on understanding the precise distribution of coal seams, their overlying and underlying rock formations, and the presence of anomalous geological bodies within these rock formations. Accurately understanding all geological attributes within and around the coal seam working face in three dimensions is a pressing need for transparent mines. Due to tunnel space constraints, conventional exploration is typically conducted in two dimensions. This mismatch between observation methods and requirements severely restricts the acquisition of 3D geological attribute parameters. Researching true 3D spatial data acquisition technologies and corresponding data processing methods is an essential path to addressing the issue of coal mine face transparency. Summary of the Invention

[0003] In view of the deficiencies in the prior art, the purpose of the present invention is to provide a three-dimensional electromagnetic perspective exploration method for the internal structure of an underground coal mine working face, so as to solve the problem of three-dimensional exploration of geological anomalies inside the coal mine working face.

[0004] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0005] A three-dimensional electromagnetic perspective method for detecting the internal structure of an underground coal mine working face is provided. The method is based on gas extraction drilling holes constructed alternately at high and low levels within the coal seam working face. The method uses a single-transmitter, multiple-receiver synchronous electromagnetic perspective method to achieve three-dimensional detection of abnormal bodies within the working face. The method comprises the following steps:

[0006] Step 1, drilling and probe layout and detection: Arrange two rows of drill holes in a broken line, determine the transmitting and receiving drill holes, and use a single-transmitter-multiple-receiver configuration to deploy the transmitting and receiving probes for detection data collection, recording the spatial positions of the transmitting and receiving points and the corresponding measurement values;

[0007] Step 2, 3D observation data processing: Divide the detection area into 3D discrete grids, convert the exponential attenuation term in the approximate formula for electromagnetic wave propagation in space into a discrete form, and perform a logarithmic equivalent transformation to solve for the absorption attenuation coefficient;

[0008] The approximate formula for electromagnetic wave propagation in space is:

[0009]

[0010] In the above formula, H is the signal value at the receiving point, H0 is the signal value at the transmitting probe, r is the distance between the transmitting point and the receiving point, and β is the absorption attenuation coefficient of the medium;

[0011] The exponential decay term is converted into discrete form:

[0012]

[0013] In the above formula, H is the signal value at the receiving point, H0 is the signal value at the transmitting probe, r is the distance between the transmitting point and the receiving point, N is the number of grids that split the ray, r i is the length of the ray in the i-th grid among the N grids it passes through, β i is the absorption attenuation coefficient of the i-th grid among the N grids that the ray passes through;

[0014] Transpose the terms in formula (2) and take logarithmic equivalent processing to transform it into the following formula (3). Combine multiple rays and solve to obtain the absorption attenuation coefficient;

[0015]

[0016] Step 3: Based on the absorption attenuation coefficient, the K-means clustering algorithm is used to perform anomaly segmentation to achieve three-dimensional detection of anomalies inside the working face.

[0017] The present invention also includes the following technical features:

[0018] Specifically, the step 1 includes:

[0019] Step 1.1: In a lane on one side of the working face or in lanes on both sides, drill two rows of holes in an alternating pattern. The holes should be oriented perpendicular to the direction of the coal mine lanes, with the trajectory kept horizontal. The spacing between holes at the same height should be 10-15m, and the drilling depth should be 100m-120m.

[0020] Step 1.2: Use the upper row of drill holes as transmitting drill holes and the lower row as receiving drill holes. Each transmitting drill hole forms a perspective combination with the 2-4 nearest receiving drill holes.

[0021] In step 1.3, one electromagnetic wave transmitting probe is placed in the transmitting borehole, and 9-11 receiving probes are placed in the receiving borehole. The transmitting and receiving probes move synchronously, exiting from the bottom of the hole toward the hole mouth point by point for detection; the perspective combination of the transmitting borehole and the receiving borehole is replaced until all boreholes are traversed; the spatial positions of the transmitting and receiving points and the corresponding measurement values ​​are recorded.

[0022] Specifically, in step 2, when the ray passes through N grids, the coordinate difference between the ray emission point coordinate (xf, yf, zf) and the receiving point coordinate (xr, yr, zr) in the three dimensions of X, Y, and Z directions is:

[0023] Δx=xr-xf

[0024] Δy=yr-yf (4)

[0025] Δz=zr-zf

[0026] Then the ray equation is established as follows:

[0027]

[0028] In formula (5), the variable t takes a value between 0 and 1.

[0029] Specifically, in step 2, the sequence of intersections of the ray with each grid in the X direction is (xf, xfr1, xfr2, ..., xfrn, xr). According to the first equation of formula (5), the equivalent transformation is performed to obtain:

[0030]

[0031] Let the right side x of formula (6) be equal to each element in the X-direction intersection sequence (xf, xfr1, xfr2, ..., xfrn, xr), and the coefficient sequence is: {0, t xfr1 ,t xfr2 ,...t xfrn ,1}, the coefficient sequence is: in the X direction, from the emission point to the receiving point, the ratio of the distance from the intersection of the ray and the grid to the emission point to the distance between the emission point and the receiving point.

[0032] Specifically, in step 2, the sequence of intersections of the ray in the Y direction with each grid is (yf, yfr1, yfr2, ..., yfrn, yr). According to the second equation of formula (5), the equivalent form is transformed to obtain:

[0033]

[0034] Let the right side y of formula (7) be equal to each element in the Y-direction intersection sequence (yf, yfr1, yfr2, ..., yfrn, yr), and the coefficient sequence is: {0, t yfr1 ,t yfr2 ,...t yfrn ,1}, the coefficient sequence is: in the Y direction, from the emission point to the receiving point, the ratio of the distance from the intersection of the ray and the grid to the emission point to the distance between the emission point and the receiving point.

[0035] Specifically, in step 2, the sequence of intersections of the ray with each grid in the Z direction is (zf, zfr1, zfr2, ..., zfrn, zr). According to the third equation of formula (5), the equivalent form is transformed to obtain:

[0036]

[0037] Let the right side z of formula (8) be equal to each element in the Z-direction intersection sequence (zf, zfr1, zfr2, ..., zfrn, zr), and the coefficient sequence is: {0, t zfr1 ,t zfr2 ,...t zfrn ,1}, the coefficient sequence is: in the Z direction, from the emission point to the receiving point, the ratio of the distance from the intersection of the ray and the grid to the emission point to the distance between the emission point and the receiving point;

[0038] Combining the intersection coefficient sequences in the above three directions, we can get the ratio of the distance from the intersection point to the transmitting point to the distance between the transmitting point and the receiving point, that is, the intersection characterization value sequence {0, t1, t2, ... t N ,1}.

[0039] Specifically, the step 3 includes:

[0040] Abnormal threshold division: Count the maximum and minimum values ​​of all absorption attenuation coefficients in the detection area, divide the maximum and minimum values ​​into 10 intervals, calculate the number of absorption attenuation coefficients falling in each interval, form an absorption attenuation frequency distribution histogram, and select the quantile value less than 10%-20% or greater than 80%-90% as the detection abnormal area; for collapse columns, faults, and water-bearing areas, select the 80% quantile value as the threshold for distinguishing abnormalities from non-abnormalities; for caves, goaf tunnels, and goaf areas, select 20% as the threshold for distinguishing abnormalities from non-abnormalities;

[0041] Classification of abnormal spatial locations: After the abnormal threshold is determined, the K-means clustering algorithm is used to spatially classify the abnormalities. All grids exceeding the threshold are classified. A series of adjacent grid cells belong to the same abnormality, and disconnected grid cells are other independent abnormalities. Clustering is performed according to the cell distance to obtain K clusters;

[0042] Anomaly attribute analysis and classification: After obtaining the elements in K classes, each class corresponds to an independent anomaly. The characteristics of the elements in the class are analyzed to obtain different attribute characteristics of the anomaly. The attribute characteristics include but are not limited to: the scale of the anomaly, the boundary shape of the anomaly, and the direction characteristics of the anomaly. The specific properties of the anomaly are inferred through the attribute characteristics.

[0043] Specifically, the K-means clustering algorithm is used to perform spatial classification of anomalies, including:

[0044] 1) The spatial distance between the i-th grid and the j-th grid in the detection area is defined as follows:

[0045]

[0046] Among them, the coordinates of the center point of the i-th grid are (x i ,y i ,z i ), the coordinates of the jth grid center point are (x j ,y j ,z j );

[0047] The distances between all abnormal grids constitute a data volume. The spatial distances between all abnormal grids are calculated to form a distance matrix D, where:

[0048]

[0049] The distance matrix D is a symmetric square matrix, d 12 is the distance between the first abnormal grid element and the second abnormal grid element. The closer the value is to 0, the closer the two are in space.

[0050] 2) Set a threshold ε to determine whether elements belong to the same class. If the difference between the elements is less than ε, they belong to the same class. If it exceeds ε, they belong to two classes.

[0051] 3) Determine the minimum value in the first column of the distance matrix D except D(1,1). If the minimum value is not greater than ε, then the minimum value and the first element in the absorption attenuation data volume constitute a class cluster1; and record the average value and element number of all elements in cluster1. At this time, there are 2 elements in cluster1, and their average is meancluster1; if the minimum value is greater than ε, then the first element and the minimum element constitute two classes, cluster1 and cluster2 respectively. Class cluster1 has the first element, and class cluster2 has the first element. Calculate the mean value of the elements in the two classes respectively;

[0052] 4) Determine whether the remaining elements belong to the existing clusters cluster1 and cluster2. If they do, place the element in that cluster. If not, construct a new cluster cluster3, update the numbers of the elements in each class, the class mean, and the numbers of the remaining elements; until all elements are classified and a total of K clusters are formed;

[0053] 5) Perform three-dimensional anomaly display on the K clusters obtained by clustering algorithm analysis.

[0054] Specifically, step 3 further includes an analysis of the scale of the anomaly: each cluster formed by clustering is regarded as an independent anomaly, and the sum of the grid volumes corresponding to all elements in the cluster is calculated to obtain the volume of the anomaly.

[0055] Specifically, step 3 also includes an analysis of the trend characteristics of the anomaly: extracting the spatial coordinates of the grid corresponding to the class element and counting the center positions of all grids; using the principal component analysis method in turn to obtain the dominant direction, sub-dominant direction, insignificant direction and corresponding length of the grid unit distribution, and analyzing to obtain the trend characteristics of the anomaly, that is, the dominant development direction and scale of the anomaly; based on the trend characteristics of the anomaly, judging the type and nature of the anomaly.

[0056] Compared with the prior art, the present invention has the following technical effects:

[0057] The present invention utilizes the height difference between adjacent drill holes to realize the positioning of the abnormal body in the height direction of the coal seam, solves the problem of spatial three-dimensional positioning of the abnormal body, and can realize the three-dimensional spatial distribution of uneven abnormal bodies inside the working face, helping coal mining to move further towards unmanned and intelligent. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 Schematic diagram of drilling arrangement;

[0059] Figure 2 This is a schematic diagram of the distribution of transmitting boreholes and receiving boreholes;

[0060] Figure 3 This is a perspective diagram of electromagnetic wave transmission and multiple reception;

[0061] Figure 4 Schematic diagram of three-dimensional space mesh division;

[0062] Figure 5 Schematic diagram of ray penetration grid;

[0063] Figure 6 is the frequency distribution diagram of attenuation coefficient;

[0064] Figure 7 This is the abnormal distribution map of the cluster analysis of the spatial distribution of the three-dimensional absorption attenuation coefficient. DETAILED DESCRIPTION

[0065] As one of the most common methods for solving the internal structure of the working face, electromagnetic perspective technology has the characteristics of non-destructive and fast operation, which is very suitable for the working needs of transparent working faces. After the working face tunnel is formed, production will not generally begin immediately. Instead, boreholes will be arranged inside the working face to extract the gas stored inside the working face to reduce the risk of later mining. The layout of the extraction boreholes is generally parallel or cross-arranged. When cross-arranged, the boreholes are located at different heights of the tunnel. The height difference between adjacent boreholes can be used to locate the anomaly in the height direction of the coal seam. The perspective itself solves the problem of planar positioning of geological anomalies, and overall solves the problem of spatial three-dimensional positioning of anomalies. Transparent working faces inevitably require understanding the three-dimensional occurrence information of strata and geological bodies rather than their planar distribution. Using electromagnetic wave perspective with high and low offset extraction boreholes can realize the three-dimensional spatial distribution of uneven anomalies inside the working face, helping coal mining to move further towards unmanned and intelligent operations.

[0066] Conventional electromagnetic wave fluoroscopy methods use two or more boreholes on the same plane to conduct inter-hole fluoroscopy to obtain the projection of the actual geological anomaly in the plane. The present invention uses a method of arranging boreholes in a 'broken line' to perform three-dimensional electromagnetic wave fluoroscopy, such as Figure 1 and 2 As shown, after the transmitting borehole position is determined, the other row of 2 or 4 adjacent boreholes are used as receiving boreholes. Figure 3 The electromagnetic wave perspective is carried out according to the rules shown in the figure to obtain perspective data. Figure 4 The scheme is divided into small hexahedral grid units. Assuming that the absorption attenuation of electromagnetic waves by each small grid unit is uniform, the electromagnetic wave rays are as follows: Figure 5 The law shown in the figure propagates in a straight line within the grid, passing through several grids from the transmitting point to the receiving point. The length of the ray in each grid is calculated, and the measurement anomaly is allocated according to the length of each ray and the length of each grid unit in the grid system. The attenuation coefficient frequency distribution diagram of all grid units is formed by weighted superposition of multiple rays, as shown in the figure below. Figure 6 , determine the classification of anomaly intensity based on the frequency distribution diagram. Use three-dimensional mapping to show the distribution of anomalies in space and calculate the anomaly scale, direction, etc.

[0067] Specific embodiments of the present invention are given below. It should be noted that the present invention is not limited to the following specific embodiments, and all equivalent modifications made on the basis of the technical solution of this application fall within the protection scope of the present invention.

[0068] Example:

[0069] This embodiment provides a method for three-dimensional electromagnetic perspective exploration of the internal structure of an underground coal mine working face. The method is based on gas extraction drilling holes constructed alternately at high and low levels within the coal seam working face. The method uses a single-transmitter, multiple-receiver synchronous electromagnetic perspective method to achieve three-dimensional exploration of abnormal bodies within the working face. The method includes the following steps:

[0070] Step 1, drilling and probe layout and detection: The zigzag arrangement of drilling is the basis for achieving three-dimensional electromagnetic perspective. Figure 1 Horizontal boreholes are arranged alternately on the top and bottom plates of the coal seam. During construction, the top plate is used as the launch borehole, and the four boreholes closest to the bottom plate are used as receiving boreholes. The relative relationship between the launch borehole and the receiving borehole is as follows: Figure 2 As shown, the data is collected using a one-to-many working method in the hole. The construction method is shown in Figure 3 The hole contains one electromagnetic wave transmitting probe and 9-11 receiving probes, each connected by a coaxial cable. The cable core serves as a common signal transmission channel for the receiving probes. During construction, the transmitting and receiving probes move synchronously, exiting point by point from the bottom of the hole toward the hole mouth to complete the measurement. Figure 1 In the figure, the zigzag drill holes are divided into two rows, upper and lower, which are staggered to ensure that the electromagnetic wave perspective method has a certain resolution in both the horizontal and vertical directions. This is consistent with the distribution of gas extraction drill holes in the actual coal seam, realizing multiple uses of one hole. Figure 2 It is a schematic diagram of the positions of the transmitting drill holes and the receiving drill holes. In the upper and lower rows of drill holes on the working face, the transmitting drill holes and the receiving drill holes are separated into the upper and lower rows. There are 2 drill holes on each side of the transmitting drill hole, a total of 4 drill holes serving as receiving drill holes, forming effective coverage of the coal seam space. Figure 3 It is a relationship diagram between the transmitting point and the receiving point in the transmitting borehole and the receiving borehole. It adopts the form of one transmit and multiple receive. Each transmitting point corresponds to 9-11 receiving points, and the transmitting position is continuously moved to achieve full coverage between the two holes.

[0071] In this embodiment, inside a lane on one side of the working face or in lanes on both sides, according to the following Figure 1 The drilling distribution shown is to construct two rows of boreholes. The borehole orientation is perpendicular to the direction of the coal mine tunnel. The trajectory is kept horizontal as much as possible. The spacing between boreholes at the same height is 10-15m, and the upper and lower rows of boreholes are staggered. The drilling depth can be 100m-120m. Optionally, the drilling can use existing gas extraction boreholes or be constructed separately. Before the three-dimensional exploration construction, the boreholes are cleaned to avoid the blockage of the boreholes by collapsed holes, rock chips, coal blocks, etc., which may cause construction interruptions. According to the borehole distribution, the transmitting borehole and the receiving borehole are designed. The relative positions of the transmitting and receiving boreholes are as follows: Figure 2Each transmitting borehole forms a perspective combination with the 2-4 nearest boreholes in another layer, and together they control a space area in the shape of a triangular prism. The receiving cable with multiple receiving probes is installed and laid flat in the receiving borehole, and the cable is connected to the receiving host. Electromagnetic wave perspective detection is carried out in the hole, using the following method: Figure 3 In the "one transmit, multiple receive" operating mode shown, each transmitting point corresponds to 9-11 receiving points. The transmitting point is continuously moved to create perspective coverage between the two holes, recording the spatial positions of the transmitting and receiving points and their corresponding measurements. By changing the combination of transmitting and receiving holes, a new perspective data volume is generated until all the designed holes have been traversed.

[0072] Step 2, 3D observation data processing: This is the process of converting the measured voltage signal into the absorption attenuation distribution of electromagnetic waves in 3D space. The approximate formula for electromagnetic wave propagation in space is:

[0073]

[0074] H is the signal value at the receiving point, H0 is the signal value at the transmitting probe, which is determined by the transmitting probe parameters, r is the distance between the transmitting point and the receiving point, and β is the attenuation coefficient of the medium. If the detection environment is a uniform space, β is a constant. If the detection environment is a non-uniform space, the electromagnetic wave absorption attenuation coefficient is different everywhere.

[0075] To deal with non-continuous problems, the discretization method is generally used in mathematics, that is, the interval to be explored is divided into small areas. For the discretization between three-dimensional detection holes, a three-dimensional discrete grid is used. The discrete grid partitioning scheme is shown in Figure 4 , the absorption attenuation coefficient in each spatial grid is regarded as a constant, and the attenuation of the signal from the emission point to the receiving point is regarded as the superposition of the attenuation on all grids through which the ray passes, transforming a continuous problem into a linear superposition problem, as follows:

[0076] The exponential decay term in formula (1) is written in discrete form, namely:

[0077]

[0078] After the three-dimensional continuous space is divided into discrete grids, any ray connecting the transmitting point and the receiving point will be divided into several discrete segments by the grid. In formula (2), N is the number of grids used to divide the ray. For different rays, the number of segments varies. By transposing the terms and performing logarithmic equivalence processing on formula (2), it is transformed into the following form:

[0079]

[0080] In formula (3), the right side of the equal sign is a specific value, which is determined by the power of the instrument itself and the data obtained by observation.i is the length of the ray in the i-th grid among the N grids it passes through, β i is the absorption attenuation coefficient of the i-th grid among the N grids that the ray passes through. By performing an equivalent mathematical transformation on formula (2), a multivariate linear equation is established between the received data value and the grid absorption attenuation value. By combining multiple rays, a combined multivariate linear equation system is formed and solved. In the traditional electromagnetic wave perspective data processing process, the method of calculating the length of each grid when the ray passes through a two-dimensional grid is relatively mature, but there is no mature method for calculating the length of each grid when the ray passes through a three-dimensional grid. The present invention proposes a method based on the spatial straight line equation to solve this problem:

[0081] Specifically, the length of the ray passing through each grid is calculated as follows:

[0082] The spatial coordinates of the transmitting point are (xf, yf, zf), and the spatial coordinates of the receiving point are (xr, yr, zr). The imaging space is divided in the X, Y, and Z directions. The distribution of the grid lines in each direction is as follows:

[0083] X direction: The starting x coordinate of the segmentation area is x0, and the ending coordinate is xn. The grid segmentation sequence in the x direction is x0, x1, x2, ... xn;

[0084] Y direction: The starting y coordinate of the segmentation area is y0, and the ending coordinate is yn. The grid segmentation sequence in the y direction is y0, y1, y2, ...yn;

[0085] Z direction: The starting z coordinate of the segmentation area is z0, and the ending coordinate is zn. The grid segmentation sequence in the z direction is z0, z1, z2, ... zn;

[0086] The ray equation is established with the coordinates of the transmitting point (xf, yf, zf) and the receiving point (xr, yr, zr) as follows:

[0087] Δx=xr-xf

[0088] Δy=yr-yf (4)

[0089] Δz=zr-zf

[0090] Δx, Δy, and Δz are the coordinate differences between the transmitting and receiving points in three dimensions. (Δx, Δy, and Δz) are the direction vectors from the transmitting point to the receiving point. The coordinates of any point on the line segment from the transmitting point to the receiving point and its extension can be determined by parametric equations. The ray equation is as follows:

[0091]

[0092] In formula (5), the range of the variable t is all real numbers greater than 0. When t = 1, the x, y, and z on the left side of the equal sign are the coordinates of the receiving point.

[0093] According to the grid characteristics, the intersection of the ray and the subdivided hexahedral unit must meet the following characteristics:

[0094] 1. There are 1 to 2 intersections between the ray and the outer surface of the hexahedron it passes through, and only the first and last hexahedrons through which the ray passes have only one intersection. The rest of the hexahedrons passed through have 2 intersections. The schematic diagram of the ray passing through the hexahedron is as follows Figure 5 As shown;

[0095] 2. If the intersection coordinates are (xc, yc.zc), they must satisfy at least one of the following conditions:

[0096] xc∈{x0,x1,…xn} or yc∈{y0,y1,…yn} or zc∈{z0,z1,…zn}; and the range can be further narrowed, that is: min(xf,xr)≤xc≤max(xf,xr) and min(yf,yr)≤yc≤max(yf,yr) and min(zf,zr)≤zc≤max(zf,zr).

[0097] min(a,b) is the smaller number between a and b, and max(a,b) is the larger number between a and b.

[0098] 3. In formula (5), the range of t varies between 0 and 1. The intersection of the ray segment and the grid surface in the X, Y, and Z directions can be calculated separately in each dimension. The calculation method is as follows:

[0099] X direction: The grid coordinates in the x direction between xf and xr are (xfr1, xfr2, …, xfrn). The sequence of all intersections of the transmit and receive lines and the grid (including the transmit and receive points) is (xf, xfr1, xfr2, …, xfrn, xr). According to the first equation of formula (5), the equivalent form is transformed to obtain:

[0100]

[0101] Let the right side x of formula (6) be equal to each element in the x-intersection sequence (xf, xfr1, xfr2, ..., xfrn, xr), and we can get the coefficient sequence: {0, t xfr1 ,t xfr2 ,...t xfrn,1}, the meaning of this coefficient sequence is the ratio of the distance between the grid intersection point and the transmitting point and the distance between the transmitting point and the receiving point in the x direction from the transmitting point (starting point) to the receiving point (end point). When Δx=0, there is no need to calculate the intersection sequence in the x direction.

[0102] Y direction: The grid coordinates in the y direction between yf and yr are (yfr1, yfr2, …, yfrn). The sequence of all intersections of the transmit and receive lines and the grid (including the transmit and receive points) is (yf, yfr1, yfr2, …, yfrn, yr). According to the second equation of formula (5), the equivalent form is transformed to obtain:

[0103]

[0104] Let the right side y of formula (7) be equal to each element in the y intersection sequence (yf, yfr1, yfr2, ..., yfrn, yr), and we can get the coefficient sequence: {0, t yfr1 ,t yfr2 ,...t yfrn ,1}, the meaning of this coefficient sequence is the ratio of the distance between the grid intersection point and the transmitting point and the distance between the transmitting point and the receiving point in the y direction from the transmitting point (starting point) to the receiving point (end point); when Δy=0, it is not necessary to calculate the intersection sequence in the y direction.

[0105] Z direction: The grid coordinates between zf and zr in the z direction are (zzfr1, zfr2, …, zfrn). The sequence of all intersections of the transmit and receive lines and the grid (including the transmit and receive points) is (zf, zfr1, zfr2, …, zfrn, zr). According to the third equation of formula (5), the equivalent form is transformed to obtain:

[0106]

[0107] Let the right side z of formula (8) be equal to each element in the z intersection sequence (zf, zfr1, zfr2, ..., zfrn, zr), and we can get the coefficient sequence: {0, t zfr1 ,t zfr2 ,...t zfrn ,1}, the meaning of this coefficient sequence is the ratio of the distance between the grid intersection point and the transmitting point and the distance between the transmitting point and the receiving point in the z direction from the transmitting point (starting point) to the receiving point (end point); when Δz=0, it is not necessary to calculate the intersection sequence in the z direction.

[0108] After obtaining the intersection coefficient sequence in the three directions, all coefficients are arranged in descending order and repeated values ​​are removed to obtain the final intersection representation value sequence: {0, t1, t2, ... t N,1}, N is the total number of intersections between the ray and the three-dimensional grid, and the characterization value refers to the ratio of the distance from the intersection point to the emission point (starting point) to the distance between the receiving and transmitting points.

[0109] The three-dimensional grid divides the transmitting-receiving line into N+1 small segments, each of which belongs to a different spatial grid. The length of the ray in each grid is calculated according to the adjacent coordinates and a linear equation is established according to formula (3). All rays are integrated to construct a linear algebraic equation system to obtain the three-dimensional distribution of the absorption attenuation coefficient β in space.

[0110] Step 3: Based on the absorption attenuation coefficient, the K-means clustering algorithm is used to perform anomaly segmentation to achieve three-dimensional detection of anomalies within the working face. Specifically, the spatial clustering analysis and attribute calculation classification of anomalies based on attenuation values ​​is the process of grouping the attenuation coefficient according to its numerical value, then confirming the anomaly demarcation threshold and determining the anomaly scale. The unsupervised K-means clustering algorithm is used to perform anomaly segmentation. This process is divided into three steps: anomaly threshold division, anomaly spatial location classification (using K-means solution), and anomaly attribute classification analysis. The specific execution process of each step is as follows:

[0111] Division of abnormal thresholds:

[0112] The purpose of abnormal threshold division is mainly to establish a segmentation standard to distinguish abnormal areas from background areas, determine a specific value as the threshold, and the absorption attenuation coefficient is less than / greater than this value to be demarcated as the normal background area, and the absorption attenuation coefficient is greater than / less than this value to be demarcated as the abnormal area. As for which division scheme to choose, it needs to be determined according to the specific detection target. For collapse columns, faults, and water-bearing areas, the scheme greater than the threshold is selected to determine the abnormal area. For karst cavities, goaf tunnels, goaf areas and other targets, the scheme less than the threshold is selected to determine the abnormal area. The determination of the threshold should be determined according to the distribution of the absorption attenuation histogram. The determination process is as follows:

[0113] Count the maximum and minimum values ​​of the attenuation coefficients of all grids in the detection area, divide the maximum and minimum values ​​into 10 intervals, calculate the number (frequency) of grid attenuation coefficients falling in each interval, and form an absorption attenuation frequency distribution histogram. The frequency distribution diagram is as follows: Figure 6 As shown in the figure, the shape of the histogram varies greatly under different detection environments. After the attenuation data are arranged from small to large, the smaller value less than the 10%-20% quantile or the larger value greater than the 80%-90% quantile is selected as the detection abnormal area. For collapse columns, faults, and water-bearing areas, the 80% quantile value is generally selected as the threshold for distinguishing abnormalities from non-abnormalities. For caves, goaf tunnels, goaf areas, etc., the 20% quantile value is selected as the threshold for distinguishing abnormalities from non-abnormalities.

[0114] Classification of abnormal spatial locations:

[0115] After the anomaly threshold is determined, the K-means clustering algorithm is used to spatially classify the anomalies. All grid cells exceeding the threshold are classified. A series of adjacent grid cells are considered to be the same anomaly. Grid cells that are not connected to each other are considered to be other independent anomalies. Clustering is performed based on cell distance. The method is as follows:

[0116] 1) For each grid in the detection area, the spatial distance between the i-th grid and the j-th grid is defined as follows:

[0117]

[0118] Among them, the coordinates of the center point of the i-th grid are (x i ,y i ,z i ), the coordinates of the jth grid center point are (x j ,y j ,z j );

[0119] The distances between all abnormal grids form a data volume. The two grids corresponding to the smaller elements in the data volume belong to the same class. The spatial distances between all abnormal grids are calculated to form a distance matrix D, where:

[0120]

[0121] The matrix D is a symmetric square matrix, d 12 The distance between the first abnormal grid element and the second abnormal grid element. The closer the value is to 0, the closer the grids where the two are located are in space, and the more likely they are to belong to the same class.

[0122] 2) Set a threshold ε to determine whether they belong to the same class. ε is a relatively small real number, generally selected from 0.1 to 0.2 times the largest element in the matrix D. If the difference between the elements is less than ε, they are considered to belong to the same class. If it exceeds ε, they belong to two classes.

[0123] 3) Determine the row i containing the minimum value other than D(1,1) in the first column of the matrix D(:,1). If D(i,1)≤ε, the i-th element and the first element in the absorption attenuation data volume form a class cluster1; and record the average value and element number of all elements in cluster1. At this time, there are 2 elements in cluster1, and their average is meancluster1; if D(i,1)>ε, the first element and the i-th element form two classes, cluster1 and cluster2 respectively. Class cluster1 has the first element, and class cluster2 has the first element. Calculate the mean value of the elements in the two classes respectively.

[0124] 4) Determine whether the remaining elements belong to the existing clusters cluster1 and cluster2. If they do, place the element in that cluster. If not, construct a new cluster3, update the numbers of the elements in each cluster, the class mean, and the numbers of the remaining elements.

[0125] Repeat step 4) until all elements have been classified and a total of K classes are formed;

[0126] 5) Perform three-dimensional anomaly display on the K clusters obtained by clustering algorithm analysis, and obtain the data processing results of three-dimensional electromagnetic wave perspective as shown in the figure below: Figure 7 .

[0127] Abnormal attribute analysis classification:

[0128] After obtaining the elements in K classes, each class corresponds to an independent anomaly. The characteristics of the elements in the class are analyzed to obtain different attribute characteristics of the anomaly. The attribute characteristics include but are not limited to: the scale of the anomaly, the boundary shape of the anomaly, the direction characteristics of the anomaly, etc. The specific properties of the anomaly are inferred through these quantitative characteristics.

[0129] Analysis of the scale of the anomaly:

[0130] Each cluster formed by cluster analysis can be considered an independent anomaly. Since the coordinate grid is a hexahedron, the volume of each grid is known and easy to calculate. Simply summing the volumes of the grids corresponding to all elements in the cluster provides a rough estimate of the anomaly's size (the grids corresponding to the actual cluster elements may not fit perfectly). To further accurately calculate the anomaly's size, the enclosing surface of the anomaly grid can be calculated before performing a precise calculation. Image boundary recognition technology offers a rich set of solutions for calculating the boundaries of 2D / 3D graphics.

[0131] Analysis of the abnormal body's trend characteristics:

[0132] The grids corresponding to the elements in the class are anomaly grids. All grids constitute an independent anomaly body. Analyzing the spatial trend characteristics of the anomaly body can facilitate further analysis of the nature of the anomaly body. The distribution of the anomaly body in three-dimensional space is generally related to its nature. Linear structures generally correspond to fault structures or goafs; structures with little difference in long and short axes on the horizontal plane are generally collapse column structures. Other geological phenomena have their own characteristics. The specific methods are as follows:

[0133] 1) Extract the spatial coordinates of the grid corresponding to the class element and count the center positions of all grids;

[0134] 2) The principal component analysis method is used to obtain the dominant direction, subdominant direction, insignificant direction and corresponding length of the grid unit distribution, and the trend characteristics of the anomaly are obtained by analysis, that is, the dominant development direction, scale and other information of the anomaly

[0135] 3) Determine the type and nature of anomalies based on their strike characteristics. Different anomaly types have different spatial distribution patterns. Fault structural anomalies have an absolute developmental advantage along a certain direction in space. The difference between the long and short axes of collapse column structural anomalies on the horizontal projection plane is not very obvious. For other human anomalies (such as goaf tunnels and goaf areas), further analysis is conducted in combination with available historical mining data.

Claims

1. A three-dimensional electromagnetic perspective exploration method for the internal structure of an underground coal mine working face, characterized in that: This method is based on gas extraction drilling holes constructed alternately at high and low levels within the coal seam working face. It uses electromagnetic perspective with simultaneous measurement of multiple receivers and one emission to achieve three-dimensional detection of abnormal bodies within the working face. The method includes the following steps: Step 1, drilling and probe layout and detection: Arrange two rows of drill holes in a broken line, determine the transmitting and receiving drill holes, and use a single-transmitter-multiple-receiver configuration to deploy the transmitting and receiving probes for detection data collection, recording the spatial positions of the transmitting and receiving points and the corresponding measurement values; Step 2, 3D observation data processing: Divide the detection area into 3D discrete grids, convert the exponential attenuation term in the approximate formula for electromagnetic wave propagation in space into a discrete form, and perform a logarithmic equivalent transformation to solve for the absorption attenuation coefficient; The approximate formula for electromagnetic wave propagation in space is: （1） In the above formula, H is the signal value at the receiving point, H 0 is the signal value at the transmitting probe, r is the distance between the transmitting point and the receiving point, and β is the absorption attenuation coefficient of the medium; The exponential decay term is converted into discrete form: （2） In the above formula, H is the signal value at the receiving point, H 0 is the signal value at the transmitting probe, r is the distance between the transmitting point and the receiving point, and N is the number of grids that split the ray. r i is the number of grids that the ray passes through i The length of the grid, β i is the number of N grids that the ray passes through i The absorption attenuation coefficient of each grid; Transpose the terms of formula (2) and take logarithmic equivalent processing to transform it into the following formula (3). Combine multiple rays and solve to obtain the absorption attenuation coefficient; （3）； Step 3: Based on the absorption attenuation coefficient, the K-means clustering algorithm is used to perform anomaly segmentation to achieve three-dimensional detection of anomalies inside the working face; The step 1 comprises: Step 1.1: In a lane on one side of the working face or in lanes on both sides, drill two rows of holes in an alternating pattern. The holes should be oriented perpendicular to the direction of the coal mine lanes, with the trajectory kept horizontal. The spacing between holes at the same height should be 10-15m, and the drilling depth should be 100m-120m. Step 1.2: Use the upper row of drill holes as transmitting drill holes and the lower row as receiving drill holes. Each transmitting drill hole forms a perspective combination with the 2-4 nearest receiving drill holes. In step 1.3, one electromagnetic wave transmitting probe is placed in the transmitting borehole, and 9-11 receiving probes are placed in the receiving borehole. The transmitting and receiving probes move synchronously, exiting from the bottom of the hole toward the hole mouth point by point for detection; the perspective combination of the transmitting borehole and the receiving borehole is replaced until all boreholes are traversed; the spatial positions of the transmitting and receiving points and the corresponding measurement values ​​are recorded.

2. The three-dimensional electromagnetic perspective exploration method for the internal structure of the working face in a coal mine according to claim 1, characterized in that: In step 2, when the ray passes through N grids, the coordinate difference between the ray emission point coordinate (xf, yf, zf) and the receiving point coordinate (xr, yr, zr) in the three dimensions of X, Y, and Z directions is: （4） The ray equation is established as follows: （5） In formula (5), the variable t The value is 0-1.

3. The three-dimensional electromagnetic perspective exploration method for the internal structure of an underground coal mine working face according to claim 2, characterized in that: In step 2, the sequence of intersections of the ray with each grid in the X direction is (xf, xfr1, xfr2, …, xfrn, xr). According to the first equation of formula (5), the equivalent form is transformed to obtain: （6） Let the x on the right side of the equation (6) be equal to each element in the X-direction intersection sequence (xf, xfr1, xfr2, …, xfrn, xr), and the coefficient sequence is: , the coefficient sequence is: in the X direction, from the emission point to the receiving point, the ratio of the distance from the intersection of the ray and the grid to the emission point to the distance between the emission point and the receiving point.

4. The three-dimensional electromagnetic perspective exploration method for the internal structure of an underground coal mine working face according to claim 3, characterized in that: In step 2, the sequence of intersections of the ray with each grid in the Y direction is (yf, yfr1, yfr2, …, yfrn, yr). According to the second equation of formula (5), the equivalent form is transformed to obtain: （7） Let the right side y of the equation (7) be equal to each element in the Y-direction intersection sequence (yf, yfr1, yfr2, …, yfrn, yr), and the coefficient sequence is: , the coefficient sequence is: in the Y direction, from the emission point to the receiving point, the ratio of the distance from the intersection of the ray and the grid to the emission point to the distance between the emission point and the receiving point.

5. The three-dimensional electromagnetic perspective exploration method for the internal structure of the working face in a coal mine according to claim 4, characterized in that: In step 2, the sequence of intersection points of the ray with each grid in the Z direction is (zf, zfr1, zfr2, ..., zfrn, zr). According to the third equation of formula (5), the equivalent form is transformed to obtain: （8） Let the right side z of formula (8) be equal to each element in the Z-direction intersection sequence (zf, zfr1, zfr2, …, zfrn, zr), and the coefficient sequence is obtained: , the coefficient sequence is: in the Z direction, from the emission point to the receiving point, the ratio of the distance from the intersection of the ray and the grid to the emission point to the distance between the emission point and the receiving point; Combining the intersection coefficient sequences in the above three directions, we can get the ratio of the distance from the intersection point to the transmitting point to the distance between the transmitting point and the receiving point, which is the intersection characterization value sequence. .

6. The three-dimensional electromagnetic perspective exploration method for the internal structure of an underground coal mine working face according to claim 1, characterized in that: The step 3 comprises: Abnormal threshold division: Count the maximum and minimum values ​​of all absorption attenuation coefficients in the detection area, divide the maximum and minimum values ​​into 10 intervals, calculate the number of absorption attenuation coefficients falling in each interval, form an absorption attenuation frequency distribution histogram, and select the quantile value less than 10%-20% or greater than 80%-90% as the detection abnormal area; for collapse columns, faults, and water-bearing areas, select the 80% quantile value as the threshold for distinguishing abnormalities from non-abnormalities; for caves, goaf tunnels, and goaf areas, select 20% as the threshold for distinguishing abnormalities from non-abnormalities; Classification of abnormal spatial locations: After the abnormal threshold is determined, the K-means clustering algorithm is used to spatially classify the abnormalities. All grids exceeding the threshold are classified. A series of adjacent grid cells belong to the same abnormality, and disconnected grid cells are other independent abnormalities. Clustering is performed according to the cell distance to obtain K clusters; Anomaly attribute analysis and classification: After obtaining the elements in K classes, each class corresponds to an independent anomaly. The characteristics of the elements in the class are analyzed to obtain different attribute characteristics of the anomaly. The attribute characteristics include but are not limited to: the scale of the anomaly, the boundary shape of the anomaly, and the direction characteristics of the anomaly. The specific properties of the anomaly are inferred through the attribute characteristics.

7. The three-dimensional electromagnetic perspective exploration method for the internal structure of an underground coal mine working face according to claim 6, characterized in that: The K-means clustering algorithm is used to perform spatial classification of anomalies, including: 1) The spatial distance between the i-th grid and the j-th grid in the detection area is defined as follows: （9） Among them, the coordinates of the center point of the i-th grid are ( x i ,y i ,z i ), the coordinates of the jth grid center point are ( x j ,y j ,z j ); The distances between all abnormal grids constitute a data body, and the spatial distances between all abnormal grids are calculated to form a distance matrix. D ,in: (10) Distance Matrix D is a symmetric square matrix, d 12 is the distance between the first abnormal grid element and the second abnormal grid element. The closer the value is to 0, the closer the two are in space. 2) Set a threshold ε to determine whether elements belong to the same class. If the difference between elements is less than ε, they belong to the same class. If it exceeds ε, they belong to two classes. 3) Determine the distance matrix D The minimum value in the first column other than D(1,1), if the minimum value is not greater than ε, then the minimum value and the first element in the absorption attenuation data volume form a class cluster1; and record the average value and element number of all elements in cluster1. At this time, there are 2 elements in cluster1, and their average is meancluster1; if the minimum value is greater than ε, then the first element and the minimum element form two classes, cluster1 and cluster2 respectively. Class cluster1 has the first element, and class cluster2 has the first element. Calculate the mean value of the elements in the two classes respectively; 4) Determine whether the remaining elements belong to the existing clusters cluster1 and cluster2. If they do, place the element in that cluster. If not, construct a new cluster cluster3, update the numbers of the elements in each class, the class mean, and the numbers of the remaining elements. This continues until all elements are classified and a total of K clusters are formed. 5) Perform three-dimensional anomaly display on the K clusters obtained by clustering algorithm analysis.

8. The three-dimensional electromagnetic perspective exploration method for the internal structure of an underground coal mine working face according to claim 6, characterized in that: The step 3 also includes an analysis of the scale of the anomaly: each cluster formed by clustering is regarded as an independent anomaly, and the sum of the grid volumes corresponding to all elements in the cluster is calculated to obtain the volume of the anomaly.

9. The three-dimensional electromagnetic perspective exploration method for the internal structure of an underground coal mine working face according to claim 6, characterized in that: Step 3 also includes an analysis of the anomaly's trend characteristics: extracting the spatial coordinates of the grids corresponding to the class elements and counting the center positions of all grids; sequentially using the principal component analysis method to obtain the dominant direction, sub-dominant direction, insignificant direction, and corresponding lengths of the grid unit distribution, and analyzing to obtain the trend characteristics of the anomaly, i.e., the dominant development direction and scale of the anomaly; and judging the type and nature of the anomaly based on the anomaly's trend characteristics.

Images