Mountain tunnel surrounding rock deformation monitoring method and system based on distributed optical fiber sensing

Abstract

The present application relates to the technical field of deformation monitoring, and more particularly to a mountain tunnel surrounding rock deformation monitoring method and system based on distributed optical fiber sensing. The method comprises the following steps: collecting Brillouin frequency shift original data and center wavelength original data, simultaneously acquiring tunnel construction progress data and surrounding rock geological survey data and performing format standardization and time sequence alignment processing to generate preprocessed optical fiber sensing basic data, construction time sequence data and geological attribute data; extracting surrounding rock deformation characteristic parameters based on the optical fiber sensing basic data and constructing a surrounding rock deformation multidimensional feature vector; performing deformation mode density clustering and surrounding rock deformation dynamic prediction based on the surrounding rock deformation multidimensional feature vector to generate a surrounding rock deformation dynamic prediction value; comparing and iteratively correcting the surrounding rock deformation dynamic prediction value and performing regional aggregation to generate a mountain tunnel surrounding rock deformation spatiotemporal distribution visualization result. The present application can realize full spatiotemporal accurate monitoring and early warning of tunnel surrounding rock deformation.

Description

Technical Field

[0001] This invention relates to the field of deformation monitoring technology, and in particular to a method and system for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing. Background Technology

[0002] In recent years, as mountain tunnels have been advanced to longer distances, greater depths, and more complex geological areas (such as karst development zones and fractured rock strata), the real-time performance, full-range coverage, and accuracy of early warning monitoring of surrounding rock deformation have become core requirements for ensuring tunnel construction safety and long-term operational stability. Traditional surrounding rock deformation monitoring relies on discrete point monitoring equipment such as total stations and multi-point displacement gauges, which suffer from problems such as numerous monitoring blind spots, poor data continuity, and significant interference from construction. It is difficult to capture the overall deformation trend of the surrounding rock along the tunnel axis, let alone predict complex deformation risks such as local abrupt changes and gradual accumulation, which can easily lead to safety accidents such as collapses and water inrushes. However, existing methods lack sufficient ability to fuse multi-source fiber optic sensor data and adapt to complex deformation scenarios. On the one hand, due to the limitations of single sensing technologies (such as low point resolution of BOTDA and limited coverage of FBG), it is easy to miss deformation at key points of the surrounding rock or to inaccurately characterize the deformation distribution along the process. On the other hand, deformation analysis is not fully combined with the tunnel construction progress and the geological properties of the surrounding rock. It relies solely on isolated judgment based on fiber optic sensor data, which is easily affected by factors such as sudden changes in geological conditions and construction disturbances, resulting in deviations in deformation calculation or delays in risk warning. Summary of the Invention

[0003] Therefore, it is necessary for the present invention to provide a method and system for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing, so as to solve at least one of the above-mentioned technical problems.

[0004] To achieve the above objectives, a method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing includes the following steps:

[0005] Step S1: Collect raw Brillouin frequency shift data along the sensing optical cable and raw center wavelength data of each sensing node in the monitoring area of ​​the mountain tunnel surrounding rock. At the same time, acquire tunnel construction progress data and surrounding rock geological survey data. Perform outlier removal and noise filtering on the raw Brillouin frequency shift data and raw center wavelength data. Perform format standardization and time sequence alignment on the tunnel construction progress data and surrounding rock geological survey data to generate preprocessed fiber optic sensing basic data, construction time sequence data and geological attribute data.

[0006] Step S2: Extract the surrounding rock deformation characteristic parameters based on the fiber optic sensing basic data, including the axial strain value of each monitoring point of the sensing cable and the strain distribution parameters along the surrounding rock; analyze the center wavelength offset of each sensing node and separate the pure strain component to generate the point strain parameters of the key points of the surrounding rock; at the same time, construct a multi-dimensional feature vector of surrounding rock deformation based on the surrounding rock deformation characteristic parameters combined with construction time sequence data and geological attribute data.

[0007] Step S3: Based on the multidimensional feature vector of surrounding rock deformation, perform deformation pattern density clustering on the surrounding rock of different monitoring sections of the tunnel to classify the deformation pattern categories corresponding to uniform deformation, local abrupt deformation, and gradual accumulation. For different deformation pattern categories, introduce the process progress weight in the construction time sequence data and the lithological influence coefficient in the geological attribute data to dynamically predict the surrounding rock deformation and generate the dynamic prediction value of surrounding rock deformation for each monitoring section.

[0008] Step S4: Compare the dynamic prediction values ​​of surrounding rock deformation with the monitoring data from the total station and multi-point displacement gauges deployed at the tunnel site, calculate the absolute error, relative error, and root mean square error between the predicted and measured values, and form a surrounding rock deformation prediction error matrix; perform iterative correction based on the surrounding rock deformation prediction error matrix, and aggregate regionally according to the tunnel mileage segment and deformation mode category to generate a visualization result of the spatiotemporal distribution of surrounding rock deformation in the mountain tunnel.

[0009] This invention also provides a mountain tunnel surrounding rock deformation monitoring system based on distributed optical fiber sensing, used to execute the mountain tunnel surrounding rock deformation monitoring method based on distributed optical fiber sensing as described above. The mountain tunnel surrounding rock deformation monitoring system based on distributed optical fiber sensing includes:

[0010] The surrounding rock area data acquisition module is used to collect raw Brillouin frequency shift data along the sensing optical cable and raw center wavelength data of each sensing node in the monitoring area of ​​the surrounding rock of the mountain tunnel. At the same time, it acquires tunnel construction progress data and surrounding rock geological survey data. The module performs outlier removal and noise filtering on the raw Brillouin frequency shift data and raw center wavelength data, and performs format standardization and time sequence alignment on the tunnel construction progress data and surrounding rock geological survey data to generate preprocessed fiber optic sensing basic data, construction time sequence data and geological attribute data.

[0011] The surrounding rock deformation feature analysis module is used to extract surrounding rock deformation feature parameters based on fiber optic sensing data, including the axial strain values ​​of each monitoring point of the sensing cable and the strain distribution parameters along the surrounding rock; it analyzes the center wavelength offset of each sensing node and separates the pure strain component to generate the point strain parameters of key points of the surrounding rock; at the same time, it constructs a multi-dimensional feature vector of surrounding rock deformation based on the surrounding rock deformation feature parameters combined with construction time series data and geological attribute data.

[0012] The surrounding rock deformation prediction module is used to perform deformation pattern density clustering on the surrounding rock of different monitoring sections of the tunnel based on the multi-dimensional feature vector of surrounding rock deformation, so as to classify the deformation pattern categories corresponding to uniform deformation, local abrupt deformation, and gradual accumulation. For different deformation pattern categories, the module introduces the process progress weight in the construction time sequence data and the lithological influence coefficient in the geological attribute data to dynamically predict the surrounding rock deformation, and generate the dynamic prediction value of surrounding rock deformation for each monitoring section.

[0013] The surrounding rock deformation visualization module is used to compare the dynamic predicted values ​​of surrounding rock deformation with the monitoring data of total station and multi-point displacement gauges deployed on-site in the tunnel, calculate the absolute error, relative error and root mean square error between the predicted and measured values, and form a surrounding rock deformation prediction error matrix. Based on the surrounding rock deformation prediction error matrix, iterative correction is performed, and regional aggregation is performed according to the tunnel mileage segment and deformation mode category, thereby generating a visualization result of the spatiotemporal distribution of surrounding rock deformation in the mountain tunnel.

[0014] The beneficial effects of this invention are:

[0015] 1. The method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing proposed in this invention has the following advantages compared with the prior art: By synchronously acquiring BOTDA along the Brillouin frequency shift data and FBG node center wavelength data, and supplementing it with construction progress and geological survey data, it achieves multi-dimensional data coverage of "large area + key points + external influences"; by removing outliers and filtering noise from optical fiber sensing data, it eliminates equipment errors and environmental interference (such as false signals caused by temperature fluctuations); by standardizing the format and aligning the time sequence of construction and geological data, it ensures the consistency of time and space dimensions of data from different sources. This multi-source data integration and preprocessing mechanism not only makes up for the coverage blind spots of single sensing technology, but also ensures data accuracy and compatibility, avoiding deformation analysis deviations caused by incomplete or distorted data from the source, and laying a solid foundation for subsequent deformation feature extraction. Secondly, axial strain is calculated using BOTDA and converted into strain distribution parameters along the tunnel, accurately presenting the continuous deformation trend of the surrounding rock throughout the tunnel (such as the overall settlement pattern of a certain mileage section). Pure strain components are separated by combining FBG with a cross-sensitive compensation model to capture minute deformations (such as millimeter-level displacements) at key points such as the tunnel arch and sidewalls. The two complement each other to resolve the contradiction between "large-scale coverage" and "precise key points". At the same time, deformation feature parameters are integrated with construction time sequence and geological attribute data to construct a multi-dimensional feature vector. This allows deformation features to not only include physical deformation data but also incorporate external influencing factors. This multi-dimensional and complementary feature extraction design completely solves the limitations of single sensing technology, ensuring that deformation characterization is both comprehensive and accurate, and deeply bound to actual construction and geological scenarios, providing rich feature support for subsequent deformation pattern recognition. Then, by using density clustering based on multidimensional feature vectors to classify deformation patterns (uniform deformation, local mutation, and gradual accumulation), deformation analysis is upgraded from "general judgment" to "classification and policy implementation". For different patterns, process progress weights and lithological influence coefficients are introduced to dynamically predict deformation, ensuring that the prediction results fit the actual scenario. This "classification and identification + scenario-based prediction" mechanism not only avoids the single-minded understanding of deformation patterns in existing methods, but also incorporates key influencing factors such as construction disturbance and lithological differences into the prediction model, which greatly improves the accuracy and timeliness of deformation prediction and effectively solves the problem of delayed risk warning caused by isolated analysis.Finally, by comparing the predicted deformation values ​​with the measured data from total stations and multi-point displacement gauges, the absolute error, relative error, and root mean square error are calculated to form an error matrix. The prediction model is then iteratively corrected and optimized (e.g., adjusting the lithological influence coefficient for cases where the error exceeds the standard in a certain lithological section) to ensure that the predicted results are highly consistent with the actual deformation. At the same time, regional aggregation is performed according to mileage sections and deformation patterns to generate a spatial and temporal distribution visualization result, allowing construction personnel to intuitively grasp the deformation risk of the entire tunnel. This not only solves the problem of the lack of credibility of the prediction results of existing methods, but also transforms abstract data into intuitive information that can directly guide construction, completely releasing the engineering value of surrounding rock monitoring data and effectively reducing the construction safety risks caused by inaccurate deformation prediction.

[0016] 2. The mountain tunnel surrounding rock deformation monitoring system based on distributed optical fiber sensing proposed in this invention consists of a surrounding rock area data acquisition module, a surrounding rock deformation feature analysis module, a surrounding rock deformation prediction module, and a surrounding rock deformation visualization module. It can realize any mountain tunnel surrounding rock deformation monitoring method based on distributed optical fiber sensing described in this invention. The system uses the combined operations of computer programs running on each module to achieve the mountain tunnel surrounding rock deformation monitoring method based on distributed optical fiber sensing. The internal structure of the system cooperates with each other, which can greatly reduce repetitive work and manpower input, and can quickly and effectively provide a more accurate and efficient mountain tunnel surrounding rock deformation monitoring process based on distributed optical fiber sensing, thereby simplifying the operation process of the mountain tunnel surrounding rock deformation monitoring system based on distributed optical fiber sensing. Attached Figure Description

[0017] Other features, objects, and advantages of the invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings:

[0018] Figure 1 This is a schematic diagram of the steps of the method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing according to the present invention.

[0019] Figure 2 for Figure 1 A detailed flowchart of step S1;

[0020] Figure 3 This is a cross-sectional schematic diagram of the strain distribution cloud map along the tunnel of the present invention;

[0021] Figure 4 This is a schematic diagram of the module of the mountain tunnel surrounding rock deformation monitoring system based on distributed optical fiber sensing of the present invention. Detailed Implementation

[0022] The technical system of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0023] To achieve the above objectives, please refer to Figures 1 to 4 This invention provides a method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing. For embodiments of this invention, please refer to... Figure 1 The diagram shown is a flowchart illustrating the steps of the mountain tunnel surrounding rock deformation monitoring method based on distributed optical fiber sensing according to the present invention. In this example, the mountain tunnel surrounding rock deformation monitoring method based on distributed optical fiber sensing includes the following steps:

[0024] Step S1: Collect raw Brillouin frequency shift data along the sensing optical cable and raw center wavelength data of each sensing node in the monitoring area of ​​the mountain tunnel surrounding rock. At the same time, acquire tunnel construction progress data and surrounding rock geological survey data. Perform outlier removal and noise filtering on the raw Brillouin frequency shift data and raw center wavelength data. Perform format standardization and time sequence alignment on the tunnel construction progress data and surrounding rock geological survey data to generate preprocessed fiber optic sensing basic data, construction time sequence data and geological attribute data.

[0025] In this embodiment of the invention, the BOTDA demodulation system is used to collect raw Brillouin frequency shift data along the optical cable once per hour in the monitoring area of ​​the mountain tunnel from K1+000 to K1+500, for a total of 720 sets of data, with a frequency shift range of 10.500-10.550GHz; the FBG demodulation system is used to simultaneously collect raw center wavelength data of 200 sensor nodes, with a wavelength range of 1540-1560nm. Construction progress data was obtained from tunnel construction management records, including the completion time of portal excavation (15 days after commencement), the start time of initial support (30 days), and the face advance rate (1.5 m / day for the K1+000-K1+100 section and 1.2 m / day for the K1+100-K1+200 section). Geological survey data of the surrounding rock was obtained through ground-penetrating radar detection and borehole sampling, including the surrounding rock integrity index (0.75 for the K1+000-K1+100 section and 0.6 for the K1+100-K1+200 section) and the uniaxial compressive strength of the rock mass (30 MPa for sandstone and 15 MPa for silty clay). An outlier removal tool was used, defining Brillouin frequency shift deviation exceeding 0.02 GHz and center wavelength deviation exceeding 0.5 nm as outliers, removing 12 sets of frequency shift anomalies and 8 nodal wavelength anomalies. Wavelet denoising was used to filter noise from the valid data, decomposing it into 5 layers and retaining the first 3 layers of low-frequency coefficients for reconstruction. For construction progress data, a time-series standardization tool is used to convert time nodes into cumulative durations relative to the start of construction, and the progress rate is generated as a continuous curve based on mileage station number; for geological exploration data, a dimensionless conversion tool is used, dividing the integrity index by 1 and the compressive strength by 50MPa to generate pre-processed fiber optic sensing basic data (including 1001 BOTDA monitoring points and 192 FBG effective node data), construction time-series data (including 50 time node parameters), and geological attribute data (including 50 cross-sectional geological parameters).

[0026] Step S2: Extract the surrounding rock deformation characteristic parameters based on the fiber optic sensing basic data, including the axial strain value of each monitoring point of the sensing cable and the strain distribution parameters along the surrounding rock; analyze the center wavelength offset of each sensing node and separate the pure strain component to generate the point strain parameters of the key points of the surrounding rock; at the same time, construct a multi-dimensional feature vector of surrounding rock deformation based on the surrounding rock deformation characteristic parameters combined with construction time sequence data and geological attribute data.

[0027] In this embodiment of the invention, based on fiber optic sensing data, the BOTDA strain calculation tool is used, combined with the system-calibrated Brillouin frequency shift-strain sensitivity coefficient of 0.5MHz / με and temperature compensation coefficient of 0.01GHz / ℃. The axial strain is calculated using the formula: Axial strain = (Measured frequency shift - Initial frequency shift) / Sensitivity coefficient - Temperature compensation. If the measured frequency shift at K1+100.0m is 10.525GHz, the initial frequency shift is 10.500GHz (25℃ reference), and the temperature is 26℃, the temperature compensation is (26-25)℃ × 0.01GHz / ℃ / 0.5MHz / με = 20με, and the axial strain is (10.525-10.500) × 10 9 Hz / 0.5×10 6 Hz / με-20με=50με-20με=30με. The three-dimensional layout trajectory of the optical cable (including the mileage of each monitoring point, circumferential angle, and radial burial depth) is retrieved from the tunnel construction BIM model. Using a line strain conversion tool, the strain average value is calculated in groups of 30° according to the circumferential angle, generating the strain distribution parameters along the surrounding rock (e.g., the strain average value of the 180°-210° group in section K1+100 is 32με). The center wavelength offset of each node was obtained using the FBG demodulation data analysis tool. For example, the offset of the right arch waist node at K1+100 is 0.19 nm. Substituting this into the calibrated FBG strain-temperature cross-sensitivity model (Δλ=1.2e+8.2T), and considering the synchronous ambient temperature of 25℃, the pure strain component was calculated using the formula (offset - 8.2T) / 1.2, yielding -12.5 με. Simultaneously, due to factors such as the coupling effect between the surrounding rock and the sensor after the FBG sensor is pasted / buried, the creep of the sensor material itself, and the aging deformation of the surrounding rock, the measured strain exhibits time-decay characteristics. This requires aging correction to eliminate systematic errors. This is achieved by using an exponential correction method. The attenuation-type aging correction model corrects the corresponding pure strain component. The correction formula is: ε_corr = ε × α(t), where ε_corr is the pure strain component after aging correction (unit: με); ε is the pure strain component before correction (here, -12.5με); α(t) is the aging attenuation coefficient (related to the number of monitoring days t, α(t)∈(0,1], the larger t is, the closer α(t) is to 1, and the weaker the attenuation effect). The specific expression for the aging attenuation coefficient α(t) is α(t) = e^(-β·t), where: t is the number of monitoring days (unit: days, for example, t=10 days); β is the attenuation coefficient factor to be fitted (unit: days). -1(Through calibration tests, the FBG sensor used in this embodiment is the same model as that used in tunnel monitoring (center wavelength 1550nm, strain sensitivity coefficient 1.2pm / με, temperature sensitivity coefficient 8.2pm / ℃). The test substrate is a specimen with the same geological properties as the surrounding rock of the tunnel (integrity index 0.75, compressive strength 0.6, cuboid specimen with dimensions of 100mm×100mm×500mm). The sensor bonding process is the same as on-site, using epoxy structural adhesive (model: 3MDP460). Before bonding, the substrate surface is sanded to remove rust and cleaned with alcohol. After bonding, it is cured at room temperature for 24 hours. At the same time, an electronic universal testing machine is used for axial static loading. The loading stress is 50% of the design working stress of the surrounding rock (corresponding to a strain of about 200με). After loading, the load is kept constant for 30 days (covering the aging period of on-site monitoring). The test environment temperature is kept constant at 25℃ (consistent with the ambient temperature of synchronous on-site monitoring) to eliminate the interference of temperature on the FBG wavelength shift. FBG is used. A demodulator (demodulation accuracy ±1pm) was used to collect the center wavelength offset of the sensor once a day at a fixed time (10:00 AM), and the loading time t (unit: day) was recorded synchronously. A total of 30 sets of data were collected (t=1,2,...,30 days). True strain benchmark acquisition: To obtain the "true strain" (the strain value unaffected by sensor aging) at each time point, two high-precision resistance strain gauges (model: BX120-3AA, sensitivity coefficient 2.10, accuracy ±0.1με) were attached to the same cross-section of the specimen as a strain value reference. The data from the resistance strain gauges were collected by a dynamic strain gauge and recorded synchronously with the FBG sensor data. The average value of the two strain gauges was taken as the true strain value ε_true(t) at each time point. FBG strain measurement value calculation: For the daily collected FBG center wavelength offset Δλ(t), temperature interference was subtracted according to the published FBG strain-temperature cross-sensitivity model (because the ambient temperature is constant, the wavelength offset caused by temperature is KT×T=8.2×10). -3 nm / ℃×25℃=0.205nm, consistent with the on-site calculation), yielding the wavelength shift caused by pure strain Δλ_ε(t)=Δλ(t)-0.205nm. Then, the FBG strain measurement value is calculated using the formula: ε_meas(t)=Δλ_ε(t) / Kε=Δλ_ε(t) / (1.2×10⁻⁶). -3The "aging attenuation ratio" k(t) is defined as the ratio of the measured strain value of the FBG to the true strain value, reflecting the measurement deviation caused by the aging effect of the sensor: k(t) = ε_meas(t) / ε_true(t). This ratio decreases as time t increases (due to creep at the sensor bonding interface, aging deformation of the surrounding rock, etc., causing the measured value to gradually approach the true value). Based on the general law of FBG sensor aging attenuation in geotechnical engineering (the attenuation ratio decreases exponentially with time), an exponential attenuation model is selected to describe the relationship between k(t) and t: k(t) = e^(-β×t), where β is the attenuation coefficient factor to be fitted (unit: days). -1 ), reflecting the decay rate; e is the natural constant (≈2.71828); the core of the least squares method is to minimize the "sum of squared residuals between the model prediction and the experimental measured value", that is, to find the optimal β such that: S(β)=Σ[k_meas(t)-e^(-β×t)] 2 →Minimum value where k_meas(t) is the measured attenuation ratio at each time point, and e^(-β×t) is the attenuation ratio predicted by the model. The corresponding fitting steps (including mathematical derivation) are as follows: Step 1) Take the natural logarithm of both sides of the model k(t)=e^(-β×t) to get: ln[k(t)]=-β×t Let y(t)=ln[k(t)], x(t)=t, then the linear model is: y(t)=-β×x(t) (the intercept is 0, because when t=0, k(0)=1, ln(1)=0, which is in line with the physical meaning). Step 2) Calculate the input data y(t) of the linear model, where the natural logarithm of each measured k_meas(t) is taken, as shown in Table 1 below:

[0028] Table 1. Examples of measured natural logarithms of k_meas(t)

[0029]

[0030] Step 3) Solve for β based on the linear model: The least squares estimation formula for the linear model y = -βx is: β = -[Σ(x(t)×y(t))] / [Σ(x(t)] 2 (Derivation basis: for the residual sum of squares S=Σ(y-(-βx)) 2 =Σ(y+βx) 2 Taking the derivative and setting it to 0, we get β = -Σ(xy) / Σ(x 2 Step 4) Substitute the complete 30-day data for calculation: For all x(t)=t and y(t)=ln[k_meas(t)] from t=1 to t=30, the specific values ​​of y(t) are shown in Table 1. Calculate the summation term: Numerator: Σ(x(t)×y(t))=1×0.0000+...+30×(-0.1665)≈-150.75; Denominator: Σ(x(t)...2 )=1 2 +2 2 +3 2 +...+30 2 =30×(30+1)×(2×30+1) / 6=9455 (calculated using the sum of squares formula); Substituting into the formula, we get: β=-[(-150.75)] / 9455≈150.75 / 9455≈0.01005 days -1 And by measuring the goodness of fit R 2 The formula used to verify the degree of matching between the model and the measured data is as follows: , where: y_meas(t) is the measured ln[k_meas(t)]; y_pred(t) is the model predicted ln[k_pred(t)]=-β×t (β=0.01005); Let y_meas(t) be the average value. Substituting the complete 30-day data, we calculate the residual sum of squares Σ(y_meas - y_pred). 2 ≈0.0008; Total sum of squares ;R 2 =1-0.0008 / 0.1617≈0.995≥0.98, which meets the goodness-of-fit requirement, proving that the model can accurately describe the aging decay law. After aging correction (substituting the aging coefficient α obtained by monitoring the 10th day), 10 =e^(-0.01005×10)=e^(-0.1005)≈0.905) becomes -12.5με×0.905≈-11.3με, generating the strain parameters of key points in the surrounding rock. Using a vector construction tool, the linear strain distribution parameters (average 32με, maximum 38με), point strain parameters (-11.3με, rate of change 0.8με / day), construction time series data, and geological attribute data (integrity index 0.75, compressive strength 0.6) are arranged in a fixed order to construct a 10-dimensional multidimensional feature vector of surrounding rock deformation, generating feature vectors for 50 cross sections.

[0031] Step S3: Based on the multidimensional feature vector of surrounding rock deformation, perform deformation pattern density clustering on the surrounding rock of different monitoring sections of the tunnel to classify the deformation pattern categories corresponding to uniform deformation, local abrupt deformation, and gradual accumulation. For different deformation pattern categories, introduce the process progress weight in the construction time sequence data and the lithological influence coefficient in the geological attribute data to dynamically predict the surrounding rock deformation and generate the dynamic prediction value of surrounding rock deformation for each monitoring section.

[0032] In this embodiment of the invention, the multidimensional feature vectors of surrounding rock deformation from 50 cross sections are clustered using a density peak clustering tool. The Euclidean distance between each sample is calculated (e.g., the feature vector distance between cross sections K1+100 and K1+110 is 3.86). A local density radius of 4, a density threshold of 5, and a distance threshold of 5 are set. Deformation modes are classified as follows: for the 20 cross sections in the K1+000-K1+100 segment, the density value ≥ 5 and the distance value < 5, it is a uniform deformation type; for the 8 cross sections in the K1+120-K1+150 segment, the density value < 5 and the distance value ≥ 5, it is a local abrupt deformation type; and for the 12 cross sections in the K1+200-K1+250 segment, the density value ≥ 5 and the distance value ≥ 5, it is a progressively cumulative deformation type. For uniformly deformable cross-sections, the excavation advance of 1.2m / day and the support strength score of 85 points were extracted from the construction sequence data. The weight of the process progress was calculated using the formula: (1.2 / 1.5)×0.4 + (85 / 100)×0.6 = 0.83. The integrity index of 0.75 and the compressive strength of 0.6 were extracted from the geological attribute data. The lithological influence coefficient was calculated using the formula: 0.75×0.5 + 0.6×0.5 = 0.675. These weights and coefficients were then substituted into the improved decision tree prediction model (training samples consisted of 200 sets of deformation data from the past 3 months). The optimized feature vector was input to generate dynamic deformation prediction values ​​for each cross-section, such as a 15-day prediction value of 12mm for cross-section K1+100 and a prediction value of 15mm for cross-section K1+120 (locally abrupt deformation).

[0033] Step S4: Compare the dynamic prediction values ​​of surrounding rock deformation with the monitoring data from the total station and multi-point displacement gauges deployed at the tunnel site, calculate the absolute error, relative error, and root mean square error between the predicted and measured values, and form a surrounding rock deformation prediction error matrix; perform iterative correction based on the surrounding rock deformation prediction error matrix, and aggregate regionally according to the tunnel mileage segment and deformation mode category to generate a visualization result of the spatiotemporal distribution of surrounding rock deformation in the mountain tunnel.

[0034] In this embodiment of the invention, displacement data from 50 cross-sections were collected using a total station at the tunnel site. For example, at cross-section K1+100, the arch crown was 11.8 mm, the right arch waist was 9.5 mm, and the left arch waist was 9.2 mm, with an average value of 10.2 mm. Deep displacement data was collected using a multi-point displacement gauge, with a depth of 1.5 m and a displacement of 7.9 mm. The actual deformation was calculated using the formula: 10.2 × 0.6 + 7.9 × 0.4 = 9.28 mm (with corresponding weights of 0.6 and 0.4). The predicted and actual values ​​were compared one-to-one, and an error calculation tool was used to calculate the following: at cross-section K1+100, the absolute error was 2.72 mm and the relative error was 29.3%; at cross-section K1+120, the predicted value was 15 mm, the actual value was 12.5 mm, the absolute error was 2.5 mm, and the relative error was 20%, forming a 50-row, 3-column error matrix. An error threshold was set (absolute error ≤ 2mm, relative error ≤ 20%), and the model parameters were iteratively corrected using the gradient descent algorithm (learning rate 0.01, 50 iterations). After correction, the predicted value for section K1+100 was 9.3mm with an error of 1.22mm, and for section K1+120 it was 13mm with an error of 0.5mm. Data was aggregated in 50m segments. The average deformation of uniform deformation in the K1+050-K1+100 segment was 9.19mm, accounting for 60%. Using a 3D geographic information system tool, the tunnel axis coordinates were imported, and the aggregated data was correlated to construct a 3D model. Blue represented uniform deformation and red represented local abrupt deformation. 1mm of deformation corresponds to a 0.1m bulge in the model. The K1+120-K1+150 segment (average deformation 12.5mm > 10mm warning value) was marked with red flashing. A visualization report containing the 3D model and error statistics table was generated.

[0035] As an embodiment of the present invention, reference Figure 2 As shown, Figure 1 A detailed flowchart of step S1 is shown below. In this embodiment, step S1 includes the following steps:

[0036] Step S11: In the monitoring area of ​​the surrounding rock of the mountain tunnel, a multi-dimensional sensing network is constructed based on the orientation of the joints and the distribution of stress concentration areas in the surrounding rock as determined by geological survey. The raw Brillouin frequency shift data along the optical cable is collected using the BOTDA demodulation system, and the amplitude attenuation rate and spectral bandwidth of the frequency shift signal are recorded simultaneously. The signal quality index is calculated based on the amplitude attenuation rate and spectral bandwidth. The raw center wavelength data of each sensing node is collected using the FBG demodulation system, and the reflected light intensity and polarization state are obtained simultaneously. The node sensing reliability coefficient is calculated based on the reflected light intensity and polarization state.

[0037] In this embodiment of the invention, BOTDA fully distributed sensing optical cables are laid in the monitoring area of ​​the surrounding rock of the mountain tunnel, based on the direction of the rock joints (30° northeast) and the stress concentration area (within 1.5m of the arch and left and right sidewalls) determined by geological survey. A spiral-encircling + radial-radial composite deployment method is used. The spiral encircling is centered on the tunnel axis, with one ring every 5m and a circumferential spacing of 0.5m. The radial radiation extends from the tunnel excavation outline into the surrounding rock, with one radiation direction every 2m, for a total of 8 directions. In the fault fracture zone (mileage K1+200-K1+210) and the lithological abrupt change interface (from sandstone to limestone at K1+150), one FBG quasi-distributed sensing node is deployed every 0.5m, forming a multi-dimensional sensing network. The BOTDA demodulation system collects raw Brillouin frequency shift data (range 10.5-11.5GHz) along the optical cable once per hour, and simultaneously records the amplitude attenuation rate (0-1) and spectral bandwidth (0.1-0.5GHz) of the frequency shift signal. The signal quality index is calculated using the formula (1-amplitude attenuation rate) × (0.5-spectral bandwidth) / 0.5. For example, when the attenuation rate is 0.3 and the bandwidth is 0.2GHz, the index is 0.7 × 0.6 = 0.42. The raw data of the center wavelength of each node (1540-1560nm) is collected synchronously by the FBG demodulation system to obtain the reflected light intensity (0-10000 counts) and polarization state (0-360°). The node sensing reliability coefficient is calculated according to the formula = (reflected light intensity / 10000) × [1-|polarization state-initial value| / 180]. For example, when the light intensity is 8000 and the polarization state deviation is 30°, the coefficient = 0.8×(1-30 / 180) = 0.67.

[0038] Step S12: Extract construction process transition time nodes, face advancement characteristic parameters, and support structure layout parameters through the tunnel BIM system to form tunnel construction progress data; obtain surrounding rock integrity index, rock mass structural shear strength, and groundwater permeability coefficient through ground-penetrating radar detection and borehole sampling analysis to form surrounding rock geological exploration data.

[0039] In this embodiment of the invention, the tunnel BIM system extracts the construction process transition time nodes, including the completion of the tunnel portal excavation (15 days after commencement), the start of initial support (30 days), and the start of secondary lining (90 days); extracts the face advancement characteristic parameters, such as the daily advancement of 1.5m for the K1+000-K1+100 section and 1.2m for the K1+100-K1+200 section; and extracts the support structure layout parameters, such as the anchor bolt length of 2.5m, spacing of 1m, and shotcrete thickness of 25cm, to form the tunnel construction progress data. Ground-penetrating radar was used to survey a 30m area in front of the tunnel face. Combined with borehole sampling analysis every 50m, the following data were obtained: surrounding rock integrity index (0.75 for K1+000-K1+100, 0.6 for K1+100-K1+200), shear strength of rock mass structural planes (30MPa for sandstone, 40MPa for limestone), and groundwater permeability coefficient (1×10⁻⁶ for K1+180-K1+210). -5 (cm / s), which constitutes the geological exploration data of the surrounding rock.

[0040] Step S13: Process the raw Brillouin frequency shift data to filter and remove abnormal data segments with a signal quality index below a threshold based on the signal quality index, and suppress noise in the effective data segments to obtain denoised Brillouin frequency shift data; calculate the initial distribution parameters of strain along the path based on the denoised Brillouin frequency shift data.

[0041] In this embodiment of the invention, by setting a signal quality index threshold of 0.3, a data filtering tool is used to process the original Brillouin frequency shift data, removing abnormal data segments with an index of 0.25 in the K1+205-K1+208 segment. Wavelet denoising is used to suppress noise in the valid data segments, decomposing them into 5 layers and retaining the first 3 layers of low-frequency coefficients for reconstruction, resulting in denoised Brillouin frequency shift data. Based on the linear relationship between Brillouin frequency shift and strain (1 GHz corresponds to 2000 με), the initial distribution parameters are calculated using the formula: strain along the path = (denoised frequency shift - initial frequency shift) × 2000. For example, at K1+100, the denoised frequency shift is 10.8 GHz, the initial frequency shift is 10.7 GHz, and the strain = 0.1 × 2000 = 200 με, generating the initial distribution parameters of strain along the path (one data point every 0.5 m).

[0042] Step S14: Process the original center wavelength data to remove failed node data with node sensing reliability coefficients below a threshold based on the node sensing reliability coefficient, and perform noise separation on the effective node data to obtain the denoised center wavelength data; calculate the initial value of node strain based on the denoised center wavelength data.

[0043] In this embodiment of the invention, by setting a node sensing reliability coefficient threshold of 0.5, a node screening tool is used to process the original center wavelength data, removing the failed node data with a coefficient of 0.45 at K1+203. A Kalman filter is used to separate noise from the valid node data, setting the process noise variance to 0.01 and the measurement noise variance to 0.02, and iteratively calculating the denoised center wavelength data. Based on the linear relationship between center wavelength and strain (1nm corresponds to 100με), the initial value is calculated according to the formula: node strain = (denoised wavelength - initial wavelength) × 100. For example, at K1+150, the denoised wavelength is 1550.2nm, the initial wavelength is 1550nm, and the strain = 0.2 × 100 = 20με, generating the initial node strain value (one data point per valid node).

[0044] Step S15: Standardize the tunnel construction progress data to convert the construction process transition time nodes into cumulative durations relative to the start of tunnel construction, and convert the face advancement characteristic parameters into advancement rate curves associated with mileage stations to generate standardized construction time sequence data; standardize the surrounding rock geological survey data to convert it into dimensionless geological feature vectors; establish a spatiotemporal reference axis based on the tunnel construction mileage stations, and spatiotemporally align the denoised fiber optic sensing data, which includes the initial distribution parameters of strain along the tunnel and the initial values ​​of nodal strain, the standardized construction time sequence data, and the geological feature vectors, to generate preprocessed fiber optic sensing basic data, construction time sequence data, and geological attribute data.

[0045] In this embodiment of the invention, the tunnel construction progress data is processed using a time-series standardization tool, converting the construction process transition time nodes into cumulative durations (15 days after tunnel entrance excavation is completed, 30 days after initial support begins); the face advancement characteristic parameters are interpolated according to mileage station to generate advancement rate curves (1.5 m / day at K1+000, 1.3 m / day at K1+150), forming standardized construction time-series data. The surrounding rock geological survey data is processed using a dimensionless transformation tool, dividing the surrounding rock integrity index by 1, the shear strength by 50 MPa, and the permeability coefficient by 1 × 10⁻⁶. -4 cm / s is converted into a geological feature vector (e.g., [0.6, 0.6, 0.1] at K1+100). A spatiotemporal reference axis is established with the tunnel construction mileage (K1+000-K1+300) as the horizontal axis and time (0-180 days) as the vertical axis. A spatiotemporal alignment tool is used to match the initial distribution parameters of strain along the tunnel (corresponding to the chainage), the initial values ​​of nodal strain (corresponding to the chainage and time), the standardized construction time series data (corresponding to the time), and the geological feature vector (corresponding to the chainage) according to the chainage and time to generate preprocessed fiber optic sensing basic data, construction time series data, and geological attribute data.

[0046] Step S2 includes the following steps:

[0047] Step S21: By retrieving the Brillouin frequency shift sequence data output by the BOTDA demodulation system from the preprocessed fiber optic sensing basic data, and combining the Brillouin frequency shift-strain sensitivity coefficient and temperature compensation coefficient calibrated by the BOTDA demodulation system, the axial strain value of each monitoring point of the sensing optical cable is calculated.

[0048] In this embodiment of the invention, a data retrieval tool is used to obtain the Brillouin frequency shift sequence data output by the BOTDA demodulation system in the K1+000-K1+500 segment from the preprocessed fiber optic sensing base data. The data acquisition interval is 0.5m, with a total of 1001 monitoring points and a frequency shift range of 10.500-10.550GHz. The Brillouin frequency shift-strain sensitivity coefficient calibrated by the BOTDA demodulation system is retrieved as 1.0GHz / 2000με (i.e., 0.5MHz / με), and the temperature compensation coefficient is retrieved as 0.01GHz / ℃. The synchronous temperature values ​​(25±2℃) of each monitoring point were obtained through environmental temperature monitoring data. The axial strain was calculated using the formula: Axial strain = (Measured frequency shift - Initial frequency shift) / Sensitivity coefficient - Temperature compensation. If the measured frequency shift at K1+100.0m is 10.525GHz, the initial frequency shift is 10.500GHz (25℃ reference), and the temperature is 26℃, then the temperature compensation is (26-25)℃×0.01GHz / ℃ / 0.5MHz / με = 20με, and the axial strain is (10.525-10.500)×10 9 Hz / 0.5×10 6 Hz / με-20με=50με-20με=30με, thus generating the axial strain values ​​at each monitoring point.

[0049] Step S22: Extract the strain gradient change rate based on the spatial distribution characteristics corresponding to the axial strain value, specifically, the strain gradient change rate = strain difference between adjacent monitoring points / optical cable laying spacing, and identify the strain abrupt change segment and smooth segment corresponding to each monitoring point along the optical cable through the strain gradient change rate;

[0050] In this embodiment of the invention, based on the spatial distribution characteristics of the axial strain value, a gradient calculation tool is used to calculate the strain gradient change rate according to the formula: |current monitoring point strain value - adjacent monitoring point strain value| / optical cable laying spacing (0.5m). For example, if the strain at K1+100.0m is 30με and the strain at K1+100.5m is 35με, the gradient change rate is |35-30| / 0.5=10με / m; if the strains at K1+100.5m and K1+101.0m are 35με and 36με, the gradient change rate is 2με / m. A strain gradient change rate threshold of 5 με / m is set. When the gradient change rate is ≥5 με / m, it is identified as a strain abrupt change segment. For example, the gradient change rates of the K1+120.0m-K1+121.0m segment are 8 με / m, 12 με / m, and 9 με / m, respectively, and are identified as abrupt change segments. When the gradient change rate is <5 με / m, it is identified as a gradual segment. For example, most gradient change rates of the K1+100.0m-K1+110.0m segment are 2-4 με / m, and are identified as gradual segments. This completes the identification of strain segments throughout the tunnel.

[0051] Step S23: Based on the strain abrupt change section and smooth section corresponding to each monitoring point along the optical cable, the three-dimensional layout trajectory of the optical cable corresponding to the sensing optical cable is obtained by combining the tunnel construction BIM model. Based on the three-dimensional layout trajectory of the optical cable, the axial strain value of each monitoring point of the sensing optical cable is statistically analyzed to generate the strain distribution parameters along the surrounding rock.

[0052] In this embodiment of the invention, the 3D deployment trajectory of the sensing optical cable is retrieved from the tunnel construction BIM model to obtain the mileage station number, cross-sectional circumferential angle, and radial burial depth parameters of each monitoring point (e.g., K1+100.0m corresponds to a circumferential angle of 180° and a radial burial depth of 0.35m). A linear strain statistical tool is used to divide the tunnel into statistical units every 10m (a total of 50 units). Within each unit, the tunnel is grouped by circumferential angle (0°-90°, 90°-180°, 180°-270°, 270°-360°), and the average, maximum, and minimum values ​​of the axial strain for each group are calculated. For example, in unit K1+100-K1+110, the 180°-270° group (sidewall area) has an average strain of 32με, a maximum of 38με, and a minimum of 28με. The statistical results of each unit group are integrated to generate a rock mass strain distribution parameter along the tunnel, including the mileage range, circumferential area, and strain statistical values.

[0053] Step S24: The center wavelength offset of each sensing node is obtained by parsing the fiber optic sensing basic data through the FBG demodulation system, and the pure strain component is separated according to the FBG strain-temperature cross-sensitivity compensation model to generate the point strain parameters of the key points of the surrounding rock; the strain distribution parameters along the surrounding rock are combined with the point strain parameters of the key points of the surrounding rock to generate the deformation characteristic parameters of the surrounding rock.

[0054] In this embodiment of the invention, the fiber optic sensing data is analyzed using an FBG demodulation system to obtain the center wavelength offset of 200 FBG sensing nodes (e.g., 0.19nm offset for the right arch waist node at K1+100 and 0.22nm offset for the arch crown node at K1+120). Based on the constructed calibrated FBG strain-temperature cross-sensitivity model (Δλ=1.2e+8.2T), the synchronous ambient temperature data (e.g., 25℃ for node K1+100) is substituted, and the pure strain component is calculated using the formula e=(Δλ-8.2T) / 1.2. For example, for node K1+100, e=(190pm-8.2×25pm) / 1.2≈(190-205) / 1.2≈-12.5με. After aging correction (attenuation coefficient of 0.905 on the 10th day of monitoring), -11.3με is obtained, generating the point strain parameters of key points in the surrounding rock. By combining the strain distribution parameters along the path of the surrounding rock with the point strain parameters, a characteristic parameter for the deformation of the surrounding rock is formed, which includes the statistical value of linear strain and the specific value of point strain.

[0055] Step S25: Based on the surrounding rock deformation characteristic parameters, combined with construction time series data and geological attribute data, construct a multidimensional feature vector of surrounding rock deformation that includes linear strain distribution, point strain parameters, construction correlation laws and geological coupling relationships.

[0056] In this embodiment of the invention, the construction process nodes corresponding to each statistical unit are extracted from the construction time sequence data (e.g., the initial support of unit K1+100-K1+110 is completed on the 30th day, and the secondary lining is completed on the 90th day), and the face advancement rate (1.2m / day) are extracted; the surrounding rock integrity index (0.75), rock mass shear strength (30MPa), and permeability coefficient (1×10⁻⁶) of each unit are extracted from the geological attribute data. -5 cm / s. Using a vector construction tool, based on each statistical unit, the following parameters were calculated: linear strain distribution parameters (mean 32με, maximum 38με, minimum 28με), point strain parameters (-11.3με, rate of change 0.8με / day), construction correlation patterns (support completion time 30 days, advancement rate 1.2m / day), and geological coupling relationships (integrity index 0.75, shear strength 30MPa, permeability coefficient 1×10⁻⁶). -5 The values ​​(cm / s) are arranged in a fixed order to form a multidimensional feature vector of surrounding rock deformation. For example, the vector of unit K1+100-K1+110 is [32,38,28,31.5,0.8,30,1.2,0.75,30,1e-5,180°,0.35] (including circumferential angle and radial burial depth), thus completing the feature vector construction of 50 units for the entire tunnel.

[0057] Step S23, which involves statistically analyzing the linear strain distribution of the axial strain values ​​at each monitoring point of the sensing optical cable based on the three-dimensional deployment trajectory of the optical cable, includes the following steps:

[0058] Step S2301: Based on the three-dimensional deployment trajectory of the optical cable and combined with the tunnel construction BIM model, obtain the circumferential curvature radius and radial burial depth parameters of the sensing optical cable at each tunnel section.

[0059] In this embodiment of the invention, a BIM model data integration tool is used to retrieve the 3D cable deployment trajectory data for each section (50 sections per 10m) from the tunnel construction BIM model from K1+000 to K1+500. Combined with previously corrected parameters, the circumferential radius of curvature and radial depth parameters of 20 monitoring points for each section are extracted. Taking section K1+100 as an example, the corrected circumferential radius of curvature for the arch crown monitoring point is 27.11m, and the radial depth is 0.396m; the circumferential radius of curvature for the right arch waist monitoring point is 26.01m, and the radial depth is 0.3476m; the circumferential radius of curvature for the sidewall monitoring point is 25.8m, and the radial depth is 0.32m. Using a coordinate association tool, the parameters of each monitoring point are bound to the section mileage station number and circumferential angle, forming a four-dimensional parameter table containing "mileage station number - circumferential angle - circumferential radius of curvature - radial depth," thus completing the acquisition of all section parameters.

[0060] Step S2302: Calculate the optical cable bending correction coefficient and the surrounding rock constraint coefficient based on the circumferential curvature radius and radial burial depth parameters, where the optical cable bending correction coefficient = 1 - circumferential reference value / circumferential curvature radius, and the surrounding rock constraint coefficient = radial burial depth parameter × rock mass elastic modulus;

[0061] In this embodiment of the invention, the circumferential reference value of the optical cable is set to 20m (a fixed value determined based on the material properties of the optical cable). A coefficient calculation tool is used to calculate the optical cable bending correction coefficient according to the formula: Optical Cable Bending Correction Coefficient = 1 - Circumferential Reference Value / Circumferential Curvature Radius. For example, the circumferential curvature radius of the monitoring point at the K1+100 section arch crown is 27.11m, and the bending correction coefficient = 1 - 20 / 27.11 ≈ 1 - 0.738 ≈ 0.262; the circumferential curvature radius of the monitoring point at the right arch waist is 26.01m, and the bending correction coefficient = 1 - 20 / 26.01 ≈ 0.231. The rock mass elastic modulus E of the corresponding monitoring point is retrieved from the BIM model rock mass mechanics parameter library (arch crown E = 30GPa = 3 × 10⁻⁶). 4 (MPa), calculated using the formula: Surrounding rock restraint coefficient = Radial burial depth parameter × Rock mass elastic modulus. For example, if the radial burial depth of the arch crown is 0.396m, the restraint coefficient = 0.396 × 3 × 10 4 ≈1.188×10 4 MPa·m; radial embedment depth of sidewall 0.32m, E=3×10 4 MPa, constraint coefficient = 0.32 × 3 × 10 4 =9.6×10 3 MPa·m, generating the bending correction coefficient and constraint coefficient for each monitoring point.

[0062] Step S2303: Couple the axial strain value of each monitoring point of the sensing optical cable with the optical cable bending correction coefficient and the surrounding rock constraint coefficient to obtain the corrected axial strain value; based on the corrected axial strain value, extract the strain accumulation and strain change rate parameters corresponding to each monitoring point.

[0063] In this embodiment of the invention, the original axial strain values ​​(in με) of each monitoring point are obtained from a distributed optical fiber sensing system. For example, the original axial strain at the arch crown of section K1+100 is 350 με, the right arch waist is 300 με, and the sidewall is 250 με. Using a strain coupling calculation tool, the corrected axial strain value is calculated as follows: Original axial strain value × Optical cable bending correction coefficient × (Rock constraint coefficient / 1 × 10⁻⁶) 4 (Divided by 1×10) 4 To unify the order of magnitude of the constraint system, the axial strain after the crown correction is calculated as follows: 350 × 0.262 × (1.188 × 10⁻⁶). 4 / 1×10 4 )≈350×0.262×1.188≈109.3με; Axial strain after correction of the right arch waist =300×0.231×(9.8×10) 3 / 1×10 4 )≈300×0.231×0.98≈67.9με. Based on the corrected axial strain value, the cumulative strain during the monitoring period (30 days) was extracted using data statistics tools (109.3με at the top of the arch and 67.9με at the right waist of the arch). The rate parameter was calculated according to the formula strain change rate = cumulative strain / number of monitoring days (109.3 / 30≈3.64με / day at the top of the arch and 67.9 / 30≈2.26με / day at the right waist of the arch).

[0064] Step S2304: Calculate the radial transformation weight using the spatial angle parameters of the optical cable layout, and calculate the circumferential transformation weight using the Poisson's ratio parameters of the rock mass. Simultaneously, based on the radial transformation weight and the circumferential transformation weight, convert the strain accumulation and strain change rate parameters corresponding to each monitoring point into the corresponding radial linear strain and circumferential linear strain.

[0065] In this embodiment of the invention, the spatial angle parameters θ of each monitoring point are extracted (arch top θ=90°, right arch waist θ=180°, sidewall θ=270°), and a weight calculation tool is used to calculate the radial transformation weight using the formula cosθ (θ is converted to radians). For example, arch top cos90°=0, right arch waist cos180°=-1, and sidewall cos270°=0. The Poisson's ratio μ of the rock mass is retrieved (arch top μ=0.25), and the circumferential transformation weight is calculated using the formula μ×sinθ. For example, arch top sin90°=1, circumferential weight=0.25×1=0.25; right arch waist sin180°=0, circumferential weight=0. Multiply the cumulative strain and rate parameters by their respective weights to obtain the radial strain (109.3×0=0με at the crown, 67.9×(-1)=-67.9με at the right arch waist), the circumferential strain (109.3×0.25≈27.3με at the crown, 67.9×0=0με at the right arch waist), the radial strain rate (3.64×0=0με / day at the crown, 2.26×(-1)=-2.26με / day at the right arch waist), and the circumferential strain rate (3.64×0.25≈0.91με / day at the crown).

[0066] Step S2305: Perform spatial interpolation on the radial and circumferential strains to generate a tunnel strain distribution cloud map; extract the spatial coordinates of strain extreme points, the density of strain contour lines, and the main vector parameters of strain direction from the tunnel strain distribution cloud map, and integrate them to form the surrounding rock strain distribution parameters that include strain amplitude, distribution pattern, and spatial gradient.

[0067] In this embodiment of the invention, a spatial interpolation tool is used, with the tunnel mileage station (K1+000-K1+500) as the x-axis, the circumferential angle of the cross-section (0°-360°) as the y-axis, and the linear strain value as the z-axis, to perform Kriging interpolation on the radial and circumferential linear strain data of all monitoring points, generating a tunnel strain distribution cloud map (the corresponding cross-section is shown in the figure). Figure 3 (As shown in the diagram) The right arch waist region in the radial cloud map is a blue negative strain zone, and the arch crown region in the circumferential cloud map is a red positive strain zone. Extreme strain points were extracted using cloud map analysis tools: the radial extreme strain point is located at the right arch waist at K1+100 (-67.9 με), with spatial coordinates (K1+100, 180°, 0.3476 m); the circumferential extreme strain point is located at the arch crown at K1+120 (32.1 με). The density of strain contour lines was statistically analyzed: the contour line spacing in the K1+080-K1+120 segment was 0.5 m (dense area), and the spacing in the K1+200-K1+250 segment was 2 m (sparse area). The principal vector of strain direction was calculated: the radial strain principal vector is along the tunnel radial direction (pointing towards the center), and the circumferential strain principal vector is along the tunnel circumferential direction (clockwise). The extreme point coordinates, contour line density, and principal vector parameters were integrated to form a table of strain distribution parameters along the surrounding rock.

[0068] Step S2301 includes the following steps:

[0069] The three-dimensional coordinate dataset corresponding to the three-dimensional cable laying trajectory is retrieved from the tunnel construction BIM model. This dataset includes the mileage station, cross-sectional circumferential angle, and initial radial depth of the cable along the tunnel axis. Based on the three-dimensional coordinate dataset, the coordinates of the spatial discrete points of the cable in each tunnel cross section are extracted, and the arc length distance and angle deviation between adjacent discrete points are calculated to generate the spatial topology parameters corresponding to the cable laying trajectory.

[0070] In this embodiment of the invention, by employing a BIM model data retrieval tool, the K1+000-K1+500 segment of the sensor optical cable deployment layer is located in the tunnel construction BIM model. At 0.5m intervals, a 3D coordinate dataset corresponding to the 3D deployment trajectory of the optical cable is extracted. Each coordinate point contains three core parameters: the mileage marker along the tunnel axis (e.g., K1+000.0, K1+000.5…K1+500.0, accurate to 0.1m), the cross-sectional circumferential angle (measured clockwise from the center of the tunnel cross-section, with 90° at the arch crown, 180° at the right arch waist, and 270° at the invert, accurate to 0.1°), and the initial radial depth (the distance extending from the tunnel excavation outline into the surrounding rock, ranging from 0.2-3.0m, 1.0m for the arch crown optical cable and 0.8m for the sidewall optical cable, accurate to 0.01m). A total of 1001 coordinate point data points are extracted. Using a coordinate filtering tool, the system divides the area into 50 sections every 10 meters based on mileage markers. For each section, 20 spatial discrete point coordinates are extracted (uniformly distributed along the circumferential direction, with an angular difference of 18° between adjacent points). A spatial distance calculation tool is used, based on the arc length formula in spherical coordinates (arc length = radial depth × angular deviation (radians)), to calculate the arc length distance between adjacent discrete points. For example, in section K1+100, two points have a radial depth of 0.9m and an angular deviation of 18° (0.314 radians), so the arc length distance is approximately 0.9 × 0.314 ≈ 0.283m. Simultaneously, the angular deviation (the difference between the actual angular difference and 18°) is calculated. For instance, if the actual angular difference between two points is 18.2°, the deviation is 0.2°. This generates the spatial topological relationship parameters for each section.

[0071] A continuous curve model of the optical cable is constructed based on spatial topological parameters. The coordinates of discrete points in space are then fitted using the least squares method to obtain the circumferential layout curve equation of the optical cable in each tunnel section. The curvature value of each point is calculated based on the circumferential layout curve equation, where the curvature value is calculated as: second derivative of the curve / (1 + square of the first derivative). 3 / 2 The initial value of the circumferential radius of curvature of the optical cable at the corresponding position is obtained by back-calculation based on the curvature value;

[0072] In this embodiment of the invention, a continuous curve model of the optical cable is constructed using a curve fitting tool based on spatial topological relationship parameters. Taking the K1+100 section as an example, the circumferential angles of 20 discrete points are set as the independent variable x (radians), and the radial depth is set as the dependent variable y (m), resulting in discrete coordinates from (x1, y1) to (x20, y20). The quadratic polynomial y=ax is then fitted using the least squares method. 2 +bx+c, construct the sum of squared errors function S=Σ(yi-(axi) 2 +bxi+c)) 2 (i=1 to 20), take the partial derivatives with respect to a, b, and c and set them to 0. Solve the system of three linear equations to obtain coefficients a=0.02, b=0.15, and c=0.7, thus obtaining the equation of the circumferential curve y=0.02x. 2 +0.15x + 0.7. Use a derivative calculation tool to find the first derivative y' = 0.04x + 0.15 and the second derivative y'' = 0.04, and substitute them into the curvature value formula (curvature value = |y''| / (1 + (y')). 2 (3 / 2)), for example, when x = 0.523 radians (30°), y' = 0.171, and the curvature value = 0.04 / (1.029)^1.5 ≈ 0.0383m -1 Based on the reciprocal relationship between curvature value and radius of curvature, the initial value of circumferential radius of curvature is calculated as 1 / 0.0383≈26.1m. Similarly, the radius of curvature calculation is completed for all discrete points of the cross-section.

[0073] The tunnel cross-section design contour data, including the excavation contour line, the initial contour line, and the thickness parameters of the support structure, are retrieved from the tunnel construction BIM model. The normal distance between the three-dimensional laying trajectory of the optical cable and the excavation contour line and the initial contour line is calculated. The difference between the normal distance and the thickness parameters of the support structure is combined to obtain the initial value of the radial burial depth of the optical cable in the surrounding rock.

[0074] In this embodiment of the invention, a BIM model data extraction tool is used to retrieve the design outline data of each section from the tunnel construction BIM model: the excavation outline is a circle with a radius of 5.5m (equation x). 2 +y 2 =5.5 2 The initial outline is a circle with a radius of 5.2m (x 2 +y 2 =5.2 2The thickness parameters of the support structure are calculated (0.15m for shotcrete, 2.5m for anchor bolts, 0.1m for steel arch, and a total thickness of 0.35m). Using a spatial distance calculation tool, based on the formula for the normal distance from a point to a circle (normal distance = |distance from point to center - radius of outline|), the normal distances between discrete points of the optical cable and the excavation and initial outlines are calculated. For example, at section K1+100, the distance from a discrete point to the center is 5.3m, the distance to the excavation outline is |5.3 - 5.5| = 0.2m, and the distance to the initial outline is |5.3 - 5.2| = 0.1m. According to the formula: initial radial burial depth = normal distance to the initial outline + support structure thickness, the initial burial depth at this point is calculated to be 0.1 + 0.35 = 0.45m. This completes the calculation of the initial burial depth for all discrete points across all sections.

[0075] Construction deviation parameters during the optical cable laying process are introduced. These parameters are obtained through comparative analysis of on-site optical cable laying records and tunnel construction BIM models, including circumferential angle deviation and radial depth deviation. The circumferential angle deviation is substituted into the initial value of the circumferential curvature radius for correction, resulting in the corrected circumferential curvature radius, specifically: corrected circumferential curvature radius = initial value of circumferential curvature radius × (1 + circumferential angle deviation coefficient).

[0076] In this embodiment of the invention, on-site construction record acquisition tools are used to obtain on-site record data of optical cable laying (such as the measured values ​​of circumferential angle and radial depth at each discrete point of section K1+100), and these values ​​are compared with the design values ​​in the BIM model to obtain construction deviation parameters. Circumferential angle deviation = measured angle - design angle. For example, if the design angle of a discrete point is 90° and the measured angle is 91.5°, the deviation is 1.5°, which is converted to a circumferential angle deviation coefficient = deviation value / 360° = 1.5 / 360 ≈ 0.0042; radial depth deviation = measured depth - design depth. For example, if the design depth of a point is 0.45m and the measured depth is 0.47m, the deviation is 0.02m. Using a parameter correction tool, the circumferential radius of curvature is corrected according to the formula: Circumferential radius of curvature = Initial value of circumferential radius of curvature × (1 + Circumferential angle deviation coefficient). For example, if the initial radius of curvature at a point is 27.0m, the corrected value is 27.0 × (1 + 0.0042) ≈ 27.11m; if the angle deviation at a point is -1.2°, the coefficient is -0.0033, and the initial radius is 26.1m, the corrected value is 26.1 × (1 - 0.0033) ≈ 26.01m. The radius of curvature correction is then completed.

[0077] By combining the convergence deformation after the surrounding rock excavation, the initial value of the radial burial depth is dynamically adjusted so as to calculate the corrected radial burial depth parameter through the coupling relationship between the convergence deformation and the radial burial depth; and the circumferential curvature radius and radial burial depth parameters of the sensing optical cable at each tunnel section are integrated and generated.

[0078] In this embodiment of the invention, the convergence deformation of each section of the surrounding rock after excavation is measured using a tunnel convergence monitoring tool (total station). For example, the arch crown converges by 5mm and the sidewall converges by 3mm at section K1+100 (converted to 0.005m and 0.003m, respectively). Based on the coupling relationship between the convergence deformation and the radial burial depth (the burial depth decreases as convergence increases, with a coupling coefficient of 0.8), the radial burial depth is calculated using the formula corrected by a dynamic adjustment tool: Radial burial depth = Initial burial depth - (Convergence deformation × Coupling coefficient). For example, if the initial burial depth of the discrete point of the arch crown is 0.45m, the corrected value is 0.45 - (0.005 × 0.8) = 0.446m; if the initial burial depth of the discrete point of the sidewall is 0.4m, the corrected value is 0.4 - (0.003 × 0.8) = 0.3976m. Using data integration tools, the corrected circumferential curvature radii (e.g., 27.11m, 26.01m) of each cross-section discrete point are associated with radial burial depth parameters (e.g., 0.446m, 0.3976m) according to mileage station and circumferential angle to generate a complete parameter dataset for each cross-section of the sensing optical cable.

[0079] Step S2304 includes the following steps:

[0080] The spatial azimuth parameter θ of each monitoring point is extracted from the three-dimensional laying trajectory of the optical cable, and combined with the rock mechanics parameter library in the tunnel construction BIM model, the surrounding rock Poisson's ratio μ and elastic modulus E of the corresponding monitoring point are retrieved; the radial foundation weight cosθ is calculated based on the spatial azimuth parameter θ, and the circumferential foundation weight μ·sinθ is calculated based on the surrounding rock Poisson's ratio μ.

[0081] In this embodiment of the invention, spatial azimuth parameters θ of 20 monitoring points at section K1+100 are obtained from the three-dimensional cable laying trajectory using a spatial parameter extraction tool (measured clockwise from the tunnel axis, ranging from 0° to 360°, e.g., θ=90° for the arch crown monitoring point, θ=180° for the right arch waist, and θ=270° for the sidewall). The surrounding rock mechanics parameters of the corresponding monitoring points are retrieved from the rock mechanics parameter library of the tunnel construction BIM model: the arch crown and sidewall areas are moderately weathered sandstone with Poisson's ratio μ=0.25 and elastic modulus E=30GPa; the invert arch area is silty clay with μ=0.35 and E=15GPa. Using a trigonometric function calculation tool, the radial foundation weight cosθ (θ converted to radians) is calculated based on θ. For example, when θ=90° (π / 2 radians), cosθ=0; when θ=180° (π radians), cosθ=-1; and when θ=0°, cosθ=1. The circumferential base weight μ·sinθ is calculated based on μ. For example, at the arch top monitoring point, μ=0.25 and sin90°=1, the circumferential base weight = 0.25×1=0.25; at the right arch waist, μ=0.25 and sin180°=0, the circumferential base weight = 0.25×0=0. The base weight calculation for all monitoring points is completed.

[0082] Obtain the measured tension and design tension during optical cable laying, and calculate the tension influence coefficient = measured tension / design tension; multiply the radial foundation weight and circumferential foundation weight by the tension influence coefficient to obtain the corrected radial conversion weight and circumferential conversion weight.

[0083] In this embodiment of the invention, the actual measured tension (in kN) during optical cable laying is collected using a tension measuring tool. For example, the designed laying tension of the optical cable at section K1+100 is 5 kN, and the actual measured tension is 5.2 kN. Using the formula Tension Influence Coefficient = Measured Tension / Design Tension, the calculated tension influence coefficient for this section is 5.2 / 5 = 1.04. If the actual measured tension at a section is 4.8 kN and the designed tension is 5 kN, the influence coefficient is 4.8 / 5 = 0.96. A weight correction tool is used to multiply the radial and circumferential foundation weights of each monitoring point by the tension influence coefficient to obtain the corrected weights: for example, the radial foundation weight of the arch crown monitoring point is 0, and it remains 0 after multiplying by 1.04; the circumferential foundation weight is 0.25, and it becomes 0.26 after multiplying by 1.04; the radial foundation weight of the right arch waist is -1, and it becomes -1.04 after multiplying by 1.04; the circumferential foundation weight is 0, and it remains 0 after multiplying by 1.04. This ensures that the weights can reflect the influence of the actual laying tension on the sensor.

[0084] The total axial strain increment at each monitoring point is extracted from the strain accumulation, and the radial strain accumulation component is calculated by combining the corrected radial transformation weight. At the same time, the circumferential strain accumulation component is calculated by combining the circumferential transformation weight.

[0085] In this embodiment of the invention, the total axial strain increment (in με) at each monitoring point of section K1+100 is extracted from the strain accumulation data of the distributed optical fiber sensing system. For example, the total axial strain increment at the arch crown monitoring point is 300 με, at the right arch waist is 250 με, and at the sidewall is 200 με. Using a strain component calculation tool and combined with the corrected radial transformation weight, the radial strain accumulation component is calculated as: total axial strain increment × radial transformation weight. For example, at the arch crown with a radial weight of 0, the radial strain accumulation component is 300 × 0 = 0 με; at the right arch waist with a radial weight of -1.04, the radial strain accumulation component is 250 × (-1.04) = -260 με. Simultaneously, the cumulative circumferential strain component is calculated by combining the circumferential transformation weight. For example, the circumferential weight of the arch crown is 0.26, and the cumulative circumferential strain component is 300 × 0.26 = 78 με; the circumferential weight of the right arch waist is 0, and the cumulative circumferential strain component is 250 × 0 = 0 με, thus generating the radial and circumferential strain cumulative components of each monitoring point.

[0086] Dynamic correction is performed based on the strain change rate parameter and a time-dimensional rock mass creep correction factor to generate the corrected rate parameter. The corrected rate parameter is then coupled with radial and circumferential transformation weights to obtain the radial strain change rate and the circumferential strain change rate.

[0087] In this embodiment of the invention, strain change rate parameters (unit με / day) for each monitoring point are extracted from strain monitoring data, such as the arch crown strain change rate = 10 με / day and the sidewall strain change rate = 8 με / day. Based on the creep characteristics of rock mass, a time-dimensional rock mass creep correction factor is introduced (calculated according to the monitoring duration; correction factor = 0.95 on the 10th day and 0.90 on the 30th day). The specific derivation of the rock mass creep correction factor is as follows: due to the creep characteristics of rock mass (the phenomenon that strain slowly increases and the rate gradually decreases under long-term load), the original strain change rate (e.g., 10 με / day) is an instantaneous monitoring value, which does not consider the rule that "the longer the time, the more obvious the rate decrease due to creep." Therefore, the creep correction factor is a decay coefficient positively correlated with the monitoring duration (value 0~1), aiming to make the rate parameters more closely match the actual deformation trend after long-term creep of the rock mass through time-dimensional correction. The potential is determined, and constraints on the value range are constructed: ① Correction factor ∈ (0,1] (since creep only causes rate decay and does not amplify, the corrected rate ≤ the original rate); ② Monotonicity constraint: the longer the monitoring time, the smaller the correction factor. After excavation disturbance, the strain change rate of the rock mass will gradually decay over time (consistent with the "decaying creep stage" characteristic of rock creep). Its decay law is described by a logarithmic decay model (verified by indoor tests and field measurements, this model can fit the creep rate decay characteristics of the surrounding rock in this example). The specific model expression is: f(t) = 1−C·ln(t+D), where the definitions and physical meanings of each parameter are clear, and f(t) is the rock mass creep correction factor on the monitoring day t (dimensionless, with a value of...). The range is (0,1], the closer to 1, the weaker the creep decay effect; t is the monitoring duration (unit: days, in this example the target duration t=10 days, t=30 days); C is the creep decay coefficient (dimensionless, determined by the rock mass mechanical properties, requires experimental calibration); D is the time offset coefficient (dimensionless, used to correct the fitting deviation in the initial stage of monitoring, requires experimental calibration); ln(·) is the natural logarithm function (base e≈2.71828). To ensure that the correction factor matches the rock mass characteristics of this embodiment, parameters C and D are determined through the following calibration test: test conditions are matched, and rock mass samples consistent with the implementation scenario of this invention are selected (geological properties: integrity index 0.75, compressive strength 0.6, and...). (The surrounding rock parameters were unified as described above), and the samples were processed into standard rock samples (size: Φ50mm×100mm). A long-term creep test was conducted using a uniaxial compression creep test, applying a constant stress equivalent to the disturbance caused by tunnel excavation (stress level was 30% of the rock mass compressive strength). Continuous monitoring was performed for 90 days, recording the daily rock mass strain change rate v_t (unit: με / day). The baseline rate was determined using the strain change rate v_1 = 25.6 με / day on the first day of monitoring as the "initial rate baseline." Based on the physical definition of the creep correction factor (correction factor = actual rate at a certain moment / initial baseline rate), the measured correction factor for key nodes during the test was calculated: the measured rate on day 10, v_10 = 24.3 με / day, measured correction factor f_test(10)=24.3 / 25.6≈0.95; measured rate on day 30 v_30=23.0 με / day, measured correction factor f_test(30)=23.0 / 25.6≈0.90), measured rate on day 90 v_90=21.8 με / day, measured correction factor f_test(90)=21.8 / 25.6≈0.85; parameter fitting optimization, substituting the three sets of key node data (t=10,f=0.95; t=30,f=0.90; t=90,f=0.85) into the logarithm The attenuation model was solved by fitting parameters C and D using the least squares method. Substituting these parameters into the calculation, we obtained D=2.0. The calculated C=0.05 / ln(12)≈0.05 / 2.4849≈0.0201 (final calibration C=0.02, simplified calculation). Therefore, the specific calculation model for the rock mass creep correction factor is f(t)=1−0.02·ln(t+2.0)). The rate parameter after correction using the rate correction tool is = original rate × creep correction factor. The calculated rate after correction of the crown is 10×0.95=9.5με / day and the rate of the sidewall is 8×0.95=7.6με / day. The corrected rate parameters are coupled with the radial and circumferential transformation weights, respectively, to obtain the radial strain rate of change = corrected rate × radial weight. For example, with a radial weight of 0 at the crown, the radial rate = 9.5 × 0 = 0 με / day; with a radial weight of -1.04 at the right arch waist, the radial rate = 9.5 × (-1.04) = -9.88 με / day. The circumferential strain rate of change = corrected rate × circumferential weight. For example, with a circumferential weight of 0.26 at the crown, the circumferential rate = 9.5 × 0.26 = 2.47 με / day.

[0088] The radial strain cumulative component and radial strain change rate are verified by time integration, and the circumferential strain cumulative component and circumferential strain change rate are also verified by time integration. The radial and circumferential transformation weights are then optimized based on the integral result deviation rate, where the integral result deviation rate = |integral value - cumulative component| / cumulative component, and the optimized weight = original weight × (1 - integral result deviation rate). The optimized radial and circumferential linear strains are then output.

[0089] In this embodiment of the invention, a time integration tool is used to integrate the radial strain change rate over time (the integration time is the monitoring period corresponding to the strain accumulation, such as 30 days), and the integral value is calculated: for example, if the radial rate of the right arch waist is -9.88 με / day, the integral value = -9.88 × 30 ≈ -296.4 με. The integral value is compared with the cumulative radial strain component (-260 με), and the deviation rate is calculated using the formula: integral value - cumulative component | / cumulative component. The deviation rate is calculated as: |-296.4 - (-260)| / |-260| ≈ 36.4 / 260 ≈ 0.14. The optimized weight is calculated using the formula: optimized weight = original radial weight × (1 - deviation rate). The optimized radial weight of the right arch waist is calculated as: -1.04 × (1 - 0.14) = -0.8944. Similarly, the circumferential strain is verified by integration. For example, if the circumferential rate at the crown is 2.47 με / day, the integral value is 2.47 × 30 ≈ 74.1 με. Compared with the cumulative circumferential component of 78 με, the deviation rate is |74.1 - 78| / 78 ≈ 3.9 / 78 ≈ 0.05. The optimized circumferential weight is 0.26 × (1 - 0.05) = 0.247. Based on the optimized weight, the optimized circumferential linear strain at the crown is calculated as 300 × 0.247 ≈ 74.1 με, and the optimized radial linear strain at the right arch waist is calculated as 250 × (-0.8944) ≈ -223.6 με, thus completing the strain optimization.

[0090] Step S24, which involves using the FBG demodulation system to analyze the fiber optic sensing data to obtain the center wavelength offset of each sensing node and separating the pure strain component based on the FBG strain-temperature cross-sensitivity compensation model, includes the following steps:

[0091] Step S2401: Extract the original spectral signal of the FBG sensing node from the fiber optic sensing basic data, and calculate the center wavelength offset by identifying the initial value of the center wavelength and combining it with the reference wavelength of the temperature reference point; extract the wavelength drift rate = the difference in offset between adjacent times / time interval based on the time series features of the center wavelength offset, and screen out the effective sensing nodes with stable signals by the wavelength drift rate.

[0092] In this embodiment of the invention, the original spectral signals (wavelength range 1540-1560nm) of 200 FBG sensing nodes in the K1+000-K1+500 segment are obtained from the fiber optic sensing basic data using a spectral signal extraction tool. The initial value of the center wavelength of each node is identified using a spectral analysis tool, such as the initial wavelength of 1550.00nm for the right arch waist node at K1+100 and 1550.50nm for the arch top node at K1+120. Three nodes in the temperature-stable zone at the tunnel entrance (annual temperature difference ≤5℃) are selected as temperature reference points, with reference wavelengths of 1549.80nm, 1549.90nm, and 1550.10nm, respectively. The center wavelength offset is calculated using the formula: Node measured wavelength - Average reference wavelength (1549.93nm). For example, the measured wavelength of node K1+100 is 1550.12nm, and the offset is 1550.12 - 1549.93 = 0.19nm. The offset time series (720 data points in total) was extracted at 1-hour intervals. The wavelength drift rate was calculated using the formula: wavelength drift rate = |offset at the next moment - offset at the previous moment| / 1 hour. If the drift rate of a node for 3 consecutive hours is 0.002nm / h, 0.003nm / h, and 0.002nm / h respectively, all of which are less than the threshold of 0.005nm / h, it is determined to be a valid sensing node with stable signal. A total of 185 valid nodes were selected.

[0093] Step S2402: Construct an FBG strain-temperature cross-sensitivity compensation model. The model input parameters include the thermal conductivity of the rock mass at the node layout location and the thermal expansion coefficient of the optical cable encapsulation material. Calculate the temperature sensitivity correction coefficient and dynamically calibrate the temperature influence term in the model based on the temperature sensitivity correction coefficient to generate a calibrated cross-sensitivity model.

[0094] In this embodiment of the invention, an FBG strain-temperature cross-sensitivity compensation model is constructed. The basic formula of the model is Δλ=Kε·e+Kt·T (Δλ is the center wavelength offset, Kε is the strain sensitivity coefficient 1.2pm / με, Kt is the temperature sensitivity coefficient 10pm / ℃, e is the strain, and T is the temperature). The thermal conductivity of the rock mass at the node layout location is retrieved from the tunnel construction BIM model rock mass parameter library (the thermal conductivity of weathered sandstone in the K1+100 region is 1.5W / (m·K)), and the thermal expansion coefficient of the encapsulation material (stainless steel pipe) is obtained from the optical cable technical parameter library (12×10). -6 / ℃. Using a correction factor calculation tool, the temperature sensitivity correction factor is calculated according to the formula: Temperature Sensitivity Correction Factor = 1 - (Rock Mass Thermal Conductivity × Coefficient of Thermal Expansion) / 10 -5 The calculation yields the correction factor as 1 - (1.5 × 12 × 10). -6 ) / 10 -5=1-0.18=0.82. Substitute the correction coefficient into the model to dynamically calibrate the temperature influence term, generating a calibrated cross-sensitivity model: Δλ=1.2e+(10×0.82)T=1.2e+8.2T, ensuring that the model can eliminate the interference of rock mass thermal conductivity and encapsulation material expansion on temperature sensitivity.

[0095] Step S2403: Substitute the center wavelength offset into the calibrated cross-sensitive model to decompose it into a mixed component that includes the coupling effect of strain and temperature; introduce synchronously acquired ambient temperature monitoring data, calculate the proportion of temperature influence, and separate the preliminary pure strain component based on the proportion threshold.

[0096] In this embodiment of the invention, by substituting the center wavelength offset of 0.19 nm (190 pm) of the effective node of the right arch waist of K1+100 into the calibrated cross-sensitivity model, we obtain 190 = 1.2e + 8.2T, forming a mixed component equation coupling strain and temperature. The ambient temperature data (25℃) synchronously collected by the node is obtained through an ambient temperature monitoring device. Substituting this data into the equation, the temperature influence term is calculated as 8.2 × 25 = 205 pm. According to the formula, the proportion of temperature influence = temperature influence term / Δλ × 100% = 205 / 190 × 100% ≈ 107.9% (due to measurement error, the proportion is allowed to fluctuate by ±10%). Setting the proportion threshold to 80%-120%, the temperature influence of the node is determined to be effective. According to the formula, the preliminary pure strain component e = (Δλ - temperature influence term) / Kε = (190 - 205) / 1.2 ≈ -12.5 με, the preliminary pure strain component is obtained. The preliminary pure strain components of 185 effective nodes are calculated.

[0097] Step S2404: Perform time-dependent correction on the initial pure strain components to determine the time-dependent reference point for strain monitoring by analyzing the support construction time nodes in the construction time sequence data; calculate the strain time-dependent attenuation coefficient, and perform time-dependent correction on the initial pure strain components based on the strain time-dependent attenuation coefficient to obtain the corrected pure strain components.

[0098] In this embodiment of the invention, the support construction time node is extracted from the construction time sequence data. For example, the initial support completion time of section K1+100 is the 30th day after the tunnel construction starts. This time point is set as the time reference point for strain monitoring (t=0). The strain time aging attenuation coefficient is calculated using the formula e^(-0.01t) (t is the monitoring duration in days). For example, on the 10th day of monitoring (t=10), the attenuation coefficient is e^(-0.1)=0.905; on the 30th day of monitoring (t=30), the attenuation coefficient is e^(-0.3)=0.741. The initial pure strain component -12.5με at node K1+100 was corrected for time dimension. The corrected pure strain component was calculated according to the formula = initial pure strain component × attenuation coefficient. After monitoring for 10 days, the corrected strain = -12.5 × 0.905 ≈ -11.3με; after monitoring for 30 days, the corrected strain = -12.5 × 0.741 ≈ -9.3με, thus eliminating the aging error caused by strain attenuation over time.

[0099] Step S2405: Construct a point strain parameter validity verification mechanism, perform spatial correlation analysis between the corrected pure strain component and the linear strain distribution parameter monitored by BOTDA in the same area, calculate the strain fit degree = 1 - absolute difference between the two / amplitude of the linear strain parameter, and perform weighted optimization on the corrected pure strain component based on the strain fit degree to generate point strain parameters that include strain amplitude, variation trend and spatial correlation.

[0100] In this embodiment of the invention, by constructing a point strain parameter validity verification mechanism, the linear strain parameter -10.8με from the BOTDA monitoring of the same area as the right arch waist of K1+100 is extracted from the previously generated surrounding rock strain distribution parameters. The correlation analysis tool is used to calculate the strain fit using the formula: Strain Fit = 1 - |Corrected Pure Strain Component - Linear Strain Parameter| / |Linear Strain Parameter|. For example, if the corrected strain on the 10th day of monitoring is -11.3με, the fit = 1 - |-11.3 - (-10.8)| / 10.8 = 1 - 0.5 / 10.8 ≈ 0.954 (95.4%). A fit weighting coefficient is set to fit (≥0.8 takes the actual value, <0.8 takes 0.8), and the optimized pure strain component is calculated using the formula: Optimized Pure Strain Component = Corrected Pure Strain Component × Weighting Coefficient, resulting in the optimized strain = -11.3 × 0.954 ≈ -10.8με. By integrating the optimized strain amplitude (-10.8με), variation trend (gradually increasing with monitoring time), and spatial correlation (95.4% agreement with BOTDA line strain) of 185 effective nodes, a point strain parameter table for key points of the surrounding rock is generated.

[0101] Step S3 includes the following steps:

[0102] Step S31: Divide the multidimensional feature vector of surrounding rock deformation according to different monitoring sections of the tunnel. Each monitoring section corresponds to a feature vector sample. Extract the linear strain distribution parameters and point strain parameters in the feature vector as the core feature dimensions, and the construction correlation law and geological coupling relationship as auxiliary feature dimensions.

[0103] In this embodiment of the invention, the 50 multidimensional feature vectors of surrounding rock deformation in the K1+000-K1+500 segment are split according to the monitoring cross sections by using a cross section segmentation tool. Each cross section (one cross section every 10m, a total of 50) corresponds to one feature vector sample. For example, the sample vector of the K1+100 cross section is [32,38,28,31.5,0.8,30,1.2,0.75,30,1e-5,180°,0.35]. Using a feature dimension extraction tool, core feature dimensions and auxiliary feature dimensions are separated from each sample vector: the core feature dimensions include linear strain distribution parameters (first 3 elements: average 32με, maximum 38με, minimum 28με) and point strain parameters (4th-5th elements: point strain -11.3με, rate of change 0.8με / day), for a total of 5 dimensions; the auxiliary feature dimensions include construction correlation patterns (6th-7th elements: support completion time 30 days, advancement rate 1.2m / day) and geological coupling relationships (8th-10th elements: integrity index 0.75, shear strength 30MPa, permeability coefficient 1e-5cm / s), for a total of 5 dimensions. The remaining circumferential angle and radial burial depth are retained as positioning parameters, thus completing the division of all cross-sectional feature dimensions.

[0104] Step S32: Using the density peak clustering algorithm, calculate the Euclidean distance between the feature vector samples of each cross section to determine the density core points in the feature space; based on the density value and distance value corresponding to the density core points, classify the deformation mode categories corresponding to uniform deformation, local mutation, and gradual accumulation.

[0105] In this embodiment of the invention, density peak clustering is used to perform cluster analysis on 50 cross-sectional feature vector samples. First, the Euclidean distance between each sample is calculated. Taking the K1+100 cross-sectional sample (core features [32,38,28,31.5,0.8]) and the K1+110 cross-sectional sample (core features [30,36,26,29.8,0.7]) as an example...

[0106] Euclidean distance = √[(32-30)] 2 +(38-36) 2 +(28-26) 2 +(31.5-29.8) 2 +(0.8-0.7) 2 ]

[0107] ≈√(4+4+4+2.89+0.01)=√14.9≈3.86. The local density calculation radius is set to 4. The number of neighboring samples within this radius for each sample is counted as the density value. For example, the density value of the sample in section K1+100 is 8 (8 neighboring samples are less than 4 apart), and the density value of the sample in section K1+120 (core feature [55,62,48,52.3,2.1]) is 3. Deformation patterns are classified according to the density value (threshold 5) and distance value (threshold 5): density value ≥ 5 and distance value < 5 indicate uniform deformation (e.g., sections K1+100 and K1+110); density value < 5 and distance value ≥ 5 indicate local mutation (e.g., section K1+120); density value ≥ 5 and distance value ≥ 5 indicate progressive accumulation (e.g., section K1+200, density value 6, distance value 6.2). This completes the deformation pattern classification for 50 sections.

[0108] Step S33: For different deformation mode categories, extract the excavation progress and support strength process progress indicators from the construction time sequence data and assign them corresponding process progress weights; extract the surrounding rock integrity coefficient and uniaxial compressive strength index from the geological attribute data and calculate the lithological influence coefficient.

[0109] In this embodiment of the invention, for a uniformly deformable cross-section (such as K1+100), the excavation advance of 1.2m / day and the support strength (25cm shotcrete thickness and 1m anchor spacing) are extracted from the construction sequence data. The weight of the process progress is calculated according to the formula: (excavation advance / 1.5m / day)×0.4+(support strength score / 100)×0.6 (1.5m / day is the standard advance, and the support strength score is quantified according to thickness and spacing, 25cm thickness gets 80 points, 1m spacing gets 90 points, and the average is 85 points). The weight is (1.2 / 1.5)×0.4+(85 / 100)×0.6=0.32+0.51=0.83. For sections with localized abrupt changes (such as K1+120), the excavation advance is 0.8m / day, the support strength score is 70 points, and the weight is (0.8 / 1.5)×0.4+(70 / 100)×0.6≈0.213+0.42=0.633. The integrity coefficient of the surrounding rock (0.75 for K1+100, 0.5 for K1+120) and the uniaxial compressive strength (30MPa for K1+100, 20MPa for K1+120) of each section were extracted from the geological attribute data. The lithological influence coefficient was calculated according to the formula: lithology influence coefficient = (integrity coefficient × 0.5) + (uniaxial compressive strength / 50MPa × 0.5) (50MPa is the standard strength). The coefficient of K1+100 = 0.75 × 0.5 + (30 / 50) × 0.5 = 0.375 + 0.3 = 0.675, and the coefficient of K1+120 = 0.5 × 0.5 + (20 / 50) × 0.5 = 0.25 + 0.2 = 0.45.

[0110] Step S34: Substitute the process progress weight and lithological influence coefficient into the multidimensional feature vector of surrounding rock deformation, and perform weighted optimization on the core feature dimension and auxiliary feature dimension; construct a deformation prediction model based on an improved decision tree, take the optimized feature vector as input, take the actual deformation of the surrounding rock monitored in history as output, optimize the decision tree nodes through a pruning algorithm, train the model and input the real-time feature vector to generate the dynamic prediction value of surrounding rock deformation corresponding to each monitoring section.

[0111] In this embodiment of the invention, a weighted optimization tool is used to substitute the process progress weight and lithological influence coefficient into the feature vector. The core feature dimension is calculated using the formula: optimized value = original value × (process progress weight + lithological influence coefficient) / 2. For example, the average core feature value of section K1+100 is 32με, and the optimized value is 32 × (0.83 + 0.675) / 2 ≈ 32 × 0.7525 ≈ 24.08με. The auxiliary feature dimension is calculated using the formula: optimized value = original value × lithological influence coefficient. For example, if the support completion time is 30 days, the optimized value is 30 × 0.675 = 20.25 days. A deformation prediction model based on an improved decision tree was constructed. The model used optimized feature vectors (10 dimensions) from 50 cross-sections as input and actual surrounding rock deformation data from the past three months (e.g., cumulative deformation of 12 mm at cross-section K1+100) as output. A pruning algorithm was used to optimize the decision tree nodes, setting a minimum number of sample splits of 5 and a minimum number of samples per leaf node of 3, and removing branches with an error rate >10%. The optimized feature vector from cross-section K1+105, collected in real-time, was input into the trained model to generate dynamic predictions of surrounding rock deformation for the next 15 days, such as a predicted deformation of 1.2 mm on day 5, 2.1 mm on day 10, and 2.8 mm on day 15, thus completing the deformation prediction.

[0112] Step S4 includes the following steps:

[0113] Step S41: Obtain the displacement data of the surrounding rock section monitored by the total station at the tunnel site and the deep displacement data of the surrounding rock monitored by the multi-point displacement gauge. Combine the two to obtain the actual deformation of the surrounding rock at each monitoring section. Match the dynamic prediction value of the surrounding rock deformation with the actual deformation of the surrounding rock according to the section. Calculate the absolute error, relative error and root mean square error of each monitoring section to form the surrounding rock deformation prediction error matrix.

[0114] In this embodiment of the invention, displacement data acquisition tools are used to obtain total station monitoring data at the tunnel site, such as the arch crown displacement of 11.8mm, right arch waist displacement of 9.5mm, and left arch waist displacement of 9.2mm at section K1+100, and the average value of 10.2mm is taken as the section displacement data; the deep displacement data of the surrounding rock monitored by multi-point displacement gauges are also obtained, with the displacement at a depth of 2m of this section being 8.5mm and at a depth of 3m being 7.3mm, and the displacement at a depth of 1.5m (near the radial burial depth of the optical cable) being 7.9mm as the deep displacement data. Using data fusion tools and by setting weight coefficients for cross-sectional displacement and deep displacement, the accuracy of the two types of monitoring data was evaluated. Total station monitoring of key cross-sectional points (arch crown, arch waist) provides strong data coverage but is easily affected by surface interference; multi-point displacement gauges monitoring deep surrounding rock provide data closer to internal deformation but are greatly affected by burial depth. The error levels (e.g., standard deviation of repeated measurements) of the two types of data were statistically analyzed, with higher-accuracy data assigned higher weights. Considering the tunnel surrounding rock grade, in hard rock, deep displacement stability is strong, and its weight can be appropriately increased; in soft rock, the apparent deformation of the cross-section better reflects risk, and the weight of cross-sectional displacement is higher. Furthermore, by searching for weight application experience in similar tunnels (same surrounding rock, cross-sectional form, construction method), a preliminary weight range of 0.5-0.7 and 0.3-0.5 was determined. A trial-and-error method was used to select multiple weight combinations within these ranges, and the error between the fused deformation and the actual deformation on site (e.g., data from pre-embedded joint gauges) was calculated, as shown in Table 2 below.

[0115] Table 2 Example of Root Mean Square Error (RMSE) Test

[0116]

[0117] With the goal of minimizing the root mean square error (RMSE), the weights that best approximate the actual deformation were determined. Table 2 shows that when w=0.6, the RMSE is 0.79mm (minimum), and the corresponding weight combination (0.6, 0.4) is the optimal solution. Calculated using the formula: Actual deformation of surrounding rock = (section displacement data × 0.6 + deep displacement data × 0.4), the actual deformation of section K1+100 is (10.2 × 0.6 + 7.9 × 0.4) = 6.12 + 3.16 = 9.28mm. The generated dynamic deformation prediction values ​​for each section (e.g., 15-day prediction of 12mm for section K1+100, 10.5mm for section K1+110) are matched one-to-one with the actual deformation. Error calculation tools are used, with the formulas: Absolute error = |predicted value - actual value|, Relative error = Absolute error / actual value × 100%, and Root mean square error = √[(predicted value - actual value)]. 2[Sampling count] Calculation shows that the absolute error of section K1+100 is |12-9.28|=2.72mm, and the relative error is 2.72 / 9.28×100%≈29.3%. The predicted value of section K1+110 is 10.5mm, and the actual value is 8.8mm, with an absolute error of 1.7mm and a relative error of 19.3%. Integrating the three types of errors from 50 sections, a 50-row, 3-column surrounding rock deformation prediction error matrix is ​​formed.

[0118] Step S42: Set an error threshold. If the prediction error matrix of the surrounding rock deformation of a certain monitoring section is greater than the error threshold, it is determined that the prediction result of the section exceeds the allowable error range. Based on the surrounding rock deformation prediction error matrix, the gradient descent algorithm is used to iteratively correct the process progress weight and lithological influence coefficient in the deformation prediction model. The partial derivatives of each error with respect to the process progress weight and lithological influence coefficient are calculated. The parameter values ​​are adjusted along the direction of error reduction until the surrounding rock deformation prediction error matrix of all sections is less than or equal to the error threshold, and the corrected accurate prediction value of the surrounding rock deformation is obtained.

[0119] In this embodiment of the invention, by setting error thresholds (absolute error ≤ 2mm, relative error ≤ 20%, root mean square error ≤ 1.5mm), and comparing the error matrix, it was found that the absolute error of section K1+100 was 2.72mm > 2mm, and the relative error was 29.3% > 20%, which was determined to exceed the allowable error range. The gradient descent algorithm was used to iteratively correct the process progress weights and lithological influence coefficients in the deformation prediction model, with a learning rate of 0.01 and 50 iterations. The partial derivatives of each error with respect to the parameters were calculated. For example, the partial derivative of the relative error of section K1+100 with respect to the process progress weight (originally 0.83) = (error change / weight change) ≈ 0.3, and the weight was adjusted along the direction of error reduction = 0.83 - 0.01 × 0.3 = 0.827; the partial derivative with respect to the lithological influence coefficient (originally 0.675) ≈ 0.2, and the adjustment coefficient = 0.675 - 0.01 × 0.2 = 0.673. Substituting the adjusted parameters into the model, the predicted values ​​were recalculated. The new predicted value for section K1+100 was 10.5 mm, with an absolute error of |10.5 - 9.28| = 1.22 mm ≤ 2 mm and a relative error of 1.22 / 9.28 × 100% ≈ 13.1% ≤ 20%. This iteration was repeated until all section errors met the threshold requirements, ultimately yielding a corrected accurate predicted value of 9.3 mm for section K1+100 and 8.5 mm for section K1+120 (where the original error exceeded the threshold).

[0120] Step S43: According to the tunnel mileage section and the surrounding rock deformation mode category, the corrected accurate prediction value of surrounding rock deformation is regionally aggregated, and the average deformation amount and deformation ratio of each type of deformation mode in each mileage section are calculated.

[0121] In this embodiment of the invention, the tunnel is divided into 10 regions (e.g., K1+000-K1+050, K1+050-K1+100, ..., K1+450-K1+500) every 50m of tunnel mileage. Each region contains 5 monitoring sections. A region aggregation tool is used to statistically analyze the deformation data of each region according to the deformation mode type (uniform deformation, local abrupt change, and gradual accumulation). For example, the K1+050-K1+100 region contains 3 uniform deformation sections (actual deformation amounts of 9.28mm, 8.8mm, and 9.5mm) and 2 gradual accumulation sections (11.2mm and 10.8mm). Calculate the average deformation for each deformation mode: Average for uniform deformation = (9.28 + 8.8 + 9.5) / 3 ≈ 9.19 mm; Average for progressive cumulative deformation = (11.2 + 10.8) / 2 = 11 mm. Calculate the deformation percentage: Percentage for uniform deformation = 3 / 5 × 100% = 60%; Percentage for progressive cumulative deformation = 40%. Similarly, perform aggregate calculations for 10 regions to generate a regionalized deformation statistics table containing region range, deformation mode, average deformation, and deformation percentage.

[0122] Step S44: Using 3D geographic information system technology, spatially correlate the tunnel axis, monitoring section location, and aggregated deformation data to generate a 3D model of the spatiotemporal distribution of surrounding rock deformation; mark areas where the deformation exceeds the warning value in red to generate a visualization result of the prediction of surrounding rock deformation in the mountain tunnel.

[0123] In this embodiment of the invention, a three-dimensional geographic information system (GIS) tool is used to import the three-dimensional coordinates of the tunnel axis (X-axis from 0-500m, Y-axis from 0-10m, Z-axis from 0-8m for the K1+000-K1+500 segment). The locations of each monitoring section (e.g., K1+100 section X=100m, Y=5m, Z=4m) are spatially correlated with the aggregated deformation data (average deformation 9.19mm, accounting for 60%). A three-dimensional model of the spatiotemporal distribution of surrounding rock deformation is constructed. The tunnel outline is presented in the model at a 1:1 scale, and different colors are used to distinguish deformation modes in different areas (uniform deformation blue, local abrupt deformation red, and gradual cumulative deformation yellow). The deformation is visually displayed through the height of the protrusion on the model surface (1mm deformation corresponds to a 0.1m protrusion). A deformation warning value of 10mm is set, and the K1+150-K1+200 area (average deformation 10.5mm > 10mm) is marked in red, with the marked area flashing in the model as a reminder. The exported model is in an interactive format, supporting zooming and rotation for viewing, and generates a visualization report of the predicted deformation results of the surrounding rock of the mountain tunnel, which includes 3D model screenshots and regional deformation statistics tables.

[0124] This invention also provides a mountain tunnel surrounding rock deformation monitoring system based on distributed optical fiber sensing, such as... Figure 4As shown, a mountain tunnel surrounding rock deformation monitoring system based on distributed optical fiber sensing is used to perform the above-described method. The system comprises:

[0125] The surrounding rock area data acquisition module is used to collect raw Brillouin frequency shift data along the sensing optical cable and raw center wavelength data of each sensing node in the monitoring area of ​​the surrounding rock of the mountain tunnel. At the same time, it acquires tunnel construction progress data and surrounding rock geological survey data. The module performs outlier removal and noise filtering on the raw Brillouin frequency shift data and raw center wavelength data, and performs format standardization and time sequence alignment on the tunnel construction progress data and surrounding rock geological survey data to generate preprocessed fiber optic sensing basic data, construction time sequence data and geological attribute data.

[0126] The surrounding rock deformation feature analysis module is used to extract surrounding rock deformation feature parameters based on fiber optic sensing data, including the axial strain values ​​of each monitoring point of the sensing cable and the strain distribution parameters along the surrounding rock; it analyzes the center wavelength offset of each sensing node and separates the pure strain component to generate the point strain parameters of key points of the surrounding rock; at the same time, it constructs a multi-dimensional feature vector of surrounding rock deformation based on the surrounding rock deformation feature parameters combined with construction time series data and geological attribute data.

[0127] The surrounding rock deformation prediction module is used to perform deformation pattern density clustering on the surrounding rock of different monitoring sections of the tunnel based on the multi-dimensional feature vector of surrounding rock deformation, so as to classify the deformation pattern categories corresponding to uniform deformation, local abrupt deformation, and gradual accumulation. For different deformation pattern categories, the module introduces the process progress weight in the construction time sequence data and the lithological influence coefficient in the geological attribute data to dynamically predict the surrounding rock deformation, and generate the dynamic prediction value of surrounding rock deformation for each monitoring section.

[0128] The surrounding rock deformation visualization module is used to compare the dynamic predicted values ​​of surrounding rock deformation with the monitoring data of total station and multi-point displacement gauges deployed on-site in the tunnel, calculate the absolute error, relative error and root mean square error between the predicted and measured values, and form a surrounding rock deformation prediction error matrix. Based on the surrounding rock deformation prediction error matrix, iterative correction is performed, and regional aggregation is performed according to the tunnel mileage segment and deformation mode category, thereby generating a visualization result of the spatiotemporal distribution of surrounding rock deformation in the mountain tunnel.

[0129] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing, characterized in that, Includes the following steps: Step S1: Collect raw Brillouin frequency shift data along the sensing optical cable and raw center wavelength data of each sensing node in the monitoring area of ​​the mountain tunnel surrounding rock. At the same time, acquire tunnel construction progress data and surrounding rock geological survey data. Perform outlier removal and noise filtering on the raw Brillouin frequency shift data and raw center wavelength data. Perform format standardization and time sequence alignment on the tunnel construction progress data and surrounding rock geological survey data to generate preprocessed fiber optic sensing basic data, construction time sequence data and geological attribute data. Step S2: Extract surrounding rock deformation characteristic parameters based on fiber optic sensing data, including axial strain values ​​at each monitoring point of the sensing cable and strain distribution parameters along the surrounding rock; analyze the center wavelength offset of each sensing node and separate the pure strain component to generate point strain parameters for key points of the surrounding rock; simultaneously, construct a multi-dimensional feature vector of surrounding rock deformation based on the surrounding rock deformation characteristic parameters combined with construction time series data and geological attribute data; Step S2 includes the following steps: Step S21: By retrieving the Brillouin frequency shift sequence data output by the BOTDA demodulation system from the preprocessed fiber optic sensing basic data, and combining the Brillouin frequency shift-strain sensitivity coefficient and temperature compensation coefficient calibrated by the BOTDA demodulation system, the axial strain value of each monitoring point of the sensing optical cable is calculated. Step S22: Extract the strain gradient change rate based on the spatial distribution characteristics corresponding to the axial strain value, specifically, the strain gradient change rate = strain difference between adjacent monitoring points / optical cable laying spacing, and identify the strain abrupt change segment and smooth segment corresponding to each monitoring point along the optical cable through the strain gradient change rate; Step S23: Based on the strain abrupt change section and smooth section corresponding to each monitoring point along the optical cable, the three-dimensional layout trajectory of the optical cable corresponding to the sensing optical cable is obtained by combining the tunnel construction BIM model. Based on the three-dimensional layout trajectory of the optical cable, the axial strain value of each monitoring point of the sensing optical cable is statistically analyzed to generate the strain distribution parameters along the surrounding rock. Step S24: The center wavelength offset of each sensing node is obtained by parsing the fiber optic sensing basic data through the FBG demodulation system, and the pure strain component is separated according to the FBG strain-temperature cross-sensitivity compensation model to generate the point strain parameters of the key points of the surrounding rock; the strain distribution parameters along the surrounding rock are combined with the point strain parameters of the key points of the surrounding rock to generate the deformation characteristic parameters of the surrounding rock. Step S25: Based on the surrounding rock deformation characteristic parameters, combined with construction time series data and geological attribute data, construct a multi-dimensional feature vector of surrounding rock deformation that includes linear strain distribution, point strain parameters, construction correlation laws and geological coupling relationships; Step S3: Based on the multidimensional feature vector of surrounding rock deformation, perform deformation pattern density clustering on the surrounding rock of different monitoring sections of the tunnel to classify the deformation pattern categories corresponding to uniform deformation, local abrupt deformation, and gradual accumulation. For different deformation pattern categories, introduce the process progress weight in the construction time sequence data and the lithological influence coefficient in the geological attribute data to dynamically predict the surrounding rock deformation and generate the dynamic prediction value of surrounding rock deformation for each monitoring section. Step S4: Compare the dynamic prediction values ​​of surrounding rock deformation with the monitoring data from the total station and multi-point displacement gauges deployed at the tunnel site, calculate the absolute error, relative error, and root mean square error between the predicted and measured values, and form a surrounding rock deformation prediction error matrix; perform iterative correction based on the surrounding rock deformation prediction error matrix, and aggregate regionally according to the tunnel mileage segment and deformation mode category to generate a visualization result of the spatiotemporal distribution of surrounding rock deformation in the mountain tunnel.

2. The method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing according to claim 1, characterized in that, Step S1 includes the following steps: Step S11: In the monitoring area of ​​the surrounding rock of the mountain tunnel, a multi-dimensional sensing network is constructed based on the orientation of the joints and the distribution of stress concentration areas in the surrounding rock as determined by geological survey. The raw Brillouin frequency shift data along the optical cable is collected using the BOTDA demodulation system, and the amplitude attenuation rate and spectral bandwidth of the frequency shift signal are recorded simultaneously. The signal quality index is calculated based on the amplitude attenuation rate and spectral bandwidth. The raw center wavelength data of each sensing node is collected using the FBG demodulation system, and the reflected light intensity and polarization state are obtained simultaneously. The node sensing reliability coefficient is calculated based on the reflected light intensity and polarization state. Step S12: Extract construction process transition time nodes, face advancement characteristic parameters, and support structure layout parameters through the tunnel BIM system to form tunnel construction progress data; obtain surrounding rock integrity index, rock mass structural shear strength, and groundwater permeability coefficient through ground-penetrating radar detection and borehole sampling analysis to form surrounding rock geological exploration data. Step S13: Process the raw Brillouin frequency shift data to filter and remove abnormal data segments with a signal quality index below a threshold based on the signal quality index, and suppress noise in the effective data segments to obtain denoised Brillouin frequency shift data; calculate the initial distribution parameters of strain along the path based on the denoised Brillouin frequency shift data. Step S14: Process the original center wavelength data to remove failed node data with node sensing reliability coefficients below a threshold based on the node sensing reliability coefficient, and perform noise separation on the effective node data to obtain the denoised center wavelength data; calculate the initial value of node strain based on the denoised center wavelength data. Step S15: Standardize the tunnel construction progress data to convert the construction process transition time nodes into cumulative durations relative to the start of tunnel construction, and convert the face advancement characteristic parameters into advancement rate curves associated with mileage stations to generate standardized construction time sequence data; standardize the surrounding rock geological survey data to convert it into dimensionless geological feature vectors; establish a spatiotemporal reference axis based on the tunnel construction mileage stations, and spatiotemporally align the denoised fiber optic sensing data, which includes the initial distribution parameters of strain along the tunnel and the initial values ​​of nodal strain, the standardized construction time sequence data, and the geological feature vectors, to generate preprocessed fiber optic sensing basic data, construction time sequence data, and geological attribute data.

3. The method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing according to claim 1, characterized in that, Step S23, which involves statistically analyzing the linear strain distribution of the axial strain values ​​at each monitoring point of the sensing optical cable based on the three-dimensional deployment trajectory of the optical cable, includes the following steps: Step S2301: Based on the three-dimensional deployment trajectory of the optical cable and combined with the tunnel construction BIM model, obtain the circumferential curvature radius and radial burial depth parameters of the sensing optical cable at each tunnel section. Step S2302: Calculate the optical cable bending correction coefficient and the surrounding rock constraint coefficient based on the circumferential curvature radius and radial burial depth parameters, where the optical cable bending correction coefficient = 1 - circumferential reference value / circumferential curvature radius, and the surrounding rock constraint coefficient = radial burial depth parameter × rock mass elastic modulus; Step S2303: Couple the axial strain value of each monitoring point of the sensing optical cable with the optical cable bending correction coefficient and the surrounding rock constraint coefficient to obtain the corrected axial strain value; based on the corrected axial strain value, extract the strain accumulation and strain change rate parameters corresponding to each monitoring point. Step S2304: Calculate the radial transformation weight using the spatial angle parameters of the optical cable layout, and calculate the circumferential transformation weight using the Poisson's ratio parameters of the rock mass. Simultaneously, based on the radial transformation weight and the circumferential transformation weight, convert the strain accumulation and strain change rate parameters corresponding to each monitoring point into the corresponding radial linear strain and circumferential linear strain. Step S2305: Perform spatial interpolation on the radial and circumferential strains to generate a tunnel strain distribution cloud map; extract the spatial coordinates of strain extreme points, the density of strain contour lines, and the main vector parameters of strain direction from the tunnel strain distribution cloud map, and integrate them to form the surrounding rock strain distribution parameters that include strain amplitude, distribution pattern, and spatial gradient.

4. The method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing according to claim 3, characterized in that, Step S2301 includes the following steps: The three-dimensional coordinate dataset corresponding to the three-dimensional cable laying trajectory is retrieved from the tunnel construction BIM model. This dataset includes the mileage station, cross-sectional circumferential angle, and initial radial depth of the cable along the tunnel axis. Based on the three-dimensional coordinate dataset, the coordinates of the spatial discrete points of the cable in each tunnel cross section are extracted, and the arc length distance and angle deviation between adjacent discrete points are calculated to generate the spatial topology parameters corresponding to the cable laying trajectory. A continuous curve model of the optical cable is constructed based on spatial topological parameters. The coordinates of discrete points in space are then fitted using the least squares method to obtain the circumferential layout curve equation of the optical cable in each tunnel section. The curvature value of each point is calculated based on the circumferential layout curve equation, where the curvature value is calculated as: second derivative of the curve / (1 + square of the first derivative). 3 / 2 The initial value of the circumferential radius of curvature of the optical cable at the corresponding position is obtained by back-calculation based on the curvature value; The tunnel cross-section design contour data, including the excavation contour line, the initial contour line, and the thickness parameters of the support structure, are retrieved from the tunnel construction BIM model. The normal distance between the three-dimensional laying trajectory of the optical cable and the excavation contour line and the initial contour line is calculated. The difference between the normal distance and the thickness parameters of the support structure is combined to obtain the initial value of the radial burial depth of the optical cable in the surrounding rock. Construction deviation parameters during the optical cable laying process are introduced. These parameters are obtained through comparative analysis of on-site optical cable laying records and tunnel construction BIM models, including circumferential angle deviation and radial depth deviation. The circumferential angle deviation is substituted into the initial value of the circumferential curvature radius for correction, resulting in the corrected circumferential curvature radius, specifically: corrected circumferential curvature radius = initial value of circumferential curvature radius × (1 + circumferential angle deviation coefficient). By combining the convergence deformation after the surrounding rock excavation, the initial value of the radial burial depth is dynamically adjusted so as to calculate the corrected radial burial depth parameter through the coupling relationship between the convergence deformation and the radial burial depth; and the circumferential curvature radius and radial burial depth parameters of the sensing optical cable at each tunnel section are integrated and generated.

5. The method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing according to claim 3, characterized in that, Step S2304 includes the following steps: The spatial azimuth parameter θ of each monitoring point is extracted from the three-dimensional laying trajectory of the optical cable, and combined with the rock mechanics parameter library in the tunnel construction BIM model, the surrounding rock Poisson's ratio μ and elastic modulus E of the corresponding monitoring point are retrieved; the radial foundation weight cosθ is calculated based on the spatial azimuth parameter θ, and the circumferential foundation weight μ·sinθ is calculated based on the surrounding rock Poisson's ratio μ. Obtain the measured tension and design tension during optical cable laying, and calculate the tension influence coefficient = measured tension / design tension; multiply the radial foundation weight and circumferential foundation weight by the tension influence coefficient to obtain the corrected radial conversion weight and circumferential conversion weight. The total axial strain increment at each monitoring point is extracted from the strain accumulation, and the radial strain accumulation component is calculated by combining the corrected radial transformation weight. At the same time, the circumferential strain accumulation component is calculated by combining the circumferential transformation weight. Dynamic correction is performed based on the strain change rate parameter and a time-dimensional rock mass creep correction factor to generate the corrected rate parameter. The corrected rate parameter is then coupled with radial and circumferential transformation weights to obtain the radial strain change rate and the circumferential strain change rate. The radial strain cumulative component and radial strain change rate are verified by time integration, and the circumferential strain cumulative component and circumferential strain change rate are also verified by time integration. The radial and circumferential transformation weights are then optimized based on the integral result deviation rate, where the integral result deviation rate = |integral value - cumulative component| / cumulative component, and the optimized weight = original weight × (1 - integral result deviation rate). The optimized radial and circumferential linear strains are then output.

6. The method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing according to claim 1, characterized in that, Step S24, which involves using the FBG demodulation system to analyze the fiber optic sensing data to obtain the center wavelength offset of each sensing node and separating the pure strain component based on the FBG strain-temperature cross-sensitivity compensation model, includes the following steps: Step S2401: Extract the original spectral signal of the FBG sensing node from the fiber optic sensing basic data, and calculate the center wavelength offset by identifying the initial value of the center wavelength and combining it with the reference wavelength of the temperature reference point; extract the wavelength drift rate = the difference in offset between adjacent times / time interval based on the time series features of the center wavelength offset, and screen out the effective sensing nodes with stable signals by the wavelength drift rate. Step S2402: Construct an FBG strain-temperature cross-sensitivity compensation model. The model input parameters include the thermal conductivity of the rock mass at the node layout location and the thermal expansion coefficient of the optical cable encapsulation material. Calculate the temperature sensitivity correction coefficient and dynamically calibrate the temperature influence term in the model based on the temperature sensitivity correction coefficient to generate a calibrated cross-sensitivity model. Step S2403: Substitute the center wavelength offset into the calibrated cross-sensitive model to decompose it into a mixed component that includes the coupling effect of strain and temperature; introduce synchronously acquired ambient temperature monitoring data, calculate the proportion of temperature influence, and separate the preliminary pure strain component based on the proportion threshold. Step S2404: Perform time-dependent correction on the initial pure strain components to determine the time-dependent reference point for strain monitoring by analyzing the support construction time nodes in the construction time sequence data; calculate the strain time-dependent attenuation coefficient, and perform time-dependent correction on the initial pure strain components based on the strain time-dependent attenuation coefficient to obtain the corrected pure strain components. Step S2405: Construct a point strain parameter validity verification mechanism, perform spatial correlation analysis between the corrected pure strain component and the linear strain distribution parameter monitored by BOTDA in the same area, calculate the strain fit degree = 1 - absolute difference between the two / amplitude of the linear strain parameter, and perform weighted optimization on the corrected pure strain component based on the strain fit degree to generate point strain parameters that include strain amplitude, variation trend and spatial correlation.

7. The method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing according to claim 1, characterized in that, Step S3 includes the following steps: Step S31: Divide the multidimensional feature vector of surrounding rock deformation according to different monitoring sections of the tunnel. Each monitoring section corresponds to a feature vector sample. Extract the linear strain distribution parameters and point strain parameters in the feature vector as the core feature dimensions, and the construction correlation law and geological coupling relationship as auxiliary feature dimensions. Step S32: Using the density peak clustering algorithm, calculate the Euclidean distance between the feature vector samples of each cross section to determine the density core points in the feature space; based on the density value and distance value corresponding to the density core points, classify the deformation mode categories corresponding to uniform deformation, local mutation, and gradual accumulation. Step S33: For different deformation mode categories, extract the excavation progress and support strength process progress indicators from the construction time sequence data and assign them corresponding process progress weights; extract the surrounding rock integrity coefficient and uniaxial compressive strength index from the geological attribute data and calculate the lithological influence coefficient. Step S34: Substitute the process progress weight and lithological influence coefficient into the multidimensional feature vector of surrounding rock deformation, and perform weighted optimization on the core feature dimension and auxiliary feature dimension; construct a deformation prediction model based on an improved decision tree, take the optimized feature vector as input, take the actual deformation of the surrounding rock monitored in history as output, optimize the decision tree nodes through a pruning algorithm, train the model and input the real-time feature vector to generate the dynamic prediction value of surrounding rock deformation corresponding to each monitoring section.

8. The method for monitoring the deformation of surrounding rock in mountain tunnels based on distributed optical fiber sensing according to claim 7, characterized in that, Step S4 includes the following steps: Step S41: Obtain the displacement data of the surrounding rock section monitored by the total station at the tunnel site and the deep displacement data of the surrounding rock monitored by the multi-point displacement gauge. Combine the two to obtain the actual deformation of the surrounding rock at each monitoring section. Match the dynamic prediction value of the surrounding rock deformation with the actual deformation of the surrounding rock according to the section. Calculate the absolute error, relative error and root mean square error of each monitoring section to form the surrounding rock deformation prediction error matrix. Step S42: Set an error threshold. If the prediction error matrix of the surrounding rock deformation of a certain monitoring section is greater than the error threshold, it is determined that the prediction result of the section exceeds the allowable error range. Based on the surrounding rock deformation prediction error matrix, the gradient descent algorithm is used to iteratively correct the process progress weight and lithological influence coefficient in the deformation prediction model. The partial derivatives of each error with respect to the process progress weight and lithological influence coefficient are calculated. The parameter values ​​are adjusted along the direction of error reduction until the surrounding rock deformation prediction error matrix of all sections is less than or equal to the error threshold, and the corrected accurate prediction value of the surrounding rock deformation is obtained. Step S43: According to the tunnel mileage section and the surrounding rock deformation mode category, the corrected accurate prediction value of surrounding rock deformation is regionally aggregated, and the average deformation amount and deformation ratio of each type of deformation mode in each mileage section are calculated. Step S44: Using 3D geographic information system technology, spatially correlate the tunnel axis, monitoring section location, and aggregated deformation data to generate a 3D model of the spatiotemporal distribution of surrounding rock deformation; mark areas where the deformation exceeds the warning value in red to generate a visualization result of the prediction of surrounding rock deformation in the mountain tunnel.

9. A mountain tunnel surrounding rock deformation monitoring system based on distributed optical fiber sensing, characterized in that, For executing the method for monitoring the deformation of surrounding rock in a mountain tunnel based on distributed optical fiber sensing as described in claim 1, the system for monitoring the deformation of surrounding rock in a mountain tunnel based on distributed optical fiber sensing comprises: The surrounding rock area data acquisition module is used to collect raw Brillouin frequency shift data along the sensing optical cable and raw center wavelength data of each sensing node in the monitoring area of ​​the surrounding rock of the mountain tunnel. At the same time, it acquires tunnel construction progress data and surrounding rock geological survey data. The module performs outlier removal and noise filtering on the raw Brillouin frequency shift data and raw center wavelength data, and performs format standardization and time sequence alignment on the tunnel construction progress data and surrounding rock geological survey data to generate preprocessed fiber optic sensing basic data, construction time sequence data and geological attribute data. The surrounding rock deformation feature analysis module is used to extract surrounding rock deformation feature parameters based on fiber optic sensing data, including the axial strain values ​​of each monitoring point of the sensing cable and the strain distribution parameters along the surrounding rock; it analyzes the center wavelength offset of each sensing node and separates the pure strain component to generate the point strain parameters of key points of the surrounding rock; at the same time, it constructs a multi-dimensional feature vector of surrounding rock deformation based on the surrounding rock deformation feature parameters combined with construction time series data and geological attribute data. The surrounding rock deformation prediction module is used to perform deformation pattern density clustering on the surrounding rock of different monitoring sections of the tunnel based on the multi-dimensional feature vector of surrounding rock deformation, so as to classify the deformation pattern categories corresponding to uniform deformation, local abrupt deformation, and gradual accumulation. For different deformation pattern categories, the module introduces the process progress weight in the construction time sequence data and the lithological influence coefficient in the geological attribute data to dynamically predict the surrounding rock deformation, and generate the dynamic prediction value of surrounding rock deformation for each monitoring section. The surrounding rock deformation visualization module is used to compare the dynamic predicted values ​​of surrounding rock deformation with the monitoring data of total station and multi-point displacement gauges deployed on-site in the tunnel, calculate the absolute error, relative error and root mean square error between the predicted and measured values, and form a surrounding rock deformation prediction error matrix. Based on the surrounding rock deformation prediction error matrix, iterative correction is performed, and regional aggregation is performed according to the tunnel mileage segment and deformation mode category, thereby generating a visualization result of the spatiotemporal distribution of surrounding rock deformation in the mountain tunnel.

Images