# Inversion Calculation Analysis of Operational Tunnel Structure Based on the Distributed Optical-Fiber Sensing System

Authors: Jie Zhao; Jiang Feng; Liang Ren
Journal: Advances in Civil Engineering
 (2017)

---

## Abstract

In combination of the actual project in Dalian Baiyun Mountain Tunnel, this paper introduces the principle of fiber optic sensor monitoring system based on optical time domain reflectometer. Then, based on the orthogonal design and even design scheme, this paper carries out a numerical experiment on the tunnel surrounding rock and establishes a regression model of the mapping relation between surrounding rock parameters of operation tunnel and the monitored displacement in order to set the difference between the monitored displacement and the calculated displacement as the fitness function. In the end, this paper carries out parameter identification based on the differential evolution algorithm. Achievements of the study proved that real-time safety warning could be realized inside the tunnel by monitoring the deformation parameters of tunnel vault at real time relying on the optical-fiber sensing system of the optical time domain reflectometer (OTDR). Parameters identification was carried out on the structure with differential evolution according to measured data and selected parameters, and great coincidence was obtained between the measured displacement and the identified parameters displacement, which proved the strong adaptability of the method.


---

## Body

## 1. Introduction

The distributed fiber optic sensing technology is a new type of sensing technology which is developed on the basis of optical time domain reflectometry (OTDR) in optical-fiber engineering [1]. Along with the development of the optical-fiber sensing technology, the distributed optical-fiber sensor is used more and more extensively in fields like civil, machinery, and aerospace structure for monitoring (Bao and Chan 2011, Feng et al., 2014). The strain monitoring system based on the distributed optical-fiber sensor is featured with long distance of monitoring and strong resistance against corrosion and interference and so forth, which solves the low-precision spatial resolution problems existing in discrete sensor based monitoring system (Glisic and Yao, 2012). Mohamad et al. (2010) analyzed the influence of newly constructed tunnels in different positions on the above old brick tunnels through BOTDR-based novel distributed optical-fiber strain measuring system [2]. Klar et al. (2014) established 2D and 3D earth shift models by utilizing the model parameters determined according to the data which is measured by the optical-fiber sensor on the spot. The earth shift model parameters have excellent anti-interference performance and accuracy according to validation [3]. When the distributed optical-fiber sensor is embedded into or surface bonded to surface of the structure, the real continuous and distributed measuring on the strain and the temperature at different points can be realized. Furthermore, different from the high cost of the discrete sensing system, the distributed optical-fiber sensing system proves effective ways to cut down the cost for structure monitoring. At present, the researches on and application of the optical time domain reflectometer (OTDR) based distributed optical-fiber sensing technology are still at the starting stage in China.As the quantum leaps in computer technology, the site safety monitoring of tunnel and method which adopts various inverse analysis calculations have been hotly focused and researched on by the engineering community. Consequently, carrying out quantitative identification analysis from the level of mechanical model is of great engineering significance and academic value [4, 5]. Scholars at home and abroad have applied particle swarm algorithm and genetic algorithm [6, 7] in the inverse calculation of parameters in rock engineering, which has achieved excellent results. However, this kind of algorithm has several disadvantages, such as slow convergence speed, long calculation time, and huge calculation amount, which are yet to be improved. Later, the differential evolution algorithm (DE) is a new global optimization algorithm brought up by Rainer Storm and Ken Price in 1995, which has no requirement for the initial values and possesses advantages, such as quick convergence speed, strong adaptability for various functions, and parallel calculation. Especially, it is applicable to the optimization of complex problems [8–11]. In order to master the stability of tunnel in which the vehicle load and environmental conditions change constantly in normal operation period and dynamically follow as well as evaluate the safety status of the tunnel, this paper developed a real-time monitoring system for tunnel vault based on OTDR distributed sensor. According to the tunnel surrounding rock structure and the vault displacement information, this paper introduces DE algorithm and then combines with orthogonal design to carry out inverse analysis on the major mechanical impact parameters of tunnel lining structure. In addition, this paper carries analysis and evaluation on the stability of tunnel structure. Furthermore, the accuracy of inverse parameters and prediction results are illustrated.


## 2. OTDR Theories and Application for Real-Time Tunnel Monitoring

As a photoelectric integrated instrument for testing optical fiber, optical Time Domain Reflectometer (OTDR) can measure the length of optical fiber, transmission loss of optical fiber, splice attenuation, and fault location. Specifically, OTDR is a photoelectric integrated instrument which uses the backward scattering of Rayleigh scattering and Fresnel reflection of light to test the optical fiber when it is transmitted through optical fiber. This backward Rayleigh scattering technology mainly measures the time it takes for the pulse light wave to transmit in the optical fiber for locating the place with scattering loss [12–14]. Widely applied in the maintenance and construction of cable circuits, it can measure the length of optical fiber, transmission loss of optical fiber, splice attenuation, and fault location.OTDR test is realized by emitting optical pulses into the optical fiber and later receiving returned information at the OTDR port. The optical pulses will scatter and reflect when transmitting in the optical fiber, due to the properties, connectors, junction points, bending, or other similar events, of the optical fiber, and some scattering and reflection will return to OTDR. The returned useful information is measured by the detector of OTDR and will be used as the time or curve slices at different points inside the optical fiber. This process repeats again and again and the final results will be averaged and shown as a track. This track indicates the strength (strong or weak) of signals in the whole optical-fiber section [15, 16]. Therefore, quick determination on the nature and location of event points in the optical-fiber link can be realized by merely observing signal graphs on the oscilloscope. This monitoring way may realize measuring of sensing volume of each point distributed along the sensor fiber and thus enjoys wide prospect for use. In addition, OTDR has advantages like easy installation, no harm to tunnel structure, low minimum resolution, strong resistance to interference, and capability for long-distance real-time online monitoring.Based on all the above merits, this paper developed a real-time monitoring system for tunnel vault based on the OTDR distributed sensor [17]. Figure 1 shows the distributed optical-fiber monitoring device for tunnel vault, with main units including the optical-fiber 1, the fiber cloth 2, the OTDR 3, and the polyimide adhesive. The optical-fiber was bonded to the tunnel vault by the polyimide adhesive. The fiber cloth was bonded, with the oxygen resin adhesive, to the external surface of the bonded optical-fiber and the polyimide adhesive, making one end of each optical-fiber connected with the OTDR. The optical fibers were equally distributed in parallel on the tunnel vault along the tunnel orientation and located at the central axis of the fiber cloth covered on the tunnel vault. Multiple optical fibers were used here according to the breadth of section of the tunnel.Figure 1
Layout of distributed optical-fiber sensors in the tunnel.The OTDR in the tunnel vault strain monitoring system emitted pulsed light into the optical fiber. Later, the pulsed light enters into the tunnel through the optical fiber and undergoes brillouin backscattering due to disturbance from the tested structure. Some scattered light and reflected light moved backward to the OTDR. The OTDR received the reflected light returned from the optical fiber and then identified loss signals in the reflected light and finally showed on the display after event study. Based on the loss conditions obtained from OTDR analysis for the full line of optical fiber, the operator could tell whether there were concrete cracks on the tunnel vault, the dislocation of cracks, and fiber cuts and locate such dislocation and fiber cuts. By arranging multiple optical fibers on the vault horizontally along the tunnel axial direction and studying the mechanism of action and location distribution, of deformations, damage, and leakages at weak points of the vault in axial direction, during the operation period of the tunnel, actual parameters can be provided to the study model and real-time safety warning can be realized inside the tunnel, so that the personal and property loss from accidents can be cut down.


## 3. Basic Principles of Differential Evolution Algorithm

Optimization of objective function is finally formed through basic principle transformation in tunnel surrounding rock structure parameter identification. The problem can be described as follows generally:(1)minEx1,x2,…,xj=1m∑i=1mYi0-Yi2xjL≤xj≤xjUi=1,2,…,n.Constraint is(2)kux=F.In the formula,Yi0 refers to measured strain value, Yi refers to calculating strain value, m refers to observation number, xj refers to to-be-inverted parameters, j refers to parameter number, xjL and xjU are, respectively, upper limit and lower limit of parameters. K and F are, respectively, finite element stiffness matrix and equivalent node force.This paper plans to carry out numerical calculation through orthogonal design and even design. In addition, range analysis and polynomial fitting are adopted to analyze the data sample generated by numerical experiment. The specific method of differential evolution algorithm realized in this paper is as follows:(1)
Three-dimensional numerical simulation model is established according to field monitoring data. Besides, L16(45) orthogonal sheet is selected to construct parameter combination scheme for numerical test. In this model, the single point restart function in ANSYS is adopted to carry out secondary loading. The strain value after secondary loading is used as the following analysis data.(2)
Sensitivity analysis is carried out on the orthogonal results. The mechanical parameters with stronger sensitivity are selected as the to-be-inversed variant while the insensitive parameters are selected as constants according to engineering experience. In addition, sample from orthogonal experiment calculation is adopted to fit in with the curve as the regression function model, thus obtaining the coefficients of the polynomial. Later, even design scheme is used to test the regression model in order to make the regression precision meet the requirements.(3)
Mean square deviation of field monitoring strain data and strain data calculated from regression model is adopted as the fitness function and substituted into the differential optimization algorithm with the to-be-inversed parameter as the optimization variant in order to achieve the surround rock parameter search.Figure2 shows displacement back analysis thinking.Figure 2
Implementation of displacement back analysis ideas.DE algorithm [18] adopts the basic principles similar to that of genetic algorithm while its evolution process is consistent with that of genetic algorithm: establishment, mutation, crossover, and selection of initial group. The only difference is that DE algorithm does not need coding and decoding operations. Instead, it directly uses the real number coding operation calculation without the need for data transformation; in addition, DE algorithm mainly adopts mutation operations of solution vector to use differential strategy in operations, which means that it makes full use of the differential vector among individuals in the group to disturb the individuals in order to achieve individual mutation operations. This mutation method can make full use of the group property, improve search capability, and compensate for the deficiency of mutation method [19].As for the problem,(3)minfx1,x2,…,xDs.txjL≤xj≤xjU,j=1,2,…,D.In above formulas,D is the dimension of solution space, xjL, xjU represent the upper limit and lower limit of xj of the number j components, respectively. The specific process of DE algorithm is as follows.The generation of initial group generally adopts random method. The initial group randomly generated should be evenly distributed in the solution space:(4)xj,i0=xj,iL+rand0,1·xj,iU-xj,iL.In above formula,xi(0)∣xj,iL≤xj,i(0)≤xj,iU,i=1,2,…,NP;j=1,2,…,D, xi0 represents the No. i chromosome in the 0th generation of the group and represents the No. j gene on the No. i chromosome in the 0th generation of the group. NP represents the group size, which is normally 5–10 times the number of variants; rand0,1 represents a random number evenly distributed among the range of (0,1).(1) Mutation Operation. In DE algorithm, the common mutation operation is to randomly select two different individuals from the group and select their vector difference and then superpose the difference onto the to-be-mutated individual after scaling it, which means (5)vig+1=xr1g+F·xr2g-xr3g.In above formula,r1,r2,r3 are mutually different integers in [1,2,…,NP] and not equal to i. F is the mutation factor, which is the main control parameter and generally within the range of 0.5~1. xig represents the No. i individual in the No. g generation group. During evolution process, each gene on the chromosome should be judged to see if they satisfy the boundary conditions in order to guarantee the effectiveness. If they do not, the method to generate initial group should be adopted to regenerate new genes.(2) Crossover Operation. Crossover operation is aimed at increasing the group diversity to avoid regional convergence phenomenon. It is divided into exponential crossover mode and binary crossover mode. In this paper, the binary crossover mode is adopted to carry out research on differential evolution algorithm for the reason that the binary crossover is more easily achieved in programming and more convenient to be used. Its specific operations are as follows:(6)uj,ig+1=vj,ig+1,if  rand0,1≤CRorj=jrandxj,ig,otherwise.In above formulas,xj,ig represents the No. g generation group, vj,ig+1 represents the intermediate mutation individual, CR represents the crossover probability and is generally within the range of [0,1], and jrand is a random value within [1,2,…,D].(3) Selection Operation. Greedy search algorithm is adopted for selection operation to substitute the experimental vector uig+1 obtained from crossover operation and the target vector xig into the target function for comparison and optimization selection. Its specific operations are as follows:(7)xig+1=uig+1,iffuig+1<fxigxig,otherwise.DE algorithm parameter search flow is shown in Figure3.Figure 3
Differences in evolutionary algorithm flow chart of the search parameters.


## 4. Engineering Case

Dalian Baiyun Tunnel was built in 1984 and completed and put into use in 1986, and it belongs to a double-hole two-way four-lane tunnel, the west hole is 374 meters long, and the east hole is 400 meters long. The tunnel belongs to a straight wall lining structure. The tunnel is designed according to the principle of mining method, and a composite structure is adopted. Existing Baiyun Tunnel has been operated for 30 years at present. Lining concrete suffers from serious aging and cracking phenomena. A new tunnel is excavated and constructed between two existing tunnels according to the requirements of engineering construction. Existing Baiyun Tunnel is reinforced through injection of lightweight concrete and grouting; the healthy and safe state monitoring still should be reinforced. Therefore, an optical-fiber health monitoring system is introduced this time. A lot of monitoring data are provided for evaluating structure status during operation period and predicating operation status and engineering service life through monitoring the tunnel structure status and other working status in order to achieve online monitoring for the safe state of existing Baiyun Tunnel arch and lining structure in real time, thereby guaranteeing safe operation of the tunnel.Since Baiyun Tunnel is constructed on rock foundation, the excavation depth is deeper and topsoil inhomogeneity is not considered. Meanwhile, since the tunnel is in the linear state during most operating period (nonlinear state is on the brink of destruction), SOILD45 unit of the linear elastic constitutive relation is adopted for simulating the mountain. The model contains a total of 14688 nodes and 12870 units. Ideal elastic-plastic constitutive model and D-P yield criterion are adopted for simulating surrounding rock. Crustal stress is not considered temporarily. Rock mass gravity stress is considered only (as shown in Figure4).Figure 4
Arrangement plan of measuring points within model grid.
(a)
(b)Firstly, we adopt prior information to select five representative mechanical parameters of surrounding rock: elastic modulus, Poisson’s ratio, cohesive force, internal friction angle, and tensile strength; then we carry out orthogonal experiment and numerical simulation according to the orthogonal sheet L16(45). In the process, we adopt single restart function in ANSYS for numerical simulation and use the load of upper covering soil and rock as the initial load. Later, we apply Hooke’s law to calculate the change of strain force by monitoring strain and assumed elastic modulus and then carry out secondary loading analysis on the tunnel based on the above. Next, we carry out a sensitivity analysis on the results of orthogonal experiment and then select parameters E and c with strong sensitivity as the to-be-inversed parameters while the rest three parameters are selected as constants according to engineering experience. Calculation effectiveness and efficiency are considered comprehensively. The particularity of tunnel belonging to a slender structure is further combined. Japan Shunsuke Sakurai is adopted as reference. Finite element is utilized for realizing back analysis displacement of plane strain problem. Formula (8) shows the three-party relational expression among displacement, elastic modulus, cohesion, and inner frictional angle, wherein elastic modulus is negatively correlated with displacement, and the cohesion is negatively correlated with displacement, the influence of inner frictional angle on displacement is complex and variable. Its regression model [20] is shown as follows:(8)y=α0+α1E0E+α2cc0+α3cc02+α4tan⁡φ+α5tan2φ.In the above formula,c0 refers to empirical estimated value of cohesion; E0 refers to empirical estimated value of rock elastic modulus; α0~α5 refer to coefficients of the regression model.The parameters of the surrounding rock have influence on the structure. The above sensitivity analysis shows that the elastic modulus and cohesive force are parameters with higher influence on surrounding rock. The two parameters are selected as inversion objects in order to simplify the numerical test. The mathematical model is shown as follows:(9)y=α0+α1E0E+α2cc0+α3cc02.In above formula, we suppose thatE0=5.0 GPa, c0=0.45 MPa. Firstly, we set the permissible range of the to-be-inversed parameters through precalculation of numerical experiment and engineering experience in order to reduce the calculation amount. Therefore, we can get that E = 3.5~6.5 GPa, c = 0.3~0.6 MPa. Then, we carry out secondary numerical experiment to obtain the result data and adopt the least square method to fit in the coefficients of the polynomial. Thus, we can get the mathematical relation between the displacement and the parameters of surrounding rock. Later, we adopt even design scheme to verify the fitting relation formula and find that the maximum error rate is smaller than 6.0%. The error values are listed in Table 1, which shows that this regression model can basically reflect the relation between the parameters and the displacement and thus can be applied into subsequent parameter identification program.Table 1
Contrast on regression model forecast with FEM calculation.
E (GPa)
C (MPa)
Error value (%)
Calculation result of finite element
Calculation result of regression model
A
3.5
0.35
-
1.788
-
5.993
-
5.886
4.2
0.50
-
4.961
-
4.987
-
4.740
4.8
0.60
4.715
-
5.755
-
6.026
5.3
0.3
1.485
-
6.856
-
6.958
6.0
0.45
-
4.827
-
4.288
-
4.081
6.5
0.55
2.38
-
3.518
-
3.602
B
3.5
0.35
5.853
-
5.141
-
5.442
4.2
0.50
-
5.797
-
5.48
-
5.162
4.8
0.60
2.961
-
6.403
-
6.593
5.3
0.3
5.117
-
6.089
-
6.401
6.0
0.45
-
4.976
-
4.255
-
4.043
6.5
0.55
-
4.398
-
4.209
-
4.024
C
3.5
0.35
-
5.785
-
5.349
-
5.040
4.2
0.50
4.278
-
5.244
-
5.468
4.8
0.60
-
4.277
-
6.832
-
6.540
5.3
0.3
4.474
-
7.388
-
7.719
6.0
0.45
1.934
-
5.426
-
5.531
6.5
0.55
-
4.275
-
4.637
-
4.439
D
3.5
0.35
-
3.426
-
4.765
-
4.602
4.2
0.50
1.178
-
5.362
-
5.425
4.8
0.60
-
1.855
-
6.891
-
6.763
5.3
0.3
-
2.585
-
6.875
-
6.697
6.0
0.45
3.448
-
5.238
-
5.419
6.5
0.55
-
1.813
-
4.359
-
4.280
E
3.5
0.35
-
5.001
-
5.623
-
5.342
4.2
0.50
5.245
-
5.561
-
5.853
4.8
0.60
-
4.53
-
5.638
-
5.383
5.3
0.3
-
5.751
-
6.457
-
6.086
6.0
0.45
3.227
-
6.168
-
6.367
6.5
0.55
1.595
-
4.564
-
4.637
F
3.5
0.35
1.709
-
4.637
-
4.716
4.2
0.50
5.67
-
5.427
-
5.735
4.8
0.60
-
1.159
-
6.824
-
6.745
5.3
0.3
1.172
-
6.835
-
6.915
6.0
0.45
-
4.231
-
5.286
-
5.062
6.5
0.55
-
3.491
-
4.324
-
4.173
G
3.5
0.35
-
3.108
-
5.335
-
5.169
4.2
0.50
4.393
-
5.267
-
5.498
4.8
0.60
-
2.384
-
6.850
-
6.687
5.3
0.3
3.193
-
7.412
-
7.649
6.0
0.45
-
4.546
-
5.427
-
5.180
6.5
0.55
5.798
-
4.627
-
4.895
H
3.5
0.35
-
5.578
-
5.216
-
4.925
4.2
0.50
4.788
-
5.452
-
5.713
4.8
0.60
1.23
-
6.473
-
6.553
5.3
0.3
-
2.907
-
6.028
-
5.853
6.0
0.45
4.773
-
4.243
-
4.446
6.5
0.55
-
2.701
-
4.242
-4.127
I
3.5
0.35
-
3.926
-
5.914
-
5.682
4.2
0.50
3.748
-
4.963
-
5.149
4.8
0.60
-
5.586
-
5.728
-
5.408
5.3
0.3
2.429
-
6.893
-
7.060
6.0
0.45
4.786
-
4.273
-
4.478
6.5
0.55
-
4.768
-
3.535
-
3.366We substitute the polynomial model obtained from the fitting process and the monitored displacement in Table1 into formula (3) to form fitness function of differential evolution. Then, we select the initial parameters for differential evolution. Besides, we select two optimization variants and set the number of groups as 20. Later, we suppose F=0.9 considering the influence different mutation factors and crossover factors have on the calculation. Figure 5 shows the relation curve between the iteration steps and the convergence value when CR value is within the range of 0.5~0.9. From the figure, we can see that when CR value is within the range of 0.7~0.9, the search effect is more stable and the convergence speed is quicker. However, when CR value is within the range of 0.5~0.6, the convergence speed becomes slower and even there might be circumstances of no convergence. When CR value is fixed as 0.7, it can be seen from Figure 6 that the optimal value can be converged when F value is within the range of 0.5~0.9. The major difference is that the search time becomes longer as F value increases.Figure 5
Relations between the iteration steps and the convergence value under different CR values whenF value is constant.Figure 6
Relations between the iteration steps and the convergence value under differentF values when CR value is constant.In this paper, we selectF=0.9, CR=0.7 in order to carry out parameter optimization and identification on this model. The optimal parameters identified by the algorithm are E=4.538, c=0.483. Then, we input these two parameters into the numerical model and calculate the comparison between the corresponding monitored displacement and the actually measured displacement as shown in Table 2. It can be known from the table that the maximum absolute error between the monitored displacement and the calculation result of to-be-inversed parameters is 5.798%, which shows that the parameter identification based on differential evolution algorithm can well identify the mechanical parameters of the surrounding rock structure of the tunnel and possess good predictability.Table 2
Comparison between tiny displacements computed via identification parameters and measured.
Measuring point
A
B
C
D
E
F
G
H
I
Monitored displacement (mm)
-
2.374
-
2.628
-
2.774
-
2.855
-
2.871
-
2.823
-
2.711
-
2.531
-
2.374
Forecasted displacement (mm)
-
2.431
-
2.512
-
2.656
-
2.803
-
2.917
-
2.724
-
2.868
-
2.462
-
2.261
Relative error (%)
1.190
-
2.199
-
2.138
-
0.907
0.798
-
1.746
2.899
-
1.351
-
2.384


## 5. Conclusion

Following conclusions can be drawn based on reverse analysis on the operational tunnel structure through the distributed optical-fiber sensing system in the paper:(1)
Real-time safety warning can be realized inside the tunnel by monitoring at real time the deformation parameters of tunnel vault with the OTDR optical-fiber sensing system.(2)
Study on the selection of parametersF and CR were made based on the differential evolution theories, which proved the influence of different F and CR levels on the convergence results.(3)
Parameters identification was done on the structure with differential evolution according to measured data and selected parameters, and great coincidence between the measured displacement and the identified parameters displacement was obtained, which proved the strong adaptability of this method.



---

*Source: 6159095-2017-08-17.xml*