基于實測松動圈的擴展卡爾曼濾波位移反分析方法

徐國文，何川，王均勇，王士民，唐銳

（1. 西南交通大學 交通隧道工程教育部重點實驗室，四川 成都，610031；2. 四川省交通運輸廳 公路規(guī)劃勘察設計研究院，四川 成都，610031）

基于實測松動圈的擴展卡爾曼濾波位移反分析方法

徐國文1，何川1，王均勇1，王士民1，唐銳2

（1. 西南交通大學 交通隧道工程教育部重點實驗室，四川 成都，610031；2. 四川省交通運輸廳 公路規(guī)劃勘察設計研究院，四川 成都，610031）

基于現(xiàn)場松動圈實測結果，采用雙重介質模型來考慮松動圈的影響，將施工中巖土體系統(tǒng)的反應看作動態(tài)隨機過程，采用擴展卡爾曼濾波與離散元數(shù)值耦合方法對鷓鴣山隧道特定斷面的現(xiàn)場監(jiān)測數(shù)據(jù)進行跟蹤反演，探討觀測誤差與模型誤差對濾波過程及參數(shù)識別結果的影響。研究結果表明：該方法得到的圍巖力學參數(shù)具有歷時性，可以反映洞室開挖后圍巖力學行為的動態(tài)變化，從而可以更加及時地了解圍巖的穩(wěn)定狀況；從誤差分析來看，模型誤差的取值對結果的準確性影響較大，而觀測誤差的取值對結果影響較小。因此，在實際操作過程中，要特別注意模型誤差取值的準確性。

松動圈；離散元；卡爾曼濾波

在巖石力學與工程中，利用現(xiàn)場量測的能反映系統(tǒng)力學行為的物理量來推算系統(tǒng)初始參數(shù)的方法被稱為反分析法［1-2］。國內外不少專家學者，通過解析法和數(shù)值法在位移反分析方面做了大量的研究工作，MIRANDA等［3］采用進化算法，對葡萄牙Venda Nova II水電站的地下洞室力學參數(shù)進行反演；王迎超等［4］以拱頂下沉和洞周收斂監(jiān)測數(shù)據(jù)為依據(jù)，采用正交試驗方案進行有限元正交反演分析，得到巖體計算參數(shù)；楊文東等［5］將流變模型中控制瞬時變形和時效變形的參數(shù)分開進行粒子群反演，得到了相應的流變參數(shù)；但這些文獻中對于松動圈的影響考慮較少。洞室開挖后，圍巖應力將經歷重分布過程，某些部位可能進入塑性狀態(tài)，隨著塑性變形的發(fā)展而導致圍巖產生松動圈。圍巖松動的同時其強度參數(shù)也會有所降低，因此松動圈的存在會導致位移測值的增大，對反分析結果有較大的影響。楊志法等［6］指出：為了進一步提高位移反分析應用研究的水平，需要認真考慮洞室開挖引起的松動圈對反演結果的影響；李鴻博［7］對是否考慮松動圈的反演結果進行了比較，發(fā)現(xiàn)不考慮松動圈時反演結果有較大的誤差。因此，本文作者基于現(xiàn)場松動圈實測結果，采用雙重介質模型來考慮松動圈的影響，并將施工中巖土體系統(tǒng)的反應看作動態(tài)隨機過程，將擴展卡爾曼濾波算法與離散元數(shù)值方法進行耦合，對鷓鴣山隧道特定斷面的現(xiàn)場監(jiān)測數(shù)據(jù)進行了跟蹤反演，驗證了此算法的可行性和合理性。

1 工程概況

圖1所示為鷓鴣山隧道隧址區(qū)域構造圖。在建鷓鴣山隧道為汶川至馬爾康高速公路的控制性工程，隧址區(qū)位于川西北高原南緣，邛崍山脈北端。隧道軸向總體走向呈東—西向，并以大角度與山體走向斜交，如圖1所示。隧道為雙向四車道分離式隧道，全長為8 784 m，隧道寬為13.42 m，高為10.49 m，最大埋深約為1 300 m。

圖1 鷓鴣山隧道隧址區(qū)域構造圖Fig.1 Regional tectonic map of Zhegu mountain tunnel

隧址區(qū)構造形跡以緊密線狀弧形褶皺為主，大中型斷裂構造不發(fā)育。洞身穿越 I3（鉆金樓倒轉背斜），地層巖性主要為三疊系變質千枚巖、板巖，巖體完整性總體呈破碎狀。隧道縱斷面圖如圖2所示。

從隧道開挖揭露的圍巖表明，掌子面以強-弱風化板巖為主，并夾有砂巖與千枚巖。圍巖受構造影響強烈，節(jié)理發(fā)育并交叉切割巖體，巖體呈破碎狀；千枚巖在地下水作用下，巖體軟化，部分呈泥塑狀。

鷓鴣山隧道區(qū)受構造影響強烈，巖層直立、倒轉、平臥、扭曲嚴重，斷層破裂帶、擠壓帶、次級斷層、褶皺普遍發(fā)育，圍巖破碎，加之地下水對圍巖軟化，圍巖穩(wěn)定性較差。隧道開挖過程中大變形現(xiàn)象普遍，且由變形演變形成的坍塌災害也十分頻繁，如圖3所示。

圖2 鷓鴣山隧道縱斷面圖Fig.2 Longitudinal profile of Zhegu mountain tunnel

圖3 鷓鴣山隧道典型大變形現(xiàn)象Fig.3 Typical phenomenon of large deformation of Zhegu mountain tunnel

2 松動圈測試

目前在松動圈測試中，常用有多點位移計、聲波測試與地質雷達測試方法。地質雷達法測試的優(yōu)點是不需鉆孔，精度、效率和分辨率高，靈活方便，剖面直觀，測試快速，因此在隧道工程中得到了廣泛應用［8-9］。

2.1 測量原理

對隧道工程松動圈探測而言，電磁波在向圍巖發(fā)射的過程中，遇到波阻抗不同的巖石界面時將產生反射波和折射波，其中反射波能量取決于反射系數(shù)R，可表示為

式中：ε1和 ε2分別為反射界面兩側介質的相對介電常數(shù)。 ε1為巖體的相對介電常數(shù)（ε1= 6.0～8.0），ε2為空氣或水的相對介電常數(shù)（ε（空氣）= 1.0，ε（水）= 81.0）。

圍巖松動圈內， 巖體破裂松弛，且空氣或水等雜質會通過裂隙滲入，而松動圈外的巖體較為完整。顯然，從式（1）可知：由于松動圈圍巖裂隙內填充物有空氣和水，其與完整巖體的相對介電常數(shù)相差較大，雷達發(fā)射的電磁波經過松動圈與完整巖石的分界面時必然發(fā)生強烈反射，且呈雜亂無章狀態(tài)傳播，無明顯同相軸，于是從收集處理的雷達探測剖面圖上即可確定圍巖松動范圍，得到松動圈厚度。

2.2 探測儀器與探測方案

采用美國GSSI公司的SIR-20探地雷達，主要工作技術參數(shù)如下：100 Hz屏蔽天線，每次掃描的采樣數(shù)為512點；掃描速度為100點/s；相對介電常數(shù)為7。

圖4所示為圍巖松動圈測線布置示意圖。選擇左洞3個未注漿斷面進行松動圈測試，測線沿上臺階環(huán)向布設，將每個待測斷面分成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附近斷面左拱肩發(fā)生了坍塌（圖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 擴展卡爾曼濾波與離散元耦合算法

大量工程實踐證明：采用確定性反分析法得出的結果與實測結果有較大出入，這主要是由于隧道開挖引起的圍巖動態(tài)不確定性對反演結果有較大的影響［10］。將巖土工程施工中巖土體系統(tǒng)的反應看作是一個動態(tài)的隨機過程，把擴展卡爾曼濾波器同離散元進行耦合，建立可反映巖土體動態(tài)隨機過程的卡爾曼濾波離散元耦合反分析方法。

3.1 擴展卡爾曼濾波算法

擴展卡爾曼濾波器是用隨機時間序列來代替隨機過程，其算法［11］如下。

1） 濾波方程為

實際計算時，假定系統(tǒng)是穩(wěn)態(tài)（steady state）的，即狀態(tài)遷移矩陣A=I；忽略系統(tǒng)噪聲的影響，即Q = 0，則相應的式（2）和式（3）可轉化為：

3.2 擴展卡爾曼濾波與離散元耦合算法

考慮松動圈的擴展卡爾曼濾波與離散元耦合反分析法算法流程如圖7所示。

圖7 耦合算法流程圖Fig.7 Flow chart of coupling algorithm

在耦合算法分析中，有2個關鍵問題：

1） 初值的選擇，初始狀態(tài)量的選擇對反演結果的精度有很大影響，本文采用智能算法，將第1次量測數(shù)據(jù)反演所得參數(shù)作為狀態(tài)量的初值。

2） 在擴展卡爾曼濾波算法中，對于反演系統(tǒng)而言， gk沒有顯式表達式，因此，觀測矩陣Ck（見式（9））不能通過直接微分求得，本文采用感度分析［12］求解，即：

式中：n為觀測向量的個數(shù)；m這待估向量的維數(shù)； xi為第i個待估參數(shù)；Δxi為第i個待估參數(shù)的增量； ei為單位陣。

4 反演分析

4.1 斷面支護參數(shù)

對K188+095斷面圍巖參數(shù)進行反演分析，該斷面所在區(qū)域初期支護參數(shù)為：厚度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 數(shù)值模型

本文采用離散元軟件UDEC［13］進行建模。模型以上臺階底面中心為原點，水平方向長度為80 m，下邊界取為 20 m，上邊界取為 30 m，隧道實際埋深為220 m，因此隧道頂部再施加相應的自重載荷。松動圈根據(jù)實際測試結果進行建模，其力學參數(shù)折減系數(shù)取為0.7。模型底面施加豎向約束，兩側施加水平約束，數(shù)值模型如圖10所示。

圖8 測線布置圖Fig.8 Layout of monitoring line

圖9 隧道拱頂下沉與水平收斂量測結果Fig.9 Measurement results of tunnel vault settlement and horizontal convergence

巖石的確定性參數(shù)不參與反演，結果如表1所示。

圖10 雙重介質離散元模型Fig.10 Discrete element model based on dual media

表1 巖石確定性參數(shù)Table 1 Deterministic parameters

表2 辨識參數(shù)初值Table 2 Initial calculation parameters

根據(jù)確定性反分析，得到初值x0如表2所示；對于誤差協(xié)方差矩陣P的確定，目前主要有2種方法：一種取/I=ΔP，另一種通過指定一個參數(shù)變化系數(shù)C1來進行近似估計。本文采用后一種方法，即矩陣P0的對角元素為（C1x0）2，非對角元素為 0；同樣地，觀測噪聲 R也采用該種方法，即矩陣 R的對角元素為（C2z）2，非對角元素為0，其中z為觀察值。

4.4 反分析計算結果

圖11所示為參數(shù)估計值隨著時間的變化曲線。由圖11可知：各參數(shù)在開挖初期變化較大，后期逐漸趨于收斂，說明濾波器穩(wěn)定工作。彈性模量的變化一直呈遞減趨勢，而內摩擦角及黏聚力的變化在初期有明顯波動。這是由于彈性模量對變形最為敏感，因此濾波器在初期就能很好工作，而內摩擦角及黏聚力的敏感程度低一些，濾波器初期工作有一個適應過程，導致上述現(xiàn)象的發(fā)生。

4.5 參數(shù)影響分析

卡爾曼濾波算法的不確定性主要表現(xiàn)在模型計算不確定性及觀測不確定性方面。其中，觀測不確定性由觀測噪音協(xié)方差矩陣 R量度，表現(xiàn)為 C2取值的不同。而模型計算不確定性由預報誤差協(xié)方差陣P量度，表現(xiàn)為C1取值的不同。

4.5.1 模型計算不確定性

為了探討模型計算的不確定性，對C1取不同的值進行數(shù)值計算（C1取值由0.1增至0.4），參數(shù)初值及方差的選擇與前述相同。濾波過程及參數(shù)識別結果如圖12所示。從圖12可以看出：不管C1觀測誤差如何取值，濾波過程的總體趨勢都相似，即各參數(shù)在開挖初期變化較大，后期逐漸趨于收斂，且彈性模量的變化一直呈遞減趨勢，而內摩擦角及黏聚力的變化在初期有明顯波動。C1越小，反映到參數(shù)識別過程為：參數(shù)前期調整更為平滑，且總體調整幅度更小。

4.5.2 觀測不確定性

為了探討觀測的不確定性，對不同C2取值進行數(shù)值計算（C2取值由0.001增至0.1），參數(shù)初值及方差的選擇與前述相同。濾波過程曲線與圖12相似，因此未列出。與圖12不同的是，同一濾波參數(shù)曲線中，不同C2對應的曲線相似性比圖12高，說明一定范圍內R的變化對濾波過程的影響很小。

圖11 參數(shù)辨識結果Fig.11 Parameters identification results

圖12 不同C1取值的濾波過程Fig.12 Filtering process of different values of C1

4.5.3 討論

卡爾曼濾波計算結果的觀測與模型不確定性規(guī)律可以有如下解釋：對于模型不確定性而言，由式（4）可知：隨著C1的增加，瞬變過程幅度增大，導致觀測矩陣Ck增大；由于觀測噪音協(xié)方差矩陣Rk較小，總體表現(xiàn)為卡爾曼增益矩陣 Gk的減小。Gk的減小將導致式（6）中位移觀測值與上一時刻的狀態(tài)量濾波估計值之差的增益加權變小，從而使參數(shù)調整具有前期更為平滑，且總體變化幅度減小的趨勢。因此，若要獲得最優(yōu)估計值就只有增加濾波次數(shù)；對于觀測不確定性而言，由式（4）可知：隨著 C2的增加，卡爾曼增益矩陣Gk將變小，但由于Rk較小，因此，其對計算結果的影響較小，導致不同C2值對應的曲線具有較高的相似性。

5 結論

1） 基于現(xiàn)場松動圈實測結果，采用雙重質模型來考慮松動圈的影響，將施工中巖土體系統(tǒng)的反應看作動態(tài)隨機過程，建立了擴展卡爾曼濾波與離散元數(shù)值耦合方法的位移反分析方法，從巖石力學眾多問題的不確定性及離散性本質來看，該方法無疑是一種較為合理的反分析方法。

2） 反演結果得到的圍巖力學參數(shù)具有歷時性，可以反映洞室開挖后圍巖力學行為的動態(tài)變化，從而可以更加及時地了解圍巖的穩(wěn)定狀況。且從反分析結果來看，由于彈性模量對變形最為敏感，因此，濾波器在初期就能很好工作，而內摩擦角及黏聚力的敏感程度低一些，導致濾波器在初期有明顯波動。

3） 卡爾曼濾波算法的不確定性主要表現(xiàn)在模型計算不確定性及觀測不確定性方面，從辨識結果來看，模型誤差的取值對結果的準確性影響較大，而觀測誤差影響較小。因此，在實際操作過程中，要特別注意模型誤差取值的準確性。

［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］ 王迎超， 尚岳全， 徐興華. 淺埋隧道巖土體參數(shù)正交反演及襯砌工作狀態(tài)評價［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］ 楊文東， 張強勇， 李樹才， 等. 粒子群算法在時效變形參數(shù)反演中的應用［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］. 中國礦業(yè)， 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.

（編輯 陳愛華）

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

XU Guowen1， HE Chuan1， WANG Junyong1， WANG Shimin1， TANG Rui2

（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）

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.

excavation damaged zone； DEM； extended Kalman filter

U451

A

1672-7207（2016）05-1722-08

10.11817/j.issn.1672-7207.2016.05.035

2015-06-15；

2015-08-25

國家自然科學基金資助項目（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