DOI: 10.11817/j.issn.1672-7207.2016.05.035

基于实测松动圈的扩展卡尔曼滤波位移反分析方法

徐国文¹，何川¹，王均勇¹，王士民¹，唐锐²

 (1. 西南交通大学 交通隧道工程教育部重点实验室，四川 成都，610031；

2. 四川省交通运输厅 公路规划勘察设计研究院，四川 成都，610031)

摘要：基于现场松动圈实测结果，采用双重介质模型来考虑松动圈的影响，将施工中岩土体系统的反应看作动态随机过程，采用扩展卡尔曼滤波与离散元数值耦合方法对鹧鸪山隧道特定断面的现场监测数据进行跟踪反演，探讨观测误差与模型误差对滤波过程及参数识别结果的影响。研究结果表明：该方法得到的围岩力学参数具有历时性，可以反映洞室开挖后围岩力学行为的动态变化，从而可以更加及时地了解围岩的稳定状况；从误差分析来看，模型误差的取值对结果的准确性影响较大，而观测误差的取值对结果影响较小。因此，在实际操作过程中，要特别注意模型误差取值的准确性。

关键词：松动圈；离散元；卡尔曼滤波

中图分类号：U451          文献标志码：A         文章编号：1672-7207(2016)05-1722-08

Displacement-based back analysis method using extend Kalman filter considering excavation damaged zone

XU Guowen¹, HE Chuan¹, WANG Junyong¹, WANG Shimin¹, TANG Rui²

 (1. Key Laboratory of Transportation Tunnel Engineering, Ministry of Education, Southwest Jiaotong University,

Chengdu 610031, China;

2. Sichuan Province Transportation Department Highway Planning, Survey, Design and Research Institute,

 Chengdu 610031, China)

Abstract: Based on the field-measured result of the excavation damaged zone (EDZ), the extended Kalman filtering & DEM coupling back analysis method ,which adopts dual media model to consider the impact of broken zone and regards the reaction of rock system as a dynamic stochastic process during construction, was established and was applied to a specified section of Zhegu mountain tunnel. The effect of observation error and model error on the filtering process and parameter identification results was discussed. The results show that the mechanical parameters obtained by this method has a diachronic characteristic, which can reflect the dynamic change of the mechanical behavior of rock after excavation, and can be more timely to understand the stability situation of surrounding rock. From the identification results, the model error has a greater impact on the outcome than the observation error, Therefore, in actual operation, particular attention must be paid to get the accurate model error.

Key words: excavation damaged zone; DEM; extended Kalman filter

在岩石力学与工程中，利用现场量测的能反映系统力学行为的物理量来推算系统初始参数的方法被称为反分析法^[1-2]。国内外不少专家学者，通过解析法和数值法在位移反分析方面做了大量的研究工作，MIRANDA等^[3]采用进化算法，对葡萄牙Venda Nova II水电站的地下洞室力学参数进行反演；王迎超等^[4]以拱顶下沉和洞周收敛监测数据为依据，采用正交试验方案进行有限元正交反演分析，得到岩体计算参数；杨文东等^[5]将流变模型中控制瞬时变形和时效变形的参数分开进行粒子群反演，得到了相应的流变参数；但这些文献中对于松动圈的影响考虑较少。洞室开挖后，围岩应力将经历重分布过程，某些部位可能进入塑性状态，随着塑性变形的发展而导致围岩产生松动圈。围岩松动的同时其强度参数也会有所降低，因此松动圈的存在会导致位移测值的增大，对反分析结果有较大的影响。杨志法等^[6]指出：为了进一步提高位移反分析应用研究的水平，需要认真考虑洞室开挖引起的松动圈对反演结果的影响；李鸿博^[7]对是否考虑松动圈的反演结果进行了比较，发现不考虑松动圈时反演结果有较大的误差。因此，本文作者基于现场松动圈实测结果，采用双重介质模型来考虑松动圈的影响，并将施工中岩土体系统的反应看作动态随机过程，将扩展卡尔曼滤波算法与离散元数值方法进行耦合，对鹧鸪山隧道特定断面的现场监测数据进行了跟踪反演，验证了此算法的可行性和合理性。

1  工程概况

图1所示为鹧鸪山隧道隧址区域构造图。在建鹧鸪山隧道为汶川至马尔康高速公路的控制性工程，隧址区位于川西北高原南缘，邛崃山脉北端。隧道轴向总体走向呈东—西向，并以大角度与山体走向斜交，如图1所示。隧道为双向四车道分离式隧道，全长为8 784 m，隧道宽为13.42 m，高为10.49 m，最大埋深约为1 300 m。

图1  鹧鸪山隧道隧址区域构造图

Fig. 1  Regional tectonic map of Zhegu mountain tunnel

隧址区构造形迹以紧密线状弧形褶皱为主，大中型断裂构造不发育。洞身穿越I3(钻金楼倒转背斜)，地层岩性主要为三叠系变质千枚岩、板岩，岩体完整性总体呈破碎状。隧道纵断面图如图2所示。

从隧道开挖揭露的围岩表明，掌子面以强-弱风化板岩为主，并夹有砂岩与千枚岩。围岩受构造影响强烈，节理发育并交叉切割岩体，岩体呈破碎状；千枚岩在地下水作用下，岩体软化，部分呈泥塑状。

鹧鸪山隧道区受构造影响强烈，岩层直立、倒转、平卧、扭曲严重，断层破裂带、挤压带、次级断层、褶皱普遍发育，围岩破碎，加之地下水对围岩软化，围岩稳定性较差。隧道开挖过程中大变形现象普遍，且由变形演变形成的坍塌灾害也十分频繁，如图3所示。

图2  鹧鸪山隧道纵断面图

Fig. 2  Longitudinal profile of Zhegu mountain tunnel

图3  鹧鸪山隧道典型大变形现象

Fig. 3  Typical phenomenon of large deformation of Zhegu mountain tunnel

2  松动圈测试

目前在松动圈测试中，常用有多点位移计、声波测试与地质雷达测试方法。地质雷达法测试的优点是不需钻孔，精度、效率和分辨率高，灵活方便，剖面直观，测试快速，因此在隧道工程中得到了广泛应   用^[8-9]。

2.1  测量原理

对隧道工程松动圈探测而言，电磁波在向围岩发射的过程中，遇到波阻抗不同的岩石界面时将产生反射波和折射波，其中反射波能量取决于反射系数R，可表示为

               (1)

式中：和分别为反射界面两侧介质的相对介电常数。为岩体的相对介电常数(= 6.0~8.0)，为空气或水的相对介电常数(= 1.0，= 81.0)。

围岩松动圈内， 岩体破裂松弛，且空气或水等杂质会通过裂隙渗入，而松动圈外的岩体较为完整。显然，从式(1)可知：由于松动圈围岩裂隙内填充物有空气和水，其与完整岩体的相对介电常数相差较大，雷达发射的电磁波经过松动圈与完整岩石的分界面时必然发生强烈反射，且呈杂乱无章状态传播，无明显同相轴，于是从收集处理的雷达探测剖面图上即可确定围岩松动范围，得到松动圈厚度。

2.2  探测仪器与探测方案

采用美国GSSI公司的SIR-20探地雷达，主要工作技术参数如下：100 Hz屏蔽天线，每次扫描的采样数为512点；扫描速度为100点/s；相对介电常数为7。

图4所示为围岩松动圈测线布置示意图。选择左洞3个未注浆断面进行松动圈测试，测线沿上台阶环向布设，将每个待测断面分成5 段(编号依次为①~⑤)进行扫描(见图4)。

图4  围岩松动圈测线布置示意图

Fig. 4  Detection line of loosing zone of surrounding rock

2.3  探测结果与分析

以K188+095断面为例，所得结果见图5和6，由于测线对称分布，故截取测线上①~③段分析。图5所示为雷达探测剖面图。

图6所示为K188+095松动圈示意图。由图6可知：该断面围岩整体松动范围较大，且呈左侧大于右侧之势，最大值在左拱肩处，达6.2 m。在开挖过程中，K188+095附近断面左拱肩发生了坍塌(图3(c))，塌方高度约为6 m，从一个侧面验证了雷达探测的可靠性。

图5  K188+095 ①~③段雷达测试剖面图

Fig. 5  GPR record of K188+095 section on lines ①-③

图6  K188+095断面松动圈示意图

Fig. 6  Excavation damage zone of section K188+095

3  扩展卡尔曼滤波与离散元耦合算法

大量工程实践证明：采用确定性反分析法得出的结果与实测结果有较大出入，这主要是由于隧道开挖引起的围岩动态不确定性对反演结果有较大的影响^[10]。将岩土工程施工中岩土体系统的反应看作是一个动态的随机过程，把扩展卡尔曼滤波器同离散元进行耦合，建立可反映岩土体动态随机过程的卡尔曼滤波离散元耦合反分析方法。

3.1  扩展卡尔曼滤波算法

扩展卡尔曼滤波器是用随机时间序列来代替随机过程，其算法^[11]如下。

1) 滤波方程为

             (2)

     (3)

式中：为时刻的状态量预报估计；为时刻的状态量滤波估计；为状态转移矩阵；为时刻的预报误差协方差阵；为时刻的滤波误差协方差矩阵；为系统噪音矩阵；为系统噪音协方差阵。

2) 更新方程

       (4)

            (5)

        (6)

式中：为卡尔曼增益矩阵；为观测矩阵；为观测噪音协方差矩阵；为时刻的滤波误差协方差矩阵；为时刻的状态量滤波估计；为时刻的观测向量。

实际计算时，假定系统是稳态(steady state)的，即状态迁移矩阵；忽略系统噪声的影响，即Q = 0，则相应的式(2)和式(3)可转化为：

               (7)

               (8)

3.2  扩展卡尔曼滤波与离散元耦合算法

考虑松动圈的扩展卡尔曼滤波与离散元耦合反分析法算法流程如图7所示。

图7  耦合算法流程图

Fig. 7  Flow chart of coupling algorithm

在耦合算法分析中，有2个关键问题：

1) 初值的选择，初始状态量的选择对反演结果的精度有很大影响，本文采用智能算法，将第1次量测数据反演所得参数作为状态量的初值。

2) 在扩展卡尔曼滤波算法中，对于反演系统而言，没有显式表达式，因此，观测矩阵(见式(9))不能通过直接微分求得，本文采用感度分析^[12]求解，即：

           (9)

      (10)

式中：n为观测向量的个数；m这待估向量的维数；为第i个待估参数；为第i个待估参数的增量；为单位阵。

4  反演分析

4.1  断面支护参数

对K188+095断面围岩参数进行反演分析，该断面所在区域初期支护参数为：厚度26 cm的C20喷射混凝土；直径8 mm的钢筋网@20×20 cm；I20b 工字钢拱@60 cm；直径22 mm的组合式锚杆@120×100 cm，长度L为3 m。

4.2  测线布置及量测结果

隧道拱顶下沉与洞周收敛采用非接触量测方式进行测量，测量精度为1 mm。测线布置如图8所示，测试结果如图9所示。

4.3  数值模型

本文采用离散元软件UDEC^[13]进行建模。模型以上台阶底面中心为原点，水平方向长度为80 m，下边界取为20 m，上边界取为30 m，隧道实际埋深为   220 m，因此隧道顶部再施加相应的自重载荷。松动圈根据实际测试结果进行建模，其力学参数折减系数取

图8  测线布置图

Fig. 8  Layout of monitoring line

图9  隧道拱顶下沉与水平收敛量测结果

Fig. 9  Measurement results of tunnel vault settlement and horizontal convergence

为0.7。模型底面施加竖向约束，两侧施加水平约束，数值模型如图10所示。

岩石的确定性参数不参与反演，结果如表1所示。

根据确定性反分析，得到初值x₀如表2所示；对于误差协方差矩阵P的确定，目前主要有2种方法：一种取，另一种通过指定一个参数变化系数C₁来进行近似估计。本文采用后一种方法，即矩阵P₀的对角元素为(C₁x₀)²，非对角元素为0；同样地，观测噪声R也采用该种方法，即矩阵R的对角元素为(C₂z)²，非对角元素为0，其中z为观察值。

图10  双重介质离散元模型

Fig. 10  Discrete element model based on dual media

表1  岩石确定性参数

Table 1  Deterministic parameters

表2  辨识参数初值

Table 2  Initial calculation parameters

4.4  反分析计算结果

图11所示为参数估计值随着时间的变化曲线。由图11可知：各参数在开挖初期变化较大，后期逐渐趋于收敛，说明滤波器稳定工作。弹性模量的变化一直呈递减趋势，而内摩擦角及黏聚力的变化在初期有明显波动。这是由于弹性模量对变形最为敏感，因此滤波器在初期就能很好工作，而内摩擦角及黏聚力的敏感程度低一些，滤波器初期工作有一个适应过程，导致上述现象的发生。

4.5  参数影响分析

卡尔曼滤波算法的不确定性主要表现在模型计算不确定性及观测不确定性方面。其中，观测不确定性由观测噪音协方差矩阵R量度，表现为C₂取值的不同。而模型计算不确定性由预报误差协方差阵P量度，表现为C₁取值的不同。

4.5.1  模型计算不确定性

为了探讨模型计算的不确定性，对C₁取不同的值进行数值计算(C₁取值由0.1增至0.4)，参数初值及方差的选择与前述相同。滤波过程及参数识别结果如图12所示。从图12可以看出：不管C ₁观测误差如何取值，滤波过程的总体趋势都相似，即各参数在开挖初期变化较大，后期逐渐趋于收敛，且弹性模量的变化一直呈递减趋势，而内摩擦角及黏聚力的变化在初期有明显波动。C₁越小，反映到参数识别过程为：参数前期调整更为平滑，且总体调整幅度更小。

4.5.2  观测不确定性

为了探讨观测的不确定性，对不同C₂取值进行数值计算(C₂取值由0.001增至0.1)，参数初值及方差的选择与前述相同。滤波过程曲线与图12相似，因此未列出。与图12不同的是，同一滤波参数曲线中，不同C₂对应的曲线相似性比图12高，说明一定范围内R的变化对滤波过程的影响很小。

图11  参数辨识结果

Fig. 11  Parameters identification results

图12  不同C₁取值的滤波过程

Fig. 12  Filtering process of different values of C₁

4.5.3  讨论

卡尔曼滤波计算结果的观测与模型不确定性规律可以有如下解释：对于模型不确定性而言，由式(4)可知：随着C₁的增加，瞬变过程幅度增大，导致观测矩阵C_k增大；由于观测噪音协方差矩阵R_k较小，总体表现为卡尔曼增益矩阵G_k的减小。G_k的减小将导致式(6)中位移观测值与上一时刻的状态量滤波估计值之差的增益加权变小，从而使参数调整具有前期更为平滑，且总体变化幅度减小的趋势。因此，若要获得最优估计值就只有增加滤波次数；对于观测不确定性而言，由式(4)可知：随着C₂的增加，卡尔曼增益矩阵G_k将变小，但由于R_k较小，因此，其对计算结果的影响较小，导致不同C₂值对应的曲线具有较高的相似性。

5  结论

1) 基于现场松动圈实测结果，采用双重质模型来考虑松动圈的影响，将施工中岩土体系统的反应看作动态随机过程，建立了扩展卡尔曼滤波与离散元数值耦合方法的位移反分析方法，从岩石力学众多问题的不确定性及离散性本质来看，该方法无疑是一种较为合理的反分析方法。

2) 反演结果得到的围岩力学参数具有历时性，可以反映洞室开挖后围岩力学行为的动态变化，从而可以更加及时地了解围岩的稳定状况。且从反分析结果来看，由于弹性模量对变形最为敏感，因此，滤波器在初期就能很好工作，而内摩擦角及黏聚力的敏感程度低一些，导致滤波器在初期有明显波动。

3) 卡尔曼滤波算法的不确定性主要表现在模型计算不确定性及观测不确定性方面，从辨识结果来看，模型误差的取值对结果的准确性影响较大，而观测误差影响较小。因此，在实际操作过程中，要特别注意模型误差取值的准确性。

参考文献：

[1] 杨林德, 冯紫良, 朱合华, 等. 岩土工程问题的反演理论与工程实践[M]. 北京: 科学出版社, 1996: 60-80.

YANG Linde, FENG Ziliang, ZHU Hehua, et al. Back analysis theory and practice in geotechnical engineering[M]. Beijing: Science Press, 1996: 60-80.

[2] 杨志法, 王思敬, 冯紫良, 等. 岩土工程反分析原理及应用[M]. 北京: 地震出版社, 2002: 12-20.

YANG Zhifa, WANG Sijing, FENG Ziliang, et al. Principle of back analysis and its application to geotechnical engineering[M]. Beijing: Earthquake Press, 2002: 12-20.

[3] MIRANDA T, ECLAIRCY-CAUDRON S, DIAS D, et al. Back analysis of geomechanical parameters by optimization of a 3D model of an underground structure[J]. Tunnelling and Underground Space Technology, 2010, 26(11): 659-673.

[4] 王迎超, 尚岳全, 徐兴华. 浅埋隧道岩土体参数正交反演及衬砌工作状态评价[J]. 中南大学学报(自然科学版), 2011, 42(6): 1764-1771.

WANG Yingchao, SHANG Yuequan, XU Xinghua. Orthogonal back analysis of geotechnical parameters and working state evaluation of lining of shallow buried tunnel[J]. Journal of Central South University (Science and Technology), 2011, 42(6): 1764-1771.

[5] 杨文东, 张强勇, 李树才, 等. 粒子群算法在时效变形参数反演中的应用[J]. 中南大学学报(自然科学版), 2013, 44(1): 282-288.

YANG Wendong, ZHANG Qiangyong, LI Shucai. Application of particle swarm optimization in time dependent parameters inversion[J]. Journal of Central South University (Science and Technology), 2013, 44(1): 282-288

[6] 杨志法, 熊顺成, 王存玉, 等. 关于位移反分析的某些考虑[J]. 岩石力学与工程学报, 1995, 14(1): 11-16.

YANG Zhifa, XIONG Shuncheng, WANG Chunyu, et al. Some consideration of the back-analysis from displacements[J]. Chinese Journal of Rock Mechanics and Engineering, 1995, 14(1): 11-16.

[7] 李鸿博. 考虑松动圈的围岩反分析在隧道工程中的应用[D]. 上海: 同济大学土木工程学院, 2006: 15-25.

LI Hongbo. Application of displacement back analysis of taking into loosing zone for tunnel engineering[D]. Shanghai: Tongji University. College of Civil Engineering, 2006: 15-25.

[8] 姜德义, 郑彦奎, 任松, 等. 地质偏压隧道松动圈探测与初衬开裂原因分析[J]. 中国矿业, 2008, 17(1): 101-104.

JIANG Deyi, ZHENG Yankui , REN Song, et al. Measuring broken rock zone around unsymmetrical-loading tunnel in geologic and analysis the reason of the liner cracks[J].China Mining Magazine, 2008, 17(1): 101-104.

[9] 伍永平, 翟锦, 解盘石, 等. 基于地质雷达探测技术的巷道围岩松动圈测定[J]. 煤炭科学技术, 2013, 41(3): 32-38.

WU Yongping, ZHAI Jin, XIE Panshi, et al. Measurement of loosing circle in surrounding rock of gateway based on technology of geological radar detection[J]. Coal Science and Technology, 2013, 41(3): 32-38.

[10] 孙均, 蒋树屏, 袁勇, 等. 岩土工程反演问题的随机理论与方法[M]. 广州: 汕头大学出版社, 1996: 150-180.

SUN Jun, JIANG Shuping, YUAN Yong, et al. Probability theory and method in geotechnical back analysis problem[M]. Guangzhou: Shantou University Press, 1996: 150-180.

[11] JAZWINSKI A H. Stochastic processes and filtering theory[M]. New York: Academic Press, 1970: 46-60.

[12] 蒋树屏. 扩张卡尔曼滤波器有限元法耦合算法及其隧道工程应用[J]. 岩土工程学报, 1996, 18(4): 11-19.

JIANG Shuping. Coupling algorithm of extended Kalman filter-FEM and its application to tunnel engineering[J]. Chinese Journal of Geotechnical Engineering, 1996, 18(4): 11-19.

[13] Itasca Consulting Group Inc. UDEC (Universal Distinct Element Code) User’s Manual, Version3.0[R]. Minnesota: Itasca Consulting Group, 2004: 36-58.

(编辑  陈爱华)

收稿日期：2015-06-15；修回日期：2015-08-25

基金项目(Foundation item)：国家自然科学基金资助项目(U1361210)；国家科技支撑计划项目(2013BAB10B04，2012BAG05B03) (Project(U1361210) supported by the National Natural Science Foundation of China; Projects(2013BAB10B04, 2012BAG05B03) supported by the National Science and Technology Pillar Program)

通信作者：何川，博士 ，教授，博士生导师，从事隧道与地下工程方向的研究；E-mail: chuanhe21@163.com