Method for dynamically predicting ground surface horizontal movement caused by underground coal mine mining based on lateral pressure coefficient

Abstract

The invention relates to the technical field of underground mining, and discloses a dynamic prediction method for ground surface horizontal movement caused by coal mine underground mining based on a lateral pressure coefficient, comprising the following steps: step S101, generating a lateral pressure coefficient time sequence curve and a goaf equivalent stiffness time sequence curve; step S102, obtaining a dynamic horizontal movement coefficient; s103, establishing a prediction model of the earth surface horizontal movement amount; step S104, positioning the horizontal movement curvature turnover window; step S105, generating a re-parameterization time axis; and S106, generating an earth surface horizontal movement amount distribution map. The method comprises the following steps: capturing horizontal movement curvature overturning characteristics generated by coupling a lateral pressure coefficient and goaf equivalent stiffness by dynamically quantifying time sequence evolution of the lateral pressure coefficient and the goaf equivalent stiffness; and the finally output dynamic horizontal movement amount distribution diagram can present a spatio-temporal evolution law, so that technical support is provided for engineering decisions such as building and structure protection and water body protection in a mining area.

Description

Technical Field

[0001] This invention relates to the field of underground mining technology, and more specifically, to a method for dynamically predicting horizontal surface movement caused by underground coal mining based on the lateral pressure coefficient. Background Technology

[0002] Underground coal mining causes horizontal surface movement, which directly affects the safety and stability of surrounding buildings, railways, and water bodies. Therefore, accurate dynamic prediction of horizontal surface movement is one of the key technologies for safe coal mining. Currently, the mainstream method in the industry is to use a probability integral method combined with Knothe-like time functions for prediction. This method calculates the amount of surface movement by constructing a coupled model of a spatial influence kernel and a time evolution function.

[0003] However, existing technologies have significant drawbacks, leading to systematic biases in prediction accuracy, particularly in determining peak times. On one hand, traditional methods treat the lateral pressure coefficient as a constant, neglecting the dynamic evolution of horizontal stress caused by vertical stress unloading during mining. The temporal changes in the lateral pressure coefficient directly amplify or weaken the surface's horizontal movement capacity, and fixed parameters cannot reflect this dynamic driving effect. On the other hand, the strength of the goaf filling increases with age, and the collapsed rock mass also compacts over time, resulting in a dynamic increase in the equivalent stiffness of the goaf. Traditional methods, using fixed stiffness parameters, struggle to reflect this dynamic inhibition of horizontal movement.

[0004] More importantly, the amplification effect of the lateral pressure coefficient and the suppression effect of the equivalent stiffness of the goaf are asynchronous on the time scale. The coupling between the two causes the second derivative of horizontal movement to exhibit a curvature reversal characteristic with a sign inversion. Traditional single time functions cannot adapt to this dynamic change, leading to system misjudgment of the peak movement time. These problems make it difficult for existing methods to accurately characterize the dynamic evolution of horizontal surface movement and meet the high-precision prediction requirements in engineering practice. Summary of the Invention

[0005] This invention provides a method for dynamically predicting horizontal surface movement caused by underground coal mining based on the lateral pressure coefficient, thereby solving the technical problems mentioned in the background.

[0006] This invention provides a method for dynamically predicting horizontal surface movement induced by underground coal mining based on lateral pressure coefficient, comprising the following steps:

[0007] Step S101: Collect the initial ground stress state and the parameters of the goaf filling medium, and generate the time series curve of the lateral pressure coefficient and the time series curve of the equivalent stiffness of the goaf, respectively.

[0008] Step S102: Using the time series curve of lateral pressure coefficient as an amplification factor and the time series curve of equivalent stiffness of goaf as an inhibition factor, the initial horizontal movement coefficient is time-varyingly corrected to obtain the dynamic horizontal movement coefficient.

[0009] Step S103: Couple the dynamic horizontal movement coefficient with the spatial direction influence kernel and the reference time function to establish a prediction model for the amount of horizontal movement of the land surface.

[0010] Step S104: Solve the second derivative of the prediction model with respect to time, identify the time period when the sign of the second derivative is reversed, and locate the horizontally moving curvature flip window in the time domain.

[0011] Step S105: Extract the second derivative components within the horizontally moving curvature flip window to construct time scaling weights, and use the time scaling weights to map the original time axis to generate a reparameterized time axis.

[0012] Step S106: Substitute the reparameterized time axis into the prediction model to replace the baseline time axis, calculate the horizontal movement vector of each point on the land surface under the reparameterized time, and synthesize and output the distribution map of the horizontal movement of the land surface.

[0013] The beneficial effects of this invention are as follows: By dynamically quantifying the temporal evolution of the lateral pressure coefficient and the equivalent stiffness of the goaf, this invention captures the horizontal movement curvature reversal characteristics generated by their coupling, overcoming the problem of insufficient dynamic evolution characterization caused by fixed parameters in traditional methods; relying on time reparameterization driven by curvature reversal window, it can adaptively correct time axis deviations without changing the calculation framework of the mainstream probability integral method and the reference time function, thereby improving the prediction accuracy of the peak time and amount of horizontal movement of the ground surface; the final output dynamic horizontal movement distribution map can present the spatiotemporal evolution law, thus providing technical support for engineering decisions such as the protection of buildings and structures in mining areas, railway safety operation, and water protection. Attached Figure Description

[0014] Figure 1 This is a flowchart of the method for dynamically predicting horizontal surface movement caused by underground coal mining based on the lateral pressure coefficient of the present invention.

[0015] Figure 2 This is a schematic diagram showing the distribution of the amount of horizontal movement of the ground surface along the working face at different time periods according to the present invention.

[0016] Figure 3 This is a schematic diagram showing the distribution of the amount of horizontal movement of the ground surface along the direction of the working face at different time periods according to the present invention. Detailed Implementation

[0017] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0018] It should be noted that, unless otherwise defined, the technical or scientific terms used in one or more embodiments of the present invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in one or more embodiments of the present invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" indicate that the element or object preceding the term encompasses the elements or objects listed following the term and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0019] like Figures 1-3 As shown, the method for dynamically predicting horizontal surface movement induced by underground coal mining based on lateral pressure coefficient includes the following steps:

[0020] Step S101: Collect the initial ground stress state and the parameters of the goaf filling medium, and generate the time series curve of the lateral pressure coefficient and the time series curve of the equivalent stiffness of the goaf, respectively.

[0021] Step S102: Using the time series curve of lateral pressure coefficient as an amplification factor and the time series curve of equivalent stiffness of goaf as an inhibition factor, the initial horizontal movement coefficient is time-varyingly corrected to obtain the dynamic horizontal movement coefficient.

[0022] Step S103: Couple the dynamic horizontal movement coefficient with the spatial direction influence kernel and the reference time function to establish a prediction model for the amount of horizontal movement of the land surface.

[0023] Step S104: Solve the second derivative of the prediction model with respect to time, identify the time period when the sign of the second derivative is reversed, and locate the horizontally moving curvature flip window in the time domain.

[0024] Step S105: Extract the second derivative components within the horizontally moving curvature flip window to construct time scaling weights, and use the time scaling weights to map the original time axis to generate a reparameterized time axis.

[0025] Step S106: Substitute the reparameterized time axis into the prediction model to replace the baseline time axis, calculate the horizontal movement vector of each point on the land surface under the reparameterized time, and synthesize and output the distribution map of the horizontal movement of the land surface.

[0026] In one embodiment of the present invention, the initial ground stress state and parameters of the goaf filling medium are collected to generate time-series curves of lateral pressure coefficient and equivalent stiffness of the goaf, respectively, including:

[0027] Based on the effective internal friction angle The initial lateral pressure coefficient is calculated according to the following formula. :

[0028] ;

[0029] Based on the average volumetric weight of the overlying strata With equivalent unloading thickness Calculate the vertical total stress increment according to the following formula. :

[0030] ;

[0031] Based on Poisson's ratio With the increase in total vertical stress Calculate the total horizontal stress increment according to the following formula. :

[0032] ;

[0033] Based on the initial horizontal total stress Initial vertical total stress Horizontal total stress increment and the increase in total vertical stress The time series curve of the lateral pressure coefficient is generated according to the following formula. :

[0034] ;

[0035] Based on the ultimate limit value of unconfined compressive strength With age-related intensity growth rate constant The unconfined compressive strength as a function of time is calculated using the following formula. :

[0036] ;

[0037] Based on conversion factor , With time-varying unconfined compressive strength The equivalent stiffness time series curve of the goaf is generated according to the following formula. :

[0038] .

[0039] It should be noted that the effective internal friction angle represents a key parameter of soil or dilatant rock mass's resistance to shear friction, reflecting the medium's ability to resist shear failure. The initial lateral pressure coefficient represents the static lateral pressure coefficient without mining disturbance. The average volumetric weight of the overburden represents the weight per unit volume of the overburden above the goaf, reflecting the influence of the overburden's own weight on the underlying stress state. The equivalent unloading thickness represents the equivalent thickness of the rock strata that actually experiences stress release during mining, reflecting the degree to which mining alters the rock strata's bearing capacity. The vertical total stress increment represents the change in vertical total stress caused by mining, reflecting the magnitude of the disturbance to the vertical stress state. Poisson's ratio represents the ratio of the medium's lateral deformation to its longitudinal deformation, reflecting the coupling characteristics of the medium's deformation. The horizontal total stress increment represents the change in horizontal total stress caused by mining. The initial horizontal total stress represents the initial value of the horizontal total stress at the observation point before mining. The instantaneous horizontal total stress represents the actual horizontal total stress value at the observation point at a certain moment. The initial vertical total stress represents the initial value of the vertical total stress at the observation point before mining. The instantaneous total vertical stress represents the actual total vertical stress value at a certain observation point at a given moment.

[0040] It should be noted that the lateral pressure coefficient time-series curve represents the curve of the lateral pressure coefficient changing over time, used to dynamically characterize the evolution of the ratio of horizontal stress to vertical stress during mining. The unconfined compressive strength limit value represents the maximum uniaxial compressive strength that the cemented backfill can withstand after complete solidification, reflecting the ultimate bearing capacity of the backfill after solidification. The age-dependent strength growth rate constant represents the rate parameter of strength increase of the cemented backfill with age, reflecting the speed of strength increase. The unconfined compressive strength represents the uniaxial compressive strength of the cemented backfill at a certain age, reflecting the actual bearing capacity of the backfill at that moment. The conversion factor represents the correlation parameter between the unconfined compressive strength and the equivalent stiffness of the goaf, used to convert the strength index of the backfill into a stiffness index. The equivalent stiffness of the goaf represents the comprehensive stiffness of the combination of the backfill and loose medium within the goaf, reflecting the goaf's ability to resist deformation. The time series curve of the equivalent stiffness of the goaf represents the curve of the equivalent stiffness of the goaf changing over time, and is used to dynamically characterize the evolution law of the goaf's resistance to deformation during mining.

[0041] It should be noted that the equivalent unloading thickness represents the equivalent thickness of the rock strata due to the loss of support in the overlying strata during mining, and is not the actual thickness of the mined coal seam. Its determination requires consideration of the mining method, overlying lithology, and strata movement patterns. For example, in longwall mining, it can be calculated comprehensively using the height of the overlying caving zone, the height of the fracture zone, and the coal seam thickness. In a certain mining area, the coal seam thickness is 3 meters, the overlying caving zone height is 15 meters, and the fracture zone height is 30 meters. Through engineering analogy and numerical simulation, the equivalent unloading thickness is determined to be 8 meters. This means that the weight of the rock strata of this thickness is no longer borne by the original coal seam due to mining, thus causing a change in stress increment. The age-dependent strength growth rate constant is obtained through laboratory tests. Test blocks with the same mix proportions as those used in the field are prepared, and the unconfined compressive strength is measured at different ages (1 day, 3 days, 7 days, 14 days, and 28 days). The test data are then fitted into the unconfined compressive strength calculation formula to obtain this constant. After testing and fitting, the strength growth rate constant of a certain cemented backfill specimen was 0.05 per day, indicating a relatively slow strength growth rate. Conversion factors, namely the constant term conversion factor and the linear term conversion factor, are used to establish a linear relationship between unconfined compressive strength and the equivalent stiffness of the goaf. This requires indoor mechanical testing of backfill specimens of different strength grades to perform compressive and stiffness tests, obtaining multiple sets of corresponding data, and then using linear regression analysis. For example, through regression calculation using 10 sets of test data, the first conversion factor is 500 MPa, and the second conversion factor is 100, meaning that for every 1 MPa increase in unconfined compressive strength, the equivalent stiffness of the goaf increases by 100 MPa, while maintaining a basic stiffness of 500 MPa.

[0042] It should be noted that by collecting samples of the overburden and surrounding soil above the goaf, direct shear tests or triaxial shear tests are conducted indoors. By applying different levels of normal stress, the shear strength of the samples at the corresponding levels is determined. The effective internal friction angle is obtained by fitting the test data according to the Mohr-Coulomb strength theory. The average volumetric weight of the overburden can be obtained in two ways: one is to select rock samples from each major layer of the overburden, dry them indoors, measure their density, calculate the weighted average density based on the thickness ratio of each layer, and then convert it into volumetric weight; the other is to consult the geological survey report of the mining area, extract the density data of each layer, and calculate the weighted average based on thickness.

[0043] It should be noted that uniaxial compression tests can be conducted indoors, using strain gauges to simultaneously measure the longitudinal and transverse strains of the rock sample under axial pressure, and calculating Poisson's ratio based on the ratio of transverse to longitudinal strain. Alternatively, in-situ acoustic testing technology can be used, measuring the longitudinal and transverse wave velocities of the rock mass and converting them to Poisson's ratio based on the rock mass density. In-situ stress testing methods can be used to collect data, commonly employing hydraulic fracturing or stress relief methods. Hydraulic fracturing involves applying water pressure within the borehole to fracture the borehole wall, and calculating the initial total horizontal stress based on parameters such as fracturing pressure and shut-off pressure. Stress relief involves drilling core samples and measuring the deformation after stress release, then inferring the initial total horizontal stress. Data collection can also be calculated based on the average volumetric weight of the overburden and the burial depth of the observation point. Specifically, the vertical distance from the observation point to the surface is determined through borehole exploration in the mining area, and the product of the previously collected average volumetric weight of the overburden is used to obtain the initial total vertical stress. According to the actual mix proportion of the cemented infill material on site, standard test blocks are prepared indoors and cured in the same curing environment as on site for a long period of time. When the strength of the test blocks no longer increases significantly (usually after 90 days or more of curing), an unconfined compressive strength test is conducted, and the test result is the unconfined compressive strength limit value.

[0044] It should be noted that this invention determines the baseline of the initial lateral pressure coefficient by using the effective internal friction angle, quantifies the vertical stress increment caused by mining by combining the average volumetric weight of the overburden and the equivalent unloading thickness, obtains the horizontal stress increment by using Poisson's ratio coupling, and generates a time-series curve of the lateral pressure coefficient, dynamically reflecting the influence of mining on the ratio of horizontal to vertical stress. At the same time, relying on the characteristic that the strength of cemented backfill increases with age, the time-dependent strength is calculated by using the unconfined compressive strength limit value and the age-dependent strength growth rate constant, and converted into a time-series curve of the equivalent stiffness of the goaf by conversion coefficients, realizing the dynamic quantification of intrinsic mechanical properties, and providing accurate time-series parameter input for subsequent calculation of dynamic horizontal movement coefficient.

[0045] In one embodiment of the present invention, the time-series curve of the lateral pressure coefficient is used as an amplification factor, and the time-series curve of the equivalent stiffness of the goaf is used as a suppression factor to perform time-varying correction on the initial horizontal movement coefficient to obtain the dynamic horizontal movement coefficient, including:

[0046] Based on the time series curve of the lateral pressure coefficient Initial lateral pressure coefficient and lateral pressure sensitivity coefficient Construct the magnification factor according to the following formula :

[0047] ;

[0048] Based on the time series curve of equivalent stiffness of goaf Stiffness sensitivity coefficient and reference modulus Construct the repressor factor according to the following formula :

[0049] ;

[0050] Based on the initial horizontal movement coefficient Amplification factor and inhibitory factors Calculate the dynamic horizontal movement coefficient according to the following formula. :

[0051] .

[0052] It should be noted that the lateral pressure sensitivity coefficient represents the degree of response of the horizontal movement coefficient to changes in the lateral pressure coefficient, reflecting the sensitivity intensity of the influence of lateral pressure evolution on horizontal movement. The unit reference value represents the reference constant for calculating the amplification factor and the suppression factor. The amplification factor represents the enhancement effect coefficient of the temporal change of the lateral pressure coefficient on horizontal movement. The reference modulus represents the reference modulus used for dimensionless goaf equivalent stiffness, reflecting the scaling standard for calculating stiffness suppression effects. The stiffness sensitivity coefficient represents the degree of response of the horizontal movement coefficient to changes in the goaf equivalent stiffness, reflecting the sensitivity intensity of the influence of goaf stiffness evolution on horizontal movement. The suppression factor represents the weakening effect coefficient of the temporal change of the goaf equivalent stiffness on horizontal movement, reflecting the degree of constraint of the evolution of goaf mechanical properties on horizontal movement. The initial horizontal movement coefficient represents the horizontal movement coefficient under no mining or baseline conditions, reflecting the basic strength of horizontal movement in the initial state. The dynamic horizontal movement coefficient represents the horizontal movement coefficient that changes over time, reflecting the dynamic evolution of horizontal movement capability under the coupling effect of lateral pressure and goaf stiffness.

[0053] It should be noted that the lateral pressure sensitivity coefficient can be determined by fitting measured data from the mining area. Surface horizontal movement observation points under the same or similar geological and mining conditions in the mining area are selected. Measured horizontal movement data and corresponding lateral pressure coefficient time-series data are obtained from these observation points. Combined with other parameters of the calculated amplification factor, these parameters are substituted into the correlation model between the amplification factor and horizontal movement to inversely deduce the lateral pressure sensitivity coefficient. For example, in a certain mining area, the lateral pressure sensitivity coefficient was inversely derived from 5 sets of observation data to be 2.5, indicating that for every 0.1 increase in the lateral pressure coefficient beyond the initial value, the amplification factor increases by 0.25, and the horizontal movement capability is correspondingly enhanced. The unit reference value is a fixed constant, set to 1. This reference value ensures that the amplification factor is 1 when the lateral pressure coefficient equals the initial value, at which point there is no additional amplification effect on the lateral pressure; and ensures that the inhibition factor is 1 when the equivalent stiffness of the goaf is 0, at which point there is no additional inhibition effect on the stiffness.

[0054] It should be noted that the reference modulus can be selected as a value of the same order of magnitude as the modulus of the surrounding strata or filling material in the goaf. This can be obtained by consulting the elastic modulus data of the strata and rock masses in the geological survey report of the mining area, or by referring to the design modulus of the filling material on site, and taking the average value as the reference modulus. For example, if the elastic modulus of the surrounding strata in a certain mining area is 2000 MPa to 3000 MPa, a reference modulus of 2500 MPa can be selected to ensure that the ratio of the equivalent stiffness of the goaf to the reference modulus is within a reasonable range, avoiding extreme values ​​in the calculation of the inhibition factor. The stiffness sensitivity coefficient can be determined by combining indoor tests and field observations. Filling material test blocks of different stiffness levels are prepared, and indoor horizontal displacement simulation tests are conducted to obtain the correlation data between stiffness and horizontal displacement. Simultaneously, combined with the horizontal movement observation data of the mining area surface, the two types of data are coupled and analyzed to obtain the stiffness sensitivity coefficient. For example, a stiffness sensitivity coefficient of 3.0, calculated through fitting, indicates that for every 0.1 increase in the ratio of the equivalent stiffness of the goaf to the reference modulus, the inhibition factor decreases by approximately 0.028, and the constraint on horizontal movement capability is strengthened. The initial horizontal movement coefficient can be calibrated using historical measured data from the mining area. Surface horizontal movement observation data from periods without mining disturbance or in the early stages of mining are selected, and combined with the geological and mining parameters at that time, the coefficient is determined using regression analysis. For example, initial stage data from 10 historical observation points in a certain mining area are selected, and regression calculation yields an initial horizontal movement coefficient of 0.005, meaning that in the initial state, the horizontal movement corresponding to the contribution of a unit geometric influence kernel and a unit time function is 0.005 meters.

[0055] In one embodiment of the present invention, a prediction model for surface horizontal movement is established by coupling the dynamic horizontal movement coefficient with the spatial direction influence kernel and the reference time function, including:

[0056] Based on the coordinates of the observation point Coordinates of the center of the sampled unit Calculate the equivalent distance using the following formula. :

[0057] ;

[0058] Based on the radius of influence and equivalent distance Construct an exponential influence function according to the following formula. :

[0059] ;

[0060] Projection of planar coordinate differences in the orientation direction Radius of influence and exponential influence function The direction of influence core is generated according to the following formula. :

[0061] ;

[0062] Based on the projection of the plane coordinate difference in the dip direction Radius of influence and exponential influence function Generate the tendency direction influence core according to the following formula. :

[0063] ;

[0064] Based on time variables Time-related parameters and shape parameters Generate the reference time function according to the following formula. :

[0065] ;

[0066] Based on dynamic horizontal movement coefficient Direction affects the nuclear and reference time function The horizontal movement of the ground surface in the direction of the strike is calculated using the following formula. :

[0067] ;

[0068] Based on dynamic horizontal movement coefficient The direction of tendency affects the nucleus and reference time function The horizontal movement of the ground surface in the dip direction is calculated according to the following formula. :

[0069] .

[0070] It should be noted that the observation point represents a specific point on the surface used to monitor horizontal movement. The center of the mined unit represents the geometric center of a single mining unit after discretizing the goaf. The equivalent distance represents the straight-line distance from the observation point to the center of the mined unit. The radius of influence represents a key geometric parameter controlling the range of surface movement influence, reflecting the spatial attenuation scale of mining impact. The exponential influence function describes the intensity of the influence of the mined unit on the movement of the observation point, used to quantify the attenuation law of spatial distance contribution to movement. The plane coordinate difference represents the coordinate difference between the observation point and the center of the mined unit in the plane coordinate system, reflecting the plane position offset between the two. The strike direction represents the extension direction of the coal seam mining face, reflecting the main extension dimension of mining impact. The dip direction represents the direction of coal seam dip, reflecting the secondary extension dimension of mining impact. The projection represents the component of the plane coordinate difference in the strike or dip direction.

[0071] It should be noted that the strike direction influence kernel represents the kernel function that quantifies the contribution of the sampled unit to the horizontal movement of the observation point in the strike direction. The dip direction influence kernel represents the kernel function that quantifies the contribution of the sampled unit to the horizontal movement of the observation point in the dip direction. The time variable represents the variable characterizing the sequence of the mining process. The time influence parameter represents the parameter controlling the convergence rate of the reference time function. The shape parameter represents the parameter controlling the shape of the reference time function. The negative exponential result represents the result of the negative exponentiation of the shape parameter power of the product of the time variable and the time influence parameter. The reference time function is used to describe the evolution of surface horizontal movement over time. The strike direction surface horizontal movement represents the magnitude of the horizontal movement of the observation point in the strike direction. The dip direction surface horizontal movement represents the magnitude of the horizontal movement of the observation point in the dip direction.

[0072] It should be noted that the radius of influence needs to be determined by inversion of measured data or engineering analogy, taking into account the geological and mining conditions of the mining area. A working face within the mining area that has been mined and has complete surface movement observation data is selected. Measured horizontal movement data of observation points and corresponding spatial location data are collected, and substituted into the correlation model of exponential influence function and directional influence kernel to inversely deduce the radius of influence. For example, in a certain mining area, a working face with a mining depth of 500 meters and overlying rock mainly composed of medium-hard sandstone, the radius of influence is found to be 150 meters through inversion of three sets of observation profile data. This indicates that the mining impact will cause significant horizontal movement of observation points within this radius, and the impact decays rapidly beyond this radius. The projection of the plane coordinate difference can be calculated by first establishing a plane coordinate system, with a corner point of the working face as the origin, the strike direction as the X-axis, and the dip direction as the Y-axis. The X-axis difference and Y-axis difference between the coordinates of the observation point and the center coordinates of the mined unit are calculated. The X-axis difference is the projection of the plane coordinate difference in the strike direction, and the Y-axis difference is the projection of the plane coordinate difference in the dip direction. If the observation point is in the positive X-axis direction of the center of the sampled unit, the projection in the strike direction is positive; if it is in the negative X-axis direction, it is negative. The same applies to determining the sign of the projection in the dip direction. For example, if the coordinates of the observation point are X100 meters and Y50 meters, and the coordinates of the center of the sampled unit are X80 meters and Y30 meters, the projection of the difference in plane coordinates in the strike direction is 20 meters, and the projection in the dip direction is also 20 meters.

[0073] In one embodiment of the present invention, solving the second derivative of the prediction model with respect to time, identifying the time period in which the sign of the second derivative reverses, and locating a horizontally moving curvature flip window in the time domain includes:

[0074] For the dynamic horizontal shift coefficient respectively With reference time function Taking the first and second derivatives with respect to time, we get , , as well as ;

[0075] The direction of the impact on the nuclear The direction of tendency affects the nucleus And the derivative result, calculate the second derivative in the direction of the direction according to the following formula. Second derivative with tendency direction :

[0076] ;

[0077] ;

[0078] Calculate the second derivative of the direction Construct a set of zero points at the moments when the time is equal to zero. :

[0079] ;

[0080] judge any point in time The product of the second derivatives of the left and right neighborhoods is determined when the following condition is satisfied. For the sign reversal moment:

[0081] ;

[0082] Selection and sign reversal time Adjacent preceding zero time points With subsequent zero point time :

[0083] ;

[0084] ;

[0085] Time period The horizontal movement curvature flip window is determined by the direction of travel;

[0086] Second derivative with respect to the direction of inclination Perform the same operation to determine the horizontal movement curvature flip window in the direction of inclination. Furthermore, the horizontal movement curvature flip windows in the direction of directional movement and the direction of inclination are merged into a time-domain window set. :

[0087] .

[0088] It should be noted that the first derivative represents the rate of change of the dynamic horizontal movement coefficient or the reference time function over time. The second derivative represents the rate of change of the first derivative of the dynamic horizontal movement coefficient or the reference time function over time. The first term represents the coupling term of the second derivative of the dynamic horizontal movement coefficient and the reference time function, reflecting the contribution of the second-order evolution of the dynamic coefficient to the second derivative of the horizontal movement. The second term represents the coupling amplification term of the first derivative of the dynamic horizontal movement coefficient and the first derivative of the reference time function, reflecting the combined contribution of their first-order evolution to the second derivative of the horizontal movement. The third term represents the coupling term of the dynamic horizontal movement coefficient and the second derivative of the reference time function, reflecting the contribution of the second-order evolution of the reference time function to the second derivative of the horizontal movement. The second derivative in the direction of strike represents the second rate of change of the horizontal movement of the surface over time in the direction of strike. The second derivative in the direction of dip represents the second rate of change of the horizontal movement of the surface over time in the direction of dip. The product of the second derivatives in the left and right neighborhoods represents the product of the second derivative values ​​in the adjacent time intervals to the left and right of the time point. The sign reversal moment represents the time point when the second derivative changes from positive to negative or vice versa, reflecting the critical time when the acceleration / deceleration state of movement changes. The preceding zero point moment represents the closest zero point of the second derivative before the sign reversal moment, i.e., the boundary time to the left of the sign reversal moment. The following zero point moment represents the closest zero point of the second derivative after the sign reversal moment, i.e., the boundary time to the right of the sign reversal moment. The horizontal movement curvature flip window in the directional direction represents the time interval during which the acceleration / deceleration state of movement changes in the directional direction. The horizontal movement curvature flip window in the dip direction represents the time interval during which the acceleration / deceleration state of movement changes in the dip direction. The time-domain window set represents the set of directional and dip curvature flip windows, reflecting all the critical periods of the second-order evolution of horizontal surface movement.

[0089] It should be noted that the left and right neighborhoods are small time intervals centered on the time point. The interval length needs to be reasonably set according to the mining time scale, usually between 0.01 days and 0.1 days. For example, if the time point is 10 days, the left neighborhood is selected as 9.99 days to 10 days, and the right neighborhood is selected as 10 days to 10.01 days. The average value of the second derivative within each interval is calculated, and then the two average values ​​are multiplied to determine the sign relationship. When there are multiple zeros in the set of zeros of the second derivative corresponding to the sign reversal time, the preceding zero time is selected as the largest among all zeros less than the sign reversal time, and the subsequent zero time is selected as the smallest among all zeros greater than the sign reversal time. For example, if the set of zeros of the second derivative in the direction of strike is 2 days, 5 days, 8 days, and 12 days, and the sign reversal time is 8 days, then the preceding zero time is 5 days, and the subsequent zero time is 12 days. The corresponding strike direction curvature flip window is 5 days to 12 days.

[0090] In one embodiment of the present invention, the second derivative components within the horizontally shifting curvature flip window are extracted to construct time-scaling weights. These time-scaling weights are then used to map the original time axis to generate a reparameterized time axis, including:

[0091] Based on dynamic horizontal movement coefficient first derivative With the second derivative and reference time function first derivative With the second derivative ;

[0092] Based on the time-strength coefficient With time domain window set Construct time-scaling weights according to the following formula :

[0093] ;

[0094] Among them, indicator function In time Belongs to the time-domain window set The value is 1 if the condition is met, and 0 otherwise.

[0095] Based on the integral independent variable Generate the reparameterized time axis according to the following formula. :

[0096] .

[0097] It should be noted that the first intermediate term represents the coupling term between the second derivative of the dynamic horizontal shift coefficient and the reference time function. The second intermediate term represents the coupling amplification term between the first derivative of the dynamic horizontal shift coefficient and the first derivative of the reference time function. The third intermediate term represents the coupling term between the dynamic horizontal shift coefficient and the second derivative of the reference time function. The numerator represents the superposition result of the first and second intermediate terms. The denominator represents the sum of the absolute values ​​of the three intermediate terms, reflecting the total amplitude of the curvature-related terms. The time stretching intensity coefficient represents the parameter controlling the stretching amplitude of the time axis. The time stretching weight represents the coefficient adjusting the evolution rate of the original time axis. Time zero represents the start of sampling or the starting node of time calculation.

[0098] It should be noted that the time stretching intensity coefficient can be determined by comparing and inverting the measured data and prediction results in the mining area. A region within the mining area with complete horizontal surface movement observation data is selected, and prediction calculations are performed using different values ​​of the time stretching intensity coefficient. The predicted peak time and movement amount of horizontal movement are compared with the measured data, and the value with the smallest prediction error is selected as the final value of the coefficient. For example, in a certain mining area, through multiple sets of tests, when the time stretching intensity coefficient is 0.8, the deviation between the predicted peak time and the measured time is only 1 day, and the movement amount deviation is less than 5 millimeters. This value is a suitable time stretching intensity coefficient. The time domain window set consists of multiple continuous time intervals. It is determined whether the current time belongs to any one of these intervals. If the value of the current time is greater than or equal to the start time of an interval and less than or equal to the end time of that interval, then the current time is considered to be within the time domain window set. If the current time is less than the start time of all intervals, or greater than the end time of all intervals, or falls between two intervals, then the current time is considered not to be within the time domain window set. For example, a time-domain window set includes intervals of 5 to 12 days and 18 to 25 days. If the current time is 8 days, it belongs to this set; if the current time is 15 days, it does not belong to this set.

[0099] In one embodiment of the present invention, the reparameterized time axis is substituted into the prediction model to replace the baseline time axis, the horizontal movement vector of each point on the land surface under the reparameterized time is calculated, and a distribution map of the horizontal movement of the land surface is synthesized and output, including:

[0100] Based on reparameterized time axis With reference time function The reparameterized time term is calculated using the following formula: ;

[0101] Based on dynamic horizontal movement coefficient Direction affects the nuclear and reparameterized time term The reparameterized directional shift component is calculated according to the following formula. :

[0102] ;

[0103] Based on dynamic horizontal movement coefficient The direction of tendency affects the nucleus and reparameterized time term The reparameterization tendency shift component is calculated according to the following formula. :

[0104] ;

[0105] Based on the discrete grid set of the ground surface Time sampling set Reparameterization towards the shift component and reparameterization tendency to shift components The horizontal movement vector of the earth's surface is generated according to the following formula. :

[0106] ;

[0107] Based on reparameterization, the direction of the shift component With reparameterization tendency to shift components The horizontal movement vector model of the earth's surface is calculated according to the following formula. :

[0108] ;

[0109] Based on the tendency to shift components by reparameterization With reparameterization towards the moving component Calculate the azimuth angle of horizontal movement of the Earth's surface using the following formula. :

[0110] .

[0111] It should be noted that the reparameterized time term represents the result of substituting the reparameterized time axis into the reference time function. The reparameterized directional movement component represents the horizontal movement in the directional direction after time correction. The reparameterized dip movement component represents the horizontal movement in the dip direction after time correction. The discrete grid set of the surface represents the set of regular points that divide the surface of the prediction area. The time sampling set represents the set of key time nodes for dividing the prediction period. The surface horizontal movement vector magnitude represents the magnitude of the horizontal movement vector, reflecting the actual intensity of the horizontal movement of the observation point. The binary arctangent operation represents the arctangent calculation method with two directional movement components as variables, used to determine the azimuth angle of the horizontal movement. The surface horizontal movement azimuth angle represents the direction angle of the horizontal movement, reflecting the specific orientation of the horizontal movement of the observation point. The surface horizontal movement distribution map represents a visual data carrier integrating the movement magnitude and azimuth angle of each grid point at different times, reflecting the spatiotemporal evolution of surface horizontal movement.

[0112] It should be noted that when dividing the surface discrete grid set, the boundary coordinates of the prediction area are determined first, and then the grid is divided according to the equal spacing rule, with the grid spacing usually ranging from 5 meters to 20 meters. For example, if the prediction area is 1000 meters long from east to west and 800 meters wide from north to south, and a grid spacing of 10 meters is selected, the surface discrete grid set will contain 101 columns and 81 rows, totaling 8181 grid points. The coordinates of each grid point are determined by accumulating the grid spacing based on the boundary starting value. The time sampling set needs to cover the initial mining stage, the main mining stage, and the stable stage, with sampling intervals being denser in the early stage and sparser in the later stage. In the initial mining stage (0 days to 30 days), sampling can be done at 1-day intervals; in the main mining stage (31 days to 180 days), sampling can be done at 3-day intervals; and in the stable stage (181 days to 360 days), sampling can be done at 7-day intervals to ensure complete capture of the entire movement and evolution process. The azimuth angle value range of the binary arctangent operation is 0 radians to 2 radians, corresponding to 0 degrees to 360 degrees. When both the reparameterized strike-shift component and the reparameterized dip-shift component are positive, the azimuth is between 0 and 90 degrees. When the strike component is negative and the dip component is positive, the azimuth is between 90 and 180 degrees. When both the strike and dip components are negative, the azimuth is between 180 and 270 degrees. When both the strike and dip components are positive, the azimuth is between 270 and 360 degrees. For example, if the reparameterized strike-shift component is 2 meters and the dip component is 2 meters, the azimuth is 45 degrees; if the strike component is -2 meters and the dip component is 2 meters, the azimuth is 135 degrees.

[0113] It should be noted that by substituting the reparameterized time axis into the reference time function, a time term adapted to the corrected time scale is generated. Combined with the dynamic horizontal movement coefficient and the directional influence kernel, the directional corrected movement components are obtained. By discretizing the surface grid and time sampling, the movement components of the entire region and all time periods are extracted. Based on the principle of vector synthesis, the movement modulus and azimuth are calculated and integrated to form distribution map data. The abstract movement parameters are transformed into intuitive visualization results, providing a clear and accurate spatiotemporal distribution basis for surface movement risk assessment and decision-making in engineering practice. This will not be elaborated further here.

[0114] It should be noted that the interval and threshold sizes are set for ease of comparison. The size of the threshold depends on the amount of sample data and the base number set by those skilled in the art for each set of sample data, as long as it does not affect the proportional relationship between the parameter and the quantized value. Furthermore, the above formulas are all dimensionless calculations, and the formulas are derived from software simulations using a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0115] The embodiments of this example have been described above. However, this example is not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms based on the guidance of this example, and all of them are within the protection scope of this example.

Claims

1. A method for dynamically predicting horizontal surface movement induced by underground coal mining based on lateral pressure coefficient, characterized in that, Includes the following steps: Step S101: Collect the initial ground stress state and the parameters of the goaf filling medium, and generate the time series curve of the lateral pressure coefficient and the time series curve of the equivalent stiffness of the goaf, respectively. Step S102: Using the time series curve of lateral pressure coefficient as an amplification factor and the time series curve of equivalent stiffness of goaf as an inhibition factor, the initial horizontal movement coefficient is time-varyingly corrected to obtain the dynamic horizontal movement coefficient. Step S103: Couple the dynamic horizontal movement coefficient with the spatial direction influence kernel and the reference time function to establish a prediction model for the amount of horizontal movement of the land surface. Step S104: Solve the second derivative of the prediction model with respect to time, identify the time period when the sign of the second derivative is reversed, and locate the horizontally moving curvature flip window in the time domain. Step S105: Extract the second derivative components within the horizontally moving curvature flip window to construct time scaling weights, and use the time scaling weights to map the original time axis to generate a reparameterized time axis. Step S106: Substitute the reparameterized time axis into the prediction model to replace the baseline time axis, calculate the horizontal movement vector of each point on the land surface under the reparameterized time, and synthesize and output the distribution map of the horizontal movement of the land surface.

2. The method for dynamically predicting horizontal surface movement caused by underground coal mining based on lateral pressure coefficient according to claim 1, characterized in that, Calculate the initial lateral pressure coefficient based on the effective internal friction angle; The vertical total stress increment is calculated based on the average volumetric weight of the overburden and the equivalent unloading thickness, and the horizontal total stress increment is calculated based on Poisson's ratio and the vertical total stress increment. The instantaneous horizontal total stress is obtained by summing the initial horizontal total stress and the increment of the horizontal total stress, and the instantaneous vertical total stress is obtained by summing the initial vertical total stress and the increment of the vertical total stress. The ratio of the instantaneous horizontal total stress to the instantaneous vertical total stress is calculated to generate the time series curve of the lateral pressure coefficient. The unconfined compressive strength as a function of time is calculated based on the unconfined compressive strength limit value and the strength growth rate constant over age; the equivalent stiffness of the goaf is calculated based on the conversion factor and the unconfined compressive strength, and a time series curve of the equivalent stiffness of the goaf is generated.

3. The method for dynamically predicting horizontal surface movement caused by underground coal mining based on lateral pressure coefficient according to claim 1, characterized in that, The difference between the time-series curve of the lateral pressure coefficient and the initial lateral pressure coefficient is calculated. The difference is multiplied by the lateral pressure sensitivity coefficient and then summed with the unit reference value to generate the amplification factor. The ratio of the time-series curve of the equivalent stiffness of the goaf to the reference modulus is calculated. The ratio is multiplied by the stiffness sensitivity coefficient and then summed with the unit reference value. The reciprocal of the summation result is taken to generate the suppression factor. The initial horizontal movement coefficient, the amplification factor, and the suppression factor are multiplied together to obtain the dynamic horizontal movement coefficient.

4. The method for dynamically predicting horizontal surface movement caused by underground coal mining based on lateral pressure coefficient according to claim 1, characterized in that, Calculate the equivalent distance between the observation point and the center of the sampled unit, and construct an exponential influence function based on the influence radius and the equivalent distance; Calculate the projection of the plane coordinate difference in the strike direction and dip direction, calculate the ratio of the projection to the square of the influence radius and take the negative sign, multiply the ratio by the influence function and sum them up for the sampled cells to generate the strike direction influence kernel and dip direction influence kernel respectively.

5. The method for dynamically predicting horizontal surface movement caused by underground coal mining based on lateral pressure coefficient according to claim 4, characterized in that, Calculate the product of the time variable and the time-affecting parameter, calculate the shape parameter power of the product, take the negative exponent of the power result, calculate the difference between the unit reference value and the negative exponent result, and generate the reference time function. The dynamic horizontal movement coefficient, the strike direction influence kernel, and the reference time function are multiplied together to obtain the strike direction surface horizontal movement. The horizontal movement coefficient, the dip direction influence kernel, and the reference time function are multiplied together to obtain the dip direction surface horizontal movement.

6. The method for dynamically predicting horizontal surface movement caused by underground coal mining based on lateral pressure coefficient according to claim 1, characterized in that, The first and second derivatives with respect to time are obtained for the dynamic horizontal movement coefficient and the reference time function, respectively. The first term is obtained by multiplying the second derivative of the dynamic horizontal shift coefficient with the reference time function. The second term is obtained by multiplying the first derivative of the dynamic horizontal shift coefficient with the first derivative of the reference time function and doubling the result. The third term is obtained by multiplying the dynamic horizontal shift coefficient with the second derivative of the reference time function. The first, second, and third terms are summed and multiplied by the directional influence kernel and the dip direction influence kernel, respectively, to generate the directional second derivative and the dip direction second derivative.

7. The method for dynamically predicting horizontal surface movement caused by underground coal mining based on lateral pressure coefficient according to claim 6, characterized in that, Calculate the moment when the second derivative of the direction of travel is equal to zero, and determine whether the product of the second derivatives in the left and right neighborhoods of the moment is less than zero. If the product is less than zero, determine that the moment is the moment of sign reversal. In the time series of sign reversal moments, the preceding zero point and the following zero point adjacent to the current sign reversal moment are selected, and the time period between the preceding zero point and the following zero point is determined as the horizontal moving curvature flip window of the direction. Perform the same operation on the second derivative of the yaw direction to determine the horizontal movement curvature flip window of the yaw direction, and merge the horizontal movement curvature flip windows of the yaw direction and the yaw direction into a time-domain window set.

8. The method for dynamically predicting horizontal surface movement caused by underground coal mining based on lateral pressure coefficient according to claim 1, characterized in that, Calculate the first derivative of the dynamic horizontal shift coefficient, the second derivative of the dynamic horizontal shift coefficient, the first derivative of the reference time function, and the second derivative of the reference time function, respectively. Calculate the product of the second derivative of the dynamic horizontal shift coefficient and the reference time function to generate the first intermediate term; Calculate the product of the first derivative of the dynamic horizontal shift coefficient and the first derivative of the reference time function, and double the product to generate the second intermediate term; Calculate the product of the dynamic horizontal shift coefficient and the second derivative of the reference time function to generate the third intermediate term; Summing the first intermediate term and the second intermediate term yields the numerator. Calculate the absolute values ​​of the first, second, and third intermediate terms, sum the three absolute values, and use the sum as the denominator.

9. The method for dynamically predicting horizontal surface movement caused by underground coal mining based on lateral pressure coefficient according to claim 8, characterized in that, Calculate the ratio of the numerator to the denominator to determine whether the current time is within the time-domain window set; When the current moment is within the time domain window set, calculate the sum of the products of the unit reference value, the time scaling strength coefficient, and the ratio to generate the time scaling weight; When the current time is not within the time domain window set, the unit baseline value is used as the time scaling weight; The time scaling weights are calculated by definite integral over the interval from time zero to the current time to generate a reparameterized time axis.

10. The method for dynamically predicting horizontal surface movement caused by underground coal mining based on lateral pressure coefficient according to claim 1, characterized in that, Substitute the reparameterized time axis as the independent variable into the baseline time function to generate the reparameterized time term; Multiply the dynamic horizontal shift coefficient, the direction influence kernel, and the reparameterized time term together to generate the reparameterized direction shift component; Multiply the dynamic horizontal shift coefficient, the directional influence kernel, and the reparameterized time term together to generate the reparameterized directional shift component; Traverse the discrete grid set and time sampling set of the ground surface to extract the reparameterized directional movement component and the reparameterized trend movement component at the current grid point and the current time. Calculate the squares of the reparameterized directional movement component and the reparameterized dip movement component, sum the two squared results and perform a square root operation to obtain the surface horizontal movement vector model. Using the reparameterized tidal movement component as the first variable and the reparameterized directional movement component as the second variable, perform a binary arctangent operation to obtain the horizontal movement azimuth of the ground surface. By combining the horizontal movement vector of the Earth's surface with the azimuth of the horizontal movement of the Earth's surface, a distribution map of the amount of horizontal movement of the Earth's surface is generated.

Images