A tunnel deformation analysis method and system based on multi-source data fusion

Abstract

This invention discloses a tunnel deformation analysis method and system based on multi-source data fusion. The method includes: calculating and sorting the loosened zone expansion index after correction by the surrounding rock integrity coefficient; selecting cross-sections in sequence for mechanical state inversion to obtain the stress release coefficient and plastic zone radius, and selecting key cross-sections based on the stress release coefficient exceeding a threshold; calculating instability characteristic values ​​using catastrophe theory based on the inversion parameters, and constructing a spatial distribution matrix of instability characteristic values ​​by combining the distance between cross-sections; performing principal component analysis on the matrix, extracting the score of the first principal component as a comprehensive instability index, and generating a longitudinal instability risk curve; identifying segments on the curve with slopes exceeding a threshold as potential instability expansion segments, and identifying weak points in the support within these segments based on the depth-to-displacement ratio acceleration. This method achieves quantitative description of the surrounding rock loosening range, mechanical parameter inversion, instability risk quantification, and precise location of weak points, significantly improving the intelligence level of tunnel stability analysis.

Description

Technical Field

[0001] This invention belongs to the field of tunnel engineering monitoring and safety assessment technology, and in particular relates to a tunnel deformation analysis method and system based on multi-source data fusion. Background Technology

[0002] During tunnel construction, surrounding rock deformation is a core indicator for evaluating tunnel stability. Currently, tunnel deformation monitoring mainly involves deploying sensors at multiple monitoring sections to acquire data such as crown settlement, perimeter convergence, internal displacement of the surrounding rock, and stress on the support structure.

[0003] However, existing analytical methods typically process data from each monitoring section independently, neglecting the correlation between the surrounding rock mechanical behavior of different sections, making it difficult to grasp the overall deformation evolution of the tunnel. Furthermore, traditional methods lack effective characterization of the expansion process of the loosened zone in the surrounding rock, failing to promptly identify the intrinsic connection between the degree of stress release and the development of the plastic zone, resulting in delayed early warnings of tunnel instability risks. Summary of the Invention

[0004] This invention aims to address the problems of insufficient data fusion, delayed identification of instability risks, and difficulty in locating weak points in support in existing tunnel deformation analysis technologies. It provides a tunnel deformation analysis method and system based on multi-source data fusion, which can screen key sections, invert stress release coefficients and plastic zone radii through loosening zone expansion index, evaluate instability characteristic values ​​using catastrophe theory, and generate longitudinal instability risk curves by combining spatial distribution matrix and principal component analysis, ultimately locating weak points in support.

[0005] In a first aspect, the present invention provides a tunnel deformation analysis method based on multi-source data fusion, comprising:

[0006] Acquire multi-source monitoring data from multiple monitoring sections within the target tunnel section. The multi-source monitoring data includes the arch settlement time series data, the displacement time series data inside the surrounding rock, and the stress time series data of the support structure for each monitoring section.

[0007] The expansion index of the loosened zone of the surrounding rock is calculated based on the proportional relationship between the internal displacement of the surrounding rock and the settlement of the arch at the monitoring section. The monitoring sections are then sorted according to the expansion index of the loosened zone of the surrounding rock to obtain the loosening trend sequence of the sections.

[0008] In the loosening trend sequence of the cross section, the cross sections to be analyzed and monitored are selected in order of the expansion index of the loosening zone of the surrounding rock from large to small to perform mechanical state inversion, so as to obtain the stress release coefficient and plastic zone radius of each cross section to be analyzed and monitored. The selection stops when the stress release coefficient of the surrounding rock obtained by the inversion exceeds the preset stress threshold. All the cross sections to be analyzed and monitored at this time are determined as the set of key monitoring sections.

[0009] Based on the stress release coefficient and plastic zone radius of each key monitoring section, the catastrophe theory model is used to calculate the instability characteristic value of the corresponding key monitoring section, and the spatial distribution matrix of the instability characteristic value is constructed according to the instability characteristic value and spatial distance of each key monitoring section.

[0010] Principal component analysis is performed on the spatial distribution matrix of the instability eigenvalues, and the score of the first principal component is extracted as the comprehensive instability index for each monitoring section. Based on the comprehensive instability index, a longitudinal instability risk curve is generated.

[0011] On the longitudinal instability risk curve, segments with slopes exceeding a preset slope threshold are identified as potential instability propagation segments, and weak points in the support are determined within these potential instability propagation segments based on the rate of change of the internal displacement ratio of the surrounding rock.

[0012] Secondly, the present invention provides a tunnel deformation analysis system based on multi-source data fusion, comprising:

[0013] The acquisition module is configured to acquire multi-source monitoring data from multiple monitoring sections within the target tunnel section. The multi-source monitoring data includes the arch settlement time series data, the displacement time series data inside the surrounding rock, and the stress time series data of the support structure for each monitoring section.

[0014] The first calculation module is configured to calculate the expansion index of the loosened zone of the surrounding rock based on the ratio of the internal displacement of the surrounding rock to the settlement of the arch at the monitoring section, and to sort each monitoring section according to the expansion index of the loosened zone of the surrounding rock to obtain the loosening trend sequence of the section.

[0015] The inversion module is configured to select the monitoring sections to be analyzed in the sequence of loosening trend of the cross section in order of the expansion index of the loosening zone of the surrounding rock from large to small, and perform mechanical state inversion to obtain the stress release coefficient and plastic zone radius of each monitoring section to be analyzed. The selection stops when the stress release coefficient of the surrounding rock obtained by the inversion exceeds the preset stress threshold, and all the monitoring sections to be analyzed at this time are determined as the set of key monitoring sections.

[0016] The second calculation module is configured to calculate the instability characteristic values ​​of the corresponding key monitoring sections based on the surrounding rock stress release coefficient and plastic zone radius of each key monitoring section using the catastrophe theory model, and construct a spatial distribution matrix of instability characteristic values ​​based on the instability characteristic values ​​and spatial distance of each key monitoring section.

[0017] The generation module is configured to perform principal component analysis on the spatial distribution matrix of the instability feature values, extract the first principal component score as the comprehensive instability index of each monitoring section, and generate a longitudinal instability risk curve based on the comprehensive instability index.

[0018] The determination module is configured to identify sections on the longitudinal instability risk curve whose slope exceeds a preset slope threshold as potential instability propagation sections, and to determine the weak points of the support within the potential instability propagation sections based on the rate of change of the internal displacement ratio of the surrounding rock.

[0019] Thirdly, an electronic device is provided, comprising: at least one processor, and a memory communicatively connected to the at least one processor, wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform the steps of the tunnel deformation analysis method based on multi-source data fusion according to any embodiment of the present invention.

[0020] Fourthly, the present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein when the program instructions are executed by a processor, the processor performs the steps of the tunnel deformation analysis method based on multi-source data fusion according to any embodiment of the present invention.

[0021] This application presents a tunnel deformation analysis method and system based on multi-source data fusion. It calculates and sorts the loosening zone expansion index, corrected for the surrounding rock integrity coefficient, based on the ratio of internal displacement of the surrounding rock to the arch settlement. Sections are selected sequentially for mechanical state inversion to obtain stress release coefficients and plastic zone radii. Key sections are selected based on stress release coefficient exceeding a threshold. Instability characteristic values ​​are calculated using catastrophe theory based on the inversion parameters, and a spatial distribution matrix of instability characteristic values ​​is constructed by combining the distance between sections. Principal component analysis is performed on the matrix, and the score of the first principal component is extracted as a comprehensive instability index to generate a longitudinal instability risk curve. Sections on the curve with slopes exceeding a threshold are identified as potential instability expansion sections. Within these sections, weak points in the support are identified based on the depth-to-displacement ratio acceleration. This achieves quantitative description of the surrounding rock loosening range, mechanical parameter inversion, instability risk quantification, and precise weak point location, significantly improving the intelligence level of tunnel stability analysis. Attached Figure Description

[0022] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0023] Figure 1 A flowchart illustrating a tunnel deformation analysis method based on multi-source data fusion, as provided in an embodiment of the present invention;

[0024] Figure 2 A structural block diagram of a tunnel deformation analysis system based on multi-source data fusion is provided in an embodiment of the present invention;

[0025] Figure 3 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention. Detailed Implementation

[0026] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0027] Please see Figure 1 The diagram shows a flowchart of a tunnel deformation analysis method based on multi-source data fusion according to this application.

[0028] like Figure 1 As shown, the tunnel deformation analysis method based on multi-source data fusion specifically includes the following steps:

[0029] Step S101: Obtain multi-source monitoring data from multiple monitoring sections within the target tunnel section. The multi-source monitoring data includes the arch settlement time series data, the displacement time series data inside the surrounding rock, and the stress time series data of the support structure for each monitoring section.

[0030] In this step, within the target tunnel section, several representative locations are selected to deploy monitoring sections based on factors such as geological conditions, tunnel depth, support type, and construction progress. The spacing between monitoring sections is generally 5-30m, with increased spacing in sections with complex geological conditions or significant depth variations. Each monitoring section should be equipped with the following three types of sensors:

[0031] Settlement monitoring points at the tunnel arch: Settlement monitoring points are buried at the center line of the tunnel arch, and non-contact measurements are taken using a precision level or total station, or automated monitoring is carried out using a static level.

[0032] Displacement monitoring points inside the surrounding rock: Drill holes at the arch or sidewall location, with a hole depth generally 1.5 to 2 times the tunnel diameter (usually 5 to 10 m). Multiple displacement gauges are installed at different depths inside the holes (such as 1 m, 2 m, 3 m, 5 m, etc.) to monitor displacement changes at different depths inside the surrounding rock.

[0033] Stress monitoring points for support structures: Pressure cells, rebar gauges or strain gauges are embedded in the initial support (such as shotcrete, steel arch) or secondary lining to monitor the contact pressure between the support structure and the surrounding rock, the axial force of the rebar or the strain of the concrete.

[0034] All sensors collect data at a preset sampling frequency (usually once a day, but this can be increased to 2-4 times a day during critical construction phases). The collected data is uploaded to the data center server in real time via wired or wireless transmission, forming the raw monitoring database. To ensure data quality, simple outlier removal logic can be set at the acquisition end, such as removing invalid data that exceeds the sensor's measurement range.

[0035] Step S102: Calculate the expansion index of the loosened zone of the surrounding rock based on the proportional relationship between the internal displacement of the surrounding rock and the settlement of the arch at the monitoring section, and sort each monitoring section according to the expansion index of the loosened zone of the surrounding rock to obtain the loosening trend sequence of the section.

[0036] In this step, for each monitoring section, at least three displacement measurement points inside the surrounding rock at different depths are selected, and the displacement values ​​of each depth measurement point at each monitoring time are extracted. The ratio of the displacement value of each depth measurement point at each monitoring time to the crown settlement value at the corresponding time is calculated to obtain the original sequence of displacement ratios of each depth measurement point.

[0037] Kalman filtering was applied to the original sequence of displacement ratios at each depth measurement point to remove measurement noise and random fluctuations, resulting in a smoothed sequence of displacement ratios at each depth measurement point.

[0038] The first-order numerical derivative of the smoothed displacement ratio sequence at each depth measuring point is performed to obtain the displacement ratio change rate sequence at each depth measuring point. The sign of the displacement ratio change rate sequence is determined, and the time period with positive change rate is marked.

[0039] The longest duration during which the displacement ratio change rate of each depth measuring point is continuously positive is statistically analyzed, and the depth measuring points whose longest duration exceeds a first preset duration threshold are identified as potential influence measuring points of the loosening ring.

[0040] From all the potential impact measurement points of the loosening zone, the continuous depth intervals where the average displacement ratio change rate exceeds the preset change rate threshold are selected, and the continuous depth intervals are defined as the loosening zone expansion zone.

[0041] The thickness of the loosened zone expansion area is measured, and the thickness is divided by the equivalent radius of the tunnel diameter of the corresponding monitoring section to obtain the initial loosened zone expansion index;

[0042] The surrounding rock integrity coefficient of the monitoring section is obtained, and the initial loosened zone expansion index is corrected according to the surrounding rock integrity coefficient to obtain the final surrounding rock loosened zone expansion index.

[0043] All monitoring sections are arranged in descending order of the surrounding rock loosening zone expansion index to generate a section loosening trend sequence.

[0044] In one specific embodiment, for each monitoring section, multiple displacement gauges are pre-deployed at different depths within the surrounding rock. In this embodiment, at least three representative depth measuring points are selected, for example, at distances of 1.0m, 2.0m, and 3.0m from the excavation face. Assume there are K depth measuring points for a given monitoring section, and the displacement value of the k-th depth measuring point at monitoring time t is d. k(t) At the same moment, the settlement value of the vault is s. (t) For each depth measuring point k, calculate the displacement ratio r at each monitoring time. k(t) :

[0045] ,

[0046] This yields the displacement ratio of the k-th depth measurement point in the original sequence. Where N is the total number of monitoring times, the displacement ratio reflects the degree of deformation coordination between the surrounding rock at that depth and the overall arch settlement. The larger the ratio, the higher the proportion of the surrounding rock participating in the deformation at that depth, which is a basic indicator for identifying the range of the loosening zone.

[0047] On-site monitoring data inevitably contains measurement noise and random fluctuations. Directly using the original sequence will lead to unstable calculation of the rate of change. In this embodiment, Kalman filtering is used to reduce noise in the displacement ratio of each depth measurement point to the original sequence.

[0048] Establish the state-space model of the Kalman filter:

[0049] Equations of state: ,

[0050] Observation equation: ,

[0051] in, The ratio of the true displacement of the k-th depth measuring point at time t. This represents the true displacement ratio of the k-th depth measurement point at time t-1. Let be the observed displacement ratio of the k-th depth measuring point at time t. Process noise represents random disturbances that occur during the evolution of the system state over time. For the observation noise, Q and R are the process noise covariance and the observation noise covariance, respectively, which can be preset according to the statistical characteristics of the monitoring data.

[0052] Through Kalman filtering prediction-update iteration, a smooth displacement ratio sequence of each depth measurement point is obtained. The filtered sequence eliminates high-frequency noise and retains the true deformation evolution trend.

[0053] The first-order numerical derivative of the displacement ratio smoothing sequence at each depth measuring point is performed to obtain the displacement ratio change rate sequence. This embodiment uses the central difference method to improve accuracy for monitoring time. :

[0054] ,

[0055] In the formula, For the first Each depth measuring point at the monitoring time The rate of change of displacement ratio For the first Each depth measuring point at the monitoring time The smooth displacement ratio, For the first Each depth measuring point at the monitoring time The smooth displacement ratio, For the first Each monitoring moment, For the first Each monitoring moment;

[0056] For the first and last moments, forward or backward differencing is used. The rate of change of displacement ratio reflects the rate of deformation of the surrounding rock at that depth. The sign of the rate of change sequence is determined, and the periods in the rate of change of displacement ratio greater than 0 are marked, that is, the periods when the displacement ratio continues to increase. This means that the deformation proportion of the surrounding rock at that depth is expanding, which is a direct indication of the expansion of the loosened zone.

[0057] For each depth measurement point, record the longest duration during which the rate of change of displacement ratio remains positive. The longer the duration, the longer the surrounding rock at that depth has been in a state of continuous loosening. A first preset duration threshold is set. (In this embodiment, T_th is set to 7 days, which can be adjusted within the range of 3 to 14 days according to the actual project situation.) > The depth measurement points were identified as potential influence points of the loosening zone. The depth range of these measurement points may be the area affected by the expansion of the loosening zone.

[0058] From all potential influence points of the loosened zone, a set of continuous measurement points along the depth direction is selected. For each continuous depth interval, the mean displacement ratio change rate of all measurement points within that interval is calculated. Set a preset rate of change threshold. (In this embodiment, V_th = 0.02 / day, meaning the daily displacement increases by 2%). > The continuous depth interval is defined as the loosening zone expansion region. This screening criterion ensures that the identified region not only has a long duration but also a significant expansion rate.

[0059] Measure the thickness H of the loosened zone expansion area, which is the depth difference between the deepest and shallowest measuring points within this zone. Obtain the equivalent radius of the tunnel diameter for this monitoring section. For circular tunnels, This refers to the tunnel radius; for non-circular tunnels, the hydraulic radius can be used for calculation. Where A is the tunnel cross-sectional area and P is the cross-sectional perimeter. Initial loosening zone expansion index. Calculate using the following formula:

[0060] ,

[0061] This index reflects the relative size of the loosened zone expansion thickness with respect to the tunnel dimensions, eliminating the influence of cross-sectional geometry and making the degree of loosening comparable across different cross-sections.

[0062] Obtain the surrounding rock integrity coefficient of the monitoring section. The surrounding rock integrity coefficient can be obtained through acoustic wave testing. It is the square of the ratio of the longitudinal wave velocity of the rock mass to the longitudinal wave velocity of the rock, and its value ranges from 0 to 1. The smaller the integrity coefficient, the more fragmented the surrounding rock, and the higher the risk of the loosened zone for the same thickness. Therefore, this embodiment uses the following modified formula to calculate the final loosened zone expansion index. :

[0063] ,

[0064] in, As a correction factor, in this embodiment, we take... It can be adjusted within the range of 1.0 to 3.0 based on engineering experience. When (integrity of surrounding rock), the correction term is 0. ;when At that time (extremely broken surrounding rock). It fully considers the amplifying effect of surrounding rock quality on the degree of harm to the loosened zone.

[0065] Repeat the above steps for all monitoring sections to calculate the final expansion index of the loosened zone of the surrounding rock for each section. All monitoring sections were arranged according to... Arrange the sections in descending order to generate a sequence of loosening trends. The sections ranked first in this sequence are those with the most significant expansion of the loosened zone and the highest degree of loosening of the surrounding rock, and will be given priority for mechanical state inversion analysis in subsequent steps.

[0066] In summary, the method in this embodiment extends the analytical perspective to different depths within the surrounding rock. By analyzing the displacement ratio and its evolution characteristics, it can accurately identify the starting position, expansion range, and dynamic changes of the loosened zone, providing deeper physical information for surrounding rock stability assessment. The dual threshold screening mechanism of "duration of positive time period" and "mean rate of change" avoids misjudgments that may occur with a single threshold, accurately distinguishing between the true expansion of the loosened zone and accidental fluctuations, ensuring the accuracy of loosened zone expansion identification. The introduction of a surrounding rock integrity coefficient to correct the initial index fully considers the amplification effect of the surrounding rock's own mass on the degree of harm to the loosened zone, ensuring that the final loosened zone expansion index not only reflects the geometric dimensions of the loosened zone but also incorporates the influence of the surrounding rock's mechanical properties, possessing a more comprehensive physical meaning.

[0067] Step S103: In the loosening trend sequence of the cross section, the cross sections to be analyzed and monitored are selected in order of decreasing expansion index of the loosening zone of the surrounding rock to be analyzed and monitored, and the mechanical state inversion is performed to obtain the stress release coefficient and plastic zone radius of each cross section to be analyzed and monitored. The selection stops when the stress release coefficient of the surrounding rock obtained by the inversion exceeds the preset stress threshold. All the cross sections to be analyzed and monitored at this time are determined as the set of key monitoring cross sections.

[0068] In this step, the current elastoplastic mechanical model of the current monitoring section to be analyzed is established;

[0069] Obtain the current crown settlement value at the latest moment in the crown settlement time series data of the current monitoring section to be analyzed, and use the current crown settlement value as the measured constraint value of the radial displacement of the tunnel wall;

[0070] Within the preset stress relief coefficient search range, the bisection method is used for iterative solution. In each iteration, the current stress relief coefficient is substituted into the current elastoplastic mechanical model to calculate the corresponding theoretical cavity wall displacement and theoretical plastic zone radius.

[0071] The relative error between the theoretical tunnel wall displacement and the measured tunnel wall displacement is calculated. If the relative error is greater than a preset accuracy threshold, the range of the stress relief coefficient is adjusted according to the sign of the error, and the iteration continues.

[0072] When the relative error between the theoretical tunnel wall displacement and the measured tunnel wall displacement is less than or equal to the preset accuracy threshold, the iteration stops, and the current stress release coefficient and the corresponding current theoretical plastic zone radius are used as the inversion result of the section to be analyzed and monitored.

[0073] The next monitoring section to be analyzed is inverted until the stress release coefficient of the surrounding rock obtained by the inversion of a certain monitoring section exceeds the preset stress threshold. Then the selection and inversion are stopped, and all the monitoring sections to be analyzed that have been inverted are identified as the set of key monitoring sections.

[0074] In one specific embodiment, for the currently selected monitoring section to be analyzed, an elastoplastic mechanical model of a circular tunnel is established based on its geometric dimensions, burial depth, and surrounding rock mechanical parameters. This embodiment adopts a modified form of the Fenner formula based on the Mohr-Coulomb yield criterion. The basic assumptions of the model are: homogeneous isotropic surrounding rock, initial geostress field under hydrostatic pressure, and stress redistribution of surrounding rock after tunnel excavation, forming elastic and plastic zones.

[0075] The model input parameters include:

[0076] tunnel radius (Unit: m), determined by the cross-sectional design dimensions;

[0077] Initial geostress of surrounding rock (Unit: MPa), estimated based on burial depth and rock mass. Where γ is the rock mass density, For burial depth;

[0078] elastic modulus of surrounding rock (Unit: GPa), provided by geological survey report or obtained from experience;

[0079] Cohesion of surrounding rock (Unit: MPa) and internal friction angle (Unit: °), provided by the geological survey report;

[0080] Support resistance provided by the support structure (Unit: MPa), calculated based on actual support parameters or taken from experience (take 0 if no support is applied).

[0081] To describe the spatial effect of the excavation face, a stress relief coefficient is introduced. (0≤) ≤1) indicates the degree of stress release in the surrounding rock at the current monitoring time. =0 corresponds to the unexcavated state. =1 corresponds to the completion of full-section excavation and complete stress release. Theoretical tunnel wall displacement. and The functional relationship can be expressed as (taking Fenner's formula as an example):

[0082] when ≤ At that time, the surrounding rock is in an elastic state, and the theoretical displacement of the tunnel wall is... The expression is: ;

[0083] when > When a plastic zone appears in the surrounding rock, the radius of the plastic zone is... and theoretical tunnel wall displacement They are respectively:

[0084] ,

[0085] ,

[0086] in, The Poisson's ratio of the surrounding rock is provided in the exploration report.

[0087] The above formulas constitute the elastoplastic mechanical model of the current cross-section, which can be based on the input... Calculate the corresponding theoretical tunnel wall displacement and radius of plastic zone .

[0088] From the time series data of the crown settlement of the current monitoring section, extract the cumulative settlement value at the latest monitoring time, and denot it as... This value represents the measured radial displacement of the tunnel wall at the current moment and will serve as the target constraint for the inversion.

[0089] Stress relief coefficient The physical meaning of determines that its value range is [0,1]. Therefore, the initial search interval is set to [ , The interval is [0,1]. A binary search method is used for iteration within this interval. The specific steps are as follows:

[0090] Take the midpoint of the current interval =( + ) / 2.

[0091] Will Substitute into the elastoplastic mechanics model and calculate the midpoint of the current interval. Corresponding theoretical tunnel wall displacement and the radius of the theoretical plastic zone .

[0092] Calculate the relative error between theoretical displacement and measured displacement. .

[0093] Set a preset accuracy threshold ε (in this embodiment, ε = 0.01, meaning the relative error does not exceed 1%). If If ≤ε, then the current value is considered to be ≤ε. Stop iteration if the accuracy requirement is met; >ε, then according to and Adjust the search range based on the size relationship:

[0094] like > ,illustrate The value is too large and should be reduced. Therefore, the search range is updated to [ , ];

[0095] like < ,illustrate The value is too small and should be increased. Therefore, the search range is updated to [ , ].

[0096] When the iteration converges, the value at this point will be... As the stress release coefficient of the surrounding rock at this cross section, the midpoint of the current interval is used. The corresponding theoretical plastic zone radius is used as the plastic zone radius of the cross section. These two parameters characterize the degree of unloading of the cross section and the range of surrounding rock disturbance, respectively.

[0097] Following the sequence of loosening trends in the cross-sections, the process iterates repeatedly for the next monitoring cross-section to be analyzed. After each cross-section is inverted, its stress release coefficient is recorded. When the stress release coefficient of a cross-section exceeds a preset stress threshold (in this embodiment, the preset stress threshold is 0.8, indicating that 80% of the surrounding rock stress has been released, entering a high-risk state), the selection and inversion process immediately stops. At this point, all monitoring cross-sections that have completed the inversion (including all cross-sections preceding this one) are identified as the set of key monitoring cross-sections. These cross-sections represent the key locations within the tunnel section with the highest degree of loosening, the most significant stress release, and the greatest likelihood of instability, and will undergo further instability risk assessment in subsequent steps.

[0098] In this embodiment, a stress release coefficient threshold is introduced as the stopping condition for screening key sections, which realizes the automatic identification of high-risk sections. Compared with the traditional method of analyzing all sections at the same depth, this method can dynamically adjust the analysis depth according to the actual state of the surrounding rock, focusing computational resources on truly dangerous parts and significantly improving analysis efficiency. At the same time, using the stress release coefficient exceeding the threshold as the judgment criterion has a clear physical meaning and avoids misjudgment that may be caused by simply relying on the magnitude of deformation.

[0099] Step S104: Based on the stress release coefficient and plastic zone radius of each key monitoring section, the catastrophe theory model is used to calculate the instability characteristic value of the corresponding key monitoring section, and the spatial distribution matrix of the instability characteristic value is constructed according to the instability characteristic value and spatial distance of each key monitoring section.

[0100] In this step, the tunnel radius of the key monitoring section, the support resistance provided by the support structure, and the surrounding rock stress release coefficient and plastic zone radius obtained through inversion are obtained;

[0101] Based on the theory of elastoplastic mechanics, a total potential energy function of the surrounding rock-support system at a key monitoring section is established. The total potential energy function includes the strain energy stored in the elastic zone of the surrounding rock, the plasticity consumed in the plastic zone, and the potential energy provided by the support structure.

[0102] Expand the total potential energy function using a Taylor series near the current equilibrium position, retaining up to the fourth-order terms, to obtain an approximate polynomial expression for the potential energy function.

[0103] By performing variable substitution and linear transformation on the approximate polynomial expression, the standard form of the cusp catastrophe model is obtained: Where p and q are control variables, and x is a state variable;

[0104] Extract the specific values ​​of the control variables from the standard form, and based on the cusp catastrophe theory, transform the bifurcation set equation... The calculation results are defined as the instability characteristic values ​​of the key monitoring sections.

[0105] It should be noted that the mileage coordinates of all key monitoring sections are obtained, and each key monitoring section is numbered in ascending order of mileage to obtain the section number sequence 1, 2, ..., n, where n is the total number of key monitoring sections.

[0106] Construct an n x n square matrix M, and set the diagonal element of the i-th row and i-th column of the square matrix M as the instability characteristic value λi of the i-th key monitoring section;

[0107] The expression for calculating the value of the non-diagonal element in the i-th row and j-th column of a square matrix M is:

[0108] ,

[0109] Where λi is the instability characteristic value of the i-th cross section, λj is the instability characteristic value of the j-th cross section, Lij is the spatial distance between the i-th key monitoring cross section and the j-th key monitoring cross section, and δ is a preset positive small quantity.

[0110] The calculated square matrix M is subjected to row normalization, that is, the elements of each row are divided by the sum of the elements of the current row to obtain the row normalized unstable eigenvalue spatial distribution matrix.

[0111] The column normalization process is performed on the row-normalized unstable eigenvalue spatial distribution matrix, that is, the element of each column is divided by the sum of the elements of the current column to obtain the final standardized unstable eigenvalue spatial distribution matrix.

[0112] In one specific embodiment, the mechanical parameters obtained from the inversion (stress release coefficient, plastic zone radius) are substituted into the catastrophe theory model to quantitatively assess the instability tendency of each key section. Furthermore, by constructing a spatial distribution matrix, the spatial correlation between sections is fused, laying the foundation for subsequent principal component analysis. This is specifically achieved through the following sub-steps:

[0113] Step S1041: Obtain basic parameters of key monitoring sections

[0114] For each key monitoring section, the following parameters are obtained from the design data and inversion results:

[0115] tunnel radius (Unit: m);

[0116] Support resistance provided by the support structure (Unit: MPa);

[0117] The stress release coefficient of the surrounding rock obtained by inversion in step S103 ;

[0118] The radius of the plastic zone obtained by inversion in step S103 (Unit: m);

[0119] Initial geostress of surrounding rock (Unit: MPa);

[0120] Cohesion of surrounding rock (Unit: MPa) and internal friction angle (Unit: °);

[0121] elastic modulus of surrounding rock (Unit: MPa) and Poisson's ratio of surrounding rock .

[0122] Step S1042: Establish the total potential energy function of the surrounding rock-support system.

[0123] Based on the theory of elastoplastic mechanics, the surrounding rock-support system is considered as an energy dissipation system, and its total potential energy is... It consists of three parts:

[0124] Strain energy stored in the elastic region The elastic strain energy stored in the elastic region surrounding the plastic region due to stress adjustment. According to elasticity mechanics, for a radius of... The strain energy of the elastic region surrounding the plastic region can be expressed as... Related functions;

[0125] Plasticity consumed in the plastic zone The irreversible energy dissipated due to yielding after the surrounding rock enters the plastic state is related to the volume of the plastic zone and the strength parameters of the surrounding rock.

[0126] Potential energy provided by the support structure Support resistance Work done on the cave wall can be expressed as: The product of the displacement of the tunnel wall.

[0127] Combining the above three factors, the total potential energy function can be written as follows: (Equation omitted for brevity) ... The expression for (considered as a state variable) is given in this embodiment using the classic cusp catastrophe theory model, approximating the system potential energy function as follows:

[0128] ,

[0129] In the formula, , , These are coefficients related to the mechanical parameters, geometric dimensions, and loads of the surrounding rock, which can be obtained at the equilibrium position. This is obtained by expanding the expression. In practical applications, the potential energy function is usually transformed into a function relating to dimensionless state variables. In the form of.

[0130] Step S1043: Taylor expansion and conversion to cusp mutation standard form

[0131] The total potential function At the current equilibrium position (i.e., the radius of the plastic zone) Performing a Taylor series expansion around the term ) and retaining up to the fourth-order terms, we obtain:

[0132] ,

[0133] In the formula, The radius of the plastic zone currently under consideration is a continuous variable. In order to be in Total potential energy at the point, For the total potential energy function pair The first derivative in The value at that location, For the total potential energy function pair The second derivative in The value at that location, For the total potential energy function pair The third derivative in The value at that location, For the total potential energy function pair The fourth derivative in The value at;

[0134] Since the first derivative of the system's potential energy is zero at the equilibrium position (i.e.) ), and ignore the constant term (Without affecting stability analysis), the above formula simplifies to:

[0135] ,

[0136] in, For state variables (increment of plastic zone radius), coefficients , , ;

[0137] By variable substitution (e.g., let...) The cubic terms can be eliminated, transforming the above expression into the standard form of a cusp mutation:

[0138] ,

[0139] in, These are all control variables, and their mapping relationship with the original mechanical parameters can be derived through substitution. Specifically:

[0140] ,

[0141] ,

[0142] In this embodiment, the calculation of these coefficients is based on the specific expression of the potential energy function in step S1042. For simplicity, the lengthy derivation process is omitted here; in actual programming, they can be directly calculated from the numerical values ​​of the mechanical parameters. .

[0143] Step S1044: Extract control variables and calculate instability characteristic values

[0144] According to the cusp catastrophe theory, the stability of a system is determined by the bifurcation set equation. The bifurcation set equation is:

[0145] ,

[0146] When discriminant At that time, the system is in a stable state;

[0147] When discriminant At this point, the system is in a critical state;

[0148] When discriminant At this time, the system is in an unstable state and may undergo sudden changes.

[0149] Therefore, the instability characteristic value of this key monitoring section is defined. for:

[0150] ,

[0151] in, The smaller the value (especially a negative value), the higher the risk of section instability. In practical engineering, it can be... It is considered a continuous indicator and used for subsequent comprehensive analysis.

[0152] Step S1045: Construct the spatial distribution matrix of unstable eigenvalues

[0153] To account for the spatial correlation among key cross-sections, the instability eigenvalues ​​and their relative distances between each cross-section need to be fused into a matrix. The specific construction process is as follows:

[0154] Section Numbering and Mileage Acquisition: Obtain the mileage coordinates of all key monitoring sections and number each section in ascending order of mileage. ,in This represents the total number of key monitoring sections;

[0155] Initialize a square matrix: Construct a OK Columns of squares ;

[0156] Fill diagonal elements: transform the square matrix The Middle Line number diagonal elements of the column Set as number Instability characteristic values ​​of each cross section ,Right now ;

[0157] Fill off-diagonal elements: for the first Line number List( The off-diagonal elements of ) are calculated using the following formula:

[0158] ,

[0159] in, and The first The instability characteristic value of the first cross section and the first Instability characteristic values ​​of each cross section, For the first The first cross section and the first The spatial distance between cross sections (i.e., the absolute value of the mileage difference, in meters). For a preset positive small quantity (e.g.) This is used to avoid the denominator being zero (when two cross-sections are located at the same mileage).

[0160] The physical meaning of this formula is that the larger the instability characteristic value of two cross sections and the closer they are, the stronger their interaction, and therefore the larger the matrix element value; conversely, the farther apart they are, the weaker their interaction, and the smaller the element value.

[0161] Row normalization: normalize the calculated square matrix Row normalization is performed by dividing each element of each row by the sum of all elements in that row, resulting in a row-normalized matrix. :

[0162] ,

[0163] In the formula, In the matrix obtained after row normalization, the position located at the th line, number Column elements, For the first Line number One off-diagonal element;

[0164] Row normalization makes the sum of the elements in each row equal to 1, eliminating the influence of the absolute magnitude of the instability characteristic values ​​of different cross sections and preserving the relative relationship between cross sections.

[0165] Column normalization: normalizing the matrix by row Column normalization is performed by dividing each element of each column by the sum of all elements in that column, resulting in the final standardized spatial distribution matrix of instability eigenvalues. :

[0166] ,

[0167] In the formula, In the final standardized matrix obtained after column normalization, the position located at the... line, number Column elements, The first normalized matrix line, number Column elements, For the row normalized matrix, the first... Liede One element;

[0168] Column normalization makes the sum of the elements in each column equal to 1, which further strengthens the symmetry and probabilistic interpretation of the matrix (it can be regarded as the spatial transfer probability matrix of the instability effects between cross sections).

[0169] In summary, the method in this embodiment integrates the instability eigenvalues ​​and the distances between cross sections through a ratio when constructing the spatial distribution matrix of instability eigenvalues. This approach considers both the inherent risk of a single cross section and the mutual influence between adjacent cross sections. This spatial correlation modeling aligns with the objective law that rock deformation in tunnel engineering exhibits continuity and transmissibility. Through dual processing of row and column normalization, the resulting standardized matrix eliminates the influence of dimensions and absolute values, making the cases of different tunnel sections and different numbers of cross sections comparable and providing standardized input data for subsequent principal component analysis.

[0170] Step S105: Perform principal component analysis on the spatial distribution matrix of the instability feature values, extract the first principal component score as the comprehensive instability index of each monitoring section, and generate a longitudinal instability risk curve based on the comprehensive instability index.

[0171] In this step, the spatial distribution matrix of the unstable eigenvalues ​​is used as the original data matrix, and the covariance matrix of the original data matrix is ​​calculated.

[0172] Solve for all eigenvalues ​​of the covariance matrix and arrange the eigenvalues ​​in descending order to obtain the eigenvalue sequence;

[0173] Calculate the variance contribution rate of each feature value in the feature value sequence, and determine whether the target variance contribution rate corresponding to the largest feature value is greater than a preset contribution rate threshold.

[0174] If the target variance contribution rate is greater than the preset contribution rate threshold, the feature vector corresponding to the largest feature value is directly extracted as the loading vector of the first principal component.

[0175] The inner product operation is performed between each row vector in the spatial distribution matrix of the instability eigenvalues ​​and the load vector of the first principal component to obtain the first principal component score of each key monitoring section, and the first principal component score is used as the comprehensive instability index of the key monitoring section.

[0176] Using the mileage coordinates of each key monitoring section as the abscissa and the comprehensive instability index of each section as the ordinate, a scatter plot is constructed. Then, cubic spline interpolation is performed on the scatter plot to obtain a continuous and smooth longitudinal instability risk curve.

[0177] In one specific embodiment, the standardized instability eigenvalue spatial distribution matrix obtained in step S104 is... The original data matrix for principal component analysis , It is The square array, in which, The total number of key monitoring sections, row index of the original data matrix Representing the A cross-section, column index Also represents the first A cross-section, element This indicates the first normalized value after double normalization. The first cross-section The standardized influence weight of each cross section.

[0178] For the original data matrix The mean of each column is subtracted from the mean of that column to obtain the centered matrix. Then calculate the covariance matrix. :

[0179] ,

[0180] In the formula, It is the transpose symbol;

[0181] covariance matrix It is Symmetric matrix, covariance matrix elements in Reflects the first Column and number The degree of linear correlation between the data (i.e. the cooperative change relationship between the instability characteristics of different sections).

[0182] For covariance matrix Perform eigenvalue decomposition to obtain all eigenvalues. and its corresponding eigenvectors Arrange the eigenvalues ​​in descending order, i.e. At the same time, the corresponding feature vectors are also arranged in the same order.

[0183] Calculate the variance contribution rate for each eigenvalue:

[0184] ,

[0185] In the formula, For the first Variance contribution rate of each eigenvalue For the first 1 eigenvalue, Let i be the i-th eigenvalue;

[0186] The variance contribution rate reflects the proportion of the total variance of the original data explained by the principal component corresponding to that eigenvalue. A preset contribution rate threshold is set. (In this embodiment, we take) (This can be adjusted according to project requirements). Determine the largest eigenvalue. Corresponding variance contribution rate Is it greater than :

[0187] like This indicates that the first principal component is sufficient to represent the main information of the original data, and there is no need to consider subsequent principal components.

[0188] The largest eigenvalue corresponding feature vector The load vector, as the first principal component, is the load vector. Each component in the table reflects the contribution weight of the original variables (i.e., the cross-sections corresponding to each column) to the first principal component;

[0189] For each key monitoring section (i.e., the first) of the original data matrix (line), its first principal component score From this row vector With load vector The inner product is obtained by:

[0190] ,

[0191] In the formula, For the standardized instability eigenvalue spatial distribution matrix, the first... line, number Column elements, Load vector The One component;

[0192] Specifically, if If the first principal component is insufficient to represent the main information of the original data, then multiple principal components need to be extracted and synthesized. The specific steps are as follows:

[0193] Determine the number of principal components: before calculation Cumulative variance contribution rate of each eigenvalue Take the satisfied The smallest value;

[0194] Extracting the load vector: Before extraction The eigenvectors corresponding to each eigenvalue ;

[0195] Calculate the weighted composite score: for each cross-section First calculate its first position. Scores on each principal component:

[0196] ,

[0197] Then, using the variance contribution rate of each principal component as a weight, a weighted composite score is calculated as the comprehensive instability index of the cross-section:

[0198] ,

[0199] The score That is, the first The comprehensive instability index of a section integrates the instability characteristics of the section itself and the spatial correlation information between it and other sections. It can more comprehensively reflect the instability risk level of the section in the whole tunnel section. The absolute value of the score does not have a direct physical unit, but its relative size can be used to rank the risks between sections: the higher the score (the larger the positive value), the higher the comprehensive instability risk of the section.

[0200] Mileage coordinates of each key monitoring section The horizontal axis represents the corresponding comprehensive instability index. Use the vertical axis to plot a scatter plot in a two-dimensional coordinate system.

[0201] Due to the limited number and potentially uneven distribution of monitoring sections, cubic spline interpolation is performed on the scatter plot to obtain a continuous risk distribution. Specifically, using... Construct cubic spline functions for nodes. So that in each interval superior It is a cubic polynomial that satisfies the condition that the function value, first derivative, and second derivative are continuous at the nodes. The smooth curve obtained after interpolation is the longitudinal instability risk curve.

[0202] The curve visually illustrates the trend of instability risk at various locations along the tunnel's longitudinal direction: the peak value corresponds to a high-risk section, the trough value corresponds to a low-risk section, and the slope of the curve reflects the point of sudden change in risk.

[0203] Step S106: Identify sections on the longitudinal instability risk curve whose slope exceeds a preset slope threshold as potential instability propagation sections, and determine the weak points of the support within the potential instability propagation sections based on the rate of change of the internal displacement ratio of the surrounding rock.

[0204] In this step, the first-order numerical derivative of the longitudinal instability risk curve is performed to obtain the instability risk slope curve;

[0205] Identify at least one continuous interval exceeding a preset slope threshold on the instability risk slope curve, and mark the corresponding segment on the longitudinal instability risk curve of the at least one continuous interval as a potential instability extension segment;

[0206] For each potential unstable expansion zone, the displacement time series data of the surrounding rock inside all monitoring sections within the potential unstable expansion zone are obtained, and the first depth measuring point closest to the excavation face and the second depth measuring point farthest from the excavation face are selected in each monitoring section.

[0207] For each monitoring section, the ratio of the displacement value of the second depth measuring point to the displacement value of the first depth measuring point is calculated at each monitoring time to obtain the depth displacement ratio time sequence of the monitoring section.

[0208] Perform a second-order difference operation on the depth displacement ratio time series to obtain the depth displacement ratio acceleration time series;

[0209] Identify the target time period in the depth displacement-to-acceleration time series where the value is continuously positive, and calculate the cumulative value of the depth displacement-to-acceleration within the target time period;

[0210] Monitoring sections whose cumulative values ​​exceed the preset cumulative threshold are identified as weak points in the support structure.

[0211] In one specific embodiment, let the longitudinal instability risk curve generated in step S105 be... ,in, For mileage coordinates, By performing a first-order numerical derivative, the instability risk slope curve is obtained. ;

[0212] In actual calculations, due to It is a smooth curve obtained through cubic spline interpolation, and its derivative can be directly obtained from the derivative expression of the spline function. The slope of the instability risk... It reflects the rate of change of the comprehensive instability index along the longitudinal direction of the tunnel: a positive slope indicates that the risk increases with the increase of mileage, a negative slope indicates that the risk decreases, and the larger the absolute value of the slope, the more drastic the risk change.

[0213] Set preset slope threshold (In this embodiment, we take) This means that the risk index increases by 0.05 per linear meter, which can be adjusted within the range of 0.02 to 0.1 depending on the actual project conditions. (This is reflected in the instability risk slope curve.) Search above for all that meet the criteria. For each continuous mileage interval, its longitudinal instability risk curve is plotted. The corresponding mileage segment is marked as a potentially unstable expansion segment.

[0214] These sections are characterized by a rapid increase in instability risk indicators, indicating that the stability of the surrounding rock is deteriorating rapidly and instability propagation may occur. Each potential instability propagation section may contain multiple monitoring sections or only one section, depending on the section length and the density of section layout.

[0215] For each potential instability propagation zone, time-series data of the internal displacement of the surrounding rock at all monitoring sections within that zone are acquired. For each monitoring section, two key depth measurement points are selected:

[0216] First depth measuring point: The depth measuring point closest to the excavation face (e.g., depth) This point mainly reflects the deformation of the shallow surrounding rock and is most directly affected by excavation disturbance.

[0217] Second depth measuring point: The depth measuring point farthest from the excavation face (e.g., depth) This point reflects the deformation of the deep surrounding rock, and its changes can reveal the expansion trend of the loosened zone;

[0218] The selection of these two depths should be based on the layout of the multi-point displacement gauges in the cross section. Usually, the two shallowest and deepest measuring points are selected to maximize the reflection of the gradient information of the deformation inside the surrounding rock.

[0219] For each monitoring section, at each monitoring time Calculate the displacement value of the second depth measuring point. Displacement value of the first depth measuring point The ratio of the depth displacement to the time series sequence is obtained. ;

[0220] Depth displacement ratio This reflects the coordination relationship between the deformation of deep and shallow surrounding rocks. When the surrounding rock is in a stable state, the deep displacement is usually smaller than the shallow displacement. Smaller and more gradual changes; as the loosened zone extends deeper, the increase in deep displacement accelerates. It shows an upward trend.

[0221] For depth displacement ratio time series Perform second-order difference operations to obtain the depth displacement-to-acceleration time series. The expression is:

[0222] ,

[0223] In the formula, In order to monitor the time The depth displacement ratio of acceleration, For the monitoring time interval, In order to monitor the time The depth displacement ratio, In order to monitor the time The depth displacement ratio, In order to monitor the time The depth displacement ratio;

[0224] Depth displacement to acceleration It characterizes the rate of increase or decrease of the deep displacement ratio change rate, that is, the "acceleration" of the displacement ratio change. When the displacement ratio increases, it indicates that the rate of increase is accelerating, meaning that the expansion of the loose zone is accelerating, which is a precursor to instability.

[0225] Depth displacement-to-acceleration time series In the middle, identify all that meet the requirements The continuous time intervals. For each consecutive positive time interval, calculate the cumulative value of the depth displacement to acceleration ratio within that time interval. :

[0226] ,

[0227] Set preset cumulative threshold (In this embodiment, we take) A point is considered a weak point when the cumulative acceleration value exceeds 0.1 (the specific value needs to be calibrated based on engineering experience). For each monitoring section, if its depth displacement is greater than the cumulative acceleration value over any positive time period... Greater than If so, then that section is identified as a weak point in the support.

[0228] Finally, output the section number, mileage location, and corresponding reinforcement suggestions for all weak points in the support (such as increasing the length of anchor bolts, grouting reinforcement, and temporary support) for on-site engineering personnel to refer to and make decisions.

[0229] In summary, the method of this application firstly, based on the proportional relationship between the internal displacement of the surrounding rock and the settlement of the arch, combines Kalman filtering for noise reduction, first-order differentiation, and positive time-period statistics to construct a loosening zone expansion index corrected by the surrounding rock integrity coefficient, thereby achieving a refined quantitative description of the loosening range of the surrounding rock and optimizing the section analysis sequence accordingly. Then, using the measured arch settlement as a constraint, the bisection method is used to iteratively solve the elastoplastic mechanical model within the physically defined stress release coefficient range, inverting the surrounding rock stress release coefficient and plastic zone radius of each section, transforming surface displacement into deep mechanical parameters with clear physical meaning, and automatically selecting key monitoring sections based on the stress release coefficient exceeding the threshold. On this basis, the cusp catastrophe theory is introduced, establishing the total potential energy function of the surrounding rock-support system based on the inversion parameters, which is then transformed into the standard form of cusp catastrophe through Taylor expansion and variable substitution, and calculated using the bifurcation set equation. The calculation results are used as instability characteristic values ​​to achieve continuous quantitative discrimination of the stability state of the surrounding rock system. Then, the instability characteristic values ​​of each section and the distance between sections are fused in the form of ratios. After double processing of row normalization and column normalization, a standardized spatial distribution matrix is ​​constructed so that the matrix elements simultaneously include the risk of the section itself and its spatial correlation. Then, principal component analysis is used to extract the first principal component with a variance contribution rate of more than 85% as the criterion. The inner product of the row vector of each section and the load vector is used to obtain the comprehensive instability index. A continuous and smooth longitudinal instability risk curve is generated by cubic spline interpolation to achieve effective dimensionality reduction and risk visualization of multidimensional information. Finally, the potential instability expansion section is automatically identified with the slope of the risk curve exceeding the threshold as the criterion. The depth displacement ratio and its second-order differential acceleration are calculated in the section. The weak points of the support are accurately located based on the continuous positive acceleration and the cumulative value exceeding the threshold.

[0230] This invention overcomes the technical bottlenecks of traditional methods, such as insufficient data utilization, difficulty in obtaining mechanical parameters, neglect of spatial correlation, difficulty in quantifying instability risk, and inaccurate weak point location. It significantly improves the automation, refinement, and intelligence of tunnel surrounding rock stability analysis, providing a scientific basis for dynamic design during tunnel construction, safety assessment during operation, and reinforcement decisions, and has important practical engineering value.

[0231] Please see Figure 2 The diagram shows a structural block diagram of a tunnel deformation analysis system based on multi-source data fusion according to this application.

[0232] like Figure 2 As shown, the tunnel deformation analysis system 200 includes an acquisition module 210, a first calculation module 220, an inversion module 230, a second calculation module 240, a generation module 250, and a determination module 260.

[0233] The acquisition module 210 is configured to acquire multi-source monitoring data from multiple monitoring sections within the target tunnel segment. This multi-source monitoring data includes time-series data on crown settlement, internal displacement of the surrounding rock, and stress on the support structure for each monitoring section. The first calculation module 220 is configured to calculate the expansion index of the loosened rock zone based on the ratio of internal displacement of the surrounding rock to crown settlement at each monitoring section, and to sort each monitoring section according to the expansion index to obtain a loosening trend sequence. The inversion module 230 is configured to sequentially select monitoring sections to be analyzed from the loosening trend sequence according to the expansion index of the loosened rock zone from largest to smallest for mechanical state inversion, obtaining the stress release coefficient and plastic zone radius of each monitoring section to be analyzed, until the stress release coefficient obtained by the inversion exceeds a preset stress threshold, at which point the selection stops. All selected monitoring sections to be analyzed at this time are determined as the set of key monitoring sections; the second calculation module 240 is configured to calculate the instability characteristic value of the corresponding key monitoring section based on the surrounding rock stress release coefficient and plastic zone radius of each key monitoring section using the catastrophe theory model, and construct the spatial distribution matrix of instability characteristic value according to the instability characteristic value and spatial distance of each key monitoring section; the generation module 250 is configured to perform principal component analysis on the spatial distribution matrix of instability characteristic value, extract the first principal component score as the comprehensive instability index of each monitoring section, and generate a longitudinal instability risk curve based on the comprehensive instability index; the determination module 260 is configured to identify the section with a slope exceeding a preset slope threshold on the longitudinal instability risk curve as a potential instability expansion section, and determine the support weak point within the potential instability expansion section according to the rate of change of the internal displacement ratio of the surrounding rock.

[0234] It should be understood that Figure 2 The modules and references described in the document Figure 1 The steps described in the text correspond to those in the method described above. Therefore, the operations, features, and corresponding technical effects described above also apply to the method described in the text. Figure 2 The various modules in the document will not be described in detail here.

[0235] In other embodiments, the present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein when the program instructions are executed by a processor, the processor performs the tunnel deformation analysis method based on multi-source data fusion in any of the above method embodiments.

[0236] In one embodiment, the computer-readable storage medium of the present invention stores computer-executable instructions, which are configured as follows:

[0237] Acquire multi-source monitoring data from multiple monitoring sections within the target tunnel section. The multi-source monitoring data includes the arch settlement time series data, the displacement time series data inside the surrounding rock, and the stress time series data of the support structure for each monitoring section.

[0238] The expansion index of the loosened zone of the surrounding rock is calculated based on the proportional relationship between the internal displacement of the surrounding rock and the settlement of the arch at the monitoring section. The monitoring sections are then sorted according to the expansion index of the loosened zone of the surrounding rock to obtain the loosening trend sequence of the sections.

[0239] In the loosening trend sequence of the cross section, the cross sections to be analyzed and monitored are selected in order of the expansion index of the loosening zone of the surrounding rock from large to small to perform mechanical state inversion, so as to obtain the stress release coefficient and plastic zone radius of each cross section to be analyzed and monitored. The selection stops when the stress release coefficient of the surrounding rock obtained by the inversion exceeds the preset stress threshold. All the cross sections to be analyzed and monitored at this time are determined as the set of key monitoring sections.

[0240] Based on the stress release coefficient and plastic zone radius of each key monitoring section, the catastrophe theory model is used to calculate the instability characteristic value of the corresponding key monitoring section, and the spatial distribution matrix of the instability characteristic value is constructed according to the instability characteristic value and spatial distance of each key monitoring section.

[0241] Principal component analysis is performed on the spatial distribution matrix of the instability eigenvalues, and the score of the first principal component is extracted as the comprehensive instability index for each monitoring section. Based on the comprehensive instability index, a longitudinal instability risk curve is generated.

[0242] On the longitudinal instability risk curve, segments with slopes exceeding a preset slope threshold are identified as potential instability propagation segments, and weak points in the support are determined within these potential instability propagation segments based on the rate of change of the internal displacement ratio of the surrounding rock.

[0243] Computer-readable storage media may include a stored program area and a stored data area, wherein the stored program area may store an operating system and an application program required for at least one function; the stored data area may store data created based on the use of the tunnel deformation analysis system based on multi-source data fusion, etc. Furthermore, the computer-readable storage medium may include high-speed random access memory, and may also include memory, such as at least one disk storage device, flash memory device, or other non-volatile solid-state storage device. In some embodiments, the computer-readable storage medium may optionally include memory remotely disposed relative to a processor, which can be connected to the tunnel deformation analysis system based on multi-source data fusion via a network. Examples of such networks include, but are not limited to, the Internet, corporate intranets, local area networks, mobile communication networks, and combinations thereof.

[0244] Figure 3 This is a schematic diagram of the structure of the electronic device provided in the embodiment of the present invention, such as... Figure 3As shown, the device includes a processor 310 and a memory 320. The electronic device may also include an input device 330 and an output device 340. The processor 310, memory 320, input device 330, and output device 340 can be connected via a bus or other means. Figure 3 Taking a bus connection as an example, the memory 320 is the computer-readable storage medium described above. The processor 310 executes various server functions and data processing by running non-volatile software programs, instructions, and modules stored in the memory 320, thereby implementing the tunnel deformation analysis method based on multi-source data fusion as described in the above embodiment. The input device 330 can receive input digital or character information and generate key signal inputs related to user settings and function control of the tunnel deformation analysis system based on multi-source data fusion. The output device 340 may include a display screen or other display device.

[0245] The aforementioned electronic device can execute the method provided in the embodiments of the present invention, and has the corresponding functional modules and beneficial effects for executing the method. Technical details not described in detail in this embodiment can be found in the method provided in the embodiments of the present invention.

[0246] In one implementation, the above-described electronic device is applied to a tunnel deformation analysis system based on multi-source data fusion, and is used as a client. It includes: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, which, when executed by the at least one processor, enable the at least one processor to:

[0247] Acquire multi-source monitoring data from multiple monitoring sections within the target tunnel section. The multi-source monitoring data includes the arch settlement time series data, the displacement time series data inside the surrounding rock, and the stress time series data of the support structure for each monitoring section.

[0248] The expansion index of the loosened zone of the surrounding rock is calculated based on the proportional relationship between the internal displacement of the surrounding rock and the settlement of the arch at the monitoring section. The monitoring sections are then sorted according to the expansion index of the loosened zone of the surrounding rock to obtain the loosening trend sequence of the sections.

[0249] In the loosening trend sequence of the cross section, the cross sections to be analyzed and monitored are selected in order of the expansion index of the loosening zone of the surrounding rock from large to small to perform mechanical state inversion, so as to obtain the stress release coefficient and plastic zone radius of each cross section to be analyzed and monitored. The selection stops when the stress release coefficient of the surrounding rock obtained by the inversion exceeds the preset stress threshold. All the cross sections to be analyzed and monitored at this time are determined as the set of key monitoring sections.

[0250] Based on the stress release coefficient and plastic zone radius of each key monitoring section, the catastrophe theory model is used to calculate the instability characteristic value of the corresponding key monitoring section, and the spatial distribution matrix of the instability characteristic value is constructed according to the instability characteristic value and spatial distance of each key monitoring section.

[0251] Principal component analysis is performed on the spatial distribution matrix of the instability eigenvalues, and the score of the first principal component is extracted as the comprehensive instability index for each monitoring section. Based on the comprehensive instability index, a longitudinal instability risk curve is generated.

[0252] On the longitudinal instability risk curve, segments with slopes exceeding a preset slope threshold are identified as potential instability propagation segments, and weak points in the support are determined within these potential instability propagation segments based on the rate of change of the internal displacement ratio of the surrounding rock.

[0253] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., including several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods of various embodiments or some parts of embodiments.

[0254] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A tunnel deformation analysis method based on multi-source data fusion, characterized in that, include: Acquire multi-source monitoring data from multiple monitoring sections within the target tunnel section. The multi-source monitoring data includes the arch settlement time series data, the displacement time series data inside the surrounding rock, and the stress time series data of the support structure for each monitoring section. The expansion index of the loosened zone of the surrounding rock is calculated based on the proportional relationship between the internal displacement of the surrounding rock and the settlement of the arch at the monitoring section. The monitoring sections are then sorted according to the expansion index of the loosened zone of the surrounding rock to obtain the loosening trend sequence of the sections. In the loosening trend sequence of the cross section, the cross sections to be analyzed and monitored are selected in order of the expansion index of the loosening zone of the surrounding rock from large to small to perform mechanical state inversion, so as to obtain the stress release coefficient and plastic zone radius of each cross section to be analyzed and monitored. The selection stops when the stress release coefficient of the surrounding rock obtained by the inversion exceeds the preset stress threshold. All the cross sections to be analyzed and monitored at this time are determined as the set of key monitoring sections. Based on the stress release coefficient and plastic zone radius of each key monitoring section, the catastrophe theory model is used to calculate the instability characteristic value of the corresponding key monitoring section, and the spatial distribution matrix of the instability characteristic value is constructed according to the instability characteristic value and spatial distance of each key monitoring section. Principal component analysis is performed on the spatial distribution matrix of the instability eigenvalues, and the score of the first principal component is extracted as the comprehensive instability index for each monitoring section. Based on the comprehensive instability index, a longitudinal instability risk curve is generated. On the longitudinal instability risk curve, segments with slopes exceeding a preset slope threshold are identified as potential instability propagation segments, and weak points in the support are determined within these potential instability propagation segments based on the rate of change of the internal displacement ratio of the surrounding rock.

2. The tunnel deformation analysis method based on multi-source data fusion according to claim 1, characterized in that, The expansion index of the loosened zone of the surrounding rock is calculated based on the proportional relationship between the internal displacement of the surrounding rock and the settlement of the arch based on the monitoring section. The monitoring sections are then sorted according to the expansion index of the loosened zone to obtain a sequence of loosening trends for each section, including: For each monitoring section, at least three displacement measuring points inside the surrounding rock at different depths are selected. The displacement values ​​of each measuring point at each depth are extracted at each monitoring time. The ratio of the displacement value of each measuring point at each depth at each monitoring time to the crown settlement value at the corresponding time is calculated to obtain the original sequence of displacement ratios of each measuring point at each depth. Kalman filtering was applied to the original sequence of displacement ratios at each depth measurement point to remove measurement noise and random fluctuations, resulting in a smoothed sequence of displacement ratios at each depth measurement point. The first-order numerical derivative of the smoothed displacement ratio sequence at each depth measuring point is performed to obtain the displacement ratio change rate sequence at each depth measuring point. The sign of the displacement ratio change rate sequence is determined, and the time period with positive change rate is marked. The longest duration during which the displacement ratio change rate of each depth measuring point is continuously positive is statistically analyzed, and the depth measuring points whose longest duration exceeds a first preset duration threshold are identified as potential influence measuring points of the loosening ring. From all the potential impact measurement points of the loosening zone, the continuous depth intervals where the average displacement ratio change rate exceeds the preset change rate threshold are selected, and the continuous depth intervals are defined as the loosening zone expansion zone. The thickness of the loosened zone expansion area is measured, and the thickness is divided by the equivalent radius of the tunnel diameter of the corresponding monitoring section to obtain the initial loosened zone expansion index; The surrounding rock integrity coefficient of the monitoring section is obtained, and the initial loosened zone expansion index is corrected according to the surrounding rock integrity coefficient to obtain the final surrounding rock loosened zone expansion index. All monitoring sections are arranged in descending order of the surrounding rock loosening zone expansion index to generate a section loosening trend sequence.

3. The tunnel deformation analysis method based on multi-source data fusion according to claim 1, characterized in that, In the loosening trend sequence of the cross-section, the cross-sections to be analyzed and monitored are selected in descending order of the expansion index of the loosened zone of the surrounding rock for mechanical state inversion, so as to obtain the stress release coefficient and plastic zone radius of each cross-section to be analyzed and monitored. The selection stops when the stress release coefficient of the surrounding rock obtained by the inversion exceeds the preset stress threshold. All the cross-sections to be analyzed and monitored at this time are determined as the set of key monitoring cross-sections, including: Establish the current elastoplastic mechanical model of the monitoring section to be analyzed; Obtain the current crown settlement value at the latest moment in the crown settlement time series data of the current monitoring section to be analyzed, and use the current crown settlement value as the measured constraint value of the radial displacement of the tunnel wall; Within the preset stress relief coefficient search range, the bisection method is used for iterative solution. In each iteration, the current stress relief coefficient is substituted into the current elastoplastic mechanical model to calculate the corresponding theoretical cavity wall displacement and theoretical plastic zone radius. The relative error between the theoretical tunnel wall displacement and the measured tunnel wall displacement is calculated. If the relative error is greater than a preset accuracy threshold, the range of the stress relief coefficient is adjusted according to the sign of the error, and the iteration continues. When the relative error between the theoretical tunnel wall displacement and the measured tunnel wall displacement is less than or equal to the preset accuracy threshold, the iteration stops, and the current stress release coefficient and the corresponding current theoretical plastic zone radius are used as the inversion result of the section to be analyzed and monitored. The next monitoring section to be analyzed is inverted until the stress release coefficient of the surrounding rock obtained by the inversion of a certain monitoring section exceeds the preset stress threshold. Then the selection and inversion are stopped, and all the monitoring sections to be analyzed that have been inverted are identified as the set of key monitoring sections.

4. The tunnel deformation analysis method based on multi-source data fusion according to claim 1, characterized in that, Based on the surrounding rock stress release coefficient and plastic zone radius of each key monitoring section, the instability characteristic values ​​of the corresponding key monitoring sections are calculated using a catastrophe theory model, including: Obtain the tunnel radius of key monitoring sections, the support resistance provided by the support structure, and the stress release coefficient and plastic zone radius of the surrounding rock obtained through inversion; Based on the theory of elastoplastic mechanics, a total potential energy function of the surrounding rock-support system at a key monitoring section is established. The total potential energy function includes the strain energy stored in the elastic zone of the surrounding rock, the plasticity consumed in the plastic zone, and the potential energy provided by the support structure. Expand the total potential energy function using a Taylor series near the current equilibrium position, retaining up to the fourth-order terms, to obtain an approximate polynomial expression for the potential energy function. By performing variable substitution and linear transformation on the approximate polynomial expression, the standard form of the cusp catastrophe model is obtained: Where p and q are control variables, and x is a state variable; Extract the specific values ​​of the control variables from the standard form, and based on the cusp catastrophe theory, transform the bifurcation set equation... The calculation results are defined as the instability characteristic values ​​of the key monitoring sections.

5. The tunnel deformation analysis method based on multi-source data fusion according to claim 1, characterized in that, The construction of the spatial distribution matrix of instability characteristic values ​​based on the instability characteristic values ​​and spatial distances of each key monitoring section includes: Obtain the mileage coordinates of all key monitoring sections, and number each key monitoring section in ascending order of mileage to obtain the section number sequence 1, 2, ..., n, where n is the total number of key monitoring sections. Construct an n x n square matrix M, and set the diagonal element of the i-th row and i-th column of the square matrix M as the instability characteristic value λi of the i-th key monitoring section; The expression for calculating the value of the non-diagonal element in the i-th row and j-th column of a square matrix M is: ， Where λi is the instability characteristic value of the i-th cross section, λj is the instability characteristic value of the j-th cross section, Lij is the spatial distance between the i-th key monitoring cross section and the j-th key monitoring cross section, and δ is a preset positive small quantity. The calculated square matrix M is subjected to row normalization, that is, the elements of each row are divided by the sum of the elements of the current row to obtain the row normalized unstable eigenvalue spatial distribution matrix. The column normalization process is performed on the row-normalized unstable eigenvalue spatial distribution matrix, that is, the element of each column is divided by the sum of the elements of the current column to obtain the final standardized unstable eigenvalue spatial distribution matrix.

6. The tunnel deformation analysis method based on multi-source data fusion according to claim 1, characterized in that, The step of performing principal component analysis on the spatial distribution matrix of the instability eigenvalues, extracting the score of the first principal component as the comprehensive instability index for each monitoring section, and generating a longitudinal instability risk curve based on the comprehensive instability index includes: Using the spatial distribution matrix of the unstable eigenvalues ​​as the original data matrix, calculate the covariance matrix of the original data matrix; Solve for all eigenvalues ​​of the covariance matrix and arrange the eigenvalues ​​in descending order to obtain the eigenvalue sequence; Calculate the variance contribution rate of each feature value in the feature value sequence, and determine whether the target variance contribution rate corresponding to the largest feature value is greater than a preset contribution rate threshold. If the target variance contribution rate is greater than the preset contribution rate threshold, the feature vector corresponding to the largest feature value is directly extracted as the loading vector of the first principal component. The inner product operation is performed between each row vector in the spatial distribution matrix of the instability eigenvalues ​​and the load vector of the first principal component to obtain the first principal component score of each key monitoring section, and the first principal component score is used as the comprehensive instability index of the key monitoring section. Using the mileage coordinates of each key monitoring section as the abscissa and the comprehensive instability index of each section as the ordinate, a scatter plot is constructed. Then, cubic spline interpolation is performed on the scatter plot to obtain a continuous and smooth longitudinal instability risk curve.

7. The tunnel deformation analysis method based on multi-source data fusion according to claim 1, characterized in that, The step of identifying sections on the longitudinal instability risk curve whose slope exceeds a preset slope threshold as potential instability propagation sections, and determining weak points in the support within these potential instability propagation sections based on the rate of change of the displacement ratio within the surrounding rock, includes: The first-order numerical derivative of the longitudinal instability risk curve is used to obtain the instability risk slope curve. Identify at least one continuous interval exceeding a preset slope threshold on the instability risk slope curve, and mark the corresponding segment on the longitudinal instability risk curve of the at least one continuous interval as a potential instability extension segment; For each potential unstable expansion zone, the displacement time series data of the surrounding rock inside all monitoring sections within the potential unstable expansion zone are obtained, and the first depth measuring point closest to the excavation face and the second depth measuring point farthest from the excavation face are selected in each monitoring section. For each monitoring section, the ratio of the displacement value of the second depth measuring point to the displacement value of the first depth measuring point is calculated at each monitoring time to obtain the depth displacement ratio time sequence of the monitoring section. Perform a second-order difference operation on the depth displacement ratio time series to obtain the depth displacement ratio acceleration time series; Identify the target time period in the depth displacement-to-acceleration time series where the value is continuously positive, and calculate the cumulative value of the depth displacement-to-acceleration within the target time period; Monitoring sections whose cumulative values ​​exceed the preset cumulative threshold are identified as weak points in the support structure.

8. A tunnel deformation analysis system based on multi-source data fusion, characterized in that, include: The acquisition module is configured to acquire multi-source monitoring data from multiple monitoring sections within the target tunnel section. The multi-source monitoring data includes the arch settlement time series data, the displacement time series data inside the surrounding rock, and the stress time series data of the support structure for each monitoring section. The first calculation module is configured to calculate the expansion index of the loosened zone of the surrounding rock based on the ratio of the internal displacement of the surrounding rock to the settlement of the arch at the monitoring section, and to sort each monitoring section according to the expansion index of the loosened zone of the surrounding rock to obtain the loosening trend sequence of the section. The inversion module is configured to select the monitoring sections to be analyzed in the sequence of loosening trend of the cross section in order of the expansion index of the loosening zone of the surrounding rock from large to small, and perform mechanical state inversion to obtain the stress release coefficient and plastic zone radius of each monitoring section to be analyzed. The selection stops when the stress release coefficient of the surrounding rock obtained by the inversion exceeds the preset stress threshold, and all the monitoring sections to be analyzed at this time are determined as the set of key monitoring sections. The second calculation module is configured to calculate the instability characteristic values ​​of the corresponding key monitoring sections based on the surrounding rock stress release coefficient and plastic zone radius of each key monitoring section using the catastrophe theory model, and construct a spatial distribution matrix of instability characteristic values ​​based on the instability characteristic values ​​and spatial distance of each key monitoring section. The generation module is configured to perform principal component analysis on the spatial distribution matrix of the instability feature values, extract the first principal component score as the comprehensive instability index of each monitoring section, and generate a longitudinal instability risk curve based on the comprehensive instability index. The determination module is configured to identify sections on the longitudinal instability risk curve whose slope exceeds a preset slope threshold as potential instability propagation sections, and to determine the weak points of the support within the potential instability propagation sections based on the rate of change of the internal displacement ratio of the surrounding rock.

9. An electronic device, characterized in that, include: At least one processor, and a memory communicatively connected to the at least one processor, wherein the memory stores instructions executable by the at least one processor to enable the at least one processor to perform the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the program is executed by the processor, it implements the method according to any one of claims 1 to 7.

Images