An irregular working face positioning method coupling surface deformation field and mining subsidence model

Abstract

The present application relates to a kind of irregular working face positioning method of coupling surface deformation field and mining subsidence model, geological parameters are inverted by mode vector method, and working face in goaf is fitted by morphological algorithm.The method is combined by setting mode vector method and morphological algorithm, and geological parameters determine the optimal spatial position solution of underground working face, and spatial position is fed back to evaluate the accuracy of geological parameter selection by evaluation index value, the two are closely related, and the overall algorithm is nested in the geological parameter inversion process Morphological algorithm inversion working face boundary, this inversion method involves multiple mining units multiple surface point movement deformation multiple cumulative calculation, and the amount of calculation is huge, and mode vector method has the characteristics of high running efficiency and fast parameter convergence, and the combination of the two realizes the integrated inversion of geological parameters and working face spatial position.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of geological exploration, and particularly relates to an irregular working face positioning method coupled with a surface deformation field and a mining subsidence model. BACKGROUND

[0002] Accurate acquisition of spatial feature information of underground working face is related to correct identification and accurate judgment of illegal mining such as underground cross-border mining. The traditional illegal mining supervision method has problems of low efficiency, poor precision and weak reliability. In recent years, many scholars have applied the related methods of positioning underground goaf by geophysical prospecting technology to the identification of underground illegal mining. This kind of exploration technology, such as seismic exploration method, transient electromagnetic method and ground penetrating radar method, mainly uses the physical characteristics (density, magnetic field, gravity field, etc.) of geological bodies to analyze stratum lithology, geological structure, ore body distribution, etc., and combines the differences of rock-soil physical characteristics of goaf to identify and position the underground mining area. This kind of method has higher requirements for measuring equipment and professional knowledge and technology, and the precision of working face positioning in complex geological environment needs to be improved.

[0003] On the basis of identifying underground mining deformation, in order to further accurately position the underground mining space, the existing technology uses geological parameters to describe the position form of underground rectangular working face, and uses genetic algorithm and annealing algorithm to invert the geological parameters of underground working face, realizes high-precision inversion of the position of underground working face, further, using InSAR monitoring data, on the basis of rectangular working face assumption, different optimization inversion algorithms are used to realize the inversion positioning of underground working face.

[0004] However, the shape of underground working face is affected by the occurrence of coal seam, ground protection zone and the spatial position of adjacent mined working face, plus the concealment characteristics of illegal mining, the working face is often not a standard rectangle, so the existing method is difficult to realize the accurate inversion positioning of such irregular working face, which limits the application of the method in the field of illegal mining identification. SUMMARY

[0005] The purpose of the present application is to provide an irregular working face positioning method coupled with a surface deformation field and a mining subsidence model to solve the above problems.

[0006] The present application achieves the above-mentioned purpose by the following technical solutions:

[0007] An irregular working face positioning method coupled with a surface deformation field and a mining subsidence model, comprising the following steps:

[0008] A spatial coordinate system is established, the measured data of the surface deformation field is obtained to initialize the geological parameters of the working face, and a predicted working face mask matrix is constructed based on the initialized working face;

[0009] The mask matrix is dilated and eroded based on an image morphological algorithm, and a residual sum of squares of the deformation prediction value and measured data of the ground deformation field is taken as a fitting effect evaluation index to obtain the mask matrix after fitting;

[0010] The geological parameters of the initial working face are iterated based on a parameter optimization method, and the iteration step of the geological parameters is updated according to a preset rule; the preset rule is specifically that the lower the index value of the evaluation index is, the better the corresponding geological parameter is;

[0011] Based on the predicted working face mask matrix of the updated geological parameters, the above steps are repeated until the iteration step is lower than a preset precision, and the optimal geological parameters and the predicted working face are output.

[0012] As a further optimization scheme of the present application, the spatial coordinate system includes H axis, X axis, Y axis, and the measured data of the ground deformation field is specifically the movement along the H axis, X axis and Y axis; the geological parameters include m, H0, n1 and n2, wherein m represents the coal seam thickness, H0 represents the coal seam buried depth of the maximum subsidence point, n1 and n2 represent the spatial normal vector based on the working face X axis and Y axis components.

[0013] As a further optimization scheme of the present application, the geological parameters of the initial working face include:

[0014] The maximum subsidence point coordinates (X0, Y0) are obtained:

[0015] The initial working face passes through the P0 (X0, Y0, H0) point, and the plane equation of the initial working face is:

[0016] n1 (X-X0) +n2 (Y-Y0) +n3 (H-H0) =0 (1).

[0017] As a further optimization scheme of the present application, the deformation prediction value is obtained through the geological parameters G, the inclination angle α of the working face and the tendency azimuth angle and is obtained based on the probability integral method;

[0018] Wherein, the inclination angle α is represented as:

[0019]

[0020] The tendency azimuth angle is represented as:

[0021]

[0022] As a further optimization scheme of the present application, the step of constructing the prediction working face mask matrix based on the initialization working face comprises: dividing the plane where the initialization working face is located into a grid, taking the grid as a grid unit to construct the mask matrix Mask, setting the initial value of all grid units as 0, and setting the value of the grid unit where the position of the maximum subsidence point in Mask as 1 as a seed unit; and setting the size of the expansion and erosion structure element as 3*3.

[0023] As a further optimization scheme of the present application, the step of expanding and eroding the mask matrix based on the image morphological algorithm specifically comprises:

[0024] S1: performing once expansion operation on all grid units with the value of 1 in the current Mask, taking the residual sum of squares of the surface deformation field measured data and the deformation prediction value after expansion and erosion as the fitting effect evaluation index, and performing fitting effect evaluation on the new units of the to-be-determined edge generated after expansion one by one, if the addition of the expansion unit reduces the evaluation index, that is, improves the fitting effect, then the new grid unit obtained by expansion is effective, and the grid unit is retained; otherwise, it is discarded;

[0025] S2: setting the value of the grid unit corresponding to all effective expansion units in Mask as 1 to realize the expansion of Mask;

[0026] S3: iteratively performing steps S1 and S2 until all new grid units generated by a certain expansion are invalid expansion units;

[0027] S4: performing once erosion operation on all grid units with the value of 1 in the current Mask, and performing fitting effect evaluation on the to-be-determined units of the edge one by one, if the discarding of the grid unit reduces the evaluation index, that is, improves the fitting effect, then the grid unit eroded is effective, that is, the grid unit is discarded; otherwise, it is retained;

[0028] S5: setting the value of the grid corresponding to all effective erosion units in Mask as 0 to realize the erosion shrinkage of Mask;

[0029] S6: iteratively performing steps S4 and S5 until the edge new to-be-determined grid units obtained by a certain erosion are all invalid erosion units, that is, any unit cannot be discarded by erosion;

[0030] S7: outputting all grid units with the value of 1 as the fitted mask matrix Mask.

[0031] As a further optimization scheme of the present application, the evaluation index is obtained in the following manner: taking the square sum of the module length of the surface three-dimensional point movement vector residual as the evaluation index value f:

[0032]

[0033] In formula (4), W H , U X , U Y are respectively the movement amount of the measured ground deformation field along the H-axis, the X-axis and the Y-axis, w H , u X , u Y are respectively the movement amount of the predicted deformation value along the H-axis, the X-axis and the Y-axis, S1 and S2 are the spatial movement vectors, and ΔS is the module length of the vector difference.

[0034] As a further optimization scheme of the present application, in step S3, the process of updating the iteration step length of the geological parameter G according to the preset rule is specifically:

[0035] An initial step length Δ and an accuracy d are set;

[0036] The initial geological parameter G is iterated one by one based on the initial step length Δ, if the evaluation index is lower after iteration, the geological parameter G is updated as a better geological parameter G and the step length Δ is increased, otherwise the geological parameter G is retained and the step length Δ is reduced.

[0037] As a further optimization scheme of the present application, the increase coefficient is 2 times when the step length Δ is increased, and the reduction coefficient is 0.5 times when the step length Δ is reduced.

[0038] The present application has the following beneficial effects:

[0039] The present application combines the modal vector method with the morphological algorithm, and the geological parameter determines the optimal spatial position solution of the underground working face, and the spatial position feeds back the accuracy of the geological parameter selection through the evaluation index value, the two are closely related, and the algorithm as a whole is nested in the morphological algorithm inversion working face boundary in the geological parameter inversion process. This inversion method involves multiple cumulative calculations of the movement and deformation of multiple mining units on multiple ground points, and the calculation amount is huge, while the modal vector method has the characteristics of high running efficiency and fast parameter convergence, and the combination of the two realizes the integrated inversion of the geological parameter and the working face spatial position. BRIEF DESCRIPTION OF DRAWINGS

[0040] Figure 1 is the overall flowchart of the present application;

[0041] Figure 2 is the iteration process flowchart of the present application;

[0042] Figure 3 is a schematic diagram of the present application for calculating the evaluation index value f based on the residual of the ground point movement vector;

[0043] Figure 4 is a working face plane description schematic diagram of the present application;

[0044] Figure 5 is a schematic diagram of verification results of the present application; DETAILED DESCRIPTION

[0045] The present application will be further described in detail below with reference to the accompanying drawings. It is necessary to point out here that the following detailed description is only used to further illustrate the present application and cannot be understood as limiting the scope of protection of the present application. Those skilled in the art can make some non-essential improvements and adjustments to the present application according to the above application content.

[0046] Example 1

[0047] As shown in Figures 1-4 , an irregular working face positioning method coupling surface deformation field and mining subsidence model adopts the modulus vector method to invert four geological parameters m, H0, n1 and n2. Each set of geological parameters calls a morphological algorithm function to invert and determine a set of optimal underground working face positions under the geological conditions, and returns a surface deformation fitting effect evaluation index. Through continuous iteration, the geological parameters G and the underground working face position gradually approach the optimal value.

[0048] In the problem of underground illegal mining identification and other problems, these geological parameters G are often difficult to obtain in advance, therefore, the underground working face positioning problem first needs to determine the related geological parameters. In the probability integral method model, the geological parameters involved mainly include coal seam thickness (mining thickness) m, coal seam inclination α, inclination azimuth n1, and coal seam buried depth H. These geological parameters must be determined as unknown quantities through algorithm inversion method, so as to improve the applicability of the algorithm.

[0049] Among these geological parameters, the inclination azimuth has the characteristics of large parameter range (generally set between 0-360°, covering all possible directions on the entire plane) and obvious periodicity (period of 360°). When the coal seam is close to horizontal, the inclination azimuth has weak correlation with the overall subsidence fitting effect, which may cause large fluctuations and failure to converge in the search optimization process, or always search for the boundary value, and even seriously affect the inversion of other geological parameters.

[0050] Therefore, according to the basic knowledge of space geometry, the problem is transformed into a spatial plane positioning problem. As shown in Figure 4 , an XYH coordinate system is constructed, in which the H axis points downward to describe the underground mining depth. The normal vector of the spatial plane and the spatial point P0(X0, Y0, H0) passing through the plane are used to determine the geometric plane where the working face is located. The normal vector n1, n2, n3 in the formula respectively are components of the normal vector along the X, Y, H axis, the degree of freedom is 2, n3 can be directly set as 1. P0(X0, Y0, H0) is a passing point of the space plane, the P0 point should be selected close to the working face, and the H0 initial value range is determined conveniently. P0(X0, Y0, H0) can be set as the coordinates of the maximum subsidence point in the ground subsidence basin and the corresponding position coal seam depth, at this time, the H0 value range can be directly determined according to the experience range of the mining depth in China. Therefore, the working face geological parameters have four, which are the coal seam thickness m, the coal seam depth H0 at the maximum subsidence point, the components n1, n2 of the normal along the X axis and the Y axis. As shown in the formula, Figure 4 The plane equation of the working face is constructed as shown in the formula.

[0051] n1(X-X0)+n2(Y-Y0)+n3(H-H0)=0 (1);

[0052] The inclination angle alpha is represented as:

[0053]

[0054] The inclination azimuth angle is represented as:

[0055]

[0056] Based on the above parameters, the positioning method of the application specifically includes the following steps:

[0057] Based on the remote sensing technology, the ground surface deformation field measured data and the maximum subsidence point coordinates (X0, Y0) are obtained, and the geological parameters G of the working face are initialized, a mask matrix Mask is constructed based on the initialized working face, and the maximum subsidence point coordinates (X0, Y0) are taken as the seed unit of the mask matrix Mask;

[0058] Based on the morphological algorithm, the seed unit of the mask matrix Mask is expanded and eroded, the residual sum of squares of the ground surface deformation field measured data and the deformation prediction value after expansion and erosion is taken as a fitting effect evaluation index, the expansion and erosion results are retained / discarded according to the evaluation index, and the fitted mask matrix Mask is obtained, wherein the ground surface deformation field measured data includes the movement along the H axis, the X axis and the Y axis, and the deformation prediction value is obtained based on the geological parameters G, the inclination angle alpha of the working face and the inclination azimuth angle

[0059] The basic idea of the method is that the entire grid plane is regarded as a binary raster image, the residual sum of squares of the ground surface deformation field measured data and the deformation prediction value is taken as a fitting effect evaluation index, and through multiple expansion and erosion of the grid unit, gradual approximation to the optimal working face position is realized, and the specific steps are as follows:​​

[0060] (1) the plane where the working face is located is divided into a grid, the grid division interval d has a certain influence on the boundary inversion accuracy, the grid division interval d should not be greater than H0 / 10, the grid is taken as a grid unit to construct a mask matrix Mask, the initial value of all grid units is set to 0; and the value of the grid unit where the position of the maximum subsidence point in Mask is set to 1 as a seed unit; the size of the expansion and erosion structure element is set to 3*3;

[0061] (2) all grid units with a value of 1 in the current Mask are subjected to an expansion operation, and the residual sum of squares of the measured data of the surface deformation field and the deformation prediction value after expansion and erosion is taken as the fitting effect evaluation index A. The new units of the to-be-determined edge generated after expansion are evaluated one by one. If the addition of the expansion unit reduces the evaluation index A, that is, the fitting effect is improved, the new grid unit obtained by expansion is effective, and the grid unit is retained. Otherwise, it is discarded;

[0062] (3) the values of the grid units corresponding to all effective expansion units in Mask are set to 1 to realize the expansion of Mask;

[0063] (4) steps (2) and (3) are iteratively executed until all new grid units generated by a certain expansion are invalid expansion units;

[0064] (5) all grid units with a value of 1 in the current Mask are subjected to an erosion operation, and the to-be-determined units of the edge after erosion are evaluated one by one. If the discarding of the grid unit reduces the evaluation index A, that is, the fitting effect is improved, the grid unit after erosion is considered effective, that is, the grid unit is discarded. Otherwise, it is retained;

[0065] (6) the values of the grid units corresponding to all effective erosion units in Mask are set to 0 to realize the erosion shrinkage of Mask;

[0066] (7) steps (5) and (6) are iteratively executed until the edge new to-be-determined grid units obtained by a certain erosion are all invalid erosion units, that is, any unit cannot be discarded by erosion;

[0067] (8) all grid units with a value of 1 are taken as the fitted mask matrix Mask and output.

[0068] In the scheme, the overall fitting effect of the deformation field adopts the square sum of the module length of the residual of the three-dimensional point movement vector of the surface as the evaluation index, as shown in formula (1), the lower the index value f, the better the fitting effect; otherwise, the worse the fitting effect: Figure 3

[0069]

[0070] In the formula: W​H , U X , U Y are respectively the measured displacement of the ground deformation field along the H-axis, X-axis, Y-axis, w H , u X , u Y are respectively the predicted displacement of the deformation along the H-axis, X-axis, Y-axis, S1, S2 are the spatial movement vectors, and ΔS is the length of the vector difference.

[0071] The initialization geological parameters G of the working face are iterated one by one based on the PS method, and the specific process is shown in FIG. Figure 2 The initialization values of the four geological parameters G are represented as G = {G[1], G[2], G[3], G[4]}. The initial step length Δ is set as Δ = {Δ[1], Δ[2], Δ[3], Δ[4]}, and the accuracy d is set as d = {d[1], d[2], d[3], d[4]}. The initialization geological parameters G[i] are iterated one by one based on the initial step length Δ. The iterated geological parameters are G[i] + Δ[i] or G[i] - Δ[i]. The fitted mask matrix Mask is reacquired based on the iterated geological parameters G, and the evaluation index value f of the Mask is obtained. If the evaluation index value f is lower after iteration, the geological parameters G are updated as the better geological parameters G, and the step length Δ is doubled. Otherwise, the geological parameters G are kept and reduced to 0.5 times of the original step length Δ. The mask matrix Mask is reconstructed based on the working face of the better geological parameters G, and the next geological parameter G[i+1] is iterated continuously. After four iterations, that is, all four geological parameters G are iterated once, it is considered that one cycle of iteration is completed.

[0072] After completing one cycle of iteration, it is determined whether the step length Δ = {Δ[1], Δ[2], Δ[3], Δ[4]} is less than the accuracy d = {d[1], d[2], d[3], d[4]}. If yes, the final geological parameters G and the mask matrix Mask are output, and the predicted working face is obtained through the final geological parameters G and the mask matrix Mask. Otherwise, the mask matrix Mask is reconstructed based on the working face of the better geological parameters G, and the iteration of the next cycle is repeated continuously.

[0073] The geological parameters G determine the optimal spatial position solution of the underground working face, and the spatial position feeds back the accuracy of the selection of the geological parameters G through the evaluation index value f, and the two are closely related. Therefore, in the inversion process of the geological parameters G, the MOR algorithm is nested to invert the working face boundary. This inversion method involves multiple cumulative calculations of the displacement deformation of multiple mining units on multiple ground points, and the calculation amount is huge. The PS method has the characteristics of high running efficiency and fast parameter convergence, and therefore the method combining the PS method and the MOR algorithm (PS&MOR) is proposed to realize the integrated inversion of the geological parameters G and the spatial position of the working face.

[0074] To verify the reliability of the algorithm, the experiment first designs a set of irregular working face under the geological parameters, through the prediction of surface subsidence and north-south, east-west horizontal movement, and takes it as the surface measured deformation field data, and then inverts the geological parameters and working face boundary, compares the designed parameters and boundary, and evaluates the inversion effect of the algorithm.

[0075] The maximum subsidence angle θ in the probability integral parameter is related to the dip angle α, and θ should change with the value of the geological parameter α. According to the mathematical relationship between the two determined by formula (7), the coefficient k is replaced as the corresponding probability integral parameter.

[0076] θ = 90 - kα (5)

[0077] An L-shaped working face is designed as a case, the ground measuring point spacing is 30m, covering the entire surface mining influence range, the unit working face grid division spacing is 20m, and other design parameters are shown in the following table:

[0078]

[0079] The surface subsidence and east-west, north-south horizontal movement predicted by the probability integral method.

[0080] In order to evaluate the robustness of the working face position inversion result of the algorithm under different parameter initial value combinations, 9 groups of orthogonal inversion experiments under 4 parameters with different initial value combinations are designed combined with the designed parameter values, and the geological parameter inversion results are shown in the following table:

[0081]

[0082] Through the analysis of 9 times of orthogonal experiments of parameter initial value, the mean error of mining height m, mining depth H0, dip angle α and dip direction azimuth ζ is 0.03m, 6.54m, 0.35° and 1.04° respectively, and the relative mean error is 1%, 2.2%, 3.5% and 0.7% respectively. It can be considered that the influence of parameter initial value on the inversion result is small, which meets the inversion accuracy requirement of underground working face geological parameters.

[0083] The three-dimensional form of the underground working face under the inversion result of the 9th group is shown in Figure 5 It can be seen that the working face boundary obtained by inversion is relatively accurate, the boundary plane error is not more than 20m, and the angle point plane position deviation is about 60m. The maximum error of the angle point mining depth is about 15m, which is located at F point. It can be considered that the algorithm is not sensitive to the selection of parameter initial value, and has certain anti-initial value interference ability.

[0084] The above embodiments only express several implementation manners of the present application, and the description is more specific and detailed, but it should not be understood as a limitation on the patent scope of the present application. It should be noted that for ordinary skilled persons in the art, without departing from the concept of the present application, several modifications and improvements can be made, which are within the protection scope of the present application.

Claims

1. An irregular working face positioning method for coupling a surface deformation field and a mining subsidence model, characterized in that: The method comprises the following steps: establishing a spatial coordinate system, obtaining measured data of a ground surface deformation field, initializing geological parameters of a working face, and constructing a prediction working face mask matrix based on the initialized working face; performing inflation and corrosion on the mask matrix based on an image morphological algorithm, taking a residual sum of squares of deformation prediction values and measured data of the ground surface deformation field as a fitting effect evaluation index, and obtaining a fitted mask matrix; iterating the geological parameters of the initialized working face based on a parameter optimization method, and updating an iteration step of the geological parameters according to a preset rule; the preset rule is specifically that the lower the index value of the evaluation index is, the better the corresponding geological parameter is; repeating the above steps based on the prediction working face mask matrix of the updated geological parameters until the iteration step is lower than a preset accuracy, and outputting optimal geological parameters and a prediction working face; The spatial coordinate system comprises H axis, X axis, Y axis, the ground surface deformation field measured data is specifically along H axis, X axis, Y axis movement; the geological parameters comprise: wherein, indicates the coal seam thickness, indicates the coal seam buried depth of the maximum subsidence point, indicates the spatial normal vector of the working face in X axis, Y axis component; the geological parameters of the initialized working face comprise: acquiring the maximum subsidence point coordinate : Initialization of the working plane The plane equation of the working plane is （1）； The deformation prediction value is obtained based on a geological parameter G, an inclination angle of the working face and a dip azimuth angle and a probability integral method. wherein the inclination angle is represented as: （2）； Tilt azimuth is represented as: （3）。 2. The method of claim 1, wherein: the step of constructing the prediction working face mask matrix based on the initialized working face comprises: dividing a plane where the initialized working face is located into grids, constructing a mask matrix Mask by taking the grids as raster cells, setting initial values of all raster cells to 0, setting a value of a raster cell where a maximum subsidence point in the Mask is located to 1 as a seed cell, and setting a size of an inflation and corrosion structure element to 3*3.

3. The method of claim 2, wherein: the step of performing inflation and corrosion on the mask matrix based on the image morphological algorithm specifically comprises: S1: performing an inflation operation on all raster cells with a value of 1 in the current Mask, taking a residual sum of squares of the measured data of the ground surface deformation field and the inflation and corrosion deformation prediction values as a fitting effect evaluation index, and performing fitting effect evaluation on new edge cells generated after the inflation one by one, if the evaluation index is reduced after adding the new edge cells, that is, the fitting effect is improved, the new raster cells generated after the inflation are effective, and the raster cells are retained; otherwise, the raster cells are discarded; S2: setting values of raster cells corresponding to all effective inflation cells in the Mask to 1, to realize inflation of the Mask; S3: iteratively performing steps S1 and S2 until all new raster cells generated after a certain inflation are invalid inflation cells; S4: performing a corrosion operation on all raster cells with a value of 1 in the current Mask, and performing fitting effect evaluation on edge new cells generated after the corrosion one by one, if the evaluation index is reduced after discarding the raster cells, that is, the fitting effect is improved, the raster cells generated after the corrosion are effective, that is, the raster cells are discarded; otherwise, the raster cells are retained; S5: setting values of raster cells corresponding to all effective corrosion cells in the Mask to 0, to realize corrosion shrinkage of the Mask; S6: iteratively performing steps S4 and S5 until all edge new raster cells generated after a certain corrosion are invalid corrosion cells, that is, no raster cell can be discarded through corrosion; S7: raster cells with a value of 1 are taken as the fitted mask matrix Mask and output.

4. The method of claim 1, wherein: the evaluation index is obtained by taking a length square sum of residual vectors of three-dimensional point movement of the ground surface as an evaluation index value: （4） In formula (4), are the measured displacements along H axis, X axis, Y axis, are the predicted displacements along H axis, X axis, Y axis, is the spatial movement vector, is the module of the vector difference.

5. The method of claim 1, wherein: in step S3, the process of updating the iteration step of the geological parameters G according to the preset rule is specifically: Setting initial step size and precision ; Based on initial step size Iterate through initial geologic parameter G, if evaluation metric is lower after iteration, update geologic parameter G as better geologic parameter G and increase step size , otherwise keep geologic parameter G and decrease step size .

6. An irregular working face positioning method of coupling ground surface deformation field and mining subsidence model according to claim 5, characterized in that: The increase step The increase coefficient is 2 times when the step is increased The decrease coefficient is 0.5 times when the step is decreased

Images