Tunnel surrounding rock stability evaluation method and system based on multi-physical field parameter inversion

Abstract

The present application relates to the technical field of construction surrounding rock stability evaluation, and discloses a tunnel surrounding rock stability evaluation method and system based on multi-physical field parameter inversion, aiming to solve the problems of insufficient data synergy, strong parameter inversion multi-solution, and poor accuracy and timeliness of surrounding rock stability evaluation of the existing method, and the scheme mainly comprises: arranging an acousto-optic-electromagnetic seismic drill multi-modal sensing array and implementing time-space synchronization; establishing an interaction entropy spectrum density model to realize multi-physical field collaborative excitation and collection; obtaining a unified feature vector through data correction, feature extraction and weighted fusion; constructing a joint inversion objective function embedded with rock physical constraints, inverting to obtain a three-dimensional physical property parameter field, and combining with acoustic emission energy to calculate a dynamic permeability field; finally, using ensemble Kalman filtering to dynamically update the model, and based on the updated parameters and seepage-uncertainty coupling correction to obtain the final risk probability. The present application improves the accuracy, real-time performance and reliability of deep-buried tunnel surrounding rock stability evaluation.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of construction risk assessment, and particularly relates to a tunnel surrounding rock stability evaluation method and system based on multi-physical field parameter inversion. BACKGROUND

[0002] As a key link of infrastructure construction, the construction safety of deep-buried tunnel engineering highly depends on the accurate understanding of the surrounding rock structure and stability. In current engineering practice, the surrounding rock state evaluation mainly relies on single physical field testing technology in stages and methods. These technologies include: arranging an array of seismic detectors on the tunnel wall, exciting and receiving seismic wave signals through an active seismic source, and inverting the wave velocity field distribution of the surrounding rock; arranging power electrodes and measurement electrodes at a specific section, and obtaining rock resistivity parameters through direct current method or induced polarization method; installing an optical fiber sensing system in the borehole to monitor local strain and temperature changes; and collecting mechanical parameters such as drilling pressure, rotation speed and drilling speed through a while-drilling recording system.

[0003] These traditional methods have significant technical limitations in the implementation process: first, the layout of various sensing systems is often carried out independently, lacking a unified space-time reference calibration. Seismic sensors, electrical method electrodes, optical fiber measuring points, etc. respectively adopt their own coordinate systems and time sampling schemes, which leads to difficulties in accurately registering different physical field data in space and strictly synchronizing in time, forming a data island effect. Secondly, the inversion interpretation based on a single physical field has inherent multi-solution problems. For example, the observed wave velocity anomaly may be caused by changes in rock mass integrity or pore pressure changes, and resistivity anomalies may reflect differences in mineral composition or water content changes, and it is difficult to make accurate geological interpretation with only single data.

[0004] To alleviate the problem of insufficient single physical field information, the existing technology proposes the idea of multi-source data fusion. The typical approach is to use the "abnormal level" or "evaluation index" of seismic, electrical, radar or while-drilling parameters as input, and use methods such as analytic hierarchy process (AHP), evidence theory (D-S), fuzzy comprehensive judgment, cloud model to construct a comprehensive evaluation index to classify the surrounding rock stability or water gushing risk. However, this kind of method generally only empirically weights or grades the multi-source data in the evaluation stage, and does not establish a unified constraint mechanism in the data acquisition scheduling and physical property parameter inversion stage. Multi-physical field data are still independent of each other in the acquisition process, and the mutual interference effect cannot be quantified and fed back to the excitation and acquisition strategy, and the data weight is also given through experience or subjective judgment, which is difficult to improve the consistency of multi-source data and the stability of inversion results from the source.

[0005] In the aspect of characterizing the seepage properties of surrounding rock, the existing researches mostly regard the permeability as a function of strain or fracture geometry parameters, adopt empirical models, or indirectly depict the evolution of seepage channels through indicators such as specific surface area and joint density. Such models can reflect the change of permeability under static loading or slow deformation to some extent, but they often do not explicitly introduce dynamic rupture indicators such as strain rate and acoustic emission energy, making it difficult to depict the sudden increase of seepage channels in the process of rapid excavation, periodic unloading and rapid penetration of micro-cracks in deep tunnel surrounding rock. Especially, when the surrounding rock enters the stage of strong disturbance, the process of micro-crack initiation, expansion and penetration is accompanied by significant changes in strain rate and acoustic emission energy release. If only static models are used, they often cannot sensitively respond to such nonlinear seepage amplification process, leading to insufficient evaluation of the formation and evolution of gushing channels.

[0006] In the aspect of evaluating the stability of surrounding rock, the existing technology mostly adopts empirical judgment methods based on single parameter threshold, or introduces fuzzy comprehensive judgment, gray evaluation, cloud model and other uncertainty analysis methods on this basis, and divides the risk of surrounding rock into several grades. Although the above methods consider parameter fluctuations and cognitive uncertainty to some extent, they usually reflect uncertainty at the evaluation index level through "confidence", "membership" or "fuzzy weight", and do not explicitly map model uncertainty to the risk probability function itself. In other words, uncertainty exists in the form of qualitative grades or additional weights, and it is difficult to dynamically correct the risk probability according to the uncertainty distribution of the physical parameter field and the evolution degree of the seepage channel, so there is a risk that major gushing events in areas with high uncertainty and high seepage may be underestimated or missed.

[0007] In addition, the existing risk evaluation method of deep tunnel surrounding rock also has shortcomings in dynamic updating mechanism. Most methods are based on one-time inversion or stage inversion results, lack an effective way to assimilate the real-time drilling parameters, acoustic emission events and optical fiber strain monitoring data collected during construction into the three-dimensional surrounding rock physical model, and are difficult to correct the time-varying deviation of the physical parameter field in time, resulting in that the risk evaluation results lag behind the actual state change of the surrounding rock.

[0008] In summary, the existing technology has obvious shortcomings in multi-source collaborative acquisition, physical constraint joint inversion, dynamic evolution of permeability and quantitative transmission of uncertainty to risk probability, and it is difficult to meet the evaluation needs of high precision, strong coupling and dynamic updating of gushing risk of surrounding rock during deep tunnel construction. SUMMARY

[0009] The present application aims to solve the problems of insufficient data collaboration, strong multi-solution of parameter inversion and poor accuracy and timeliness of surrounding rock stability evaluation in the existing surrounding rock state evaluation method, and proposes a tunnel surrounding rock stability evaluation method and system based on multi-physical field parameter inversion.

[0010] The technical solution adopted by the present application to solve the above technical problems is:

[0011] In a first aspect, a tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion is provided, and the method comprises:

[0012] Step 1, in a deep-buried tunnel, a multi-modal sensor array including six types of sensors of sound, light, electricity, magnetism, shock and drilling is arranged, and the multi-modal sensor array is subjected to spatial coordinate calibration and time synchronization to form original multi-physical field signals with unified time and space reference;

[0013] Step 2, a multi-field mutual interference entropy spectrum density model for quantifying the interference degree between different physical field signals is established, the mutual interference degree of each physical field combination is quantitatively characterized based on the multi-field mutual interference entropy spectrum density model and an adaptive scheduling strategy is generated, the excitation source of each physical field is controlled to perform time-sharing or parallel cooperative excitation according to the adaptive scheduling strategy, and the response signals of the multi-modal sensor array are synchronously collected to obtain a multi-source cooperative collection data set;

[0014] Step 3, the multi-source cooperative collection data set is subjected to data quality verification and space-time consistency correction, physical field characteristics including seismic wave first arrival time, apparent resistivity, Cole-Cole polarization parameters, GPR reflection interface depth, optical fiber strain, acoustic emission energy and normalized drilling speed index are extracted from the corrected data, a multi-modal characteristic vector including the characteristic vectors corresponding to the physical field characteristics is constructed, the reliability weight of each physical field characteristic is calculated based on the mutual interference entropy and observation consistency residual, the reliability weight corresponding to each physical field characteristic is formed into a data weight matrix, and each multi-modal characteristic vector is weighted and fused to obtain a unified fusion characteristic vector;

[0015] Step 4, according to the unified fusion characteristic vector and the data weight matrix, a joint inversion objective function embedded with Gassmann equation and Archie formula is constructed, the data weight matrix is applied to the weighted residual term of the joint inversion objective function, the joint inversion objective function is solved through Gauss-Newton iterative algorithm, a three-dimensional surrounding rock physical parameter field including wave velocity, resistivity, porosity and strain field is obtained, a permeability dynamic evolution function with strain, strain rate and acoustic emission energy as independent variables is constructed based on the strain field in the three-dimensional surrounding rock physical parameter field and the acoustic emission energy extracted in step 3, the initial permeability in the three-dimensional surrounding rock physical parameter field is subjected to space-time dynamic updating, a dynamic permeability field is calculated, and finally a three-dimensional surrounding rock physical model including wave velocity, resistivity, porosity and permeability is formed;

[0016] Step 5: Using the three-dimensional surrounding rock physical property model as the initial background field, the ensemble Kalman filter algorithm is used to assimilate the drilling parameters, acoustic emission events, and fiber optic strain monitoring data collected in real time during construction into the three-dimensional surrounding rock physical property model, and the physical property parameter field in the three-dimensional surrounding rock physical property model is dynamically updated. Based on the updated physical property parameters, Logistic regression is used to calculate the basic risk probability at each spatial location and time. Based on the updated physical property parameter field, a physical property parameter uncertainty tensor characterizing the uncertainty of the resistivity, porosity, and permeability models is constructed. The Frobenius norm of the physical property parameter uncertainty tensor is calculated, and the seepage-uncertainty coupling correction coefficient is determined by combining the relative increment between the dynamic permeability and the initial permeability. The seepage-uncertainty coupling correction coefficient is used as a multiplicative factor on the basic risk probability to perform spatiotemporal dynamic correction on the basic risk probability, and the final risk probability is obtained. The surrounding rock stability is evaluated based on the final risk probability value.

[0017] Furthermore, in step 2, the multi-field interfering entropy spectral density model quantifies the degree of interference between different physical field signals using the following formula:

[0018] ;

[0019] in, Indicates the first The physical field and the first A physical field signal in time ,frequency The normalized cross-interference energy ratio, , Indicates the first The physical field and the first A physical field signal in time ,frequency The cross-correlation spectral density, Indicates the first A physical field in time ,frequency The mutual perturbation entropy is such that the larger the value of the mutual perturbation entropy, the greater the degree of interference from other physical fields.

[0020] Furthermore, in step 3, physical field features are extracted from the corrected data, which is achieved through the following formula:

[0021] The first arrival time of the seismic waves was automatically identified using the AIC criterion.

[0022] ;

[0023] in, This indicates the possible first arrival time of the seismic wave, i.e., the first arrival time in the data sequence. One sampling point, AIC criterion value calculated when assuming the first arrival time as the

[0024] The calculation formula of apparent resistivity is as follows:

[0025]

[0026] Wherein, represents apparent resistivity, represents device coefficient, represents measured potential difference, represents injected current;

[0027] The Cole-Cole polarization parameters include time constant, which is obtained by fitting the Cole-Cole model to the complex resistivity dispersion curve:

[0028]

[0029] Wherein, represents complex resistivity, represents zero frequency resistivity, represents charging rate, represents time constant, represents frequency correlation coefficient, represents imaginary unit, represents angular frequency;

[0030] The GPR reflection interface depth is calculated by the time-depth conversion model:

[0031] ;

[0032] Wherein, represents GPR reflection interface depth, represents propagation speed of electromagnetic wave in rock mass, represents two-way propagation time of electromagnetic wave, represents speed of light in vacuum, represents relative dielectric constant of rock mass;

[0033] The optical fiber strain is calculated by the phase-strain conversion formula:

[0034] ​​​​​​​​​​​

[0035] in, This indicates the amount of phase change measured. Indicates the effective refractive index of the optical fiber. Indicates the sensing length of the optical fiber. Indicates the wavelength of the laser source. Indicates fiber strain;

[0036] The acoustic emission energy is calculated by energy integration:

[0037] ;

[0038] in, This represents the acoustic emission energy, that is, the cumulative energy of an acoustic emission event. Representing discrete time points The amplitude of the acoustic emission signal at that location. Indicates the sampling time interval;

[0039] The acoustic emission events are detected using a short-window to long-window ratio algorithm.

[0040] ;

[0041] in, express The ratio of the amplitude of the short time window to that of the long time window at a given moment is used to determine a valid acoustic emission event when this ratio exceeds a preset threshold. This represents the number of sample points within the short time window. This represents the number of sample points within a long time window. This represents the index of the sample points within the short time window. This represents the index of the sample points within the long time window. Indicates time The amplitude of the acoustic emission signal, Indicates time The amplitude of the acoustic emission signal;

[0042] The standardized drilling rate index is calculated using drilling parameters:

[0043] ;

[0044] in, Indicates the standardized drilling rate index. Indicates drilling speed. Indicates drilling pressure. This indicates the drill bit rotation speed.

[0045] Furthermore, in step 3, the formula for calculating the credibility weight is as follows:

[0046] ;

[0047] wherein, represents the credibility weight of the physical field characteristic of the th physical field at the spatial position , time , represents the mutual interference entropy of the th physical field at time , frequency , represents the consistency residual of the observation data of the th physical field, and represents the empirical adjustment factor, represents a minimum constant;

[0048] The multi-modal feature vectors are weighted and fused, and the formula is as follows:

[0049] ;

[0050] wherein, represents the unified fusion feature vector at the spatial position , time , represents the feature vector of the th physical field at the spatial position , time , represents a feature normalization operator.

[0051] Further, in step 4, the mathematical expression of the joint inversion objective function is as follows:

[0052] ;

[0053] wherein, represents the joint inversion objective function, and represents the surrounding rock physical parameter field vector to be inverted, represents the forward response vector of the th physical field, that is, the theoretical response calculated from the surrounding rock physical parameter field vector , represents a data weight matrix composed of credibility weights, represents a smoothing constraint coefficient, represents a model smoothing degree difference operator, represents a structure coupling constraint coefficient, represents a mutual information weight, wherein, represents a mutual information calculation operator for quantifying the spatial structure dependence between and , represents a spatial position , time a set of petrophysical parameters of the surrounding rock, and represent spatial gradients of different petrophysical parameter fields, represents the L2 norm.

[0054] Further, in step 4, the Gassmann equation is as follows:

[0055] ;

[0056] wherein, represents the bulk modulus of fluid-saturated rock, represents the bulk modulus of dry rock skeleton, represents the bulk modulus of rock matrix minerals, represents the bulk modulus of pore fluid, represents the rock porosity;

[0057] The Archie formula is as follows:

[0058] ;

[0059] wherein, represents the rock conductivity, represents the pore water conductivity, represents the water saturation, represents the lithology empirical coefficient.

[0060] Further, in step 4, the dynamic permeability field is calculated, and the corresponding formula is as follows:

[0061] ;

[0062] wherein, represents the permeability at spatial position , time , represents the initial permeability in the petrophysical parameter field, represents the strain at spatial position , time obtained by joint inversion in step 4, represents the strain rate at spatial position , time , represents the acoustic emission energy at spatial position , time , , and represent the permeability coupling coefficients of strain, strain rate and acoustic emission energy, respectively, denotes the natural exponential function.

[0063] Further, in step 5, the physical property parameter field in the three-dimensional surrounding rock physical property model is dynamically updated, including the following processes:

[0064] Taking the physical property parameter field in the three-dimensional surrounding rock physical property model as an initial value, an initial set is generated through random disturbance:

[0065] ;

[0066] ;

[0067] wherein, denotes the initial set, denotes the i-th set member in the initial set, denotes the number of set members, denotes the physical property parameter field in the three-dimensional surrounding rock physical property model, and the corresponding physical property parameters include wave velocity, resistivity, porosity and permeability, denotes the random disturbance vector of the i-th set member, which is subject to a multi-dimensional normal distribution with a mean of 0 and a covariance of ; ; In each assimilation cycle, based on the geology-mechanics-seepage coupling model, each set member is subjected to forward evolution prediction, and the evolution equation is as follows:

[0068]

[0069] ;

[0070] wherein, denotes the predicted value of the i-th set member at time t, denotes the geology-mechanics-seepage coupling model, denotes the assimilated value of the i-th set member at time t, denotes the process noise of the i-th set member at time t; After new observation data including while-drilling parameters, acoustic emission events and optical fiber strain monitoring data are obtained during the construction process, the state of the model is mapped to the observation space through the observation operator, the Kalman gain matrix is calculated, and the new observation data is assimilated into the model by using the Kalman gain matrix: ;

[0071] wherein,

[0072] ;

[0073] wherein, ​​​​​​represents the analysis value of the assimilated first ensemble member, represents the forecast value of the first ensemble member, represents the perturbation observation data of the first ensemble member, represents the theoretical observation value corresponding to the forecast value of the first ensemble member, represents the Kalman gain matrix, , represents the forecast error covariance matrix, represents the observation operator, represents the observation error covariance matrix, represents the matrix transpose.

[0074] Further, in step 5, the calculation formula of the basic risk probability is as follows:

[0075] ;

[0076] wherein, represents the basic risk probability at spatial position and time , represents the Logistic regression model, and respectively represent the resistivity and porosity of the first ensemble member in the updated physical property parameter field at spatial position , and respectively represent the permeability and pore pressure change of the first ensemble member in the updated physical property parameter field at spatial position and time , represents the coefficient of the Logistic regression model.

[0077] Further, in step 5, the seepage-uncertainty coupling correction coefficient is introduced to dynamically correct the basic risk probability, which is realized by the following formula:

[0078] ;

[0079] ;

[0080] wherein, represents the basic risk probability at spatial position and time , represents the seepage-uncertainty coupling correction coefficient at spatial position and time The basic risk probability, Indicates the probability of basic risk The final risk probability after correction This represents the correction factor for the seepage-uncertainty coupling. The Frobenius norm represents the uncertainty tensor. ,in, , and These represent the model uncertainties for resistivity, porosity, and permeability, respectively. Indicates spatial location ,time penetration rate Indicates the initial penetration rate. It represents the average uncertainty of all physical property parameters.

[0081] Furthermore, in step 5, the surrounding rock stability evaluation includes calculating the confidence interval for the final risk probability, using the following formula:

[0082] ;

[0083] in, Indicates spatial location ,time The 95% confidence interval.

[0084] Secondly, a tunnel surrounding rock stability evaluation system based on multi-physics parameter inversion is provided to implement the tunnel surrounding rock stability evaluation method based on multi-physics parameter inversion as described in the first aspect. The system includes:

[0085] A multimodal sensor array, deployed in a deep-buried tunnel, includes six types of sensors: acoustic, optical, electrical, magnetic, vibration, and drilling sensors. It is calibrated with spatial coordinates and synchronized with time to form a unified spatiotemporal reference for the original multiphysics field signal.

[0086] The scheduling and control module is used to establish a multi-field mutual interference entropy spectral density model for quantifying the interference level between signals from different physical fields. Based on the multi-field mutual interference entropy spectral density model, the mutual interference level of each combination of physical fields is quantitatively characterized and an adaptive scheduling strategy is generated. According to the adaptive scheduling strategy, the excitation sources of each physical field are controlled to perform time-division or parallel collaborative excitation, and the response signals of the multi-modal sensing array are collected synchronously to obtain a multi-source collaborative acquisition dataset.

[0087] a data processing module, configured to perform data quality verification and space-time consistency correction on the multi-source collaborative acquisition dataset, extract physical field characteristics including seismic wave first arrival time, apparent resistivity, Cole-Cole polarization parameters, GPR reflection interface depth, optical fiber strain, acoustic emission energy, and normalized drilling speed index from the corrected data, construct multi-modal feature vectors containing feature vectors corresponding to the physical field characteristics, calculate the reliability weights of each physical field characteristic based on mutual interference entropy and observation consistency residuals, form a data weight matrix with the reliability weights corresponding to each physical field characteristic, and perform weighted fusion on each multi-modal feature vector to obtain a unified fusion feature vector;

[0088] an inversion module, configured to construct a joint inversion objective function embedded with Gassmann equation and Archie formula according to the unified fusion feature vector and the data weight matrix, apply the data weight matrix to a weighted residual term of the joint inversion objective function, solve the joint inversion objective function through a Gauss-Newton iterative algorithm, and obtain a three-dimensional surrounding rock physical parameter field including wave velocity, resistivity, porosity, and strain field; construct a permeability dynamic evolution function with strain, strain rate, and acoustic emission energy as independent variables based on the strain field in the three-dimensional surrounding rock physical parameter field and the acoustic emission energy extracted in step 3, perform space-time dynamic updating on the initial permeability in the three-dimensional surrounding rock physical parameter field, calculate a dynamic permeability field, and finally form a three-dimensional surrounding rock physical model including wave velocity, resistivity, porosity, and permeability.

[0089] a dynamic updating and surrounding rock stability evaluation module, configured to take the three-dimensional surrounding rock physical model as an initial background field, use an ensemble Kalman filter algorithm to assimilate real-time acquisition of while-drilling parameters, acoustic emission events, and optical fiber strain monitoring data into the three-dimensional surrounding rock physical model, and perform dynamic updating on the physical parameter field in the three-dimensional surrounding rock physical model; based on the updated physical parameters, calculate the basic risk probability of each spatial position and time using Logistic regression, construct a physical parameter uncertainty tensor representing the uncertainty of the resistivity, porosity, and permeability model based on the updated physical parameter field, calculate the Frobenius norm of the physical parameter uncertainty tensor, determine a seepage-uncertainty coupling correction coefficient in combination with the relative increment between the dynamic permeability and the initial permeability, apply the seepage-uncertainty coupling correction coefficient as a multiplicative factor to the basic risk probability, perform space-time dynamic correction on the basic risk probability, obtain a final risk probability, and perform surrounding rock stability evaluation according to the final risk probability value.

[0090] In a third aspect, a computer readable storage medium is provided, and the computer readable storage medium stores a computer program, and when the computer program is executed, the steps of the tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion as described in the first aspect are implemented.

[0091] The tunnel surrounding rock stability evaluation method and system based on multi-physical field parameter inversion provided by the present application has the following beneficial effects: a unified acquisition reference is established by constructing a multi-modal sensor array including six types of sensors of sound, light, electricity, magnetism, shock and drilling and implementing strict space-time synchronous calibration, and the synergy and comparability of multi-source data are guaranteed from the source; the multi-field mutual interference entropy spectrum density model is established, and the reliability weight is calculated based on the mutual interference entropy and the observation consistency residual, so that the multi-field mutual interference entropy is promoted from a simple data quality index to a key variable for controlling the data weight of the joint inversion target function, a strong coupling chain of "mutual interference entropy-reliability weight-multi-field joint inversion" is formed, the interference of strong mutual interference and low reliable data on the physical property inversion is effectively inhibited, and the multi-solution and overall error of single field inversion are reduced; the seepage rate dynamic evolution function driven by strain, strain rate and acoustic emission energy is constructed, the opening and expansion rate of the micro-cracks of the surrounding rock and the energy release process are explicitly mapped to the space-time change of the seepage rate, the cross-physical field coupling representation of "sound-strain-acoustic emission-seepage rate evolution-seepage channel development" is realized, and the key input reflecting the evolution degree of the seepage channel is provided for the water gushing risk assessment; the seepage-uncertainty coupling correction coefficient is constructed by the norm of the physical property parameter uncertainty tensor and the relative increment of the seepage rate, and directly acts on the basic risk probability calculation, a multiplicative correction mechanism of "uncertainty-correction coefficient-risk probability" is formed, the risk probability of the high-uncertainty and high-seepage area is adaptively amplified, and the reliability of the water gushing risk identification is significantly improved. Through the above technical design, the whole-process closed-loop control of the deep-buried tunnel surrounding rock from multi-field collaborative perception, physical property constraint joint inversion, dynamic model assimilation to risk probability quantitative assessment is realized, and the accuracy, real-time performance and engineering usability of the construction period surrounding rock risk assessment are significantly improved. BRIEF DESCRIPTION OF DRAWINGS

[0092] Figure 1 A flowchart of the tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion provided for the embodiment is shown in the figure;

[0093] Figure 2 A flowchart of the layout and calibration of the multi-modal sensor network provided for the embodiment is shown in the figure;

[0094] Figure 3 A flowchart of the collaborative excitation and data synchronous acquisition based on the mutual interference matrix provided for the embodiment is shown in the figure;

[0095] Figure 4A flowchart of multi-source data preprocessing and physical feature extraction provided for the embodiment is shown in the figure;

[0096] Figure 5 A flowchart of multi-field parameter joint inversion under the constraint of physical properties provided for the embodiment is shown in the figure;

[0097] Figure 6 A flowchart of dynamic updating and risk quantitative partitioning based on set Kalman filtering provided for the embodiment is shown in the figure;

[0098] Figure 7 A structural diagram of a tunnel surrounding rock stability evaluation system based on multi-physical field parameter inversion provided for the embodiment is shown in the figure. DETAILED DESCRIPTION

[0099] In the existing deep-buried tunnel surrounding rock stability evaluation method, the data is isolated due to the non-uniformity of the time and space reference in the multi-physical field data collection, the single field parameter inversion has strong multi-solution property, leading to insufficient accuracy of geological interpretation, and the surrounding rock stability evaluation lags behind the actual state change of the surrounding rock due to the lack of effective dynamic updating mechanism, which is difficult to meet the precise early warning demand during the construction period.

[0100] Based on this, the technical scheme of the present application is proposed, and the present application constructs three key technical chains that are mutually coupled around the "multi-source collaborative collection - physical property constraint inversion - dynamic risk assessment": first, the mutual interference entropy spectrum density model is used to measure the interference degree between various types of monitoring signals, the mutual interference entropy and the observation consistency residual are jointly converted into the reliability weight of each physical field feature, and the weight is introduced into the multi-field joint inversion objective function as the data weight matrix, so that the physical field with smaller mutual interference entropy and more reliable data automatically obtains higher weight in the inversion, thereby weakening the influence of strong mutual interference data and reducing the multi-solution property and error of physical property inversion; second, the permeability of the surrounding rock is explicitly constructed as a function of strain, strain rate and acoustic emission energy, and the dynamic evolution function of the permeability coupled by sound-strain-acoustic emission is used to describe the influence of micro-crack opening, expansion and energy release on the development of seepage channel, thereby providing a key parameter reflecting the evolution degree of seepage channel for subsequent risk assessment; third, the seepage-uncertainty coupling correction coefficient is constructed based on the norm of the physical property parameter uncertainty tensor and combined with the relative increment of permeability, the model uncertainty and the development degree of seepage channel are jointly mapped into multiplicative correction of the basic risk probability, the risk probability is dynamically amplified or contracted in time and space, and the reliability of water inrush risk identification in high uncertainty area is improved.

[0101] Specifically, the application first establishes a unified acquisition reference by laying out six types of sensors including sound, light, electricity, magnetism, shock and drilling to form a multi-modal sensor array and performing time and space synchronous calibration; generates an adaptive scheduling strategy based on a multi-field mutual interference entropy spectrum density model to realize collaborative excitation and synchronous acquisition of multi-physical field signals; after quality checking and space-time consistency correction of the obtained multi-source collaborative acquisition data set, extracts various physical field features and constructs a multi-modal feature vector, calculates a reliability weight based on mutual interference entropy and observation consistency residual, and obtains a unified fusion feature vector through weighted fusion; based on this, a joint inversion objective function embedded with Gassmann equation and Archie formula is constructed, the data weight matrix is applied to the weighted residual term of the joint inversion objective function, and a three-dimensional surrounding rock physical parameter field containing wave velocity, resistivity, porosity and strain field is obtained by solving through Gauss-Newton iteration, and based on the strain field in the three-dimensional surrounding rock physical parameter field and the acoustic emission energy extracted in step 3, a permeability dynamic evolution function with strain, strain rate and acoustic emission energy as independent variables is constructed, the initial permeability in the three-dimensional surrounding rock physical parameter field is updated dynamically in space and time, and a dynamic permeability field is calculated; finally, taking the three-dimensional physical model as the initial background field, the real-time monitoring data is assimilated into the model to realize dynamic updating by using the ensemble Kalman filter algorithm, and based on the updated physical parameters, the final risk probability is obtained by Logistic regression and seepage-uncertainty coupling correction, and the whole process closed loop from multi-field data acquisition to dynamic surrounding rock stability evaluation is completed.

[0102] The technical solutions in the embodiments will be clearly and completely described below with reference to the drawings in the embodiments. Obviously, the described embodiments are only part of the embodiments of the present application, not all embodiments.

[0103] Figure 1 A flowchart of a tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion is shown, please refer to Figure 1 , which comprises the following steps:

[0104] S1, layout and calibration of multi-modal sensor network:

[0105] In a deep-buried tunnel, a multi-modal sensor array containing six types of sensors including sound, light, electricity, magnetism, shock and drilling is laid out, and its spatial coordinate calibration and time synchronization are performed to form a time and space reference unified original multi-physical field signal.

[0106] This step aims to lay the foundation for the stability and accuracy of the whole surrounding rock stability evaluation system, ensuring that the sensor network can work efficiently under multi-physical field conditions, and can simultaneously consider anti-seismic performance and space utilization efficiency. The core is to optimize the collaborative work of "sound-optical-electric-magnetic-seismic-drilling" multi-physical field through reasonable sensor layout, flexible circuit board design and system calibration.

[0107] Please refer to Figure 2 In practical applications, S1 specifically includes S101 to S105:

[0108] S101, spatial layout and coordinate calibration:

[0109] A hierarchical layout structure is constructed inside the tunnel: the outer layer is the rock wall sensor array layer, the inner layer is the signal acquisition and communication backbone layer, and the central layer is the flexible interconnection core layer; Fixed sensors (including acoustic emission sensors, electrodes, optical fiber end points, and UWB anchor points) are accurately coordinate calibrated by total station to ensure unified docking with the BIM model; Mobile platforms (including TBM and track trolley) integrate UWB tags and IMU, and use tight coupling filtering algorithm to solve six degrees of freedom pose in real time, and the state update formula is:

[0110] ;

[0111] Where, represents the state vector at time , including the position, attitude and speed of the mobile platform, represents the state transition matrix, represents the state vector at time , represents the control input matrix, represents the control input vector, represents the process noise vector.

[0112] To avoid cross interference between different physical field signals, a differentiated layout strategy is adopted: acoustic sensors are placed away from electric method electrodes and optical fiber end points, and acoustic and seismic sensors are preferentially placed in the rock wall close area; Optical fiber sensors are arranged in the flexible guide groove along the ring of the tunnel; Electric method electrodes and electromagnetic sensing units are placed in the contact area between the middle of the tunnel and the supporting structure; The while-drilling signal acquisition module is placed at the tail end of the drill pipe and above the drill bit.

[0113] S102, time synchronization calibration:

[0114] Adopting distributed hybrid synchronization architecture, the master control end is configured with high-stability rubidium atomic clock, and time signals are distributed through fiber ring network, combined with White Rabbit or enhanced precision time protocol to realize nanosecond-level synchronization; high-stability local oscillators are embedded in each acquisition node end, and UWB timing module and wireless synchronization signal compensation mechanism are introduced in wireless or relay section to realize hybrid synchronization of optical fiber and wireless signals; two-way delay measurement algorithm and temperature drift real-time compensation model are adopted to dynamically correct fiber transmission delay error; dual master clock redundancy structure and adaptive switching strategy are designed, and when any master clock is abnormal, the standby node automatically takes over the timing function; through field alignment test and hardware trigger verification, the system time synchronization error is kept within ±1 μs in the long term;

[0115] S103, flexible integrated system design and implementation:

[0116] Adopting rigid-flex modular structure: the rigid module undertakes signal conditioning, power distribution, data acquisition and communication functions; the flexible connection area realizes sensor interconnection and deformation buffering; the rigid part adopts multi-layer shielding PCB board structure and is configured with anti-interference grounding layer; the flexible part selects high-strength polyimide substrate, and the cable layout adopts serpentine wiring and shielding braid layer design; A30-40 hardness silicon-based elastic potting glue layer is filled between the circuit board and the probe rod inner wall to form a buffering and energy-absorbing and heat conduction channel; the system key nodes adopt waterproof quick connector and redundant signal channel design; the embedded multi-channel digital interface bus realizes unified communication and power supply of various sensors; the power supply part adopts distributed DC-DC isolation power supply architecture, and each node has overvoltage, short circuit and self-recovery protection function;

[0117] S104, sensor interconnection and layout:

[0118] According to the physical field type and interference sensitivity, differential layout is adopted, and physical isolation and protocol layering mechanism is adopted: electromagnetic shielding layer and optical fiber transmission isolation are set between acoustic and electromagnetic sensors; acoustic and seismic sensor signals are transmitted through LVDS high-speed link, optical fiber sensor data is transmitted through distributed optical fiber bus, and electric and electromagnetic signals are collected through RS485 isolated bus; all data are converged in the flexible integrated mainboard, and are sent to the main controller after unified packaging and time stamp correction by the FPGA total control module; the communication link adopts star-ring hybrid topology structure, and TSN industrial switching node is set in the backbone layer; the power supply adopts distributed DC bus and node local voltage stabilization scheme, and each sensor unit has overcurrent and short circuit self-recovery function;

[0119] S105, comprehensive layout scheme:

[0120] The system adopts a hierarchical integration and regional control architecture: the top layer is the central control layer, responsible for clock synchronization, data aggregation, and remote scheduling; the middle layer is the zone acquisition layer, containing excitation and acquisition modules for each physical field; the bottom layer is the sensing and execution layer, composed of multi-type sensor arrays; the system is deployed according to a three-level layout principle of "main measurement area - auxiliary area - reference area". The main measurement area is equipped with multi-physical field intersection nodes, the auxiliary area is equipped with fiber optic strain, temperature, and electrode arrays, and the reference area is used for long-term monitoring and noise correction; a combination of time slot allocation and spectrum isolation is used for acquisition and scheduling. Only physical fields with low mutual interference levels are allowed to run in parallel within the same time slot, while systems with high mutual interference levels are started in a time-sharing manner; scheduling commands are issued in real time by the central control layer through the TSN industrial network, and all acquisition nodes automatically execute synchronization triggering and data caching to form a distributed collaborative acquisition system;

[0121] Through the above technical solutions, a comprehensive deployment system with high robustness, high synchronization and multi-field coordination capabilities has been constructed, providing stable, controllable and high-precision basic support for subsequent collaborative excitation and synchronous acquisition.

[0122] S2. Cooperative excitation and synchronous data acquisition based on mutual interference matrix:

[0123] A multi-field mutual interference entropy spectral density model is established to quantify the interference level between signals from different physical fields. An adaptive scheduling strategy is generated based on the multi-field mutual interference entropy spectral density model. According to the adaptive scheduling strategy, the excitation sources of each physical field are controlled to perform time-division or parallel collaborative excitation, and the response signals of the multi-modal sensing array are collected synchronously to obtain a multi-source collaborative acquisition dataset.

[0124] This step aims to resolve the mutual interference problem when multiple physics fields operate in parallel, ensuring that the multi-modal signals of "acoustic-optical-electrical-magnetic-vibration-drilling" remain consistent in time and space, and achieving high-precision multi-source collaborative acquisition. Logically, this step follows the integrated deployment scheme of S1 and is a key link in the system's transition from "static deployment" to "dynamic collaborative operation".

[0125] Please see Figure 3 In practical applications, S2 specifically includes S201 to S203:

[0126] S201, Interference Matrix Modeling and Adaptive Scheduling:

[0127] A multi-field interfering entropy spectral density model is constructed based on the power spectral density of the acquired sound field, light field, electric field, magnetic field, seismic field, and drilling signal. , define the first The mutual perturbation entropy of the physical fields is:

[0128] ;

[0129] in, Indicates the first The physical field and the first A physical field signal in time ,frequency The normalized cross-interference energy ratio, , Indicates the first The physical field and the first A physical field signal in time ,frequency The cross-correlation spectral density, Indicates the first A physical field in time ,frequency The cross-interference entropy is used to quantify the degree of interference between signals from different physical fields. The larger the cross-interference entropy value, the greater the degree of interference from other physical fields, and the lower its reliability.

[0130] A standard mutual interference template library is established, with template parameters including physical field type, frequency band distribution, power level, duty cycle, and sampling period characteristics. The FPGA main control unit dynamically selects the optimal mutual interference template based on the environmental noise spectrum, signal-to-noise ratio, real-time calculated mutual interference entropy, and node operating status. An embedded online learning module dynamically corrects mutual interference template deviations and fine-tunes parameters to achieve adaptive optimization of time slot allocation and power control. Interference suppression is achieved using frequency domain isolation, time-division multiplexing, and pseudo-random coding.

[0131] S202, Collaborative Excitation Control and Acquisition Management:

[0132] The system adopts a modular independent excitation unit and a centralized timing control architecture. Each type of physical field corresponds to an independent excitation module, including a seismic exciter, an electrical current source, a GPR pulse transmitter, an acoustic emission exciter, and an optical interference excitation device. All excitation modules are connected to the unified timing bus of the FPGA master control, and nanosecond-level synchronous triggering is achieved through a high-precision clock signal.

[0133] Taking the seismic excitation signal as an example, its time-varying waveform is expressed as follows:

[0134] ;

[0135] in, This represents the time-varying amplitude envelope, that is, the amplitude envelope that changes with time. Indicates the starting frequency. Frequency modulation (FM) indicates the fundamental frequency of the signal. This represents the time variable. The formula reflects the linear frequency modulation (LFM) characteristic of the excitation signal, which can improve the penetration and spectral resolution of the excitation signal under limited energy conditions.

[0136] According to the mutual interference template distribution result, the excitation timing is automatically scheduled, the excitation power, waveform parameters and duration are dynamically adjusted; for high-power systems, an interlock protection mechanism and hardware isolation control are set, and for low-power systems, a short-time high-frequency pulse mode is used for operation.

[0137] S203, synchronous work of cooperative excitation and collection:

[0138] A centralized control-distributed collection-unified synchronization operation mechanism is established, each collection node automatically completes time alignment and signal normalization, and the collected data is attached with a high-precision time stamp and uploaded in real time to a central data processing unit through a TSN industrial network or a fiber bus.

[0139] Time partitioning and signal priority management strategies are adopted: real-time collection of transient response sensitive signals (such as seismic and acoustic emission data) is prioritized, and low-frequency steady-state signals (such as resistivity and electromagnetic field data) use a cache and delayed upload mechanism; after collection is completed, synchronization verification and quality checking are performed, delayed offset data is automatically corrected through a time stamp backtracking algorithm, and a multi-physical field original data set with strict spatiotemporal consistency and traceability is output.

[0140] S3, multi-source data preprocessing and physical feature extraction:

[0141] The multi-source cooperative collection data set is subjected to data quality checking and spatiotemporal consistency correction, physical field features including seismic wave first arrival time, apparent resistivity, Cole-Cole polarization parameters, GPR reflection interface depth, optical fiber strain, acoustic emission energy and normalized drilling speed index are extracted from the corrected data, a multi-modal feature vector containing the corresponding feature vectors of the physical field features is constructed, the reliability weights of each physical field feature are calculated based on mutual interference entropy and observation consistency residual, the reliability weights corresponding to each physical field feature are formed into a data weight matrix, and each multi-modal feature vector is weighted and fused to obtain a unified fusion feature vector.

[0142] This step is based on the multi-source cooperative collection data set obtained by S2 synchronous collection, and carries out data quality checking, noise suppression and feature extraction, aiming to convert complex multi-field observations into a set of physically meaningful invertible parameters. Through unified correction of the original data in the time domain, spatial domain and frequency domain, the data between different physical fields are ensured to have spatiotemporal consistency and feature comparability. This stage is a key link for the transition from "synchronous collection" to "multi-field fusion modeling", and provides an input basis for subsequent physical property constraint joint inversion.

[0143] Please refer to Figure 4 In actual application, S3 specifically includes S301 to S308:

[0144] S301, data quality checking and spatiotemporal consistency correction:

[0145] This step aims to perform quality inspection and space-time registration on various "sound-optical-electric-magnetic-seismic-drilling" original data after multi-source collaborative collection is completed, to ensure that the data entering the feature extraction stage have unified space-time reference and reliable physical consistency, which is the key link of the entire multi-field test chain. Specifically, it includes:

[0146] Time consistency verification is performed by comparing the PTP (Precision Time Protocol) synchronization signal of the master control end with the local time of each collection node to calculate the time deviation:

[0147] ;

[0148] Where, represents the time deviation, represents the local timestamp of the sensor node, represents the system master clock time, and when ( usually 1 μs) meets the synchronization accuracy requirements. This verification process not only verifies the stability of the hardware time link, but also ensures the strict correspondence of different physical field signals in the collection time, thereby laying a unified time reference for subsequent signal feature extraction and inversion.

[0149] Space consistency registration is performed to correct the spatial offset of fixed and mobile sensing units (such as probes installed on TBM shield heads and track cars) through rigid body coordinate transformation:

[0150] ;

[0151] Where, represents the original coordinate point, represents the rotation matrix, represents the translation vector, is the calibrated coordinate. This transformation ensures that the data of each sensing node is expressed in a unified BIM / GIS coordinate system, achieving spatial co-reference between different collection platforms and different sensing modes.

[0152] This step logically connects S2 and is the transition link from "spatial-temporal collaborative collection" to "physical feature extraction"; at the same time, it provides a unified data geometry and time framework for the multi-field joint inversion of S4, ensuring the physical comparability and precision stability of the model input.

[0153] S302, seismic data preprocessing and feature extraction:

[0154] Seismic signals reflect the wave velocity structure and stress state of surrounding rock, and its preprocessing aims to improve the signal-to-noise ratio and extract the propagation characteristics of seismic waves.

[0155] The Wiener filter method is used in this embodiment to optimize seismic data:

[0156] ;

[0157] wherein, represents the response characteristic of the filter at frequency , is the signal power spectrum, is the noise power spectrum. The filter improves the signal-to-noise ratio by suppressing the noise frequency band and enhancing the effective frequency band, and provides clearer waveforms for subsequent first arrival identification. This formula connects the results of the previous stage of synchronous acquisition, making the signal transition from multiple fields to pure and identifiable.

[0158] The AIC (Akaike Information Criterion) criterion is used to automatically identify the seismic wave first arrival time:

[0159] ;

[0160] wherein, represents the possible seismic wave first arrival time, i.e., the th sampling point in the data sequence, represents the AIC criterion value calculated when the first arrival time is assumed to be the th sampling point, represents the variance of the first samples, represents the variance of the th to th sample, represents the total number of samples in the time window. By minimizing AIC, the wave arrival time can be accurately determined, and the velocity structure information can be extracted. This formula connects the results of the filtering step in the time domain and provides high-precision arrival time data for the next step of velocity inversion.

[0161] S303, electrical data preprocessing and feature extraction:

[0162] The electrical signal reveals the conductivity and polarization characteristics of the surrounding rock, and the data preprocessing aims to remove noise and extract electrical properties.

[0163] The apparent resistivity is calculated by the electrode potential and the injected current:

[0164] ;

[0165] wherein, represents the apparent resistivity, represents the device coefficient, represents the measured potential difference, The injected current is represented. The formula converts the original electrical signal into a physical quantity, reflecting the electrical conductivity characteristics of the rock mass, and is the direct link between the electrical field data and the geological properties.

[0166] The Cole-Cole polarization parameters include time constants, which are obtained by fitting the Cole-Cole model to the complex resistivity frequency dispersion curve:

[0167] ;

[0168] wherein, represents the complex resistivity, represents the zero-frequency resistivity, represents the charge rate, represents the time constant, represents the frequency-dependent coefficient, represents the imaginary unit, represents the angular frequency. The model parameterizes the frequency domain data into physical quantities that characterize the polarization behavior of the medium, providing electrical constraints for multi-field coupling inversion, converting time-domain electrical signals into frequency-domain physical parameters, and providing a unified parameter space for subsequent multi-modal fusion.

[0169] S304, GPR / EM data preprocessing and feature extraction:

[0170] GPR (Geological Radar) and electromagnetic subsystems are mainly used to detect the internal structure characteristics, joint fissures and water-bearing abnormal layers of surrounding rocks. Since the signal is a high-frequency pulsed electromagnetic wave, its time-domain data needs to be converted into depth-domain information to facilitate fusion with other physical field data.

[0171] Based on this, the original GPR data is time-zero offset corrected and band-pass filtered (commonly 10-1000MHz) in this embodiment. Time-zero correction is used to eliminate excitation delay errors, and band-pass filtering removes high-frequency noise and system noise, providing a clean signal basis for subsequent time-depth conversion.

[0172] The GPR reflection interface depth is calculated by the time-depth conversion model:

[0173] , ;

[0174] wherein, represents the GPR reflection interface depth, represents the propagation speed of electromagnetic waves in the rock mass, represents the two-way propagation time of electromagnetic waves, represents the speed of light in a vacuum, represents the relative permittivity of the rock mass. This formula completes the conversion from the time domain to the spatial domain, enabling high-frequency electromagnetic data to be fused with seismic, acoustic and electrical data in the same three-dimensional coordinate system.

[0175] The reflection coefficient is calculated for the time-depth converted signal to reflect the electrical property difference of layer boundary:

[0176]

[0177] wherein, represents the reflection coefficient, , , respectively, represent the electromagnetic wave impedance of the upper and lower layers, and respectively represent the dielectric density of the upper and lower layers, and respectively represent the electromagnetic wave propagation speed of the upper and lower layers.

[0178] S305, optical sensing data preprocessing and feature extraction:

[0179] The fiber strain is calculated by the DAS phase-strain conversion formula:

[0180]

[0181] wherein, represents the measured phase change amount, represents the effective refractive index of the optical fiber, represents the sensing length of the optical fiber, represents the wavelength of the laser light source, represents the fiber strain. This formula maps the optical signal from the phase domain to the strain domain, realizes a physically interpretable quantitative expression, and makes the optical signal have spatial correspondence.

[0182] FBG temperature-strain separation: the fiber Bragg grating (FBG) sensor is sensitive to both temperature and strain, and its wavelength drift formula is:

[0183]

[0184] wherein, is the center wavelength change, is the initial center wavelength, is the photoelastic constant, is the thermal expansion coefficient, is the thermo-optic coefficient, is the temperature change. By decoupling the double channels, the strain and temperature contributions can be distinguished, providing accurate thermal-strain input for subsequent structure health inversion.

[0185] S306, acoustic emission data preprocessing and feature extraction:

[0186] The acoustic emission (AE) system is used to identify rock mass microcrack activity; the ultrasonic system is used to obtain the elastic modulus and crack closure characteristics of the surrounding rock. ​​​

[0187] The embodiment uses short-time window and long-time window ratio algorithm (STA / LTA) to detect events:

[0188] ;

[0189] wherein, represents the short-time window and long-time window amplitude ratio at time t, and the ratio is greater than a preset threshold to determine an effective acoustic emission event, represents the number of sample points in the short-time window, represents the number of sample points in the long-time window, represents the sample point index in the short-time window, represents the sample point index in the long-time window, represents the acoustic emission signal amplitude at time . represents the acoustic emission signal amplitude at time . The algorithm identifies the micro-fracture rupture time in the rock mass through local energy mutation, and realizes the extraction of event-level features from continuous waveforms.

[0190] Acoustic emission energy is calculated by energy integration:

[0191] ;

[0192] wherein, represents the acoustic emission energy, i.e., the cumulative energy of the acoustic emission event, represents the acoustic emission signal amplitude at discrete time point , represents the sampling time interval.

[0193] S307, pre-processing and feature extraction of while-drilling data:

[0194] The normalized drilling speed index is calculated by drilling parameters:

[0195] ;

[0196] wherein, represents the normalized drilling speed index, represents the drilling speed, represents the drilling pressure, represents the drilling speed. The formula reflects the drillability of the formation through energy normalization, making the parameters of different drilling rigs comparable, and providing a unified index for the lithology strength constraint in joint inversion.

[0197] S308, multi-source feature fusion and uncertainty quantification:

[0198] The step unifies the characteristics of seismic, electrical, electromagnetic, optical, acoustic emission and while drilling and quantifies the uncertainty, realizes the fusion from multi-signal domain to unified physical parameter domain, provides weighted and reliable input data set for multi-field parameter joint inversion under physical property constraint in subsequent steps, and the core is to construct multi-modal fusion characteristics with dynamic confidence weight, and ensure that different physical field data contribute effective information as needed in inversion.

[0199] Specifically, after the pre-step processing, a series of quantifiable key physical characteristics are obtained, including seismic wave first arrival time, apparent resistivity, Cole-Cole polarization parameter, GPR reflection interface depth, optical fiber strain, normalized drilling speed index, etc. In order to ensure the space-time consistency and mathematical operability of these characteristics, they are constructed into feature vectors under the space-time reference , wherein, is the unified three-dimensional space coordinate, is the time stamp, which is consistent with the space-time reference of S1 "space-time calibration" and S2 "synchronous acquisition", and the subscript corresponds to six types of physical fields: corresponds to seismic, corresponds to electrical, corresponds to electromagnetic, corresponds to optical, corresponds to acoustic emission, corresponds to while drilling. The feature vector is aligned with the BIM / GIS unified coordinate through coordinate mapping, realizing the space co-reference and time synchronization of different sensors and different sampling frequency data, and providing a unified data carrier for subsequent weighted fusion.

[0200] The multi-modal adaptive weight function driven by mutual interference entropy is used to ensure the stability and physical reasonableness of subsequent joint inversion, and to avoid the weight deviation caused by single dependence on instrument accuracy. The embodiment proposes a multi-modal adaptive weight function driven by mutual interference entropy, which dynamically calculates the confidence weight of each physical field feature corresponding to the feature vector . The function core integrates two types of key parameters: multi-field mutual interference entropy and observation consistency residual , and the calculation formula is as follows:

[0201] ;

[0202] , wherein, represents the confidence weight of the th physical field feature in the spatial position , time , and represents the confidence weight of the th physical field in time​ , frequency of mutual interference entropy, represents the consistency residual of the first physical field observation data, and represents an empirical adjustment factor, which is dynamically adjusted according to the complexity of the tunnel geology (for example, the interference of high stratum is increased , and the model uncertainty in the high area is increased ), represents a very small constant, which is used to avoid zero denominator. The output range is normalized to [0, 1], which directly reflects the effective contribution of the corresponding feature vector .

[0203] Based on the feature vector calculated above and the corresponding dynamic confidence weight , in order to realize the unified mapping of different physical field heterogeneous characteristics and eliminate the fusion deviation caused by the dimension difference, the multi-modal feature vectors are weighted and fused, and the formula is as follows:

[0204] ;

[0205] wherein, represents the unified fusion feature vector at the spatial position , time , represents the feature vector of the first physical field at the spatial position , time , represents a feature normalization operator, which adopts a min-max normalization method, and is used to unify the numerical range of different physical field features and ensure the mathematical operability of the fused features. The unified fusion feature vector integrates the effective information of each physical field, and gives higher contribution to the data with high confidence by the confidence weight , avoiding the interference of low confidence data on the inversion result.

[0206] After the unified fusion feature vector is constructed, the spatial interpolation and smoothing constraint of the fused feature field are performed to ensure the continuity and physical rationality of each physical parameter in the spatial distribution. Through the Kriging interpolation or inverse distance weighting (IDW) method, the sparse point data is filled in space to generate a high-resolution multi-modal feature body, which establishes the physical continuous boundary condition for the subsequent “three-dimensional joint inversion”.

[0207] Through the multi-source feature extraction and fusion processing of this step, the final output is a unified fusion feature vector and a corresponding confidence weight matrix. This dataset not only reflects the physical response characteristics of the surrounding rock structure, but also quantifies the data reliability through dynamic weights, directly serving as the core input for subsequent multi-field parameter joint inversion, realizing the key transition from "multi-source heterogeneous data" to "unified inversion input", and ensuring the process loop from multi-source observation to model inversion.

[0208] Further, by introducing a confidence weight mechanism jointly driven by mutual interference entropy and observation consistency residual, the embodiment completes the quantitative discrimination of the reliability of different physical field data in the feature fusion stage: the smaller the mutual interference entropy value and the lower the observation residual of the physical field, the greater the confidence weight, and the higher the contribution in the unified fusion feature vector; the weight of data with large mutual interference entropy or large residual is adaptively weakened. Since the weight matrix is directly used as data weight in the subsequent multi-field joint inversion objective function, the mutual interference entropy and the consistency residual are not only used for data quality evaluation, but also change the sensitivity and dependence of the inversion process on different field data by affecting the weighted residual term of the objective function, realizing a strong coupling chain of "mutual interference entropy→confidence weight→joint inversion". Compared with the traditional equal weight or empirical weight method, this scheme can effectively suppress the interference of abnormal data on the physical property inversion result and reduce the multi-solution and overall error level of the physical property parameter field under the condition of significant multi-source data interference and uneven local observation quality.

[0209] S4, multi-field parameter joint inversion under physical property constraints:

[0210] According to the unified fusion feature vector and the data weight matrix, a joint inversion objective function embedded with Gassmann equation and Archie formula is constructed, the data weight matrix is applied to the weighted residual term of the joint inversion objective function, the joint inversion objective function is solved by Gauss-Newton iterative algorithm, and a three-dimensional surrounding rock physical property parameter field including wave velocity, resistivity, porosity and strain field is obtained; based on the strain field in the three-dimensional surrounding rock physical property parameter field and the acoustic emission energy extracted in step 3, a permeability dynamic evolution function with strain, strain rate and acoustic emission energy as independent variables is constructed, the initial permeability in the three-dimensional surrounding rock physical property parameter field is updated in space and time, and a dynamic permeability field is calculated, and finally a three-dimensional surrounding rock physical property model including wave velocity, resistivity, porosity and permeability is formed.

[0211] This step solves the parameter field model consistent in statistics and physics by establishing an objective function with weight adjustment, rock physics embedding and structure coupling term, thereby providing a stable and reliable model basis for subsequent dynamic assimilation and risk prediction.

[0212] Please refer toFigure 5 In practical applications, S4 specifically includes S401 to S405:

[0213] S401. Construct the joint inversion objective function:

[0214] In this stage, the unified fused feature vector output from the previous steps Based on its credibility weight matrix, a coupled inversion objective function capable of simultaneously processing different physical field data is established. The objective function is defined as follows: Simultaneously, it minimizes the residuals of various observation data and introduces a mutually information-weighted structural coupling constraint term to ensure the boundary structure and unified fusion feature vector between physical fields such as acoustics, electrical methods, seismicity, and electromagnetics. Maintaining consistency. The mathematical expression for the joint inversion objective function is as follows:

[0215] ;

[0216] in, Let represent the joint inversion objective function, and represent the field vector of surrounding rock physical parameters to be inverted. Indicates the first The forward response vector of a physical field, i.e., the vector of the surrounding rock physical property parameters. The theoretical response obtained from the calculation Represents credibility weight The constructed data weight matrix, Represents the smoothing constraint coefficient. This represents the model smoothness difference operator. Represents the structural coupling constraint coefficient. Represents mutual information weights, ,in, Represents the mutual information computation operator, used for quantization. and The degree of spatial structural dependence between them; a larger value indicates a stronger demand for matching the fusion features and model structure in that region, and a more stringent structural constraint. Indicates spatial location ,time The set of surrounding rock physical properties, and Represents the spatial gradient of fields with different physical property parameters. This represents the L2 norm.

[0217] The above objective function is passed through By dynamically adjusting the structural constraint strength in different regions to replace the original fixed coefficients, the inversion process can both follow the inherent coupling relationship of the physical field and adapt to the spatial distribution differences of the unified fusion feature vectors, thus avoiding the constraint deviation of "one-size-fits-all".

[0218] S402, embedding of rock physics and electrical equations:

[0219] After the mathematical framework of the objective function is established, this stage explicitly embeds the physical constraint relationship between different physical fields into the inversion equation by introducing the rock physics model (including the pore fluid coupling model and the resistivity and porosity coupling model), so that the model solution not only satisfies mathematical optimality, but also satisfies physical rationality.

[0220] The pore fluid coupling model (Gassmann equation) is as follows:

[0221] ;

[0222] wherein, represents the bulk modulus of the fluid-saturated rock, represents the bulk modulus of the dry rock skeleton, represents the bulk modulus of the rock matrix mineral, represents the bulk modulus of the pore fluid, represents the porosity of the rock. This equation is used as a soft constraint term in the inversion to maintain a consistent relationship between seismic velocity and porosity, fluid type. Through this constraint, physical conflicts between seismic inversion results and electrical inversion results can be prevented, and collaborative interpretation of acoustic wave propagation characteristics and fluid state can be achieved.

[0223] The resistivity and porosity coupling model (Archie formula) is as follows:

[0224] ;

[0225] wherein, represents the rock conductivity, represents the pore water conductivity, represents the water saturation, represents the lithology empirical coefficient. This model combines electrical parameters with pore fluid characteristics to ensure that electrical changes can match rock pore structure. Its role in multi-field constraints is to connect “electrical response” and “rock properties”, thereby achieving dual constraints on fluid content and pore structure in joint inversion.

[0226] S403, numerical solution and model iteration:

[0227] After obtaining the unified objective function and physical constraints, this stage solves the parameter model through numerical optimization algorithms. The linearized form of the Gauss-Newton iteration strategy is as follows:

[0228] ;

[0229] wherein, , represents the sensitivity (Jacobian) matrix, represents the partial derivative of the forward response with respect to the model parameters, , represents the observation data covariance matrix, representing the uncertainty and correlation of each physical field observation data, , represents the smoothing constraint coefficient, controlling the model smoothness weight, , represents the model smoothness difference operator, , represents the structure coupling constraint coefficient, controlling the structure consistency weight between different parameter fields, , represents the structure coupling operator, realizing cross-gradient and other structure constraints, , represents the model parameter increment, representing the parameter correction amount of the current iteration step, , represents the observation data vector, coming from the unified fusion feature vector, , represents the forward response vector, calculated by the current model parameters, , represents the matrix transpose.

[0230] The model update formula is:

[0231] ;

[0232] , where, , represents the model parameters of the th iteration, , represents the iteration step length coefficient, which is adaptively determined by line search or trust region algorithm to ensure stable convergence. In each iteration, the system will calculate the predicted data , the observation residual and the target function change , so as to dynamically evaluate the convergence trend and ensure that the iteration process converges stably to a physically reasonable solution.

[0233] Through each round of update, the model continuously approaches the real geological structure while maintaining consistency with multi-physical field observations. The intermediate model output each time not only provides the basis for convergence judgment in this stage, but also provides the spatial gradient basis for the next step of structure consistency analysis.

[0234] S404, structure coupling and uncertainty evaluation:

[0235] After iterative solution, the model obtained in this stage is subjected to structure consistency control and uncertainty evaluation to ensure that the inversion result is not only credible but also quantifiable. Specifically, it includes:

[0236] The cross-gradient structure constraint method is adopted to minimize the cross-gradient term , to ensure the consistency of the structural interface direction of different physical field parameters (such as seismic velocity field and resistivity field); when the cross gradient is zero, it indicates that different parameter fields have the same boundary trend, realizing the matching and coordination of multi-field results in space, effectively reducing the "false interface" between different modes, and enhancing the geological interpretability of the model.

[0237] Posterior uncertainty evaluation. According to the linear approximation, the model covariance matrix can be evaluated as follows:

[0238] ;

[0239] wherein, represents the model covariance matrix, representing the uncertainty distribution of the inversion result, and the diagonal elements represent the variance distribution of the model parameters, represents the sensitivity matrix, represents the observation data covariance matrix, represents the smoothing constraint coefficient, represents the structural coupling constraint coefficient, represents the model smoothing degree difference operator, represents the structural coupling operator, represents the matrix transpose. The area with large variance usually corresponds to the high-risk or information-deficient area, which can be used as the key monitoring object for subsequent dynamic assimilation.

[0240] This step realizes the strengthening of model structural consistency through cross gradient constraint, and completes statistical precision evaluation combined with posterior uncertainty quantification, forming a double constraint closed loop in physical and statistical aspects, providing initial conditions for subsequent real-time assimilation and risk analysis, and ensuring the coherent transmission of parameters from inversion to surrounding rock stability evaluation.

[0241] S405, constructing a sound-strain driven permeability dynamic evolution function:

[0242] This step aims to quantify the real-time variation law of permeability of surrounding rock caused by microfracture development during deep tunnel excavation, and to build a cross-physical field correlation mechanism of "mechanical deformation-acoustic response-seepage characteristics", providing a key parameter reflecting the real-time seepage state of surrounding rock for subsequent surrounding rock stability evaluation.

[0243] Specifically, based on the input parameters obtained in the previous steps, combined with the internal physical correlation between microfracture and seepage channel expansion of deep tunnel surrounding rock, a sound-strain driven permeability dynamic evolution function is constructed, and the formula is as follows:

[0244] ;

[0245] wherein, represents the permeability at spatial position , time , Represents the initial permeability in the physical property field. This indicates the spatial location obtained from the joint inversion in S4. ,time The response, Indicates spatial location ,time strain rate, Indicates spatial location ,time Acoustic emission energy, , and The permeability coupling coefficients, representing strain, strain rate, and acoustic emission energy respectively, need to be obtained through indoor core mechanics-seepage coupling tests or on-site multi-field synchronous calibration based on the actual surrounding rock lithology of the tunnel (such as granite, sandstone, shale, etc.). Typical values ​​range from [value range missing in original text]. , , , This represents the natural exponential function.

[0246] In the above formula, the exponent term is obtained through... It reflects the degree of crack opening and closing caused by static deformation of the surrounding rock, through Reflecting the effect of deformation rate on crack propagation, through The energy release intensity reflects the development of microfractures, and the synergistic effect of these three factors increases the initial permeability. To dynamic penetration rate Evolution—When the surrounding rock strain increases, the strain rate accelerates, or the acoustic emission energy density increases, The exponential growth directly reflects the physical process by which the expansion of microfractures leads to an increase in seepage channels and enhanced seepage capacity.

[0247] Therefore, this embodiment does not merely use acoustic emission as an independent early warning indicator, but rather incorporates strain, strain rate, and acoustic emission energy into the permeability evolution function. This allows the opening and closing degree, propagation rate, and energy release of microfractures in the surrounding rock to be uniformly characterized within the seepage parameters. When the surrounding rock enters a stage of strong disturbance, a sudden increase in strain, an accelerated strain rate, and frequent acoustic emission events will jointly drive an exponential increase in permeability, reflecting the rapid evolution of seepage channels from localized dispersion to interconnected penetration. This dynamic permeability field serves as the input for subsequent seepage analysis and risk assessment, making the probability of water inrush risk more sensitive to the "degree of development of seepage channels," thus achieving a physical chain characterization of "sound-strain-acoustic emission → permeability evolution → seepage channel development → amplified water inrush risk."

[0248] S5. Dynamic updating and quantitative risk partitioning based on ensemble Kalman filtering:

[0249] Using the three-dimensional surrounding rock physical property model as the initial background field, the ensemble Kalman filter algorithm is used to assimilate the drilling parameters, acoustic emission events, and fiber optic strain monitoring data collected in real time during construction into the three-dimensional surrounding rock physical property model, dynamically updating the physical property parameter field in the three-dimensional surrounding rock physical property model. Based on the updated physical property parameters, Logistic regression is used to calculate the basic risk probability at each spatial location and time. Based on the updated physical property parameter field, a physical property parameter uncertainty tensor characterizing the uncertainty of the resistivity, porosity, and permeability models is constructed. The Frobenius norm of the physical property parameter uncertainty tensor is calculated, and the seepage-uncertainty coupling correction coefficient is determined by combining the relative increment between the dynamic permeability and the initial permeability. The seepage-uncertainty coupling correction coefficient is used as a multiplicative factor on the basic risk probability to perform spatiotemporal dynamic correction on the basic risk probability, obtaining the final risk probability. The surrounding rock stability is evaluated based on the final risk probability value.

[0250] This step uses the three-dimensional surrounding rock physical property model and its model covariance matrix output from the previous steps. By introducing an ensemble Kalman filter algorithm, dynamic self-updating and quantitative risk probability partitioning of the three-dimensional surrounding rock physical property model are achieved. The core of this step lies in using real-time construction monitoring data (including drilling parameters, acoustic emission events, and fiber optic strain monitoring data) to assimilate the three-dimensional surrounding rock physical property model, enabling the model to adjust in real time with time and tunneling progress, forming a "dynamic geological cognitive system" with self-learning and feedback capabilities.

[0251] Please see Figure 6 In practical applications, S5 specifically includes S501 to S505:

[0252] S501, Set initialization and background field generation:

[0253] Using the physical property parameter field in the three-dimensional surrounding rock physical property model as initial values, an initial set is generated through random perturbation:

[0254] ;

[0255] ;

[0256] in, Represents the initial set. Represents the first element in the initial set. A set of members, Indicates the number of members in the set (usually 50–200). This represents the physical property parameter field in the three-dimensional surrounding rock physical property model. The corresponding physical property parameters include wave velocity, resistivity, porosity, and permeability. Indicates the first The random perturbation vectors of *n* set members follow a pattern with mean 0 and covariance of 0. The multidimensional normal distribution, i.e. .

[0257] S502, Model Prediction and Evolution:

[0258] Within each assimilation cycle (e.g., every 1m of tunneling or every 1min of sampling), based on the geological-mechanical-seepage coupling model, a forward evolution prediction is performed for each ensemble member, and the evolution equation is as follows:

[0259] ;

[0260] in, express Time of the first Predicted values ​​for each set member This represents a geological-mechanical-seepage coupling model. express Time of the first The assimilation value of each set member. express Time of the first Process noise of a set member.

[0261] Through this evolution equation, the system can automatically consider dynamic factors such as formation disturbance, drilling pressure, and stress migration. The periodic update mechanism is reflected here: whenever the tunneling progress or monitoring data meets the set conditions (such as time intervals or abnormal triggering events), it automatically enters the next assimilation cycle, realizing real-time self-learning and rolling prediction of the model.

[0262] S503, Observation Assimilation and Model Update:

[0263] After obtaining new observational data (such as drilling NDX, acoustic emission event rate, fiber strain increment, annular pressure, etc.) during the construction process, including drilling parameters, acoustic emission events, and fiber strain monitoring data, the Kalman gain matrix is ​​calculated by mapping the state of the observation operator model to the observation space and then using the Kalman gain matrix to assimilate the new observational data into the model.

[0264] ;

[0265] in, Indicates the first The analytical values ​​after assimilation of each set member, Indicates the first Predicted values ​​for each set member Indicates the first Perturbation observation data of each set member, Indicates the first The predicted values ​​corresponding to the theoretical observations of each set member Represents the Kalman gain matrix. , This represents the forecast error covariance matrix. Represents the observation operator. Represents the observation error covariance matrix. This indicates the matrix transpose.

[0266] Forecast error covariance matrix The calculation formula is as follows:

[0267] ;

[0268] ;

[0269] in, This represents the mean of the set.

[0270] To avoid excessively rapid ensemble shrinkage, this embodiment also applies localization and inflation techniques to adjust the forecast error covariance matrix:

[0271] ;

[0272] in, This represents the adjusted forecast error covariance matrix. Indicates the inflation coefficient. Represents the spatial correlation coefficient matrix (localization matrix). This represents the Schur product (element-by-element product).

[0273] This step integrates real-time observations into the model evolution, enabling the system to continuously correct prediction errors based on the latest data, thus achieving "adaptive learning" of the model. The periodic assimilation mechanism allows the model to be continuously corrected as tunneling progresses, while inflation and localization strategies ensure long-term model stability and avoid overfitting.

[0274] S504. Risk Probability Calculation and Dynamic Quantification:

[0275] This step is based on the three-dimensional update model output by S503. It achieves accurate quantification of water inrush risk through "basic probability calculation + original formula coupling correction". It retains the core logic of the original text and adds a coupling mechanism between dynamic seepage and model uncertainty.

[0276] Specifically, the basic risk probability is calculated using a Logistic regression model combined with a ensemble statistical average method:

[0277] ;

[0278] in, Indicates spatial location ,time The basic risk probability, This represents a Logistic regression model. and They represent the first in the updated physical property parameter field. The spatial location of a set member Resistivity and porosity at that location and They represent the first in the updated physical property parameter field. The spatial location of a set member ,time Changes in permeability and pore pressure, The coefficients of the Logistic regression model are dynamically updated using historical water inrush case data or an online least squares algorithm for the assimilation process. This formula, through ensemble statistical averaging, directly transfers the uncertainty of the physical model to the risk probability. If abnormal signals such as a sudden increase in acoustic emission frequency or abrupt changes in fiber strain are detected during the assimilation process, the system will automatically adjust the regression coefficients or update the observation error covariance matrix, forming an "online risk self-learning" function to ensure the timeliness and adaptability of the basic risk probability.

[0279] To overcome the limitations of traditional static parameter risk calculations and quantify the synergistic effect of "dynamic seepage characteristics - model uncertainty," this embodiment introduces a seepage-uncertainty coupling correction coefficient to dynamically correct the basic risk probability, achieved through the following formula:

[0280] ;

[0281] ;

[0282] in, Indicates spatial location ,time The basic risk probability, Indicates the probability of basic risk The final risk probability after correction Indicates spatial location ,time penetration rate This represents the initial penetration rate.

[0283] This represents the seepage-uncertainty coupling correction coefficient, with a value range of [1, 2.2]. Its core function is to dynamically adjust the risk probability weights—the higher the permeability and the greater the model uncertainty, the larger the correction coefficient and the stronger the adjustment of the risk probability, which is in line with the risk evolution law of "seepage channel expansion + model bias" in deep-buried tunnels.

[0284] Frobenius norm of uncertainty tensor, where, , and represent the model uncertainty of resistivity, porosity and permeability, respectively, taken from S404 the square root of all diagonal elements, e.g., , denotes the model covariance matrix at spatial location , denotes the diagonal elements of the extraction matrix and picks the element with index , here specifically for permeability, the same applies to resistivity and porosity.

[0285] denotes the average uncertainty of all petrophysical parameters, i.e. where, denotes the index of petrophysical parameter type, used to iterate over the three petrophysical parameters: resistivity, porosity and permeability, denotes the calculation of the model uncertainty of the th petrophysical parameter at spatial location .

[0286] With the above seepage-uncertainty coupling correction factor, the model uncertainty and the degree of seepage channel development are jointly applied to the risk probability calculation process in this embodiment: when the uncertainty tensor norm of a petrophysical parameter at a certain spatial location is large, and the permeability is significantly increased relative to the initial value, the correction factor is amplified, and the final risk probability obtained by multiplicatively correcting the basic risk probability is significantly higher than that in the uncorrected case, so that a more conservative risk assessment is automatically given in the high-uncertainty and high-seepage area; on the contrary, when the model uncertainty is low or the permeability change is not obvious, the correction factor is close to 1, and the influence on the basic risk probability is limited, avoiding over-amplification of the risk in the stable area. In this way, uncertainty no longer stays at the qualitative weight or grade division level, but directly enters the water inrush risk probability function, realizing an explicit mapping mechanism of "uncertainty → correction factor → risk probability", which is conducive to reducing the false negative risk of major water inrush events in high-uncertainty areas.

[0287] Through the above correction, the final risk probability is output, which is normalized to [0, 1]. In order to further distinguish the risk level and the uncertainty level, the confidence interval of the corrected final risk probability is calculated:

[0288] ;

[0289] in, Indicates spatial location ,time The 95% confidence interval, measured by the median, reflects the statistical average level of the risk probability, thus representing central tendency. The width of the confidence interval represents the statistical volatility of the risk probability, quantifying uncertainty. The confidence interval directly reflects the reliability of the risk probability estimate—a wider interval indicates greater volatility of the risk probability among different set members, poorer statistical stability, and lower reliability. This result provides a quantitative basis for subsequent identification of high uncertainty areas, supporting refined decision-making in risk classification.

[0290] S505. Uncertainty Assessment and Risk Classification:

[0291] This step is based on the final risk probability output by S504. seepage-uncertainty coupling correction coefficient and 95% confidence interval It completes the quantitative classification of surrounding rock water inrush risk, uncertainty labeling, and output of construction treatment suggestions, and constructs a closed-loop system of "risk calculation - classification assessment - decision feedback".

[0292] Specifically, based on the final risk probability Using this as the core indicator and considering the actual needs of the project, the monitoring area is divided into three risk levels to ensure that the surrounding rock stability assessment directly corresponds to the construction measures:

[0293] High risk (red alert): When At that time, it was identified as a high-risk area. The dynamic permeability of the surrounding rock in this area was significantly higher than expected, and the model uncertainty had been corrected by coupling coefficients. Integrating risk quantification, there is an immediate risk of water inrush. Immediate measures such as sealing and grouting reinforcement must be taken to reduce the tunneling speed to below 50% of the normal rate, and the sampling interval for monitoring such as acoustic emission and fiber optic strain should be shortened from 5 minutes to 1 minute to achieve real-time risk tracking.

[0294] Medium risk (yellow alert): When At that time, it was determined to be a medium-risk area. The seepage characteristics or model uncertainties of the surrounding rock in this area show an abnormal tendency, requiring strengthened monitoring and verification—one preliminary verification borehole should be drilled every 5 meters to track the pore pressure change trend in real time; if... If the value continues to rise above 0.6, pre-grouting preparation work should be started immediately.

[0295] Low risk (green normal): When At this time, it is determined to be a low-risk area. The surrounding rock seepage state is stable, the model uncertainty level is low, the normal excavation rhythm can be maintained, and monitoring can be carried out according to the conventional frequency (recording 1 multi-field parameter every 10m).

[0296] In practical application, a high uncertainty area identification mechanism can also be constructed: combined with the coupling correction coefficient output by S504 and the confidence interval , a double high uncertainty area identification logic is established to make up for the limitations of traditional risk probability evaluation:

[0297] High seepage-high uncertainty area: when , no matter what risk level it is, it is additionally marked with the label "high seepage-high uncertainty". This area needs to be prioritized for monitoring data due to the dynamic and drastic changes in the permeability coefficient and the large model bias. The electrode spacing of the electrical method is reduced from 2m to 1m, the number of seismic detectors is increased, and the model is updated in the next assimilation period to reduce the uncertainty level;

[0298] High probability fluctuation area: when the width of the confidence interval (interval upper limit-interval lower limit) exceeds 0.4, it is marked as "high probability fluctuation area". This area has poor statistical stability of risk probability, so it needs to optimize the Logistic regression model coefficients by supplementing historical water inrush case data, and at the same time, the data collection time is extended from 1min to 3min to reduce the interference of random noise on risk calculation and improve the reliability of probability results.

[0299] In practical application, dynamic feedback and decision are output, which not only serves the site construction decision, but also realizes feedback optimization of the dynamic assimilation process described above:

[0300] Three-dimensional dynamic risk zoning map: marked with red, yellow and green corresponding to the risk level, superimposed with special labels of "high seepage-high uncertainty area" and "high probability fluctuation area", to intuitively present the risk spatial distribution and uncertainty characteristics;

[0301] Uncertainty statistical report: contains the mean, width statistical data of each risk area, quantifying the spatial distribution of model uncertainty, providing basis for monitoring scheme optimization;

[0302] Construction disposal suggestion list: for different risk levels and uncertainty types, specific measures such as monitoring frequency adjustment, grouting reinforcement parameters, drilling scheme design are clear, which are directly fed back to the site construction control system to guide real-time construction decision.

[0303] Meanwhile, the risk grading results of this step are returned to the S501 "set initialization and background field generation" link together with the uncertainty statistical data, to optimize the initial distribution of set members in the next assimilation cycle (such as appropriately increasing the disturbance amplitude of set members in high-uncertainty areas), realize the self-learning cycle of "surrounding rock stability evaluation-model updating", and continuously improve the accuracy and adaptability of the entire dynamic assimilation and surrounding rock stability evaluation system.

[0304] The tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion provided by the embodiment establishes a unified acquisition reference by constructing a multi-modal sensor array containing six types of sensors of sound, light, electricity, magnetism, shock, and drilling and implementing strict space-time synchronous calibration; generates an adaptive scheduling strategy based on a multi-field mutual interference entropy spectrum density model, to realize the coordinated excitation and synchronous acquisition of multi-physical field signals; extracts multiple types of physical field features from the corrected data (and calculates a credibility weight based on mutual interference entropy and observation consistency residual, to construct a unified fusion feature vector with weight; constructs a joint inversion objective function embedded with Gassmann equation and Archie formula based on this, and obtains a three-dimensional surrounding rock physical parameter field with physical consistency through Gauss-Newton iteration solution; then calculates a dynamic permeability field through a sound-strain driven permeability dynamic evolution function; finally, uses the set Kalman filter algorithm to assimilate real-time monitoring data into the model to realize dynamic updating based on the three-dimensional physical model as the initial background field, and obtains the final risk probability through Logistic regression and seepage-uncertainty coupling correction based on the updated physical parameters, to form a full-process closed-loop system from multi-field data acquisition, dynamic modeling to accurate surrounding rock stability evaluation.

[0305] Compared with the prior art, the embodiment has the following remarkable effects:

[0306] In the aspect of the inversion of the physical parameters of the surrounding rock, the embodiment is configured to construct a multi-physical-field mutual interference entropy spectrum density model, to quantitatively characterize the mutual interference level between the physical fields under different acquisition schemes, to preferentially select the excitation and observation combination with a smaller mutual interference entropy in the acquisition stage, and to improve the consistency of the multi-source observation data from the source. On this basis, the embodiment combines the mutual interference entropy and the observation consistency residual to calculate the reliability weight of each physical field, and constructs a multi-field joint inversion objective function containing a physical property constraint term with the reliability weight, so that the observation data that is less disturbed and more consistent with the prior information of the physical properties obtains a higher weight in the objective function. Compared with the traditional multi-source inversion method that only relies on the linear superposition of the empirical weight, the embodiment uses the multi-field mutual interference entropy to constrain the acquisition scheduling, reduces the systematic deviation between the observation data of different physical fields, and improves the overall synergy of the multi-source data. In the weight construction process, the mutual interference entropy and the consistency residual are introduced, so that the data weight can be adaptively adjusted according to the observation quality, and the interference of abnormal or low reliability data on the inversion result is weakened. Through the physical property constraint term, the unreasonable physical property combination solution is suppressed, and the multi-solution and local oscillation phenomena of the inversion result are weakened. Therefore, under the same observation conditions, the multi-field mutual interference entropy constrained joint inversion of the embodiment is beneficial to obtain a spatially more smoothly distributed, consistent and more consistent with the actual geological conditions of the surrounding rock physical property parameter field, to improve the problems of the traditional multi-source empirical weighted inversion, such as the sensitivity to abnormal data, the unstable solution space and the like, and to improve the reliability of the inversion of the physical properties of the surrounding rock.

[0307] In the aspect of the characterization of the seepage characteristics of the surrounding rock and the evaluation of the water gushing risk, the embodiment explicitly constructs the permeability as a function of the strain, the strain rate and the acoustic emission energy, so that the permeability can make a more sensitive response to the cracking, expansion and penetration of the micro-cracks of the surrounding rock under strong disturbance conditions. When the surrounding rock enters the stage of rapid loading and unloading or significant cracking activity, the significant changes in the strain rate and the acoustic emission energy will be reflected as the nonlinear amplification of the permeability through the model, which is more conducive to depicting the sudden increase behavior of the formation and evolution of the water gushing channel.

[0308] In terms of risk assessment, the embodiment is based on the Frobenius norm of the uncertainty tensor of the physical property parameters and the relative increment of the permeability to construct a seepage-uncertainty coupling correction coefficient, directly acts the model uncertainty and the evolution degree of the seepage channel on the Logistic-type water gushing risk probability function in the form of a multiplicative coefficient, and realizes the dynamic correction of the basic risk probability in space and time. Compared with the traditional method of introducing fuzzy weight or grade division only at the evaluation index level, the embodiment amplifies the risk probability conservatively in the area where the permeability increases significantly and the parameter uncertainty is large, reduces the possibility of underestimating or missing the major water gushing event; in the area where the permeability changes little or the uncertainty is relatively low, it avoids over-amplification of the risk, and improves the discrimination and pertinence of the risk assessment results; through seepage-uncertainty coupling correction, the risk assessment results are dynamically adjusted with the update of multi-source monitoring information in the construction process, and are closer to the actual evolution process of the surrounding rock state.

[0309] In summary, the embodiment introduces the seepage channel evolution and the model uncertainty into the quantitative characterization system of the risk probability through the acoustic-strain-acoustic emission coupling permeability evolution function and the seepage-uncertainty coupling risk correction coefficient, which is beneficial to improve the sensitivity and reliability of the water gushing risk assessment of the deep-buried tunnel, and to improve the problem of missing the risk in the high-uncertainty area in the traditional method, so as to further improve the safety margin of the water gushing disaster prevention and control during the tunnel construction period.

[0310] Based on the above technical solution, the embodiment further provides a tunnel surrounding rock stability evaluation system based on multi-physical field parameter inversion, which is used to realize the tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion as described in the embodiment, please refer to Figure 7 , the system comprises:

[0311] A multi-modal sensor array is arranged in a deep-buried tunnel, which contains six types of sensors including acoustic, optical, electrical, magnetic, seismic and drilling, and is calibrated by spatial coordinates and synchronized by time, and is used to form original multi-physical field signals with unified time and space reference;

[0312] A scheduling control module is used to establish a multi-field mutual interference entropy spectrum density model for quantifying the interference degree between different physical field signals, to quantitatively characterize the mutual interference degree of each physical field combination based on the multi-field mutual interference entropy spectrum density model and generate an adaptive scheduling strategy; according to the adaptive scheduling strategy, the excitation source of each physical field is controlled to be excited in time-sharing or parallel cooperation, and the response signals of the multi-modal sensor array are synchronously collected to obtain a multi-source cooperative acquisition data set;

[0313] a data processing module, configured to perform data quality verification and space-time consistency correction on the multi-source collaborative acquisition dataset, extract physical field characteristics including seismic wave first arrival time, apparent resistivity, Cole-Cole polarization parameters, GPR reflection interface depth, optical fiber strain, acoustic emission energy, and normalized drilling speed index from the corrected data, construct multi-modal feature vectors containing feature vectors corresponding to the physical field characteristics, calculate the reliability weights of each physical field characteristic based on mutual interference entropy and observation consistency residuals, form a data weight matrix with the reliability weights corresponding to each physical field characteristic, and perform weighted fusion on each multi-modal feature vector to obtain a unified fusion feature vector;

[0314] an inversion module, configured to construct a joint inversion objective function embedded with Gassmann equation and Archie formula according to the unified fusion feature vector and the data weight matrix, apply the data weight matrix to a weighted residual term of the joint inversion objective function, solve the joint inversion objective function through a Gauss-Newton iterative algorithm, and obtain a three-dimensional surrounding rock physical parameter field including wave velocity, resistivity, porosity, and strain field; construct a permeability dynamic evolution function with strain, strain rate, and acoustic emission energy as independent variables based on the strain field in the three-dimensional surrounding rock physical parameter field and the acoustic emission energy extracted in step 3, perform space-time dynamic updating on the initial permeability in the three-dimensional surrounding rock physical parameter field, calculate a dynamic permeability field, and finally form a three-dimensional surrounding rock physical model including wave velocity, resistivity, porosity, and permeability.

[0315] a dynamic updating and surrounding rock stability evaluation module, configured to take the three-dimensional surrounding rock physical model as an initial background field, use an ensemble Kalman filter algorithm to assimilate real-time acquisition of while-drilling parameters, acoustic emission events, and optical fiber strain monitoring data into the three-dimensional surrounding rock physical model, and perform dynamic updating on the physical parameter field in the three-dimensional surrounding rock physical model; based on the updated physical parameters, calculate the basic risk probability of each spatial position and time using Logistic regression, construct a physical parameter uncertainty tensor representing the uncertainty of the resistivity, porosity, and permeability model based on the updated physical parameter field, calculate the Frobenius norm of the physical parameter uncertainty tensor, determine a seepage-uncertainty coupling correction coefficient in combination with the relative increment between the dynamic permeability and the initial permeability, apply the seepage-uncertainty coupling correction coefficient as a multiplicative factor to the basic risk probability, perform space-time dynamic correction on the basic risk probability, obtain a final risk probability, and perform surrounding rock stability evaluation according to the final risk probability value.

[0316] It can be understood that the tunnel surrounding rock stability evaluation system based on the multi-physical field parameter inversion described in the embodiment is a system for implementing the tunnel surrounding rock stability evaluation method based on the multi-physical field parameter inversion described in the embodiment. For the system disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant part can be referred to the part of the method. Therefore, it will not be described here.

[0317] The system embodiments described above are only schematic, wherein the units described as separate components can or can not be physically separate, that is, they can be located in one place, or can be distributed to multiple network units. Part or all of the modules can be selected according to actual needs to achieve the purpose of the embodiment. Those skilled in the art can understand and implement without creative labor.

[0318] Through the description of the above embodiments, those skilled in the art can clearly understand that the embodiments can be realized by means of software plus necessary universal hardware platforms, and of course, can also be realized by hardware. Based on such understanding, the above technical solutions can be embodied in the form of software products, and the computer software products can be stored in a computer readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and include a plurality of instructions to make a computer device (which can be a personal computer, a server, or a network device, etc.) execute the methods described in each embodiment or some parts of the embodiments.

Claims

1. A tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion, characterized in that, The method comprises: Step 1, in a deep buried tunnel, a multi-modal sensor array comprising six types of sensors of sound, light, electricity, magnetism, shock and drilling is laid out, and the multi-modal sensor array is spatially coordinate calibrated and time synchronized to form original multi-physical field signals with unified time and space reference; Step 2, a multi-field mutual interference entropy spectrum density model for quantifying the interference degree between different physical field signals is established, based on the model, the mutual interference degree of each physical field combination is quantitatively characterized and an adaptive scheduling strategy is generated; according to the adaptive scheduling strategy, the excitation sources of each physical field are controlled to be excited in time sharing or in parallel cooperation, and the response signals of the multi-modal sensor array are synchronously collected to obtain a multi-source cooperative acquisition data set; Step 3, the multi-source cooperative acquisition data set is subjected to data quality verification and space-time consistency correction, the physical field characteristics including seismic wave first arrival time, apparent resistivity, Cole-Cole polarization parameters, GPR reflection interface depth, optical fiber strain, acoustic emission energy and normalized drilling speed index are extracted from the corrected data, a multi-modal feature vector comprising the feature vectors corresponding to the physical field characteristics is constructed, the reliability weight of each physical field characteristic is calculated based on the mutual interference entropy and the observation consistency residual, the reliability weight corresponding to each physical field characteristic is formed into a data weight matrix, and each multi-modal feature vector is weighted and fused to obtain a unified fusion feature vector; Step 4, according to the unified fusion feature vector and the data weight matrix, a joint inversion objective function embedded with Gassmann equation and Archie formula is constructed, the data weight matrix is applied to the weighted residual term of the joint inversion objective function, the joint inversion objective function is solved by Gauss-Newton iterative algorithm to obtain a three-dimensional surrounding rock physical parameter field comprising wave velocity, resistivity, porosity and strain field; based on the strain field in the three-dimensional surrounding rock physical parameter field and the acoustic emission energy extracted in step 3, a permeability dynamic evolution function with strain, strain rate and acoustic emission energy as independent variables is constructed, the initial permeability in the three-dimensional surrounding rock physical parameter field is spatio-temporally dynamically updated, a dynamic permeability field is calculated to finally form a three-dimensional surrounding rock physical model comprising wave velocity, resistivity, porosity and permeability; Step 4, according to the unified fusion feature vector and the data weight matrix, a joint inversion objective function embedded with Gassmann equation and Archie formula is constructed, the data weight matrix is applied to the weighted residual term of the joint inversion objective function, the joint inversion objective function is solved by Gauss-Newton iterative algorithm to obtain a three-dimensional surrounding rock physical parameter field comprising wave velocity, resistivity, porosity and strain field; based on the strain field in the three-dimensional surrounding rock physical parameter field and the acoustic emission energy extracted in step 3, a permeability dynamic evolution function with strain, strain rate and acoustic emission energy as independent variables is constructed, the initial permeability in the three-dimensional surrounding rock physical parameter field is spatio-temporally dynamically updated, a dynamic permeability field is calculated to finally form a three-dimensional surrounding rock physical model comprising wave velocity, resistivity, porosity and permeability; Step 5, using the three-dimensional surrounding rock physical property model as the initial background field, the real-time collected while-drilling parameters, acoustic emission events and optical fiber strain monitoring data during construction are assimilated into the three-dimensional surrounding rock physical property model by using the ensemble Kalman filter algorithm, and the physical property parameter field in the three-dimensional surrounding rock physical property model is dynamically updated; based on the updated physical property parameters, the basic risk probability of each spatial position and time is calculated by using Logistic regression, the physical property parameter uncertainty tensor representing the uncertainty of the resistivity, porosity and permeability model is constructed based on the updated physical property parameter field, the Frobenius norm of the physical property parameter uncertainty tensor is calculated, and the relative increment between the dynamic permeability and the initial permeability is combined to determine the seepage-uncertainty coupling correction coefficient, the seepage-uncertainty coupling correction coefficient is used as a multiplicative factor to act on the basic risk probability, and the basic risk probability is dynamically corrected in space and time to obtain the final risk probability, and the surrounding rock stability is evaluated according to the final risk probability value.

2. The tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion according to claim 1, characterized in that, In step 2, the multi-field mutual interference entropy spectrum density model quantifies the interference degree between different physical field signals by the following formula: ； wherein, represents the normalized mutual interference energy ratio of the first physical field and the first physical field signal in time , frequency , , represents the cross-correlation spectral density of the first physical field and the first physical field signal in time , frequency , represents the mutual interference entropy of the first physical field in time , frequency , the greater the value of the mutual interference entropy represents the greater the degree of interference from other physical fields.

3. The tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion according to claim 1, characterized in that, In step 3, the physical field characteristics are extracted from the corrected data, and the calculation is realized by the following formula: The first arrival time of the seismic wave is automatically identified by using the AIC criterion: ； wherein denotes the possible seismic wave first arrival time, i.e. the sample point in the data sequence, denotes the AIC criterion value calculated when assuming the first arrival time to be the sample point, denotes the variance of the first samples, denotes the variance of the th to the th sample, denotes the total number of samples within the time window; The calculation formula of the apparent resistivity is as follows: ； wherein, represents the apparent resistivity, represents the device coefficient, represents the measured potential difference, represents the injected current; The Cole-Cole polarization parameters include a time constant, which is obtained by fitting the complex resistivity dispersion curve by using the Cole-Cole model: ； wherein, represents complex resistivity, represents zero frequency resistivity, represents chargeability, represents time constant, represents frequency dependent coefficient, represents imaginary unit, represents angular frequency; The GPR reflection interface depth is calculated by using the time-depth conversion model: ， ； wherein, represents the GPR reflection interface depth, represents the propagation velocity of the electromagnetic wave in the rock mass, represents the two-way propagation time of the electromagnetic wave, represents the speed of light in a vacuum, represents the relative dielectric constant of the rock mass; The optical fiber strain is calculated by using the phase-strain conversion formula: ； wherein, represents a measured phase change amount, represents an effective refractive index of the optical fiber, represents an induced length of the optical fiber, represents a wavelength of the laser light source, represents a strain of the optical fiber; The acoustic emission energy is calculated by using energy integration: ； wherein represents the acoustic emission energy, i.e. the cumulative energy of the acoustic emission event, represents the acoustic emission signal amplitude at the discrete time point represents the sampling time interval;​ The acoustic emission event is detected by using the short-time window and long-time window ratio algorithm: ； wherein, represents a ratio of the short-time window and the long-time window amplitude at the time instant, the ratio being greater than a predetermined threshold value, the ratio being determined as an effective acoustic emission event, represents the number of sample points within the short-time window, represents the number of sample points within the long-time window, represents the sample point index within the short-time window, represents the sample point index within the long-time window, represents the time of the acoustic emission signal amplitude, represents the time of the acoustic emission signal amplitude; The normalized drilling speed index is calculated by using drilling parameters: ； wherein, RPI represents a normalized rate of penetration index, ROP represents a rate of penetration, WOB represents a weight on bit, RPM represents a bit rotation speed.

4. The tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion according to claim 1, characterized in that, In step 3, the calculation formula of the credibility weight is as follows: ； wherein, represents the confidence weight of the th physical field characteristic in spatial position , time , represents the mutual interference entropy of the th physical field in time , frequency , represents the consistency residual of the th physical field observation data, and represents the empirical adjustment factor, represents a minimum constant; The multi-modal feature vectors are weighted and fused, and the formula is as follows: ； wherein represents a unified fusion feature vector at a spatial position , time , represents a feature vector of the th physical field at a spatial position , time , represents a feature normalization operator.

5. The tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion according to claim 4, characterized in that, In step 4, the mathematical expression of the joint inversion objective function is as follows: ； in, Let represent the joint inversion objective function, and represent the field vector of surrounding rock physical parameters to be inverted. Indicates the first The forward response vector of a physical field, i.e., the vector of the surrounding rock physical property parameters. The theoretical response obtained from the calculation This represents a data weight matrix composed of credibility weights. Represents the smoothing constraint coefficient. This represents the model smoothness difference operator. Represents the structural coupling constraint coefficient. Represents mutual information weights, ,in, Represents the mutual information computation operator, used for quantization. and The degree of spatial structural dependence between them Indicates spatial location ,time The set of surrounding rock physical properties, and Represents the spatial gradient of fields with different physical property parameters. This represents the L2 norm.

6. The tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion according to claim 1, characterized in that, In step 4, the Gassmann equation is as follows: ； wherein, B represents the bulk modulus of the fluid-saturated rock, Bdry represents the bulk modulus of the dry rock matrix, Bm represents the bulk modulus of the rock matrix mineral, Bf represents the bulk modulus of the pore fluid, B represents the bulk modulus of the rock porosity; The Archie formula is as follows: ； wherein, represents the rock conductivity, represents the pore water conductivity, represents the water saturation, represents the lithology empirical coefficient.

7. The tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion according to claim 1, characterized in that, In step 4, the dynamic permeability field is calculated, and the corresponding formula is as follows: ； wherein represents the permeability at spatial location , time , represents the initial permeability in the petrophysical parameter field, represents the strain at spatial location , time obtained in step 4 of the joint inversion, represents the strain rate at spatial location , time , represents the acoustic emission energy at spatial location , time , , and represent the permeability coupling coefficients for strain, strain rate and acoustic emission energy, respectively, represents the natural exponential function.

8. The tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion according to claim 1, characterized in that, In step 5, the physical property parameter field in the three-dimensional surrounding rock physical property model is dynamically updated, including the following processes: The physical property parameter field in the three-dimensional surrounding rock physical property model is used as the initial value, and the initial ensemble is generated by random disturbance: ； ； wherein, denotes an initial set, denotes the i-th set member in the initial set, denotes the i-th set member in the initial set, denotes the number of set members, denotes a property parameter field in the three-dimensional surrounding rock property model, and corresponding property parameters include wave velocity, resistivity, porosity, and permeability, denotes a random perturbation vector of the i-th set member, which is subject to a multi-dimensional normal distribution with a mean of 0 and a covariance of ;​​ In each assimilation period, based on the geology-mechanics-seepage coupling model, each ensemble member is forward evolution predicted, and the evolution equation is as follows: ； wherein, denotes the forecast value of the ensemble member at time step , denotes a geomechanical-seepage coupling model, denotes the assimilated value of the ensemble member at time step , denotes the process noise of the ensemble member at time step , After obtaining new observation data including while-drilling parameters, acoustic emission events and optical fiber strain monitoring data during construction, the model state is mapped to the observation space through the observation operator, the Kalman gain matrix is calculated, and the new observation data is assimilated into the model by using the Kalman gain matrix: ； wherein, represents the analysis value after assimilation of the th ensemble member, represents the forecast value of the th ensemble member, represents the perturbed observation data of the th ensemble member, represents the theoretical observation value corresponding to the forecast value of the th ensemble member, represents the Kalman gain matrix, , represents the forecast error covariance matrix, represents the observation operator, represents the observation error covariance matrix, represents the matrix transpose.

9. The tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion according to claim 8, characterized in that, In step 5, the calculation formula of the basic risk probability is as follows: ； wherein, represents the base risk probability at spatial location , time , represents a logistic regression model, and represent the resistivity and porosity, respectively, of the th ensemble member in the updated petrophysical parameter field at spatial location , and represent the permeability and pore pressure change, respectively, of the th ensemble member in the updated petrophysical parameter field at spatial location , time , represents the coefficients of the logistic regression model.

10. The tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion according to claim 8, characterized in that, In step 5, the seepage-uncertainty coupling correction coefficient is introduced to dynamically correct the basic risk probability, and the correction is realized by the following formula: ； ； wherein, denotes the basis risk probability at spatial location , time , denotes the final risk probability after correction of the basis risk probability , denotes the seepage-uncertainty coupling correction factor, denotes the Frobenius norm of the uncertainty tensor, wherein, , and denote the model uncertainty of resistivity, porosity and permeability, respectively, denotes the permeability at spatial location , time , denotes the initial permeability, denotes the average uncertainty of the full set of petrophysical parameters. 11.The tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion according to claim 10, characterized in that, In step 5, the surrounding rock stability evaluation includes confidence interval calculation of the final risk probability, and the calculation formula is as follows: ； wherein, indicates the 95% confidence interval for the spatial location , time .

12. A tunnel surrounding rock stability evaluation system based on multi-physical field parameter inversion, characterized in that, The system for implementing the tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion according to any one of claims 1 to 11 comprises: A multi-modal sensor array is arranged in a deep-buried tunnel and contains six types of sensors, namely, acoustic, optical, electrical, magnetic, seismic, and drilling sensors, and is calibrated in space coordinates and synchronized in time to form original multi-physical field signals with unified time-space reference; A scheduling control module is configured to establish a multi-field mutual interference entropy spectral density model for quantifying the interference degree between different physical field signals, to quantitatively characterize the mutual interference degree of each physical field combination based on the multi-field mutual interference entropy spectral density model and generate an adaptive scheduling strategy, and to control the excitation sources of each physical field to perform time-sharing or parallel cooperative excitation according to the adaptive scheduling strategy and synchronously collect the response signals of the multi-modal sensor array to obtain a multi-source cooperative acquisition data set; A data processing module is configured to perform data quality verification and space-time consistency correction on the multi-source cooperative acquisition data set, to extract physical field features including seismic wave first arrival time, apparent resistivity, Cole-Cole polarization parameters, GPR reflection interface depth, optical fiber strain, acoustic emission energy, and normalized drilling speed index from the corrected data, to construct multi-modal feature vectors containing the corresponding feature vectors of the physical field features, to calculate the credibility weight of each physical field feature based on mutual interference entropy and observation consistency residual, to form a data weight matrix with the credibility weight corresponding to each physical field feature, and to perform weighted fusion on each multi-modal feature vector to obtain a unified fusion feature vector; An inversion module is configured to construct a joint inversion objective function embedded with Gassmann equation and Archie formula based on the unified fusion feature vector and the data weight matrix, to apply the data weight matrix to the weighted residual term of the joint inversion objective function, to solve the joint inversion objective function through a Gauss-Newton iterative algorithm, and to obtain a three-dimensional surrounding rock physical property parameter field including wave velocity, resistivity, porosity, and strain field; to construct a permeability dynamic evolution function with strain, strain rate, and acoustic emission energy as independent variables based on the strain field in the three-dimensional surrounding rock physical property parameter field and the extracted acoustic emission energy, to perform space-time dynamic updating on the initial permeability in the three-dimensional surrounding rock physical property parameter field, to calculate a dynamic permeability field, and to finally form a three-dimensional surrounding rock physical property model including wave velocity, resistivity, porosity, and permeability. The dynamic updating and surrounding rock stability evaluation module is used for taking the three-dimensional surrounding rock physical property model as an initial background field, using a set Kalman filtering algorithm, assimilating the real-time collected while-drilling parameters, acoustic emission events and optical fiber strain monitoring data in the construction process into the three-dimensional surrounding rock physical property model, and dynamically updating the physical property parameter field in the three-dimensional surrounding rock physical property model; based on the updated physical property parameters, a Logistic regression is used to calculate the basic risk probability of each spatial position and time, a physical property parameter uncertainty tensor is constructed based on the updated physical property parameter field to represent the uncertainty of the resistivity, porosity and permeability models, the Frobenius norm of the physical property parameter uncertainty tensor is calculated, and in combination with the relative increment between the dynamic permeability and the initial permeability, a seepage-uncertainty coupling correction coefficient is determined, the seepage-uncertainty coupling correction coefficient is taken as a multiplicative factor to act on the basic risk probability, the basic risk probability is dynamically corrected in time and space, the final risk probability is obtained, and the surrounding rock stability is evaluated according to the final risk probability value.

13. A computer-readable storage medium, characterized in that, The computer readable storage medium stores a computer program, when the computer program is run, the steps of the tunnel surrounding rock stability evaluation method based on multi-physical field parameter inversion as claimed in any one of claims 1 to 11 are realized.

Images