Multi-rock mass fracture parameter prediction method and rock mass fracture model construction method based on Markov chain and Kriging method

Abstract

The invention discloses a rock mass fracture parameter prediction method and a rock mass fracture model construction method. The prediction method comprises the following steps: determining the scale of a fracture unit by using a first-order Markov test theory, and dividing the fracture unit of a known parameter region according to the scale of the fracture unit; performing optimal transfer matrix calculation on the fracture density of the divided fracture units by using a first-order Markov chain theory; performing fracture density prediction on the unknown parameter region by using the optimal transfer matrix obtained by calculation, and completing height data distribution of a trace midpoint and a fracture surface in a uniform distribution form; using a Kriging interpolation estimation method to complete occurrence prediction of the fracture at the distributed position; and using a Kriging interpolation estimation method to complete trace length prediction of the fracture at the distributed position. According to the method, the Markov chain and Kriging interpolation estimation methods are fused to predict and estimate the density, occurrence, trace length and other parameters of rock mass fractures, and statistical errors are greatly reduced.

Description

Technical Field

[0001] The invention belongs to the technical field of geological survey, in particular to a method for predicting rock mass crack parameters for foundation engineering, and a rock mass crack grid simulation method constructed by using the method. Background Art

[0002] The existing rock mass fracture network modeling method is mainly based on the discrete fracture network method. However, the correlation between parameters is not considered when using the discrete fracture network modeling method to estimate the parameters of unknown areas. In addition, there are studies on fracture parameter modeling using geostatistical methods, but this application is mostly used for the prediction of fracture density, while other parameters, such as occurrence, size and other information, still need to use mathematical statistics methods to achieve fracture network modeling. Summary of the invention

[0003] 1. Technical issues to be resolved

[0004] The technical problem to be solved by the present invention is that the existing method for predicting discrete rock fracture network parameters cannot take into account the correlation between rock fracture parameters. In addition, although the fracture parameter prediction method combined with the geostatistical method takes into account the parameter correlation, it is only used for parameter prediction of fracture density, and has the problem of low accuracy and reliability compared with the geological fracture network distribution model established by using a scientific correlation prediction method for multiple fracture parameters.

[0005] (II) Technical solution

[0006] In order to solve the above technical problems, the present invention proposes a rock mass fracture parameter prediction method, which is used to predict the rock mass fracture parameters of a target area, wherein the target area includes a known parameter area and an unknown parameter area, and the method comprises the following steps:

[0007] S1. Determine the scale of the fracture unit using the first-order Markov property test theory, and divide the fracture units in the known parameter region according to the scale of the fracture unit;

[0008] S2, using the first-order Markov chain theory to calculate the optimal transfer matrix for the crack density of the crack unit in the known parameter area divided in step S1 above;

[0009] S3, using the optimal transfer matrix calculated above to predict the crack density in the unknown parameter area, and using a uniform distribution form to complete the distribution of the trace midpoint and crack surface height data;

[0010] S4. Use Kriging interpolation estimation method to complete the prediction of the fracture occurrence at the assigned location;

[0011] S5. Use the Kriging interpolation estimation method to complete the prediction of the trace length of the crack at the assigned location.

[0012] According to a preferred embodiment of the present invention, step S1 comprises:

[0013] Step S1.1: Calculate the adjacent order likelihood ratio test statistic η t-1,t , the calculation formula is as follows:

[0014] Where T ij,…,kl and T j,…,kl The order of continuous t-step state transition is S i →S j →…→S k →S l The number of occurrences and the order of continuous t-1-step state transitions are S j →…→S k →S l The number of occurrences of

[0015] Step S1.2: Calculate the likelihood ratio test statistic η of any order k,q , the calculation formula is as follows:

[0016] η k,q =η k,k+1 +η k+1,k+2 +…+η q-1,q ;

[0017] Step S1.3:

[0018] Calculate the loss function, R(k), that is

[0019] R(k)=η k,q -(m M -m k )(m-1);

[0020] Step S1.4: Query degrees of freedom (m q -m k The chi-square critical value corresponding to the 95% significance level under (m-1) is compared with the calculated R value. The k corresponding to the minimum R value is the most likely order of the Markov chain.

[0021] Step S1.5: Determine whether the best order is 1st order. If not, change the size of unit L and repeat the above steps S1.1 to S1.5.

[0022] According to a preferred embodiment of the present invention, step S5 comprises:

[0023] Step 5.1:

[0024] Calculate the semivariance function value between known data sample points. The calculation formula is as follows:

[0025]

[0026] Where H is the step size between sample points, γ(H) represents the experimental semivariogram function, N(H) represents the logarithm of two sample points with a step size of H, s represents the sth pair of sample points, and x(s) represents the sample value of one of the sth pair of sample points;

[0027] Step 5.2:

[0028] The spherical model is used to fit the theoretical variogram. The spherical model formula is as follows:

[0029]

[0030] Where γ(h) represents the theoretical semivariogram function value between sample points with a step size of h, c0 represents the nugget value, c represents the sill value, a represents the range, and h represents the step size (h is a variable);

[0031] Step 5.3:

[0032] Construct the Kriging equation system, the formula is as follows:

[0033]

[0034] where γ ij represents the theoretical semivariance function value between the i-th sample point and the j-th sample point, γ i0 represents the theoretical semivariance function value between the i-th sample point and the predicted point x0, λ i is the weight coefficient, φ is the deviation value;

[0035] Step 5.4:

[0036] The formulas for calculating the Kriging estimate and Kriging variance are as follows:

[0037]

[0038]

[0039] Among them, Z * (x0) is the Kriging estimate, σ E 2 is the Kriging variance.

[0040] A second aspect of the present invention provides a method for constructing a rock mass fracture model, comprising a rock mass fracture parameter prediction method according to any one of claims 1 to 3.

[0041] According to a preferred embodiment of the present invention, the rock mass fracture model construction method further comprises the following steps:

[0042] S6. Use the coordinate inverse calculation method to complete the xOy coordinate inverse calculation of the fracture center and the calculation of the diameter through the existing coordinates, occurrence and trace length information;

[0043] S7. Establish a three-dimensional fracture network model of the target area based on the predicted rock fracture parameters.

[0044] According to a preferred embodiment of the present invention, step S6 comprises:

[0045] Step 6.1: Determine the inclination of the fracture surface and π / 2 relationship;

[0046] Step 6.2: Based on the above relationship, make the following calculation: When the xOy coordinates of the trace are equal to π / 2, the coordinates of the center point of the fracture surface are consistent. When π / 2 is not equal to π / 2, the following formula is used for calculation:

[0047] Where s2 is the inverse vector of the projection vector of the normal vector to the xOy plane, V UnitX and V UnitY represents the unit vectors of the x-axis and y-axis, H is the distance from the center of the crack surface to the xOy plane, is the inclination angle of the crack surface, and x0 and y0 represent the coordinates of the midpoint of the crack trace.

[0048] (III) Beneficial effects

[0049] The present invention combines the Markov chain and Kriging interpolation estimation methods to predict and estimate multiple parameters such as the density, occurrence and trace length of rock fractures, and realizes the calculation of fracture surface diameter parameters in combination with the above parameters. Only the distribution hypothesis method of mathematical statistics is used to determine the position. The distribution of the fracture center position points in such a small unit is arranged relative to the position distribution in a large space, which greatly reduces the statistical error.

[0050] The present invention adopts a Markov chain method for predicting rock mass fracture density, which is determined by the previous and next states between adjacent units and supports estimation using a small amount of data, which has higher practicability in practical applications.

[0051] The rock mass fissure parameter prediction method and rock mass fissure simulation method of the present invention can be applied to water conservancy and hydropower engineering, tunnel excavation engineering, etc., and have broad application prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1It is a schematic diagram of simulation units of an embodiment of a rock mass fracture simulation method constructed using the rock mass fracture parameter prediction method of the present invention.

[0053] Figure 2 It is a schematic diagram of the principle of inverse calculation of the fracture surface diameter and the xOy coordinates of the center point using the fracture trace length and its midpoint coordinates, the height coordinates of the fracture surface center point and the occurrence data.

[0054] Figure 3 It is a flow chart of the rock mass fracture parameter prediction method based on the combination of Markov chain and Kriging method of the present invention.

[0055] Figure 4 It is the best classification graph determined by using the elbow rule combined with the K-means++ algorithm in one embodiment of the present invention.

[0056] Figure 5 FIG. 1 is a schematic diagram of classification poles according to an embodiment of the present invention. Figure 5 As shown, the grouping results are displayed as a polar diagram obtained using the stereographic projection method.

[0057] Figure 6 1 is a stability probability diagram corresponding to four groups of cracks under different K according to an embodiment of the present invention.

[0058] Figure 7 It is a comparison chart between the prediction results and the actual results of the prediction interval of an embodiment of the present invention.

[0059] Figure 8 It is a prediction result occurrence distribution pole diagram of an embodiment of the present invention.

[0060] Fig. 9 It is a three-dimensional fracture network distribution map of the prediction area according to an embodiment of the present invention. DETAILED DESCRIPTION

[0061] In order to solve the above-mentioned technical problems, the present invention adopts the Markov chain model and Kriging interpolation estimation method in the geostatistical method, and in specific cases involves using the above methods in combination to predict multiple rock fracture parameters, and using the predicted parameters to establish a more reliable fracture network model.

[0062] Specifically, the present invention proposes a method integrating Markov chain and Kriging interpolation estimation to predict and estimate multiple parameters such as density, occurrence and trace length of rock fractures, and only uses the distribution hypothesis method of mathematical statistics to determine the position.

[0063] The present invention uses a Markov chain to grid the known parameter area in the target area to obtain a number of small-scale units. Arranging the distribution of crack center position points in such small-scale units will greatly reduce the statistical error relative to the position distribution in a large space.

[0064] In linear distribution areas, it is difficult to achieve a good prediction of the fracture density in linear unknown areas by simply relying on the Kriging interpolation estimation method, because the better application scenario of the Kriging method is the prediction of spatial distribution. A big reason for this feature is that there are abundant fracture data in the geological space, so the parameters of other unknown parameter areas in the target space can be statistically analyzed. However, when the data is concentrated on a line, its position information variable is single, so it is difficult to have a good correlation with the geological attributes. Therefore, for the prediction of fracture density, we adopted the Markov method, which is determined by the previous and next states between adjacent units and supports the use of a small amount of data for estimation, which has higher practicality in practical applications and is therefore very suitable for the prediction of this parameter. On this basis, after determining the position of the trace, the input-output mapping relationship between the fracture trace position and the occurrence and trace length can be established in the target area. After that, we can use the Kriging interpolation method to predict the occurrence and trace length.

[0065] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with specific embodiments and with reference to the accompanying drawings.

[0066] Figure 1 It is a schematic diagram of simulation units of an embodiment of a rock mass fracture simulation method constructed using the rock mass fracture parameter prediction method of the present invention.

[0067] like Figure 1 As shown in the figure, in this embodiment, the blue area in the simulation unit represents the location of the known parameter area distribution. The data of the parameters such as fracture density, trace length and occurrence used for prediction are all from this area. The gray area represents the location where the fracture parameters to be predicted are distributed, that is, the unknown parameter area. The red and green circular planes represent the fracture surface. The numbers in the table in the figure represent the number of fractures in each unit after the geological simulation unit division.

[0068] Figure 2 It is a principle diagram for inverse calculation of fracture surface diameter and center point xOy coordinates using fracture trace length and its midpoint coordinates, fracture surface center point height coordinates and occurrence data. As shown in the figure, n1 is the normal vector of the fracture surface, is the inclination angle, θ is the dip angle, H is the distance from the center of the fracture surface to the xOy plane, L is the vertical distance from the projection point of the fracture surface center point O1 to xOy to the trace, s1 is the projection vector of the normal vector to the xOy plane, s2 is the inverse vector of s1, and O2 is the coordinate of the midpoint of the trace.

[0069] Figure 3 This is a flow chart of the rock mass fracture parameter prediction method based on the combination of Markov chain and Kriging method of the present invention. Figure 3 As shown, the method includes:

[0070] S1. Use the first-order Markov property test theory to determine the scale of the fracture unit, and divide the fracture units in the known parameter area according to the fracture unit scale.

[0071] First, the first-order Markov property test theory is used to determine the scale of the fracture unit along the known data direction, and the fracture units in the known parameter area in the target area are divided according to the fracture unit scale. Moreover, after the fracture unit division is completed, the fracture density distribution of the fracture units in all known parameter areas is obtained, and a reasonable upper limit value is set for these data.

[0072] Generally, the setting of the upper limit is determined by the proportion of the maximum density value. When the proportion is small, the density value can be set to a level lower than the density value, so that the transition property can better meet the first-order requirements. Markov Chain is a random process used to describe the transition of a system between a series of states. Its core feature is "memorylessness". For a first-order Markov chain, the probability of the system transitioning from state A to state B is only affected by the current state A, not by earlier states. For multi-order Markov chains, the current state depends on the previous states (for example, a second-order Markov chain depends on the first two states). Higher-order models can better capture complex dependencies, but they also increase computational complexity.

[0073] The process and corresponding formula of the first-order test are as follows:

[0074] First, select the initial fracture element size L and the highest test order M

[0075] Step S1.1: Calculate the adjacent order likelihood ratio test statistic η t-1,t , the calculation formula is as follows:

[0076]

[0077] Where T ij,…,kl and T j,…,kl The order of continuous t-step state transition is S i →S j →…→S k →S l The number of occurrences and the order of continuous t-1-step state transitions are S j →…→S k →S lThe number of occurrences of i, j, k, l, and m respectively represent the current fracture density values ​​in the geological unit, and m represents the highest fracture density value in the unit.

[0078] Step S1.2:

[0079] Compute the likelihood ratio test statistic η of any order k,q , the calculation formula is as follows:

[0080] η k,q =η k,k+1 +η k+1,k+2 +…+η q-1,q (2)

[0081] Where k and q represent the Markov transition order as k and q respectively.

[0082] Step S1.3:

[0083] Calculate the loss function, R(k), that is

[0084] R(k)=η k,q -(m M -m k )(m-1) (3)

[0085] Where M represents the highest test order, and m represents the highest unit fracture density in each group of simulated geological units.

[0086] Step S1.4:

[0087] The query freedom is (m q -m k The chi-square critical value corresponding to the 95% significance level under )(m-1) is compared with the calculated R values. The k corresponding to the minimum R value is the most likely order of the Markov chain.

[0088] Step S1.5:

[0089] Determine whether the best order obtained is 1st order. If not, change the size of unit L and repeat the above steps S1.1 to S1.5.

[0090] According to the above first-order Markov property test theory, the code for testing the first-order Markov property of drilling is edited using Matlab language. The pseudo code steps are as follows:

[0091] Input: fracture coordinate data D along the predicted direction, initial unit size L and highest order M; Output: k / / optimal order calculated according to L;

[0092] step:

[0093]

[0094]

[0095] The crack coordinate data D along the length direction is input into the UnitDensit algorithm. The algorithm counts the density of cracks in each unit along the length direction according to the coordinates and the set unit scale L, and finally obtains J NUM The cell density data is obtained, and then the TransferSt algorithm is used to count the transfer attribute categories and times between adjacent cells of the above-divided cell data, and the result is O Order , using the adjacent order parameter method AdjacentMethod to calculate the parameter N Adj For example, the 1st and 2nd orders are defined as adjacent orders, and the 1st and 3rd orders are defined as non-adjacent orders. The adjacent order statistics method AdjacentSta is used to calculate the result η Adj , using the non-adjacent order method NotAdjacentSta to calculate η NotAdj , use the loss function calculation method LossFuncMethod to calculate the loss function value R LossFunc , use the optimal order judgment method OrderJudgement method to determine whether the first-order property is met. If it is not satisfied, change the unit size L and recalculate.

[0096] S2. Using the first-order Markov chain theory, the optimal transfer matrix is ​​calculated for the fracture density in the known parameter region of the fracture unit divided in step S1.

[0097] This step determines the optimal transfer probability matrix of the known parameter area of ​​the target area. The theory used in this process is Walther's phase sequence law. This law can be summarized as the rock and soil body has the same deposition order in different directions, but the deposition scale is different, and the deposition scale in the horizontal direction is much larger than the deposition scale in the vertical direction. The conditions applied in the present invention are based on the fact that the frequency of transfer between the same attributes has a high proportion in geology, and the transfer between different attributes can be simulated according to the results obtained from the known unit transfer data, so as to only correct the transfer of these same attributes to improve the accuracy of the prediction. The theory is as follows:

[0098] When the correction coefficient K is taken, the principle of getting a better transfer matrix is ​​as follows: First, select an observation interval in the unknown parameter region, in which the attributes of the data are known, and each K value corresponds to its own prediction result. The correction value with the highest degree of match between the actual value and the predicted value is recorded as the optimal K. The degree of match here is expressed as the stable probability P l To judge, the formula is as follows:

[0099]

[0100] Where P(Z i =S ki ) is the attribute S in unit i ki The probability of occurrence is expressed as N(Z i =S ki ) / N sim Calculate its probability value, where N(Z i =S ki ) is N sim Simulation times Z i The attribute in the unit is S ki The number of times.

[0101] This section includes two theoretical contents, namely the first-order Markov chain model and the determination of stability probability. The pseudo code flow is as follows:

[0102] Input: The density data D of the fracture unit determined according to the above L Density , correction K, simulation times N s , the target area has known unit location X LocKnown and the number of parameters N XLocKnown ; The number of K is N K ; Output: optimal K / / optimal correction coefficient of transfer matrix;

[0103] step:

[0104]

[0105]

[0106] The D calculated in the first step Density As input data, it is input into the TransferMatrix algorithm, which is used to calculate the transfer probability of different densities within the unit, and finally obtain the transfer probability matrix M TransferProb , using the transfer probability matrix correction algorithm CorrectTransProb, the transfer probability matrix is ​​corrected according to the set K value, and the result M is obtained CTP Then, the one-dimensional Markov chain method SimProbrobMethod is used to calculate the transition probability from the known attribute to the specified attribute of the prediction unit, and then the roulette method RouletteMethod is used to determine the attribute of the final prediction unit. The above process is repeated until the attributes in all units are predicted. The stable probability calculation method StationaryPM is combined with the selected verification interval data V VarificationSec Calculate the stability probability P l Repeat the above process N sFinally, MeanMethod is used to calculate the average stability probability, and the ComparedKMethod method is used to compare and determine the K value corresponding to the maximum stability probability, which is recorded as the optimal correction coefficient K of the fracture density transfer in this group.

[0107] S3. Use the optimal transfer matrix calculated above to predict the crack density in the unknown parameter area, and use a uniform distribution to complete the distribution of trace midpoint and crack surface height data.

[0108] The implementation steps of step S3 are the same as those of step S2, except that the calculation of the stability probability is reduced.

[0109] The uniform distribution is used to simulate the crack position, and the three-dimensional position data is obtained based on the Monte Carlo method. The pseudo code steps are as follows:

[0110] Input: The crack density data D of the prediction unit obtained from step S3 DensityPre , length L, width W, height H parameters of each unit, D DensityPre The number of data N DensitPre ;

[0111] Output: D Location / / Fissure point location information;

[0112] step:

[0113] <1> for i=1…N Density do

[0114] <2> D Location ←FractureLocMethod(L,W,H,D Density ) / / Calculate the three-dimensional crack position point

[0115] <3> end for

[0116] D DensityPre , L, W and H are input into the FractureLocMethod algorithm as input data. The algorithm retrieves the density value in each prediction unit one by one, determines the number of randomly generated fracture 3D coordinate points according to uniform distribution, and finally outputs the 3D coordinate data D of all predicted fracture points. Location .

[0117] S4. Use the Kriging interpolation estimation method to complete the prediction of the fracture occurrence at the assigned location.

[0118] The process of this step is as follows:

[0119] Step 5.1:

[0120] Calculate the semivariance function value between known data sample points. The calculation formula is as follows:

[0121]

[0122] Where H is the step size, which is determined according to the number of groups set by oneself, and the above semivariance function value is the average value calculated from the data in the H interval. That is, the distance value between the two sample points selected to be calculated is a constant.

[0123] Step 5.2:

[0124] The spherical model is used to fit the theoretical variogram. The spherical model formula is as follows:

[0125]

[0126] Where c0 is the nugget value, c is the base value, a is the range, and h is the step size (h is a variable).

[0127] Step 5.3:

[0128] Construct the Kriging equation system, the formula is as follows:

[0129]

[0130] where λ i is the weight coefficient, and φ is the deviation value.

[0131] Step 5.4:

[0132] The formulas for calculating the Kriging estimate and Kriging variance are as follows:

[0133]

[0134] The pseudo code steps for using the Kriging interpolation estimation method to predict the occurrence are as follows:

[0135] Input: known unit crack coordinate data D and direction data O, predicted unit two-dimensional coordinate information D Location ;

[0136] Output: D PredOrientation / / Crack prediction location point occurrence information;

[0137] step:

[0138]

[0139] The fracture coordinate data D and the orientation data O are input into this program as input data. The normalization transformation method Boxcox is used to normalize D. Then the experimental semivariogram function calculation function in the Kriging interpolation toolkit in Matlab is used to calculate γ Real , using the theoretical semivariance fitting function combined with the spherical model, the fitting model type can also be changed according to actual needs to calculate S Struct , Kriging interpolation function Kriging is used to treat D in the prediction area Location The fracture orientation parameters at the position are predicted to obtain D PredONormal Finally, the formula in the above pseudo code algorithm is used to perform inverse transformation to obtain the occurrence direction data D of the final prediction point. PredOrientation .

[0140] S5. Use the Kriging interpolation estimation method to complete the density prediction of the cracks at the assigned locations.

[0141] The implementation steps of step S5 are the same as those of step S4, except that the output result of the prediction process is changed from occurrence angle to trace length data.

[0142] The above is the rock mass fracture parameter prediction method to be invented. Based on the rock mass fractures predicted by the above method, the present invention also proposes a rock mass fracture model construction method. In addition to the steps of the rock mass fracture parameter prediction method, the rock mass fracture model construction method also includes the following steps.

[0143] S6. Use the coordinate inverse calculation method to complete the xOy coordinate inverse calculation of the fracture center and the diameter calculation through the existing coordinates, occurrence and trace length information.

[0144] The sixth step is to use the coordinate inversion method to determine the coordinates of the center point of the crack surface on the xOy plane and its diameter D. The specific theoretical process is as follows:

[0145] Step 6.1: Determine the inclination angle and π / 2 relationship;

[0146] Step 6.2: Based on the above relationship, make the following calculation. When the two are equal, the coordinates of the trace xOy plane and the coordinates of the center point of the crack plane are consistent. When the two are not equal, use the following formula to calculate:

[0147] Where s2 is the inverse vector of the projection vector of the normal vector to the xOy plane, V UnitX and V UnitY represents the unit vectors of the x-axis and y-axis, H is the distance from the center of the crack surface to the xOy plane, is the inclination angle of the fracture surface, x0 and y0 represent the midpoint coordinates of the fracture trace length, and the above parameters are explained in detail in Figure 2 .

[0148] The calculation formulas for other parameters used in the calculation process are as follows:

[0149]

[0150] The pseudo code steps for back-calculating the coordinates of the midpoint of the original fracture disk from the midpoint of the trace using the above parameters are as follows:

[0151] Input: height coordinate H of the center point of the fracture surface, inclination angle θ and dip angle of the fracture surface Coordinates of the center point of the trace x0 and y0;

[0152] Output: plane coordinates x and y of the center point xOy of the crack surface, and the vertical distance L from the projection point of the center point O1 of the crack surface to xOy to the trace line;

[0153] step:

[0154]

[0155] The height coordinate H of the center point of the fracture surface, the dip angle θ and the inclination angle The coordinates x0 and y0 of the center point of the trace are input into the above program as known data. First, a judgment is performed. For vertical cracks, the midpoint of the trace coincides with the midpoint of the disk model, so the two are equal. For cracks at other angles, the coordinate inversion formula in the above pseudo code is used for calculation.

[0156] The method for inversely calculating the diameter length using the existing trace length data is as follows: The distance between O1 and O2 is The calculation formula of the disk model diameter D is as follows:

[0157]

[0158] The pseudo code steps for calculating the diameter D of the disk model are as follows:

[0159] Input: fracture surface trace length data L Trace and the height coordinate H of the center point, the vertical distance L from the projection point of the fracture surface center point O1 to xOy to the trace line;

[0160] Output: fracture surface diameter D;

[0161] step:

[0162] <1> for i=1…N LTrace do / / Calculate the diameter of the crack surface

[0163] <2>

[0164] <3> end for

[0165] The trace length data LTrace The center point height coordinate H and the vertical distance L from the projection point of the fracture surface center point O1 to xOy to the trace are input into the algorithm, and the above formula is used to calculate the diameter D corresponding to all trace lengths.

[0166] S7. Establish a three-dimensional fracture network model of the target area based on the predicted rock fracture parameters.

[0167] Finally, a relevant drawing method is used to combine the parameter information obtained from the above predictions and calculations to obtain a three-dimensional fracture network model.

[0168] An example is used below to illustrate the above process of the present invention, based on the tunnel inner wall crack data of the curtain grouting gallery.

[0169] First, the first part is to divide the simulation unit scale. In terms of data selection, for small-scale units such as tunnels, the crack data of the tunnel can only be counted from the crack data of the bottom plate trace. At the same time, in the first part of data processing, data can be deleted according to the actual data distribution. For example, when the density value of a certain crack in the unit is relatively small, the upper limit of the crack density can be set in the form of reduction, so that the first-order requirements can be better met during the simulation.

[0170] Figure 4 It is the best classification graph determined by using the elbow rule combined with the K-means++ algorithm in one embodiment of the present invention.

[0171] From the clustering results obtained by the K-means++ method and combined with the elbow rule, it can be seen that the best grouping result of these crack data is four groups.

[0172] Figure 5 FIG. 1 is a schematic diagram of classification poles according to an embodiment of the present invention. Figure 5 As shown, the grouping results are displayed as a polar diagram obtained using the stereographic projection method.

[0173] Next, the first-order Markov property test method is used to determine the size of the simulation unit for each set of results. The calculation results are as follows:

[0174] Table 1 Calculation results of statistics for each group

[0175]

[0176]

[0177] Table 2 Calculation results of loss functions for each group

[0178]

[0179] The above results show that the unit length of the four groups is 2 meters. At the same time, considering the requirement of data change stability, the upper limit of the density in the unit is restricted after the above unit division is completed. The upper limits of the four groups are set to 4, 1, 2 and 3 respectively. The first-order Markov chain theory will be used later to determine the optimal transfer matrix.

[0180] Figure 6 : is a stability probability diagram corresponding to four groups of cracks under different K in one embodiment of the present invention. Figure 6 As shown in the figure, by analyzing the stability probability of the known data interval, it can be seen that when K of the four groups is 1, 0.25, 0.85 and 1, the maximum stability probability Pl is 4.456×10-5, 0.0065, 0.0398 and 0.03 respectively. The optimal transfer matrix of each group obtained by the above correction coefficient will be used to simulate the fracture density parameters of the prediction area.

[0181] The four groups of transfer probability matrices were modified using the above parameters to simulate the prediction interval. Due to the randomness of this method, we set the number of simulations to 10,000. The average accuracy of the prediction results was 39.8%. The comparison histogram of the 51st prediction results and the actual results is shown in the figure below. In order to facilitate the comparison, the density value is uniformly increased by 1.

[0182] Figure 7 : is a comparison chart of the prediction results and actual results of the prediction interval of an embodiment of the present invention. Figure 7 As shown, based on the simulation results, the center position point of the fracture surface in the unit is simulated.

[0183] Figure 8 is the prediction result occurrence distribution extreme point diagram of this embodiment. Figure 8 As shown, the three-dimensional coordinate information of the fracture surface obtained above is combined with the Kriging interpolation method to generate the occurrence parameter information at the corresponding position.

[0184] Afterwards, the Kriging interpolation method was used to predict the trace length data of the unknown area, and the geometric parameters of the fracture surface in the unknown area were obtained by inverse calculation of the center point coordinates and the fracture surface diameter calculation method as shown in the following table.

[0185] Table 3 Prediction results of crack geometric parameters in unknown areas of tunnel

[0186]

[0187] Finally, the three-dimensional fracture network distribution model of the predicted interval is drawn by combining the unit fracture position, occurrence and diameter information of the unknown area obtained through the above prediction.

[0188] Fig. 9This is a network distribution diagram of the fracture structure surface in the prediction area of ​​this embodiment. In this diagram, Group 1 and Group 3 respectively represent the fracture distribution results obtained by using the K-means grouping method to divide the known data set into four groups, and then using the known data of the above groups to predict the prediction area. The reason why there are only Group 1 and Group 3 is that in this prediction, Group 2 and Group 4 are used for prediction, which shows that there are no fractures belonging to the above groups in the area. This prediction result reflects the development law of aggregation of the same occurrence categories in the actual fracture distribution.

[0189] The specific embodiments described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A rock mass fracture parameter prediction method for predicting rock mass fracture parameters in a target area, characterized in that: The target area includes a known parameter area and an unknown parameter area, and the method includes the following steps: S1. Determine the scale of the fracture unit using the first-order Markov property test theory, and divide the fracture units in the known parameter region according to the scale of the fracture unit; S2, using the first-order Markov chain theory to calculate the optimal transfer matrix for the crack density of the crack unit in the known parameter area divided in step S1; S3, using the optimal transfer matrix to predict the crack density of the unknown parameter area, and using a uniform distribution form to complete the distribution of trace midpoints and crack surface height data; S4. Use Kriging interpolation estimation method to complete the prediction of the fracture occurrence at the assigned location; S5. Use the Kriging interpolation estimation method to complete the prediction of the trace length of the crack at the assigned location.

2. The rock mass fracture parameter prediction method according to claim 1, characterized in that: Step S1 includes: Step S1.1: Calculate the adjacent order likelihood ratio test statistic η t-1,t , the calculation formula is as follows: Where T ij,…,kl and T j,…,kl The order of continuous t-step state transition is S i →S j →…→S k →S l The number of occurrences and the order of continuous t-1-step state transitions are S j →…→S k →S l The number of occurrences of Step S1.2: Calculate the likelihood ratio test statistic η of any order k,q , the calculation formula is as follows: or k,q =the k,k+1 +n k+1,k+2 +…+h q-1,q ; Step S1.3: Calculate the loss function, R(k), that is R(k)=η k,q -(m M -m k )(m-1); Step S1.4: Query the degrees of freedom (m q -m k The chi-square critical value corresponding to the 95% significance level under (m-1) is compared with the calculated R value. The k corresponding to the minimum R value is the most likely order of the Markov chain. Step S1.5: Determine whether the best order is 1st order. If not, change the size of unit L and repeat the above steps S1.1 to S1.

5.

3. The rock mass fracture parameter prediction method according to claim 1, characterized in that: The step S5 comprises: Step 5.1: Calculate the semivariance function value between known data sample points. The calculation formula is as follows: Where H is the step size between sample points, γ(H) represents the experimental semivariogram function, N(H) represents the logarithm of two sample points with a step size of H, s represents the sth pair of sample points, and x(s) represents the sample value of one of the sth pair of sample points; Step 5.2: The spherical model is used to fit the theoretical variogram. The spherical model formula is as follows: Where γ(h) represents the theoretical semivariogram function value between sample points with a step size of h, c0 represents the nugget value, c represents the sill value, a represents the range, and h represents the step size (h is a variable); Step 5.3: Construct the Kriging equation system, the formula is as follows: where γ ij represents the theoretical semivariance function value between the i-th sample point and the j-th sample point, γ i0 represents the theoretical semivariance function value between the i-th sample point and the predicted point x0, λ i is the weight coefficient, φ is the deviation value; Step 5.4: The formulas for calculating the Kriging estimate and Kriging variance are as follows: Among them, Z * (x0) is the Kriging estimate, σ E 2 is the Kriging variance.

4. A method for constructing a rock mass fracture model, characterized in that: The method comprises a rock mass fracture parameter prediction method as claimed in any one of claims 1 to 3.

5. The rock mass fracture model construction method according to claim 4, characterized in that: The following steps are also included: S6. Use the coordinate inverse calculation method to complete the xOy coordinate inverse calculation of the fracture center and the calculation of the diameter through the existing coordinates, occurrence and trace length information; S7. Establish a three-dimensional fracture network model of the target area based on the predicted rock fracture parameters.

6. The rock mass fracture model construction method according to claim 5, characterized in that: Step S6 includes: Step 6.1: Determine the inclination of the fracture surface and π / 2 relationship; Step 6.2: When the inclination When the inclination angle is equal to π / 2, the coordinates of the xOy plane of the trace are consistent with the coordinates of the center point of the crack surface; When π / 2 is not equal to π / 2, use the following formula to calculate: Where s2 is the inverse vector of the projection vector of the normal vector to the xOy plane, V UnitX and V UnitY represents the unit vectors of the x-axis and y-axis, H is the distance from the center of the crack surface to the xOy plane, is the inclination angle of the crack surface, and x0 and y0 represent the coordinates of the midpoint of the crack trace.

Images