魯光銀，羅帥，朱自強，石克亮，夏成志

基于柱坐標系的隧道空間全波場數(shù)值模擬與分析

魯光銀，羅帥，朱自強，石克亮，夏成志

(中南大學 地球科學與信息物理學院，湖南 長沙 410083)

不同于地面地震，隧道空間的特殊結(jié)構(gòu)使得其波場十分復雜，基于實際隧道模型的全波場數(shù)值模擬是有效認識隧道空間內(nèi)波場傳播規(guī)律的重要手段。關(guān)于這方面的研究以往多是基于直角坐標系，而隧道作為類似空氣柱的特殊地質(zhì)模型，采用柱坐標系下的全波場數(shù)值模擬對認識隧道空間的波場特征更有意義?；谥鴺讼迪碌囊浑A彈性波速度-應(yīng)力方程，推導任意偶數(shù)階交錯網(wǎng)格差分格式，設(shè)計了人工截斷邊界和隧道自由邊界，對巖性分界面、斷層破碎帶、溶洞等隧道典型不良地質(zhì)模型進行數(shù)值模擬與波場特征分析。研究結(jié)果表明：基于柱坐標系下的一階彈性波速度-應(yīng)力方程，采用交錯網(wǎng)格有限差分算法結(jié)合FCT技術(shù)，可實現(xiàn)隧道空間全波場高精度數(shù)值模擬，波場特征符合波的運動學與動力學特征，為隧道空間復雜波場特征的認識和資料解譯提供有效的依據(jù)。

超前預報；隧道空間；柱坐標系；全波場；數(shù)值模擬

隨著交通運輸?shù)目焖侔l(fā)展，越來越多的隧道被建造。然而，隧道作為一項隱蔽性地下工程，在勘察設(shè)計階段很難精確的發(fā)現(xiàn)不良地質(zhì)體。在隧道開挖過程中，有可能遇到斷層、溶洞，暗河等不良地質(zhì)現(xiàn)象，極易誘發(fā)突水、涌泥、巖爆等地質(zhì)災(zāi)害，進而造成工期延誤，經(jīng)濟損失，甚至人員傷亡等重大安全事故[1?2]。地震法作為隧道超前預報的主要方法之一，包括TSP，TST，HSP，TRT和陸地聲納法等[3?8]，以其長距離、相對高精度而得以廣泛應(yīng)用。但是隧道空間內(nèi)的波場是十分復雜的，接收記錄上往往存在直達波、反射波、面波、縱波、橫波轉(zhuǎn)換波[9]，這無疑為數(shù)據(jù)處理與成果解譯增加了難度?；谒淼揽臻g的全波場數(shù)值模擬是認識隧道空間內(nèi)復雜波場傳播規(guī)律的基礎(chǔ)，同時也是觀測系統(tǒng)設(shè)計，數(shù)據(jù)采集、處理、解釋的重要依據(jù)。近年來有不少學者進行了相關(guān)研究，其中絕大部分是基于直角坐標系，魯光銀等[10]采用直接法處理隧道自由邊界條件，并進行了地質(zhì)界面模型的正演模擬與偏移成像；劉江平，程飛等[9, 11]采用吸收邊界條件處理隧道自由邊界，并進行了不同觀測系統(tǒng)和不同傾角地質(zhì)界面對波場特征和偏移成像的影響研究；凌飛等[12]基于黏彈性介質(zhì)并在隧道周圍引入一開挖損傷帶對隧道典型不良地質(zhì)現(xiàn)象進行了全波場數(shù)值模擬；宋杰[13]針對隧道工程建設(shè)中可能遇到的典型不良地質(zhì)體，建立了巖性分界面、斷層以及巖溶體等多個三維地質(zhì)模型，進行了三維地震波探測正演模擬，分析總結(jié)了隧道典型不良地質(zhì)體三維波場響應(yīng)和地震成像特征；王京[14]基于二階彈性波位移波動方程采用有限元方法對隧道典型不良地質(zhì)體進行了數(shù)值模擬，并重點探討了有限元方法實現(xiàn)自由邊界條件的處理方法；Jetschny等[15]考慮了隧道空間的影響，利用面波和轉(zhuǎn)換橫波對隧道前方的地質(zhì)構(gòu)造進行預測，并且討論了這2種波的轉(zhuǎn)換機理；YANG等[16]基于三維彈性波方程交錯網(wǎng)格有限差分算法，實現(xiàn)了煤礦巷道模型的三維槽波勘探數(shù)值模擬；Harmankaya等[17]實現(xiàn)了利用隧道掘進過程中產(chǎn)生的噪聲作為震源進行超前探測的數(shù)值模擬，并驗證了該方法對散射體的定位能力；LIU等[18]進行了極坐標系下的全波場數(shù)值模擬，并處理了復雜的邊界條件；查欣杰等[19]通過構(gòu)建二維含低速異常的隧道介質(zhì)模型，研究了隧道彈性波場傳播規(guī)律和異常體邊界成像準確性。模擬隧道空間三維波場最經(jīng)濟的方法是把構(gòu)造模型包括震源近似為軸對稱模型。而相較于直角坐標系采用柱坐標系進行數(shù)值模擬具有以下2點優(yōu)勢[20]：1) 直角坐標系中的有限差分算法的公式比較簡單, 但是處理軸對稱地質(zhì)構(gòu)造需要比較多的計算成本；2) 對于具有軸對稱性的波場問題, 可以轉(zhuǎn)化為二維問題, 大大簡化計算。論文基于柱坐標系下的一階彈性波速度?應(yīng)力方程，將隧道與圍巖邊界作自由邊界處理，同時引入柱坐標系下的分裂式完全匹配層，采用時間二階，空間十階的高階交錯網(wǎng)格有限差分法并結(jié)合FCT技術(shù)，實現(xiàn)了柱坐標系下隧道空間的高精度全波場數(shù)值模擬，通過計算巖性分界面、斷層破碎帶、溶洞等復雜地質(zhì)模型，分析了隧道典型不良地質(zhì)現(xiàn)象的波場特征，為隧道空間復雜波場特征的認識和資料解譯提供了依據(jù)。

1 波動方程及其有限差分離散

考慮到隧道空間的近似軸對稱結(jié)構(gòu)以及三維數(shù)值模擬的計算成本，本文將隧道空間作軸對稱處理。對于具有軸對稱的完全彈性介質(zhì)，介質(zhì)參數(shù)僅是和的函數(shù)，三維模型可以簡化為二維。此時柱坐標系下的一階彈性波速度?應(yīng)力方程為：

式中：為介質(zhì)密度；和為拉梅常數(shù)；v和v分別為軸向(方向)和徑向(方向)的速度；和分別為軸向和徑向的正應(yīng)力；為切應(yīng)力。

對上述方程采用交錯網(wǎng)格[21]進行離散，得到任意偶數(shù)階差分格式：

2 穩(wěn)定性條件與邊界條件

正演模擬過程中為保證離散數(shù)值計算的穩(wěn)定性需要設(shè)置穩(wěn)定性條件，與直角坐標系不同, 柱坐標系要求時間步長更小,穩(wěn)定性條件為[22]：

式中：min為最小波長；min(Δ,Δ)為空間最小網(wǎng)格步長；max為最大波速。

若考慮為深埋隧道，軸對稱柱坐標系下隧道空間全波場數(shù)值模擬需要處理2種邊界條件，即計算區(qū)域外的人工截斷邊界和隧道自由邊界條件(圖1)。對于人工截斷邊界采用LIU[23]提出的柱坐標系下的分裂式完全匹配層處理；而隧道自由邊界條件則采用橫向各項同性介質(zhì)替換法處理[24]，即用橫向各向同性介質(zhì)近似代替自由邊界，直接令自由邊界處的正應(yīng)力為0，切應(yīng)力為0的條件則通過對自由邊界上物性參數(shù)的設(shè)定，在波動方程交錯網(wǎng)格有限差分迭代求解中實現(xiàn)。二維情況下的處理方式為：

式中：τxx和τrr分別為隧道掌子面和側(cè)壁的正應(yīng)力；ρ0，λ0和μ0為自由邊界上的密度和拉梅常數(shù)。ρ，λ和μ為彈性介質(zhì)的密度和拉梅常數(shù)。

3 觀測系統(tǒng)與地質(zhì)模型

為了獲取較好的波場記錄，采用如圖2所示的觀測系統(tǒng)。隧道長70 m，寬10 m，與實際隧道相仿。炮點置于隧道掌子面與右邊墻的交點處，埋深1 m；在隧道左右邊墻各布置了47檢波器，偏移距15 m，道間距0.5 m，埋深1.5 m。

圖2 觀測系統(tǒng)示意圖

為驗證本文算法的正確性，以及分析隧道典型不良地質(zhì)現(xiàn)象的波場特征設(shè)計了如圖3所示的3個地質(zhì)模型。模型Ⅰ為垂直巖性分界面，界面距掌子面40 m；模型Ⅱ為斷層破碎帶，破碎帶寬10 m，傾角45°；模型Ⅲ為溶洞，溶洞半徑5 m，各模型的彈性參數(shù)見表1。

(a) ModelⅠ；(b) ModelⅡ；(c) Model Ⅲ

表1 模型彈性參數(shù)表

4 數(shù)值模擬結(jié)果及分析

針對上述理論模型采用高階差分格式進行數(shù)值模擬，其中模型大小為，200 m×200 m，Δ=Δ=0.25 m，Δ=0.025 ms，震源選擇中心頻率為100 Hz的雷克子波。

4.1 巖性分界面模型

如圖4為垂直界面模型徑向分量20，35和50 ms的波場快照。由圖4可知，在=20 ms時，由于沒有遇到巖性分界面，在波場快照中可以看到直達縱波(P波)、直達橫波(S波)、隧道自由表面上產(chǎn)生的面波(R波)；同時可以看到，面波延隧道側(cè)壁傳播并在掌子面與側(cè)壁交點處轉(zhuǎn)換為繞射橫波(RS波)[25]，繞射橫波傳至側(cè)壁時又重新轉(zhuǎn)換成面波；在=35 ms時，在巖性分界面處處產(chǎn)生的反射縱波(PP波)和透射縱波(TP波)，以及由P波產(chǎn)生的反射轉(zhuǎn)換橫波(PS波)以及透射轉(zhuǎn)換橫波(TPS波)；在=50 ms時SS波出現(xiàn)，但由橫波產(chǎn)生的轉(zhuǎn)換縱波(SP波)能量較弱，不太明顯。從數(shù)值模擬結(jié)果可以看出，波場傳播特征符合波的運動學和動力學規(guī)律，證明了本文算法的正確性。

圖5為垂直界面模型模擬記錄。由圖5可知，徑向分量中的SS波、PS波能量強于軸向分量，而軸向分量中的PP波能量則強于徑向分量；軸向分量中炮檢同側(cè)接收的PP波能量弱于炮檢異側(cè)，而徑向分量中炮檢同側(cè)接收的PS波、SS波弱炮檢異側(cè)；界面一次反射縱波到達掌子面時發(fā)生反射，該反射波向前傳至界面處再次反射形成多次反射波(MV波)，當其到達檢波器的時間和界面一次反射波相近時就會與界面一次反射波相互疊加，若該多次反射波能量較強，就會嚴重影響有效波的識別，這在軸向分量的記錄中體現(xiàn)的十分明顯；同時，由單炮記錄可估算出P波和R波的視速度分別約為1 879 m/s和1 000 m/s，而PP波和SS波的旅行時分別約為48.5 ms和83.5 ms這與模型設(shè)計相符，進一步證明了數(shù)值模擬的正確性。

(a) t=20 ms；(b) t=35 ms；(c) t=50 ms

(a) X分量炮檢同側(cè)接收；(b) X分量炮檢異側(cè)接收；(c) R分量炮檢同側(cè)接收；(d) R分量炮檢異側(cè)接收

4.2 斷層破碎帶模型

圖6為斷層破碎帶模型模擬記錄。如圖6所示，與巖性分界面模型相比斷層破碎帶模型仍然具有軸向分量PP波能量較強，徑向分量PS波、SS波能量較強的特征；炮檢同側(cè)接收的記錄中頂、底界面的反射縱波(PP1，PP2)能量較強且同相軸較為清晰，可較好的區(qū)分破碎帶頂、底界面，但在一定程度上受波場疊加干擾，而炮檢波異側(cè)接收記錄波場疊加干擾嚴重，影響了有效波的識別；同時頂、底界面的反射橫波(SS1，SS2)受波場疊加干擾較小，利于頂、底界面的的識別與橫波成像；不同于單反射界面模型，斷層破碎帶模型在其頂?shù)捉缑婢鶗a(chǎn)生反射及透射，并且在頂?shù)捉缑嬷g還會發(fā)生多次反射，因此模擬記錄上波場比較復雜。

4.3 溶洞模型

圖7為溶洞模型模擬記錄。如圖7所示，對于上述溶洞作模型，溶洞頂部的反射縱波(PP1)在軸向分量及徑向分量的炮檢同側(cè)的接收記錄中有比較清晰的體現(xiàn)，但溶洞底部的反射縱波(PP2)由于受波場疊加影響，對溶洞底部的識別會有一定程度的干擾，徑向分量炮檢異側(cè)接收的反射縱波能量十分微弱，在記錄中幾乎無法體現(xiàn)；溶洞頂、底的反射橫波(SS1，SS2)在徑向分量的記錄中有較好的體現(xiàn)，利于橫波成像；由于溶洞內(nèi)部的低速填充物對地震波能量具有聚集作用[26]，所以溶洞底部的反射波PP2能量較溶洞頂部反射波PP1能量更強；不同于界面反射，溶洞作為局部地質(zhì)體，產(chǎn)生的繞射波會嚴重干擾有效波的識別。

(a) X分量炮檢同側(cè)接收；(b) X分量炮檢異側(cè)接收；(c) R分量炮檢同側(cè)接收；(d) R分量炮檢異側(cè)接收

(a) X分量炮檢同側(cè)接收；(b) X分量炮檢異側(cè)接收；(c) R分量炮檢同側(cè)接收；(d) R分量炮檢異側(cè)接收

5 結(jié)論

1) 考慮隧道空間的情況下，基于柱坐標系下的一階彈性波速度?應(yīng)力方程，采用高階交錯網(wǎng)格算法結(jié)合FCT技術(shù)，實現(xiàn)了柱坐標系下隧道空間典型不良地質(zhì)模型的高精度全波場數(shù)值模擬，系統(tǒng)分析了其傳播機理與波場特征，為隧道空間復雜波場的認識以及數(shù)據(jù)處理和解譯提供了有效的依據(jù)。

2) 巖性分界面模型，當界面傾角為90°時，掌子面及地質(zhì)界面間的多次反射波對有效波的識別會產(chǎn)生干擾，在數(shù)據(jù)處理和解譯時應(yīng)當加以區(qū)分。

3) 斷層破碎帶模型，炮檢同側(cè)接收的記錄中，頂、底界面的反射縱波能量較強且同相軸較為連續(xù)，可較好的區(qū)分頂、底界面，但多次波以及波場疊加干擾會影響有效波的識別。

4) 溶洞頂?shù)走吔绲姆瓷淇v波以及反射橫波在記錄中均有體現(xiàn)，但是溶洞底部反射縱波受波場疊加的影響，會對溶洞頂?shù)走吔绲淖R別造成一定程度的影響，而反射橫波能夠有效的區(qū)分溶洞頂?shù)走吔?，利于橫波成像。

[1] 張楊, 楊君, 周黎明, 等. TSP在隧道工程施工中的常見干擾和對巖體裂隙水及軟弱夾層等的預報研究[J]. 地球物理學進展, 2018, 33(2): 892?899. ZHANG Yang, YANG Jun, ZHOU Liming, et al. Common interference and the prediction of rock fissure water and weak interlayer in tunnel construction using TSP[J]. Progress in Geophysics, 2018, 33(2): 892?899.

[2] 周輪, 李術(shù)才, 許振浩, 等. 隧道綜合超前地質(zhì)預報技術(shù)及其工程應(yīng)用[J]. 山東大學學報(工學版), 2017, 47(2): 55?62. ZHOU Lun, LI Shucai, XU Zhenhao, et al. Integrated advanced geological prediction technology of tunnel and its engineering application[J]. Journal of Shandong University (Engineering Science), 2017, 47(2): 55?62.

[3] 李術(shù)才, 劉斌, 孫懷鳳, 等. 隧道施工超前地質(zhì)預報研究現(xiàn)狀及發(fā)展趨勢[J]. 巖石力學與工程學報, 2014, 33(6): 1090–1113. LI Shucai, LIU Bin, SUN Huaifeng, et al. State of art and trends of advanced geological prediction in tunnel construction[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(6): 1090?1113.

[4] 趙永貴. 國內(nèi)外隧道超前預報技術(shù)評析與推介[J]. 地球物理學進展, 2007, 22(4): 1344?1352. ZHAO Yonggui. Analysis and recommendation of tunnel prediction techniques at home and abroad[J]. Progress in Geophysics, 2007, 22(4): 1344?1352.

[5] ZHAO Yonggui, JIANG Hui, ZHAO Xiaopeng. Tunnel seismic tomography method for geological prediction and its application[J]. Applied Geophysics, 2006, 3(2): 69? 74.

[6] Dickmann T, Sander B. Drivage concurrent tunnel seismic prediction[J]. Felsbau-Rock and Soil Engineering, 1996(14): 406?411.

[7] Inazaki T, Isahai H, Kawamura S, et al. Stepwise application of horizontal seismic profiling for tunnel prediction ahead of the face[J]. The Leading Edge, 1999, 18(12): 1429?1431.

[8] Neil D W, Haramy K Y, Hanson D H, et al. Tomography to evaluate site condition during tunneling[J]. Geotechnical Speccial Publicaion, 1999(90): 71?281.

[9] 劉江平, 程飛, 范承余, 等. 基于隧道空間全波場二維數(shù)值模擬與特征分析[J]. 巖土工程學報, 2012, 34(9): 1705?1711. LIU Jiangping, CHENG Fei, FAN Chengyu, et al. Two- dimensional numerical simulation of tunnel-based seismic full-wave fields[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(9): 1705?1711.

[10] 魯光銀, 熊瑛, 朱自強. 隧道反射波超前探測有限差分正演模擬與偏移處理[J]. 中南大學學報(自然科學版), 2011, 42(1): 136?141. LU Guangyin, XIONG Ying, ZHU Ziqiang. Detection simulation ahead of tunnel face and reverse-time migration with reflection wave method[J]. Journal of Central South University (Science and Technology), 2011, 42(1): 136?141.

[11] CHENG Fei, LIU Jiangping, QU Niannian, et al. Two-dimensional pre-stack reverse time imaging based on tunnel space[J]. Journal of Applied Geophysics, 2014 (104): 106?113. DOI:10.1016/j.jappgeo.2014.02. 013.

[12] 凌飛, 肖宏躍, 朱夏樂, 等. 基于黏彈性介質(zhì)的隧道地震二維正演模擬[J]. 長江科學院院報, 2015, 32(5): 121?126. LING Fei, XIAO Hongyue, ZHU Xiale, et al. Two- dimensional forward seismic modeling for tunnel based on viscoelastic medium[J]. Journal of Yangtze River Scientific Research Institute, 2015, 32(5): 121?126.

[13] 宋杰. 隧道施工不良地質(zhì)三維地震波超前探測方法及其工程應(yīng)用[D]. 濟南: 山東大學, 2016. SONG Jie. The three-dimensional seismic ahead prospecting method and its application for adverse geology in tunnel construction[D]. Jinan: Shandong University, 2016.

[14] 王京. 隧道空間有限元地震波場模擬與逆時偏移成像[D]. 武漢: 中國地質(zhì)大學, 2017. WANG Jing. Finite element seismic wave field simulation and reverse time migration based on tunnel space[D]. Wuhan: China University of Geosciences, 2017.

[15] Jetschny S, Bohlen T, Kurzmann A. Seismic prediction of geological structures ahead of the tunnel using tunnel surface waves[J]. Geophysical Prospecting, 2011, 59(5): 934?946.

[16] YANG Sitong, WEI Jiuchuan, CHENG Jiulong, et al. Numerical simulations of full-wave fields and analysis of channel wave characteristics in 3-D coal mine roadway models[J]. Applied Geophysics, 2016, 13(4): 621?630.

[17] Harmankaya U, Kaslilar A, Wapenaar K, et al. Locating scatterers while drilling using seismic noise due to tunnel boring machine[J]. Journal of Applied Geophysics, 2018, 152: 86?99.

[18] LIU Yutao, LIU Jiangping, CHENG Fei, et al. Full wave field finite difference modeling of tunnel space in polar coordinates[C]// International Geophysical Conference, Qingdao, China, 17-20 April 2017, Qingdao, China. Society of Exploration Geophysicists and Chinese Petroleum Society, 2017.

[19] 查欣潔, 高星, 王偉, 等. 隧道工程勘察中的超前預報成像方法研究[J]. 地球物理學報, 2018, 61(3): 1151? 1156. ZHA Xinjie, GAG Xing, WANG Wei, et al. Advanced prediction migration method research in tunnel engineering investigation[J]. Chinese J Geophys, 2018, 61(3): 1151?1156.

[20] Hiroshi Takenaka, Hiroki Tanaka, Taro Okamoto, et al. Quasi-cylindrical 2.5D wave modeling for large-scale seismic surveys[J]. Geophysical Research Letters, 2003, 30(21):2086.DOI:10.1029/2003GL018068.

[21] Graves R W. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite difference[J]. Bulletin of the Seismological Society of America, 1996, 86(4): 1091?1106.

[22] LIU Qinghuo, Sinha B K. A 3D cylindrical PML/FDTD method for elastic waves in fluid-filled pressurized boreholes in triaxially stressed formations[J]. Geophysics, 2003, 68(5): 1731?1743.

[23] LIU Q H. Perfectly matched layers for elastic waves in cylindrical and spherical coordinates[J]. The Journal of the Acoustical Society of America, 1999, 105(4): 2075? 2084.

[24] Mittet R. Free-surface boundary conditions for elastic staggered-grid modeling schemes[J]. Geophysics, 2002, 67(5): 1616?1623.

[25] Bohlen T, Lorang U, Rabbel W, et al. Rayleigh-to-shear wave conversion at the tunnel face: From 3D-FD modeling to ahead-of-drill exploration[J]. Geophysics, 2007, 72(6): T67?T79.

[26] 董良國, 黃超, 劉玉柱, 等. 溶洞地震反射波特征數(shù)值模擬研究[J]. 石油物探, 2010, 49(2): 121?124, 15. DONG Liangguo, HUANG Chao, LIU Yuzhu, et al. Numerical simulation of seismic wave propagation in cave carbonate reservoir[J]. Geophysical Prospecting for Petroleum, 2010, 49(2): 121?124, 15.

Full wave filed numerical simulation and analysis of tunnel space based on cylindrical coordinate system

LU Guangyin, LUO Shuai, ZHU Ziqiang, SHI Keliang, XIA Chengzhi

(School of Geosciences and Info-physic, Central South University, Changsha 410083, China)

Differing from ground seismic exploration, the special structure of tunnel space makes its wave field very complex. Full-wave field numerical simulation based on actual tunnel model is an important method to effectively understand the propagation law of wave field in tunnel space. Previous studies on this field are mostly based on Cartesian coordinate system. As a special geological model similar to air column, numerical simulation of full wave field in cylindrical coordinate system is more meaningful for understanding the characteristics of wave field in tunnel space. An arbitrary even-order staggered grid difference scheme was deduced, artificial truncated boundary and tunnel free boundary were designed, and numerical simulation and wave field characteristics analysis were performed for typical problematic geological models of tunnels such as lithologic interface, fault fracture zone and karst cave based on the first-order elastic wave velocity-stress equation in cylindrical coordinates. The results show that the staggered grid finite difference method based on the first-order elastic wave velocity-stress equation in cylindrical coordinates combined with FCT technology can achieve high-precision numerical simulation of the full-wave field in tunnel space. The wave field characteristics are in accordance with the kinematic and dynamic characteristics of the wave, which provides an effective basis for understanding the complex wave field characteristics in tunnel space and interpreting the data.

advance prediction; tunnel space; cylindrical coordinates; full wave field; numerical simulation

P631.4

A

1672 ? 7029(2020)02 ? 0388 ? 08

10.19713/j.cnki.43?1423/u.T20190267

2019?04?08

國家自然科學基金資助項目(41974148)；湖南省安全生產(chǎn)監(jiān)督管理局資助項目(201907)

魯光銀(1976?)，男，湖北宜昌人，教授，博士，從事工程地球物理勘探研究；E?mail：13975894898@139.com

(編輯 蔣學東)