Improved unconstrained nonlinear programming probability integral parameter test method

Abstract

The invention relates to the technical field of coal mining, in particular to an improved unconstrained nonlinear programming probability integral parameter test method, which comprises the following steps of: inputting working face position parameters, actual measurement data and geological mining parameters; according to the method, robust estimation is applied to unconstrained nonlinear programming inversion probability integral parameters, and a weight adjustment mechanism is introduced, so that the sensitivity to abnormal values is reduced. The method is beneficial to weakening adverse effects of abnormal values on unknown number estimation to a certain extent, so that the estimation stability is improved, the method has the advantages of being high in parameter obtaining precision and good in curve fitting effect, meanwhile, the method has good parameter estimation performance and stability and strong anti-interference capability, and the method has good application prospects to a certain extent. According to the method, data of incomplete and nonstandard observation stations can be effectively utilized, the accuracy and reliability of parameter inversion can be ensured, incomplete data can be fully utilized, and resource waste is avoided.

Description

Technical Field

[0001] The invention relates to the technical field of coal mining, in particular to an improved probability integral parameter testing method for unconstrained nonlinear programming. Background Art

[0002] The disasters associated with coal mining have brought severe challenges to the ecological and social stability of mining areas. Therefore, it is particularly important to study the laws and characteristics of surface deformation caused by coal mining and to accurately predict mining subsidence. The most commonly used prediction method is the probability integral method. By establishing a probability integral nonlinear model, the formula is derived and then the surface deformation law is calculated and analyzed. The formula contains parameters that need to be inverted in advance. Only after they are obtained can the specific formula be obtained.

[0003] The existing parameter inversion methods include characteristic point method, linear approximation method, orthogonal experiment method, intelligent optimization algorithm, nonlinear programming method, etc. The linear approximation method has the disadvantages of high requirements for the layout of working surface and measuring station, and sensitivity to initial solution; the orthogonal experiment design method selects representative experiments from a large number of experiments to determine the optimal solution, which has the disadvantages of large workload and low efficiency, and is difficult to promote universally; when the optimization algorithm is used for different initial solutions, the parameters estimated by the algorithm are quite different, resulting in unsatisfactory fitting effect; in recent years, the intelligent optimization algorithm has become a widely used research method in the field of probability integral parameter inversion by scholars. Compared with the previous three prediction methods, the reliability and stability have been significantly improved, but there is a problem of unstable prediction results for the same group of data each time. In recent years, the nonlinear programming method has been gradually applied to the field of parameter prediction due to its superlinear convergence and strong stability, and has achieved good results.

[0004] Regardless of the parameter prediction method used, measured data is required for calculation. However, due to the complexity of geological mining factors in various regions during the actual mining process, such as the complex structure of the rock and soil above the goaf, the changing structure and properties of the rock strata during the mining process, etc., as well as the errors caused by on-site monitoring and the model deviations of the probability integral model itself, the observed data often contain errors. In order to obtain more accurate prediction results, we need to introduce robust estimation theory. Robust estimation methods can effectively identify observations that are obviously inconsistent with other data or abnormal, thereby improving the quality and reliability of the data. By eliminating gross error data and limiting the use of available erroneous data, we can reduce the impact of data errors and make the final prediction results more accurate and reliable. Therefore, it is necessary to develop a new robust parameter prediction method to improve the accuracy and stability of the prediction. Summary of the invention

[0005] The purpose of the present invention is to provide an improved probability integral parameter testing method for unconstrained nonlinear programming, which has good parameter estimation performance, stability and strong anti-interference ability, and can effectively utilize the data of incomplete and non-standard observation stations to a certain extent.

[0006] To achieve the above object, the present invention provides the following technical solution: an improved probability integral parameter testing method for unconstrained nonlinear programming, the method comprising the following steps:

[0007] Step 1: Input the working face location parameters, measured data and geological mining parameters, select the initial probability integral parameters of the parameters to be determined based on the empirical values ​​of the probability integral parameters distinguished by the properties of the overburden, and set the iterative convergence accuracy and initial Hessian matrix;

[0008] Step 2: Use the unconstrained linear programming method to calculate the mining impact propagation angle θ0;

[0009] Step 3: Substitute the obtained θ0 into the probability integral prediction formula, establish a new parameter vector to be determined, calculate the predicted subsidence formula and predicted subsidence value of the observation point, combine the difference between the measured subsidence value and the predicted subsidence value, and update the weight matrix according to the IGGIII weight function;

[0010] Step 4: Establish the objective function f2 of the sinking coefficient q, the main influencing angle tangent tanβ, and the inflection point offsets S1, S2, S3, and S4;

[0011] Step 5: Use the unconstrained nonlinear programming method to invert the parameters q, tanβ, S1, S2, S3, and S4;

[0012] Step 6: Determine whether the desired probability integral parameter meets the convergence condition. If not, loop through steps 3-6 until the convergence condition is met and output the optimal parameter solution.

[0013] Step 7: At this time, θ0, q, tanβ, S1, S2, S3, and S4 are known. The measured values ​​of horizontal movement of the surface points and the predicted values ​​of horizontal movement are used to update the weight matrix according to the IGGIII weight function to establish the objective function f3;

[0014] Step 8: Use the unconstrained nonlinear programming method to invert the horizontal shift coefficient b;

[0015] Step 9: Determine whether the required parameter b meets the convergence condition. If not, loop through steps 7-9 until the convergence condition is met and output the solution for the optimal horizontal shift coefficient b.

[0016] Preferably, in step 1, the input geological mining parameters include normal mining thickness m, coal seam inclination α, dip length D1, strike length D3 and mining depth H, and the parameter vector A=[q, tanβ, θ0, S1, S2, S3, S4] is established. The initial probability integral parameters A0=[q0, tanβ0, θ00, S10, S20, S30, S40] and b0 of the parameters to be determined are randomly selected within a certain range according to the empirical values ​​of the probability integral parameters distinguished by the properties of the overlying rocks, and the iterative convergence precision ε1, ε2, ε3 and the initial Hessian matrix H0 are set.

[0017] Preferably, in step 2, the specific calculation steps are as follows:

[0018] Step 2.1: Substitute the x, y coordinates of the observation point at the working face and the geological mining parameters into the probability integral formula to calculate the predicted subsidence formula W(A) of the observation point;

[0019] The calculation model of the probability integral formula of the subsidence at any point on the surface is:

[0020] W(x,y)=W0*C xm m*C ym

[0021] Where m is the normal mining thickness; W0 is the maximum subsidence value, Cxm and Cym are the mining coefficients in the x and y directions respectively, and the calculation formula is:

[0022] W0=m*q*cosα

[0023] Where q is the subsidence coefficient; α is the coal seam inclination;

[0024]

[0025] In the formula, dλ is the unit thickness of the coal seam; H is the mining depth of the mining boundary; tanβ is the tangent of the main influencing angle; r is the main influencing radius, l is the calculated length of the strike working face, and L is the calculated length of the dip working face. Their calculation formulas are:

[0026]

[0027] l=D3-s3-s4

[0028]

[0029] Among them, D1 is the dip length of the working face, D3 is the strike length of the working face, θ0 is the propagation angle of mining influence, s1, s2, s3, s4 are the inflection point offsets of the downhill boundary, uphill boundary, left boundary and right boundary respectively;

[0030] Step 2.2: Establish the objective function f1 with parameter vector A as follows:

[0031]

[0032] Step 2.3: Use the unconstrained nonlinear programming method to invert the unknown parameter A. This time, the BFGS method is used, as shown below:

[0033] Step 2.3.1: Calculate the gradient of the objective function f1, that is, the partial derivative G k ;

[0034] Step 2.3.2: Set BFGS to obtain the convergence criterion ε1, and stop the BFGS calculation iteration when the norm of the gradient is less than ε1;

[0035] Step 2.3.3: Using the BFGS update rule, in each iteration, calculate the search direction d, where k = 0 and the initial Hessian matrix H0 is the identity matrix;

[0036] d=-H k *G k

[0037] Step 2.3.4: Use Armjio line search to determine the step size α so as to obtain feasible and significantly improved points in the search direction. The Armjio line search rule is as follows:

[0038] f(X k +αd k )≦f(X k )+cαG k T d k

[0039] Where c∈(0,1), this time c=10 -3 ;

[0040] Step 2.3.5: Update the current iteration point X k+1 , gradient G k+1 , and update the Hessian matrix using the BFGS update rule, the formula is as follows:

[0041] X k+1 =X k +α*d

[0042]

[0043] Among them, Y k =G k+1 -G k , S k =X k+1 -X k , if S k T *Yk >0, then v=Y k T *S k , otherwise H k+1 =H k ;

[0044] Step 2.3.6: Return to step 2.3.2 until the convergence condition is met and output the optimal parameter solution A k+1 =[q k+1 , tanβ k+1 ,θ0 k+1 、S1 k+1 、S2 k+1 、S3 k+1 、S4 k+1 ].

[0045] Preferably, in step 3, the specific steps of updating the weight matrix according to the IGGIII weight function are as follows:

[0046] Step 3.1: Create a new parameter vector Z to be determined, and set the initial probability integral parameter Z0 = [q k+1 , tanβ k+1 ,θ0 k+1 、S1 k+1 、S2 k+1 、S3 k+1 、S4 k+1 ], θ0 k+1 Substitute the probability integral subsidence calculation formula to calculate the predicted subsidence formula W(Z) of the observation point;

[0047] Step 3.2: Substitute Z0 into W(Z), calculate the predicted subsidence value W' of the observation point, and calculate the residual Vw and unit weighted mean error σ0 between the measured subsidence value W and the predicted subsidence value W':

[0048] Vw=W′-W

[0049]

[0050] Among them, P is the weight matrix, this time we take the unit matrix, n is the number of observations, and m is the number of unknowns;

[0051] Step 3.3: Update the weight matrix Pw according to the IGGIII weight function variable weight rule. The IGGIII weight variable rule is:

[0052]

[0053] The initial weight matrix PW0 is set to the unit matrix, where k0 and k1 are harmonic coefficients, k0 is set to 1.1 and k1 is set to 1.7; is the standardized residual, and

[0054]

[0055] Among them, V Wi is the residual of the ith observation.

[0056] Preferably, in step S4, the objective function f2 is as follows:

[0057]

[0058] Preferably, in step 5, the BFGS method is used to calculate the parameters q, tanβ, S1, S2, S3, and S4, and the specific steps are the same as step 2.3.

[0059] Preferably, in step 6, the convergence conditions are set as follows:

[0060] Z k+1 -Z k <ε2

[0061] Determine whether the desired probability integral parameter meets the convergence condition. If not, repeat steps 3-6 until the convergence condition is met and output the optimal parameter solution Z. k+1 .

[0062] Preferably, in step 7, the specific steps of establishing the objective function f3 are as follows:

[0063] Step 7.1: At this time, θ0, q, tanβ, S1, S2, S3, and S4 are known. Substituting the above seven parameters into the probability integral surface arbitrary point horizontal movement prediction model can generate the function U(b) about b:

[0064] The probability integral horizontal movement formula for any point on the surface is:

[0065]

[0066] Where φ is the angle from the positive direction of the x-axis to the specified direction counterclockwise; U ax and U ay are the horizontal displacement values ​​of the projection point positions on the strike and dip main sections respectively;

[0067]

[0068] Where b is the horizontal movement coefficient, b1 represents the horizontal movement coefficient in the downhill direction, b2 represents the horizontal movement coefficient in the uphill direction, H0, H1, H2 are the average mining depth, lower boundary mining depth and upper boundary mining depth respectively; tanβ1, tanβ2 are the main influencing angle tangent values ​​in the downhill direction and uphill direction respectively:

[0069] Step 7.2: Substitute the value of b (b0 for the first time) into U(b), calculate the predicted horizontal shift value U' of the observation point, calculate the predicted residual VU and the unit weight error σ0, and update the weight matrix PU according to the IGGIII weight function, where the initial weight matrix PU is the unit matrix, then the following formula is obtained:

[0070] V U =U′-U

[0071]

[0072] Among them, P is the weight matrix, this time we take the unit matrix, n is the number of observations, and m is the number of unknowns; Step 7.3: Establish the objective function f2 of the horizontal shift coefficient b:

[0073]

[0074] Preferably, in step 8, the specific steps of using the unconstrained nonlinear programming method to invert the horizontal shift coefficient b are the same as step 2.3.

[0075] Preferably, in step 9, the convergence condition ε3 is set as follows:

[0076] b k+1 -b k <ε3

[0077] Repeat steps 7-9 until the convergence condition is met and output the optimal horizontal shift coefficient b k+1 The solution.

[0078] Compared with the prior art, the present invention has the following beneficial effects:

[0079] The present invention reduces the sensitivity to outliers by applying robust estimation to unconstrained nonlinear programming inversion probability integral parameters and introducing a weight adjustment mechanism. It helps to weaken the adverse effects of outliers on unknown estimation to a certain extent, thereby improving the stability of the estimation, and has the advantages of high parameter accuracy and good curve fitting effect. At the same time, the method has good parameter estimation performance, stability and strong anti-interference ability, and to a certain extent, can effectively utilize the data of incomplete and non-standard observation stations, which can not only ensure the accuracy and reliability of parameter inversion, but also make full use of incomplete data and avoid waste of resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0080] Figure 1 is a flow chart of the method of the present invention;

[0081] Figure 2 is the sinking fitting diagram;

[0082] Figure 3 It is the horizontal shift fitting diagram. DETAILED DESCRIPTION

[0083] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0084] An improved probability integral parameter test method for unconstrained nonlinear programming is proposed. Taking a certain working face as an example, the real probability integral parameters of the working face are: [q, b, tanβ, θ0, S1, S2, S3, S4] = [0.89, 0.3, 2.37, 85, 50, 50, 50, 50]. The method comprises the following steps:

[0085] Step 1: Input the working face location parameters, measured data and geological mining parameters, select the initial probability integral parameters of the parameters to be determined based on the empirical values ​​of the probability integral parameters distinguished by the properties of the overburden, and set the iterative convergence accuracy and initial Hessian matrix;

[0086] The input geological mining parameters are: normal mining thickness m = 4m, coal seam dip angle α = 8°, dip length D1 = 300m, strike length D3 = 1000m and mining depth H = 300m. The parameter vector A = [q, tanβ, θ0, S1, S2, S3, S4] to be determined is established. The initial probability integral parameters A0 = [q0, tanβ0, θ00, S10, S20, S30, S40] and b0 of the parameters to be determined are randomly selected within a certain range according to the empirical values ​​of the probability integral parameters distinguished by the properties of the overburden. The iterative convergence precision ε1, ε2, ε3 and the initial Hessian matrix H0 are set.

[0087] Step 2: Use the unconstrained linear programming method to calculate the mining impact propagation angle θ0. The specific calculation steps are as follows:

[0088] Step 2.1: Substitute the x, y coordinates of the observation point at the working face and the geological mining parameters into the probability integral formula to calculate the predicted subsidence formula W(A) of the observation point;

[0089] The calculation model of the probability integral formula of the subsidence at any point on the surface is:

[0090] W(x,y)=W0*C xm *C ym

[0091] Where m is the normal mining thickness; W0 is the maximum subsidence value, Cxm and Cym are the mining coefficients in the x (strike) and y (dip) directions respectively. The calculation formula is:

[0092] W0=m*q*cosα

[0093] Where q is the subsidence coefficient; α is the coal seam inclination;

[0094]

[0095] In the formula, dλ is the unit thickness of the coal seam; H is the mining depth of the mining boundary; tanβ is the tangent of the main influencing angle; r is the main influencing radius, l is the calculated length of the strike working face, and L is the calculated length of the dip working face. Their calculation formulas are:

[0096]

[0097] l=D3-s3-s4

[0098]

[0099] Among them, D1 is the dip length of the working face, D3 is the strike length of the working face, θ0 is the propagation angle of mining influence, s1, s2, s3, s4 are the inflection point offsets of the downhill boundary, uphill boundary, left boundary and right boundary respectively;

[0100] Step 2.2: Establish the objective function f1 with parameter vector A as follows:

[0101]

[0102] Step 2.3: Use the unconstrained nonlinear programming method to invert the unknown parameter A. This time, the BFGS method is used. The specific steps are as follows:

[0103] Step 2.3.1: Calculate the gradient of the objective function f1, that is, the partial derivative G k ;

[0104] Step 2.3.2: Set BFGS to obtain the convergence criterion ε1, and stop the BFGS calculation iteration when the norm of the gradient is less than ε1;

[0105] Step 2.3.3: Using the BFGS update rule, in each iteration, calculate the search direction d, where k = 0 and the initial Hessian matrix H0 is the identity matrix;

[0106] d=-H k *G k

[0107] Step 2.3.4: Use Armjio line search to determine the step size α so as to obtain feasible and significantly improved points in the search direction. The Armjio line search rule is as follows:

[0108] f(X k +αdk )≦f(X k )+cαG k T d k

[0109] Where c∈(0,1), this time c=10 -3 ;

[0110] Step 2.3.5: Update the current iteration point X k+1 , gradient G k+1 , and update the Hessian matrix using the BFGS update rule, the formula is as follows:

[0111] X k+1 =X k +α*d

[0112]

[0113] Among them, Y k =G k+1 -G k , S k =X k+1 -X k , if S k T *Y k >0, then v=Y k T *S k , otherwise H k+1 =H k ;

[0114] Step 2.3.6: Return to step 2.3.2 until the convergence condition is met and output the optimal parameter solution A k+1 =[q k+1 , tanβ k+1 ,θ0 k+1 、S1 k+1 、S2 k+1 、S3 k+1 、S4 k+1 ];

[0115] Step 3: Substitute the obtained θ0 into the probability integral prediction formula, establish a new parameter vector to be determined, calculate the predicted subsidence formula and predicted subsidence value of the observation point, combine the difference between the measured subsidence value and the predicted subsidence value, and update the weight matrix according to the IGGIII weight function. The specific steps are as follows:

[0116] Step 3.1: Create a new parameter vector Z to be determined, and set the initial probability integral parameter Z0 = [q k+1 , tanβ k+1 ,θ0 k+1 、S1k+1 、S2 k+1 、S3 k+1 、S4 k+1 ], θ0 k+1 Substitute the probability integral subsidence calculation formula to calculate the predicted subsidence formula W(Z) of the observation point;

[0117] Step 3.2: Substitute Z0 into W(Z), calculate the predicted subsidence value W' of the observation point, and calculate the residual Vw and unit weighted mean error σ0 between the measured subsidence value W and the predicted subsidence value W':

[0118] Vw=W'-W

[0119]

[0120] Among them, P is the weight matrix, this time we take the unit matrix, n is the number of observations, this time we take 112, m is the number of unknowns, the value is 6;

[0121] Step 3.3: Update the weight matrix Pw according to the IGGIII weight function variable weight rule. The IGGIII weight variable rule is:

[0122]

[0123] The initial weight matrix PW0 is set to the unit matrix, where k0 and k1 are harmonic coefficients, k0 is set to 1.1 and k1 is set to 1.7; is the standardized residual, and

[0124]

[0125] Among them, V Wi is the residual of the ith observation;

[0126] Step 4: Establish the objective function f2 of the sinking coefficient q, the main influencing angle tangent tanβ, and the inflection point offsets S1, S2, S3, and S4 as follows:

[0127]

[0128] Step 5: Use the unconstrained nonlinear programming method to invert the parameters q, tanβ, S1, S2, S3, and S4, the same as step 2.3;

[0129] Step 6: Set the convergence conditions as follows:

[0130] Z k+1 -Z k <ε2

[0131] Determine whether the desired probability integral parameter meets the convergence condition. If not, repeat steps 3-6 until the convergence condition is met and output the optimal parameter solution Z.k+1 ;

[0132] Step 7: At this time, θ0, q, tanβ, S1, S2, S3, and S4 are known. The measured values ​​of horizontal movement of the surface points and the predicted values ​​of horizontal movement are used to update the weight matrix according to the IGGIII weight function and establish the objective function f3. The specific steps are as follows:

[0133] Step 7.1: At this time, θ0, q, tanβ, S1, S2, S3, and S4 are known. Substituting the above seven parameters into the probability integral surface arbitrary point horizontal movement prediction model can generate the function U(b) about b:

[0134] The probability integral horizontal movement formula for any point on the surface is:

[0135]

[0136] Where φ is the angle from the positive direction of the x-axis to the specified direction counterclockwise; U ax and U ay are the horizontal displacement values ​​of the projection point positions on the strike and dip main sections respectively;

[0137]

[0138] Where b is the horizontal movement coefficient, b1 represents the horizontal movement coefficient in the downhill direction, b2 represents the horizontal movement coefficient in the uphill direction, H0, H1, H2 are the average mining depth, lower boundary mining depth and upper boundary mining depth respectively; tanβ1, tanβ2 are the main influencing angle tangent values ​​in the downhill direction and uphill direction respectively:

[0139] Step 7.2: Substitute the value of b (b0 for the first time) into U(b), calculate the predicted horizontal shift value U' of the observation point, calculate the predicted residual VU and the unit weight error σ0, and update the weight matrix PU according to the IGGIII weight function, where the initial weight matrix PU is the unit matrix, then the following formula is obtained:

[0140] V U =U'-U

[0141]

[0142] Among them, P is the weight matrix, this time we take the unit matrix, n is the number of observations, 112, m is the number of unknowns, 1;

[0143] Step 7.3: Establish the objective function f2 of the horizontal shift coefficient b:

[0144]

[0145] Step 8: Use the unconstrained nonlinear programming method to invert the horizontal movement coefficient b, the same as step 2.3;

[0146] Step 9: Set the convergence condition ε3 as follows:

[0147] b k+1 -b k <ε3

[0148] Determine whether the desired parameter b meets the convergence condition. If not, loop through steps 7-9 until the convergence condition is met and output the optimal horizontal shift coefficient b. k+1 The solution.

[0149] Substitute the probability integral parameters obtained by inversion into the probability integral sinking amount and horizontal movement calculation formula to obtain the predicted sinking value and horizontal movement value. Compare them with the measured values ​​and draw a fitting curve. The mean error of sinking value fitting is 3.89mm, and the mean error of horizontal movement value fitting is 4.23mm. For details, please refer to Figure 2 and Figure 3 .

[0150] In order to better illustrate the effectiveness and superiority of the technical solution of the present invention, the following comparative analysis is made on the embodiments of the present invention after applying the above solution:

[0151] Some outliers are artificially added to the original data, and the probability integral parameters obtained by the conventional unconstrained nonlinear programming method and the method of the present invention are compared with the true probability integral parameters. Each experiment selects 5% and 10% of the observation points with the maximum sinking value, the edge and the inflection point, and adds 200mm of gross error as new test data. In order to prevent the accidental occurrence of the experiment, each experiment uses two methods to perform 30 inversions respectively. Table 1 gives the relative errors of the mean values ​​of the parameters inverted by the two algorithms. It can be seen that the method of the present invention shows better robustness: when the contamination rate is 5%, the relative errors of the three main parameters of q, b, and tanβ are controlled within 0.3%, 0.4%, and 0.2%, respectively, and the relative errors of the original method are all above 2%; when the contamination rate is 10%, the relative errors of the three main parameters of q, b, and tanβ are controlled within 0.6%, 0.4%, and 0.1%, respectively, and the relative errors of the original method are within 2.9%, 6.8%, and 10.3%, respectively.

[0152] Table 1 Comparison of parameter inversion results

[0153]

[0154]

[0155] The above description of the embodiments is for ordinary researchers in the field to understand and apply the present invention. Relevant personnel familiar with the technology in the field can easily modify the above embodiments and apply the relevant principles in the invention to other embodiments without creative work. Therefore, the present invention is not limited to the above embodiments. Those skilled in the art, based on the disclosure of the present invention, all technical solutions belonging to the ideas of the present invention belong to the protection scope of the present invention, and improvements and modifications made on the basis of the principles of the present invention should be within the protection scope of the present invention.

Claims

1. An improved probability integral parameter testing method for unconstrained nonlinear programming, characterized in that: The method comprises the following steps: Step 1: Input the working face location, measured data and geological mining parameters, select the initial probability integral parameters, set the iterative convergence accuracy and the initial Hessian matrix; Step 2: Use the unconstrained linear programming method to calculate the mining impact propagation angle θ0; Step 3: Substitute θ0 into the probability integral formula, calculate the predicted subsidence value, compare it with the measured value, and update the weight matrix through the IGGIII weight function; Step 4: Establish the objective function f2, including the subsidence coefficient q, the main influence angle tangent tanβ and the inflection point offsets S1, S2, S3 and S4, and use the unconstrained nonlinear programming method to invert the above parameters; Step 5: Use the unconstrained nonlinear programming method to invert the parameters q, tanβ, S1, S2, S3, and S4; Step 6: Determine whether the inversion parameters meet the convergence conditions. If not, repeat steps 3-5 until convergence; Step 7: Use the measured and predicted values ​​of horizontal movement to update the weight matrix and establish the objective function f3; Step 8: Invert the horizontal shift coefficient b; Step 9: Determine whether parameter b meets the convergence condition. If not, loop through steps 7-9 until the convergence condition is met and output the optimal solution.

2. The improved probability integral parameter testing method for unconstrained nonlinear programming according to claim 1, characterized in that: In the step 1, the input geological mining parameters include normal mining thickness m, coal seam inclination α, dip length D1, strike length D3 and mining depth H, establish the parameter vector A to be determined, select the initial probability integral parameters A0 and b0 according to the empirical values ​​of the probability integral parameters distinguished by the properties of the overburden, set the iterative convergence precision ε1, ε2, ε3, and the initial Hessian matrix H0.

3. The improved probability integral parameter testing method for unconstrained nonlinear programming according to claim 1, characterized in that: In step 2, the specific calculation steps are as follows: Step 2.1: Substitute the coordinates x, y of the observation point at the working face and the geological mining parameters into the probability integral formula to calculate the predicted subsidence formula W(A) of the observation point; The calculation model of the probability integral formula of the subsidence at any point on the surface is: W(x,y)=W0*C xm **C ym Where m is the normal mining thickness; W0 is the maximum subsidence value, Cxm and Cym are the mining coefficients in the x and y directions respectively, and the calculation formula is: W0=m*q*cosα Where q is the subsidence coefficient; α is the coal seam inclination; Where dλ is the unit thickness of the coal seam; H is the mining depth of the mining boundary; tanβ is the tangent of the main influencing angle; r is the main influencing radius, 1 is the calculated length of the strike working face, and L is the calculated length of the dip working face. The calculation formula is: l=D3-s3-s4 Among them, D1 is the dip length of the working face, D3 is the strike length of the working face, θ0 is the propagation angle of mining influence, s1, s2, s3, s4 are the inflection point offsets of the downhill boundary, uphill boundary, left boundary and right boundary respectively; Step 2.2: Establish the objective function f1 with parameter vector A as follows: Step 2.3: Use the BFGS method to invert the parameter A, as shown below: Step 2.3.1: Calculate the gradient G of the objective function f1 k ; Step 2.3.2: Set BFGS to obtain the convergence criterion ε1; Step 2.3.3: Use the BFGS update rule to calculate the search direction d in each iteration, where k = 0 and the initial Hessian matrix H0 is the identity matrix; d=-H k *G k Step 2.3.4: Use Armjio line search to determine the step size α. The Armjio line search rule is as follows: f(X k +αd k )≤≦f(X k )+cαG k T d k Where c∈(0,1), this time c=10 -3 ; Step 2.3.5: Update the current iteration point X k+1 , gradient G k+1 , and update the Hessian matrix using the BFGS update rule, the formula is as follows: X k+1 =X k +a*d Among them, Y k =G k+1 -G k , S k =X k+1 -X k , if S k T *Y k >0, then v=Y k T *S k , otherwise H k+1 =H k ; Step 2.3.6: Return to step 2.3.2 until the convergence condition is met and output the optimal parameter solution A k+1 .

4. The improved probability integral parameter testing method for unconstrained nonlinear programming according to claim 1, characterized in that: In step 3, the specific steps of updating the weight matrix according to the IGGIII weight function are as follows: Step 3.1: Create a new parameter vector Z to be determined, set the initial probability integral parameter Z0, and replace θ0 k+1 Substitute the probability integral subsidence calculation formula to calculate the predicted subsidence formula W(Z) of the observation point; Step 3.2: Substitute Z0 into W(Z), calculate the predicted subsidence value W' of the observation point, and calculate the residual Vw and unit weighted mean error σ0 between the measured subsidence value W and the predicted subsidence value W': Vw=W'-W Among them, P is the weight matrix, n is the number of observations, and m is the number of unknowns; Step 3.3: Update the weight matrix Pw according to the IGGIII weight function variable weight rule. The IGGIII weight variable rule is: The initial weight matrix PW0 is set to the unit matrix, where k0 and k1 are the harmonic coefficients. is the standardized residual, and: Among them, V Wi is the residual of the ith observation.

5. The improved probability integral parameter testing method for unconstrained nonlinear programming according to claim 1, characterized in that: In step S4, the objective function f2 is as follows:

6. The improved probability integral parameter testing method for unconstrained nonlinear programming according to claim 3, characterized in that: In step 5, the specific steps of using the BFGS method to calculate parameters are the same as step 2.

3.

7. The improved probability integral parameter testing method for unconstrained nonlinear programming according to claim 1, characterized in that: In step 6, the convergence conditions are set as follows: WITH k+1 -WITH k <ε2 Z k+1 Output the optimal parameter solution.

8. The improved probability integral parameter testing method for unconstrained nonlinear programming according to claim 1, characterized in that: In step 7, the specific steps of establishing the objective function f3 are as follows: Step 7.1: Substitute θ0, q, tanβ, S1, S2, S3, and S4 into the probability integral surface arbitrary point horizontal movement prediction model to generate the function U(b) about b: The probability integral horizontal movement formula for any point on the surface is: Where φ is the angle from the positive direction of the x-axis to the specified direction counterclockwise; U ax and U ay are the horizontal displacement values ​​of the projection point positions on the strike and dip main sections respectively; Where b is the horizontal movement coefficient, b1 represents the horizontal movement coefficient in the downhill direction, b2 represents the horizontal movement coefficient in the uphill direction, H0, H1, H2 are the average mining depth, lower boundary mining depth and upper boundary mining depth respectively; tanβ1, tanβ2 are the main influencing angle tangent values ​​in the downhill direction and uphill direction respectively: Step 7.2: Substitute the value of b into U(b), calculate the predicted horizontal shift value U' of the observation point, calculate the predicted residual VU and the unit weight error σ0, update the weight matrix PU according to the IGGIII weight function and use it as the unit matrix, then we have the following formula: V U =U'-U Among them, P is the weight matrix, n is the number of observations, and m is the number of unknowns; Step 7.3: Establish the objective function f2 of the horizontal shift coefficient b:

9. The improved probability integral parameter testing method for unconstrained nonlinear programming according to claim 3, characterized in that: In step 8, the specific steps for inverting the horizontal shift coefficient b are the same as step 2.

3.

10. The improved probability integral parameter testing method for unconstrained nonlinear programming according to claim 1, characterized in that: In step 9, the convergence condition ε3 is set as follows: b k+1 -b k <e3 b k+1 is the optimal horizontal shift coefficient solution for the output.

Images