CN104965969A - Inversion method for surrounding rock mechanical parameters of large cavern group - Google Patents

Abstract

The invention discloses an inversion method for surrounding rock mechanical parameters of a large cavern group. The method aims at solving a bottleneck problem that the calculation in the inversion for surrounding rock mechanical parameters of the large cavern group through a stochastic global optimization algorithm consumes too much time. In the method, minimum errors of monitoring displacement and calculation displacement of surrounding rocks are employed as object functions, the FLAC 3D numerical simulation technology is employed, and a backtracking and searching - information vector machine cooperation optimization algorithm provided by the invention is also employed, so that the inversion, which consumes too much time on numerical calculation by one time, on the optimal value for the surrounding rock mechanical parameters of the large cavern group can be rapidly carried out. In practical engineering applications, compared with the parameter inversion directly performed through the stochastic global optimization algorithm, the value simulation reanalysis frequency in a parameter inversion process can be significantly reduced on the premise that a global optimal solution is obtained through the inversion method, so that the time for the parameter inversion can be effectively reduced. The method is quite suitable for the inversion of the surrounding rock mechanical parameters of the large cavern group, is efficient and rapid, and is easy and useful.

Description

A kind of Large-scale Underground Tunnels And Chambers rock reaction force inversion method

Technical field

The invention belongs to Geotechnical Engineering field, relate to a kind of Large-scale Underground Tunnels And Chambers rock reaction force inversion method, refer to particularly, relate to a kind of Large-scale Underground Tunnels And Chambers rock reaction force inversion method based on backtrack search-information vector machine Cooperative Optimization Algorithm.

Background technology

The exploitation of current China hydraulic power potentials will mainly be distributed in western high-mountain gorge areas, along with the enforcement of Western Development Strategy, tens of large-scale hydroelectric projects enter the construction phase, wherein quite a lot of power station is all using crisscross Large-scale Underground Tunnels And Chambers as underground power house buildings, as Yalongjiang River Jinping hydropower station underground power house, adjust hole (222.6m × 12.0m × 67.5m) three large caverns to be the cavity group of main body with main building (335.20m × 25.20m × 61.85m, long × wide × high), transformer chamber (266.10m × 18.40m × 38.75m) and tail; Yellow River Laxiwa Hydroelectric Station main building excavation size reaches 312 × 30 × 74m; The major-minor factory building size of Xiluodu Hydropower Station on Jinsha River underground power house reaches 430.3m × 31.9m × 75.1m, and cavity group excavated volume is up to 1,700 ten thousand m3.The cavity group of underground power house is made up of multiple caverns such as main building, auxiliary power house, transformer chamber, down stream surge-chamber, diversion tunnel, tailrace tunnel, bus tunnel, access tunnels, and stressed very complicated, Stability Analysis of The Surrounding Rock must be undertaken by numerical value emulation method.

Current, there is the outstanding problem calculating large, rock reaction force consuming time and be difficult to determine in the numerical simulation of Large-scale Underground Tunnels And Chambers.The reliability of the selection logarithm value emulation of rock reaction force has material impact, due to the discreteness in space and scale effect remarkable, the accurate mechanics parameter of Large-scale Underground Tunnels And Chambers country rock is difficult to be determined by shop experiment.The rock reaction force inversion method of Large-scale Underground Tunnels And Chambers a kind ofly utilizes the effective ways of field monitoring information determination rock reaction force: parameter inversion problems is converted to mathematical unconstrained optimization problem, using the difference between displacement calculating value and displacement monitoring value as objective function, using rock reaction force as optimized variable, by certain Optimization Method optimization problem, obtain optimum solution.

Large-scale Underground Tunnels And Chambers has the advantages that cavern is many, build is complicated, adds the nonlinear mechanics character feature of country rock, makes to be optimized Back Analysis Problem when solving, the non-linear multi peak value of implicit expression of an objective function high complexity often.Adopt traditional gradient optimizing method to solve and often can only obtain locally optimal solution, in recent years, adopt the stochastic global optimization Algorithm for Solving Optimized Back-analysis problems such as genetic algorithm, particle swarm optimization algorithm, differential evolution algorithm to become a kind of trend.But, when adopting Stochastic Global Optimization Algorithms, often need invocation target function in large quantities, but due to the high complexity of huge underground cavity group structure, be difficult to the explicit expression set up between rock reaction force and displacement, often need by numerical evaluation to determine this implication relation, a large amount of invocation target functions just means and repeats a large amount of numerical evaluation, thus cause calculating too great high calculation cost problem consuming time, and then cause stochastic global optimization algorithm in the application of Large-scale Underground Tunnels And Chambers rock reaction force inverting, be subject to great restriction.The present invention is directed to the outstanding problem that stochastic global optimization algorithm application runs into when the inverting of Large-scale Underground Tunnels And Chambers rock reaction force, propose backtrack search-information vector machine Cooperative Optimization Algorithm (BSA-IVM), and in conjunction with FLAC 3D numerical value emulation method, propose the Large-scale Underground Tunnels And Chambers rock reaction force inversion method based on backtrack search-information vector machine Cooperative Optimization Algorithm, by reducing the numerical value weight analysis number of times of rock reaction force inverting as much as possible, thus reach the object that reduction calculating is consuming time, improve inverting efficiency.

FLAC 3D is a kind of numerical simulation software being applicable to simulate rock mass non-linear mechanical behavior, it adopts Explicit Lagrangian algorithm and mixing-discrete partition technology, can simulate plastic failure and the flowing of country rock material very exactly, be one of significant in value emulation technology of cavern's stability analysis at present.

Backtrack search optimized algorithm (Backtracking Search Optimization Algorithm, BSA) is a kind of new evolution algorithm, is proposed in 2013 by Pinar Civicioglu.Research shows, compares with the stochastic global optimization algorithm such as genetic algorithm, population, differential evolution algorithm, and BSA algorithm global optimizing ability is stronger, speed of convergence is faster, input parameter is less.The searching process of BSA algorithm is made up of global optimizing process and local optimal searching process, and wherein, the number of times of local optimal searching invocation of procedure objective function accounts for absolute majority ratio.If the objective function call number of local optimal searching process can be effectively reduced, solve calculating larger Large-scale Underground Tunnels And Chambers rock reaction force inverse problem consuming time by making this algorithm to be applicable to undoubtedly.

Machine learning is a new and developing branch of artificial intelligence, and it automatically finds rule from known embodiment, sets up the forecast model to unknown example, compared with traditional regression method, is more suitable for the regression problem of complicated nonlinearity.Artificial neural network and support vector machine are current representational machine learning methods.But artificial neural network and support vector machine all exist some open problems, such as, artificial neural network also exists optimal network topological structure and optimum hyper parameter is not easily determined, there is (owing) learns the problem such as risk, small sample Generalization Ability difference; Kernel function and the reasonable hyper parameter of support vector machine do not have feasible theoretical method for solving, are difficult to ensure forecasting reliability.

Information vector machine (Informative Vector Machine, IVM) is a kind of new machine learning method, is proposed in 2002 by Neil Lawrence.The method adopts the method based on information entropy theory, the most informational sample composition active set of part is optimized from a large amount of training samples, by the results of learning identical with former training sample set can be reached to the study of active set, represent in conjunction with rarefaction nuclear matrix simultaneously, thus greatly simplify time complexity and the space complexity of study.In addition, IVM by the screening of supposition density with minimize KL divergence (Relative Entropy) and achieve close approximation to non-gaussian distribution noise model Posterior distrbutionp.IVM has excellent recurrence performance, and its hyper parameter can self-adaptation obtain, and has stronger applicability to nonlinearity regression problem.In the inventive method, after BSA enters local optimal searching state, adopt the real goal function near IVM matching locally optimal solution, utilize IVM regression proxy model to have the key link that objective function returns:

(1) learning process of IVM regression proxy model

IVM regression proxy model is in learning process, and maintain two sample index collection I and J, wherein I is active set, and J treats selected works, time initial, j={1,2 ..., N}, and at any time, i ∪ J={1,2 ..., N} (assuming that from initial N training sample, will screen d information vector), information vector obtains in the mode of a kind of continuous print, similar on-line study: first, and application ADF is approximate has i information vector, i.e. I _itime Posterior distrbutionp and likelihood distribution: (for the situation of Gaussian distribution, approximate solution is consistent with Exact Solutions)

$\left\{\begin{array}{c}{q}_{I_i}\left(f\right)=N(f;{\mathrm{\μ}}_{I_i},{\mathrm{\Σ}}_{I_i})\≈p\left(f\right|X_{I_i:},y_{I_i},\mathrm{\θ})\\ p\left(y_{n_{I_i}}\right|f_{n_{I_i}})\≈N(m_{n_{I_i}};f_{n_{I_i}},{\mathrm{\β}}_{n_{I_i}}^{-1})\end{array}\right.---\left(1\right)$

In formula: p represents probability distribution, q represents APPROXIMATE DISTRIBUTION, and μ represents Gaussian distribution average, and Σ is covariance matrix, and m represents likelihood substitute variable, and β represents noise profile variance, for the input vector of effective training sample, θ represents covariance function hyper parameter.Afterwards, according to a following method choice i+1 information vector

$\left\{\begin{array}{c}{\mathrm{\ΔH}}_{I_i,j}=-\frac{1}{2}\log \left|{\mathrm{\Σ}}_{I_{i+1}}\right|+\frac{1}{2}\log \left|{\mathrm{\Σ}}_{I_i}\right|=-\frac{1}{2}\log \left|{\mathrm{\Σ}}_{I_{i+1}}{\mathrm{\Σ}}_{I_i}^{-1}\right|\\ \underset{j\∈J}{\arg }\max {\mathrm{\ΔH}}_{I_i,j}\end{array}\right.---\left(2\right)$

Above formula represents: select current treating can maximize the sample j reducing Posterior distrbutionp information entropy in selected works J, as the i-th+1 information vector.Circulation performs said process, until complete selection (the i.e. I=I of d information vector _d).Now, can obtain

$\left\{\begin{array}{c}p\left(y_I\right|X_{I,:},\mathrm{\θ})\≈N(m_I;0,K_I+B_I^{-1})\\ p\left(f\right|y_I,X_{I,:},\mathrm{\θ})\≈N(f;{\mathrm{\μ}}_I,{\mathrm{\Σ}}_I)\\ {\mathrm{\μ}}_I={\mathrm{\Σ}}_I{y}_I,{\mathrm{\Σ}}_I={(B_I+K_I^{-1})}^{-1}\end{array}\right.---\left(3\right)$

In formula, B represents noise profile variance, K or Σ represents Gaussian distribution covariance matrix.

In IVM regression proxy model, the optimum solution of covariance function hyper parameter θ is just by maximizing edge likelihood p (y _i| X _i:, θ) and self-adaptation obtains.Concrete, by getting negative logarithm-log (p (y _i| X _i:, θ)), maximization problems is converted into minimization problem, and then utilizes conjugate gradient decent to realize optimum hyper parameter self-adaptation obtain.

(2) regression process of IVM regression proxy model

Said process achieves the study substituting raw sample data collection with active set I, and regression process is afterwards consistent with the way in Bayesian regression learning process, substitutes into corresponding variable, vector or matrix, obtains IVM and returns Posterior distrbutionp:

$\left\{\begin{array}{c}p\left(f_{*}\right|y,X,x_{*},{\mathrm{\θ}}_l)\≈p\left(f_{*}\right|y_I,X_{I,:},x_{*},{\mathrm{\θ}}_{l_I})=N(f_{*};{\mathrm{\μ}}_{*},{{\mathrm{\σ}}_{*}}^2)\\ {\mathrm{\μ}}_{*}=K_{*I}^T{K}_I^{-1}{\mathrm{\Σ}}_I{B}_I{y}_I,{\mathrm{\σ}}_{*}^2=k_{*}+K_{*I}^T{K}_I^{-1}({\mathrm{\Σ}}_I-K_I)K_I^{-1}{K}_{*I}\end{array}\right.---\left(4\right)$

Formula (4) implies individual position coordinates x _*with fitness f _*corresponding relation, in order to replace real fitness function curve.

Newton method is a kind of efficient classical Local Optimization Algorithm, be applicable to the local extremum problem solving nonlinear function, its ultimate principle utilizes quafric curve to carry out pointwise to be similar to former objective function, and the minimum point approaching gradually carrying out approximate former objective function with the minimum point of quafric curve is changed the time.The inventive method for objective function, adopts Newton method to carry out high efficiency local optimal searching with IVM agent model.

The ultimate principle of BSA-IVM algorithm is: when after the certain step number of BSA algorithm iteration, IVM is utilized to carry out performance matching to objective function in the local neighborhood of current optimum individual, obtain the approximate expression of former objective function, thus by former objective function explicitization, obtain single order and second derivative information thus, and then utilize Newton method to inquire into locally optimal solution fast, thus realize accelerating BSA algorithm local optimal searching process, to reduce the object of numerical value weight analysis.

Summary of the invention

In order to overcome the defect existed in prior art, the present invention proposes a kind of Large-scale Underground Tunnels And Chambers rock reaction force inversion method.

Its technical scheme is as follows:

A kind of Large-scale Underground Tunnels And Chambers rock reaction force inversion method, comprises the following steps:

(1) according to the constitutive model of geology survey report or indoor and outdoor rock mechanics experiment result determination cavity group country rock, the FLAC 3D mathematical calculation model of cavity group country rock is set up.

(2) parameter inversion problems is converted into optimization problem, sets up the objective function of inverting optimization problem: wherein, x is one group of Mechanics Parameters of Rock Mass, d _ix () is the measured displacements d of the monitoring point for displacement _i(x), it is the displacement calculating of the monitoring point for displacement.Objective function is less, and displacement calculating is more close to measured displacements, and the confidence level of the result of calculation of corresponding FLAC3D numerical model is higher.

(3) adopt backtrack search-information vector machine Cooperative Optimization Algorithm (BSA-IVM) to solve objective function, Optimization Steps is as follows:

1. algorithm parameter is arranged: according to the number determination population number NP treating inverting rock mass parameter, the condition of convergence of set algorithm, comprises objective function minimum value ε and maximum permission iterative steps T _max;

2. stochastic generation experimental population P _ijand oldP _ij, wherein i is population scale, and j is the number of the Rock And Soil parameter treating inverting, and the individuality of two populations is all randomly distributed in optimizing region;

3. to experimental population P _ijcarry out fitness evaluation, obtain target function value E (i) of all individualities, then individual optimal value E _f(i)=E (i), current P _ijthe particle P that middle target function value is minimum _gjfor current global optimum's particle, target function value E (g) of its correspondence is current globally optimal solution, now iterations t=1;

4. as E (g) < ε, t<T _max, then step is below performed; Otherwise, the parameter value of output recover;

5. enter circulation and carry out global optimizing state, and record all individual informations.

6. P is used _ijrandom replacement oldP _ij: generate the random number a between (0,1) and b, as a<b, use P _ijreplace oldP _ij, otherwise do not replace;

7. by oldP _ijin the order random alignment again of individuality, generate new individual population oldP1 _ij;

8. to initial population P _ijmake a variation, generate the population T after variation _ij, the formula that wherein makes a variation is T _ij=P _ij+ F. (oldP1 _ij-P _ij), F is constant, is used for gating matrix (oldP1 _ij-P _ij) amplitude;

9. to the population T after variation _ijcarry out hybridization to calculate, obtain the matrix T 1 after hybridizing _ij: generate the i*j be made up of " 0 " and " 1 " and tie up matrix maP _ij.MaP _ijbe used for controlling population T _ijin individuality will by original seed group P _ijin correspondence individuality replace position, i.e. maP _ijmiddle all values is the position of " 0 ", T _ijin the individuality of these positions will by initial population P _ijthe individuality of middle correspondence position is replaced.Generate maP _ijmethod be introduce parameter composite rate mixrate, control the individual number that will be replaced by composite rate;

10. T1 is obtained by evaluation _ijin target function value H (i) of all individualities;

as H (k) <E (k), namely individual k evolves rear than original more excellent, then upgrade individual optimal value E _fk ()=H (k), obtains new ideal adaptation angle value set E _fi (), upgrades personal best particle P simultaneously _kj=T1 _kj;

upgrade global optimum individual: with the minimum value E of the fitness value when former generation individuality _f(d)=min (E _j(i)) compare with previous generation global optimum individuality E (g), work as E _j(d) <E (g), then E (g)=E _f(d);

after the certain number of times of global optimizing iteration, from all individualities evolutionary history, choose Euclidean distance individual as learning sample from current optimum individual g nearest m, set up new learning sample collection (X _sample, Y _sample);

log-on message vector machine (IVM) carries out the IVM agent model learning to obtain true fitness function to learning sample;

differentiate is carried out to IVM agent model, obtains the first order derivative for current optimum sample point and second derivative;

newton method is adopted to utilize derivative information to carry out effective search to local optimum individual; Local optimum individuality is substituted into FLAC3D program once just calculate, obtain the real function fitness value of optimum individual; The prediction true fitness function value of optimum individual and the real function fitness value of current optimum individual are compared, if be better than current individual, then replace current optimum individual with prediction optimum individual, and judge whether to meet convergence of algorithm criterion, if meet, then termination of iterations; Otherwise, return step 5.;

(4) if objective function reaches the aimed at precision requirement of setting, then stop calculating, the parameter of output recover; Otherwise, continue to get back to step (1), carry out new round calculating, constantly repeatedly, until objective function arrives aimed at precision.

Preferably, the amplitude controlling constant F=3*randn that described Mutation Strategy is introduced;

Preferably, the hybrid rate mixrate=1 that described Hybridization Strategy is introduced;

Preferably, described ideal adaptation angle value method for solving is: adopt business mathematics software MATLAB and business numerical evaluation software FLAC 3D to carry out joint inversion, in MATLAB environment, the individual position coordinates (one group of rock reaction force) of BSA algorithm is saved in data-interface file A, start FLAC 3D numerical evaluation software by self-defined software transfer order and enter duty, and call one group of rock reaction force that the embedded FISH program of FLAC 3D reads data-interface file A, substitute into the FLAC 3D numerical model set up and obtain displacement calculating and then acquisition target function value; Then, target function value is stored in data file B, by the target function value in MATLAB program file reading B, obtains the fitness value of this individuality thus.

Beneficial effect of the present invention is:

(1) FLAC 3D program of numerical calculation embeds in BSA global optimization approach as a separate modular by the present invention, the efficient kernel solver of business numerical simulation software auxiliary under, significantly can accelerate the carrying out of back analysis, provide very big convenience to user of the present invention.

(2) the present invention is by global optimization approach, Local Optimization Algorithm, machine learning method calculates with numerical simulation and combines, give full play to the characteristic and advantage of various method, achieve and significantly reduce numerical value weight analysis number of times in rock reaction force refutation process under the precondition obtaining globally optimal solution, compared with traditional inversion method, can the more efficient reasonable value obtaining rock reaction force in Large-scale Underground Tunnels And Chambers simulation calculation rapidly, calculate long problem consuming time provide an effective way for effectively solving in the inverting of current Large-scale Underground Tunnels And Chambers rock reaction force.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of the inventive method;

Fig. 2 be BSA-IVM algorithm centered by current optimum individual, select learning sample schematic diagram;

Fig. 3 is the FLAC 3D illustraton of model of hydraulic tunnel;

Fig. 4 is the point layout figure of hydraulic tunnel;

Fig. 5 is underground power house excavation figure by stages;

Fig. 6 is the M7 point layout figure of 00+15.00 section;

Fig. 7 is the displacement monitoring curve of M7 measuring point;

Fig. 8 is the arrangenent diagram of M19 and the M24 measuring point of 00+54.00 section;

Fig. 9 is the displacement monitoring curve of M19 measuring point;

Figure 10 is the displacement monitoring curve of M24 measuring point;

Figure 11 is the comparison diagram of displacement calculating and measured displacements.

Embodiment

Below in conjunction with the drawings and specific embodiments, technical scheme of the present invention is described in more detail.

With reference to Fig. 1-Fig. 2, the hydraulic tunnel diameter of certain engineering is 3.0m, average lithology, and level and vertical initial field stress are 10MP.Assuming that the constitutive model of this model is Morh-Coulomb ideal elastoplastic model, it treats that the Mechanics Parameters of Rock Mass of inverting is elastic modulus E, cohesive strength c and angle of internal friction to choose in excavation face 1,2,3,4,5 measuring points as displacement monitoring point.In order to the feasibility of verification algorithm, assuming that the true mechanics parameter of country rock is E=2.1GPa, c=1.25MPa, substituted into displacement that FLAC 3D program just obtaining as " measured displacements ".

As shown in Figure 3, Figure 4, the constitutive model of this hydraulic tunnel is Morh-Coulomb ideal elastoplastic model, sets up the mathematical calculation model of tunnel with FLAC3D, sets up the objective function that inverting is optimized

Be E=2.1GPa, c=1.25MPa by true for the country rock of supposition mechanics parameter, the displacement that substitution FLAC 3D program is just obtaining is as " displacement monitoring " d _i(x);

Start MATLAB software, adopt particle cluster algorithm (PSO), backtrack search algorithm (BSA), backtrack search-information vector machine Cooperative Optimization Algorithm (BSA-IVM) to solve objective function respectively, Optimization Steps is as follows:

1. the optimum configurations of initialization algorithm: PSO algorithm parameter: NP=80, c ₁=c ₂=2.0, Vmax=[1,1,1]; BSA algorithm parameter: NP=80; BSA-IVM algorithm parameter: NP=80; IVM machine learning is started after BSA algorithm iteration 7 step, learning sample in the individual neighborhood of the local optimum chosen adds up to 3 × 80=240, IVM machine learning selects 100 as information vector at every turn in 240 learning samples, and continuous circulation selects 8 times, 100 information vectors obtaining " useful " are carried out study and are returned.The convergence criterion of the objective function of three kinds of algorithms is ε=1 × 10 ^-3;

2. arrange the scope of the individual search volume of algorithm, this scope is exactly the span (see table 1) of the rock reaction force treating inverting;

The region of search of table 1 rock reaction force

Rock mass parameter	Play mould E (GPa)	Cohesive strength c (MPa)	Angle of internal friction Φ (°)
				Interval	2.0-2.5	1.0-1.4	27-32

3. starting guide algorithm carries out evolutionary operation, after individuality evolution strategy is algorithmically evolved, the reposition (80 groups of rock reaction force) of individuality is deposited data-interface file A;

4. start FLAC 3D numerical evaluation software by self-defined software transfer order, by the rock reaction force in the FISH program fetch interface file A that FLAC 3D is embedded, and substituted in set up FLAC 3D numerical model, obtain displacement calculating thus will " displacement monitoring " d _ix () substitutes into objective function the target function value that middle acquisition 80 groups of rock reaction force are corresponding, and target function value is stored in data file B;

5. the target function value in MATLAB program read data files B is used, obtain the fitness value of all individualities, and select minimum fitness value by comparing, the corresponding individual position coordinates of minimum fitness value i.e. generation individual corresponding optimum Analysis of Field Geotechnical Parameters combination for this reason;

6. by these minimum fitness value and convergence criterion ε=1 × 10 ^-3compare, if be less than convergence criterion, then stop calculating, export individual position coordinates (i.e. optimum rock reaction force combination); Otherwise continue the optimizing of a new round, until reach convergence criterion;

The parametric inversion of three kinds of optimized algorithms the results are shown in Table 2, as can be seen from Table 2, compares with PSO with BSA, and the result of BSA-IVM algorithm parameter inverting is more close to actual value;

The inversion result of table 2 rock reaction force

The cost that algorithms of different carries out parametric inversion is listed in table 3 with calculating is consuming time, the calculation cost of BSA-IVM is consuming time far below BSA and PSO with calculating as can be seen from Table 3, shows that BSA-IVM Cooperative Optimization Algorithm has significant superiority when carrying out the inverting of Large-scale Underground Tunnels And Chambers rock reaction force.

Table 3 three kinds of algorithms calculate comparison consuming time

Algorithm	PSO	BSA	BSA-IVM
				FLAC 3D call number	47200	4960	1004
Calculate (s) consuming time	2.86×10 4	2.18×10 3	4.12×10 2

Engineering example

The underground powerhouse of some hydropower station is formed primarily of a few part such as main building, auxiliary power house, main transformer hole, bus tunnel, lock chamber hole.Underground power house main building arrangement is in massif, and the thick about 210 ~ 240m of top covering rockmass, it is the vein invaded in fresh migmatitic granite and a small amount of later stage that factory building excavation discloses country rock, and rock is hard, and representational vein is N π 24 Diabase dyke and δ π 201 diorite arteries and veins.According to geological record and topographic data, carry out classification by " Standard for classification of engineering rock masses " (GB 50218-94) to factory building rock mass, its gross index (BQ) is much in 450, and rank mostly is I ~ II grade.Carry out engineering of surrounding rock geological classes by " code for geological investigation of water resources and hydropower engineering " (GB 50287-99), country rock mostly is II ~ III class.Main transformer cavity parallel is arranged in factory building downstream, the thick 35m of the dike between factory building and main transformer hole, top covering rockmass thickness 188 ~ 220m.Country rock is the vein invaded in fresh migmatitic granite and a small amount of later stage, and rock is hard, and rock mass is much more complete, and representational vein is δ π 202 diorite arteries and veins.Main transformer hole carries out classification by " Standard for classification of engineering rock masses " (GB 50218-94) to rock mass, and its gross index (BQ) is much in 450, and rank mostly is I ~ II grade, and only local fault shatter belt and joint density area are III ~ IV grade.Carry out engineering of surrounding rock geological classes by " code for geological investigation of water resources and hydropower engineering " (GB-50287-99), country rock mostly is II ~ III class.Through statistics, II class surrounding rock accounts for 39.4%, III class surrounding rock and accounts for 60.6%.Bus tunnel totally four, long 35m, adopts large duck eye to be connected scheme.Bus tunnel country rock is that fresh mix grouan 1# bus tunnel Rock Mass Integrality is poor, and country rock is III class rock, and 2# ~ 4# bus tunnel rock mass is more complete, and country rock is II class rock.The sectional view of cavity group stage excavation is shown in Fig. 6.

The measured displacements choosing M19, M24 tri-measuring points of M7 and the 0+54.00 section of 0+15.00 section carries out the parametric inversion of second phase excavation (three layers with four layers), and the concrete layout of three measuring points and actual displacement monitoring situation are shown in Fig. 7-Figure 11.The excavation part of FLAC 3D numerical simulation model is shown in Fig. 5, country rock constitutive model is set to Morh-Coulomb ideal elastic-plastic constitutive model, needs the rock reaction force of inverting to have: elastic modulus of surrounding rocks E1,1# unit section main building downstream, main building upstream elastic modulus of surrounding rocks E2,2#-4# unit section main building downstream elastic modulus of surrounding rocks E3, transformer chamber downstream elastic modulus of surrounding rocks E4, main building upstream country rock cohesive strength C1, main building downstream country rock cohesive strength C2, country rock angle of internal friction parameter value interval is in table 4.

The region of search of table 4 underground rock cavern mechanics parameter inverting

Parameter	E1(GPa)	E2(GPa)	E3(GPa)	E4(GPa)	C1(GPa)	C2(GPa)	Φ(°)
								The region of search	60-80	10-30	20-40	100-300	5-7	5-7	40-60

The optimum configurations of BSA algorithm and BSA – IVM algorithm: population scale NP=10, meets f<1 × 10 with objective function ^-2as iteration convergence criterion, other optimum configurations of algorithm is with the example in " method validation ".

The step of parametric inversion by specification is implemented.The optimal value of each parameter is in table 5, and the inverting contrast consuming time of two kinds of inversion algorithms is in table 6.

Table 5 is based on the optimal value of each parameter of BSA – IVM algorithm inverting

Parameter	E1(GPa)	E2(GPa)	E3(GPa)	E4(GPa)	C1(GPa)	C2(GPa)	Φ(°)
								The region of search	65.3	21.6	28.9	215.1	6.7	6.8	56.0

The calculating comparison consuming time of table 6 two kinds of inversion methods

Method	BSA	BSA-IVM
			FLAC 3D call number	560	84
Calculate (s) consuming time	1.4×10 4	1.9×10 3

In order to verify the accuracy of the optimized parameter of inverting, the optimized parameter above-mentioned inverting obtained substitutes into FLAC 3D program of numerical calculation and just calculates, and Figure 10 is shown in the displacement calculating value of each monitoring point and the contrast of measured displacements value that participate in inverting.

As seen from Figure 11, M7 measuring point displacement calculating value is slightly smaller than displacement monitoring value, and relative error is 2.39%; The displacement calculating value of M19 measuring point is a bit larger tham displacement monitoring value, and relative error is 4.53%; The displacement calculating value of M24 measuring point and displacement monitoring value are substantially suitable, and relative error is only 0.35%; The average relative error of three measuring points is 2.42%.Relative error and the average relative error of each measuring point are all less than 5%, indicate the rationality of inversion result thus.

Engineering example result shows, the inventive method is feasible, there is the advantage that global optimizing is strong, inverting efficiency is high, can overcome and when adopting stochastic global optimization algorithm to carry out inverting, calculate excessive, the deficiency that adopts traditional gradient optimal method inversion accuracy poor consuming time, logarithm value calculates larger Large-scale Underground Tunnels And Chambers rock reaction force back analysis consuming time and has stronger applicability, has good future in engineering applications.

The above; be only the present invention's preferably embodiment; protection scope of the present invention is not limited thereto; anyly be familiar with those skilled in the art in the technical scope that the present invention discloses, the simple change of the technical scheme that can obtain apparently or equivalence are replaced and are all fallen within the scope of protection of the present invention.

Claims (4)

1. a Large-scale Underground Tunnels And Chambers rock reaction force inversion method, is characterized in that, comprises the following steps:

(1) according to the constitutive model of geology survey report or indoor and outdoor rock mechanics experiment result determination cavity group country rock, the FLAC 3D mathematical calculation model of cavity group country rock is set up;

(2) parameter inversion problems is converted into optimization problem, sets up the objective function of inverting optimization problem: wherein, x is one group of Mechanics Parameters of Rock Mass, d _ix () is the measured displacements of i-th monitoring point for displacement, it is the displacement calculating of i-th monitoring point for displacement; Objective function is less, and displacement calculating is more close to measured displacements, and the confidence level of results of numerical model calculation is higher;

(3) adopt backtrack search-information vector machine Cooperative Optimization Algorithm (BSA-IVM) to solve objective function, Optimization Steps is as follows:

1. algorithm parameter is arranged: according to the number determination population number NP treating inverting rock mass parameter, the condition of convergence of set algorithm;

2. stochastic generation two experimental population P _ijand oldP _ij, wherein i is population scale, and j is the number of the rock mass parameter treating inverting, and the individuality of two populations is all randomly distributed in optimizing region;

3. to experimental population P _ijcarry out fitness evaluation, obtain target function value E (i) of all individualities, then individual optimal value E _f(i)=E (i), current P _ijthe particle P that middle target function value is minimum _gjfor current global optimum's particle, target function value E (g) of its correspondence is current globally optimal solution, now iterations t=1;

4. as E (g) < ε, t<T _max, then step is below performed; Otherwise, the parameter value of output recover;

5. enter circulation and carry out global optimizing state, and record all individual informations;

6. P is used _ijrandom replacement oldP _ij: generate the random number a between (0,1) and b, as a<b, use P _ijreplace oldP _ij, otherwise do not replace;

7. by oldP _ijin the order random alignment again of individuality, generate new individual population oldP1 _ij;

8. to initial population P _ijmake a variation, generate the population T after variation, the formula that wherein makes a variation is T _ij=P _ij+ F. (oldP1 _ij-P _ij), F is constant, is used for gating matrix (oldP1 _ij-P _ij) amplitude;

9. to the population T after variation _ijcarry out hybridization to calculate, obtain the matrix T 1 after hybridizing _ij: generate the i*j be made up of " 0 " and " 1 " and tie up matrix maP _ij; maP _ijbe used for controlling population T _ijin individuality will by original seed group P _ijin correspondence individuality replace position, i.e. maP _ijmiddle all values is the position of " 0 ", T _ijin the individuality of these positions will by initial population P _ijthe individuality of middle correspondence position is replaced; Generate maP _ijmethod be introduce parameter composite rate mixrate, control the individual number that will be replaced by composite rate;

10. T1 is obtained by evaluation _ijin target function value H (i) of all individualities;

as H (k) <E (k), namely individual k evolves rear than original more excellent, then upgrade individual optimal value E _fk ()=H (k), obtains new ideal adaptation angle value set E _fi (), upgrades personal best particle P simultaneously _kj=T1 _kj;

upgrade global optimum individual: with the minimum value E of the fitness value when former generation individuality _f(d)=min (E _j(i)) compare with previous generation global optimum individuality E (g), work as E _j(d) <E (g), then E (g)=E _f(d);

after the certain number of times of global optimizing iteration, enter local optimal searching state; From all individualities evolutionary history, choose Euclidean distance individual as learning sample from current optimum individual g nearest m, set up new learning sample collection (X _sample, Y _sample);

log-on message vector machine (IVM) learns learning sample, obtains the IVM agent model of true fitness function;

differentiate is carried out to IVM agent model, obtains the first order derivative for current optimum sample point and second derivative;

adopt Newton method to utilize derivative information to carry out displacement to sample point, obtain forecast sample array, choose the minimum point of fitness function value in forecast sample array as prediction optimum individual; Prediction optimum individual substitutes into FLAC 3D program and once just calculates, and obtains the real function fitness value of prediction optimum individual; The prediction true fitness function value of optimum individual and the real function fitness value of current optimum individual are compared, if be better than current individual, then replace current optimum individual with prediction optimum individual, and judge whether to meet convergence of algorithm criterion, if meet, then termination of iterations; Otherwise, return step 5.;

(4) if objective function reaches the aimed at precision requirement of setting, then stop calculating, the parameter of output recover; Otherwise, continue to get back to step (1), carry out new round calculating, constantly repeatedly, until objective function arrives aimed at precision.

2. huge underground cavity rock reaction force inversion method according to claim 1, is characterized in that, preferably, and the amplitude controlling constant F=3*randn that described Mutation Strategy is introduced.

3. huge underground cavity rock reaction force inversion method according to claim 1, is characterized in that, the hybrid rate mixrate=1 that described Hybridization Strategy is introduced.

4. Large-scale Underground Tunnels And Chambers rock reaction force inversion method according to claim 1, it is characterized in that, described ideal adaptation angle value method for solving is: adopt business mathematics software MATLAB and business numerical evaluation software FLAC 3D to carry out joint inversion, in MATLAB environment, the individual position coordinates of BSA algorithm is saved in data-interface file A, start FLAC 3D numerical evaluation software by self-defined software transfer order and enter duty, and call one group of rock reaction force that the embedded FISH program of FLAC 3D reads data-interface file A, substitute into the FLAC 3D numerical model set up and obtain displacement calculating and then acquisition target function value, then, target function value is stored in data file B, by the target function value in MATLAB program file reading B, obtains the fitness value of this individuality thus.