路基下伏地質(zhì)缺陷彈性波數(shù)值模擬

李慶春，吳 華，周學明

1.長安大學地質(zhì)工程與測繪學院，西安 710054

2.長安大學理學院，西安 710064

3.鐵道第三勘察設計院集團有限公司，天津 300251

0 引言

我國西部地區(qū)，尤其是西北地區(qū)，存在特殊的自然地理環(huán)境和惡劣的地質(zhì)環(huán)境，地下地質(zhì)結(jié)構(gòu)十分復雜，地質(zhì)缺陷廣泛發(fā)育，主要有黃土洞穴（濕陷性暗穴、陷穴，地下采空區(qū)）、地面塌陷、活斷層、地裂縫、滑坡面、破碎帶及軟弱結(jié)構(gòu)層等。這些隱伏的地質(zhì)缺陷常常引發(fā)大量的地質(zhì)災害，給公路及鐵路的施工和安全運營造成極大威脅。

在路基下伏地質(zhì)缺陷探測方面，美國等西方發(fā)達國家以物探方法為主，而我國目前則以鉆探為主、物探為輔［1－2］。國內(nèi)近年來在利用地球物理勘探技術(shù)查明路基下伏地質(zhì)缺陷方面做了大量的工作，提出了多種方法：蔣正紅等［2］為了能較準確地揭示出路基病害的位置、規(guī)模及延展范圍等，提出應用地球物理層析成像技術(shù)對鐵路路基病害進行探測；孫忠弟［3］針對高等級公路下伏空洞勘探，利用電法、微重力、地震、地質(zhì)雷達等方法進行了較為系統(tǒng)的研究；童立元等［1］從采空區(qū)探測、穩(wěn)定性評價、治理及質(zhì)量監(jiān)控技術(shù)等方面對高速公路下伏采空區(qū)問題進行了綜合探討和分析，指出在采空區(qū)探測技術(shù)方面，有必要從理論上和實踐上來優(yōu)選最佳的探測技術(shù)組合，提高勘測精度；李慶春等［4］利用多道瞬態(tài)面波分析、淺層轉(zhuǎn)換波地震技術(shù)探測穿過公路的濕陷性黃土暗穴，取得了良好的應用效果；李志輝等［5］采用CT速度成像技術(shù)，對巖溶、洞穴、軟弱土層等病害存在部位進行定位，取得了良好的應用效果；韓永琦［6］、韓金河［7］采用地震折射層析法分別對青藏鐵路路基開裂和石膏巖溶進行了探測，取得了較好的應用效果；羅維斌等［8］利用高密度電法對黃土路基下伏地質(zhì)缺陷（地裂縫、黃土暗穴、斷層破碎帶、滑坡以及地下空洞等）進行模擬及實際探測，總結(jié)出針對不同缺陷的探測技術(shù)組合；張華等［9］將地震層析成像（CT）應用于公路采空區(qū)調(diào)查。

總之，在公路及鐵路建設中，地球物理探測已顯示并發(fā)揮了重要作用，在路面、路基、橋隧等探測檢測方面做過相應工作，使用了包括瑞利面波、聲波無損檢測、淺層反射地震、隧道地震預測、探地雷達、直流電法、自然電場法等諸多技術(shù)。但是針對西部黃土地區(qū)路基下伏地質(zhì)缺陷開展有效的綜合地球物理探測，目前尚缺乏實用成熟的整套技術(shù)。

典型路基下伏地質(zhì)缺陷的地球物理數(shù)值模擬是節(jié)約地球物理探測成本和驗證各種地球物理方法探測路基下伏地質(zhì)缺陷可行性和有效性的重要技術(shù)手段，但是，在路基下伏地質(zhì)缺陷數(shù)值模擬方面，前人的研究［10－11］主要集中在缺陷的變形、沉降、穩(wěn)定性和承載力模擬方面，關(guān)于模擬路基下伏地質(zhì)缺陷對彈性波響應方面的研究為數(shù)不多。周學明等［12］采用交錯網(wǎng)格高階有限差分數(shù)值模擬技術(shù)對一個包含有多種地質(zhì)缺陷的綜合地質(zhì)模型進行了模擬計算，通過模擬表明工程地震勘探在探測淺層斷層、裂縫、溶洞等地質(zhì)缺陷方面具有有效性。

筆者擬對西部黃土地區(qū)路基下伏地質(zhì)缺陷進行地球物理數(shù)值模擬，以驗證各種地球物理方法探測路基下伏地質(zhì)缺陷的可行性和有效性，以期為實際施工中地球物理方法技術(shù)組合以及觀測系統(tǒng)的優(yōu)化設計提供技術(shù)支持。

1 地震波場記錄有限差分數(shù)值模擬

數(shù)值模擬可模擬研究地震波在地下各種介質(zhì)中的傳播規(guī)律，是研究地下構(gòu)造和和巖性等各種地質(zhì)因素與地震波運動學和動力學響應特征之間關(guān)系的重要技術(shù)手段，對于認識和研究地震波在介質(zhì)中的傳播具有重要作用。地震數(shù)值模擬已廣泛應用在地震勘探的資料采集、處理和解釋的各個階段。筆者根據(jù)一階速度－應力彈性波方程，采用交錯網(wǎng)格高階差分方法，引入高效的完全匹配層（PML）技術(shù)解決邊界吸收問題，可進行黏彈介質(zhì)、雙相介質(zhì)和各向異性介質(zhì)的彈性波數(shù)值計算。

1.1 波動方程

根據(jù)彈性介質(zhì)位移、應力和應變之間的關(guān)系，可以推導出二維均勻各向同性介質(zhì)彈性波一階速度－應力方程［13－16］：

式中：vx、vz分別為質(zhì)點位移速度的水平和垂直分量；σx、σz分別為x和z方向的正應力；τzx為切應力分量；ρ為密度；t為時間；λ和μ為拉梅系數(shù)。

1.2 差分格式

1976 年，Madariaga［17］提出了一種較為先進的波動方程交錯網(wǎng)格有限差分解法，并首次利用交錯網(wǎng)格來模擬彈性介質(zhì)內(nèi)圓形擴展破裂產(chǎn)生的波動。該方法與常規(guī)差分格式相比，具有局部精度高、收斂速度快等優(yōu)點，且計算量和占用內(nèi)存與常規(guī)方法相當。交錯網(wǎng)格差分格式不僅要求空間網(wǎng)格交錯，而且要求在時間上也要交錯。交錯網(wǎng)格中速度和應力節(jié)點如圖1所示。

圖1 交錯網(wǎng)格中速度和應力節(jié)點示意圖Fig.1 Sketch map of velocity and stress node in staggered grid

時間上的高階差分會極大地增加計算量，在地震勘探所需的時間采樣上時間頻散并不明顯，增大時間步長可能會導致算法不穩(wěn)定。而空間頻散是數(shù)值頻散的主要部分，對空間導數(shù)需采用高階差分格式來有效抑制數(shù)值頻散。考慮到運算效率和精度，筆者選擇了時間二階、空間十階的差分格式［12］：

其中：

將式（2）、（3）代入式（1）即可得到彈性波方程的高階差分格式。

1.3 邊界條件

1）自由邊界條件

在水平自由表面上（z＝0），應力分量必須滿足σzz＝0，τzx＝0。由于采用的是交錯網(wǎng)格，所有的速度及應力分量不可能同時位于自由表面，設定vx，σxx，σzz位于自由表面。當采用高階差分格式時，還要求τzx，σzz在高階差分算子上一定是關(guān)于z＝0反對稱的，即虛像原理，并在自由界面上設置N層真空層，即縱、橫波速度分別為0（vP＝0，vS＝0）。則在計算速度分量時滿足：

其中：i表示x方向上的離散值；k表示時間上的離散值。

2）完全匹配層吸收邊界條件

采用完全匹配層吸收邊界條件吸收來自人工邊界的反射，如圖2所示。

將vx分解為x方向部分和z方向部分，則相應的匹配層方程［12］可寫為

其中，d（x）、d（z）采用 Coffino等［18］導出的如下衰減模型：

圖2 完全匹配層吸收邊界示意圖Fig.2 Sketch map of perfect matched layer absorbing boundary

式中：vPmax為最大縱波速度；δ為匹配層寬度；R為理想的邊界反射系數(shù)（這里取為10－6）；系數(shù)a＝0.25，b＝0.75。當d（x），d（z）不等于0時表示衰減；d（x），d（z）為0時表示不衰減。

利用完全匹配層作為吸收邊界的基本做法是在研究區(qū)域的四周引入完全匹配層，左右邊界和上下邊界采用完全匹配層吸收邊界。左右匹配層在x方向衰減，在z方向不衰減，即d（x）≠0、d（z）＝0。上下匹配層z方向衰減，x方向不衰減，即d（x）＝0，d（z）≠0。角點區(qū)域x方向和z方向都衰減，即d（x）≠0，d（z）≠0。模擬瑞利波場時，研究區(qū)域的上邊界取為自由邊界。

3）穩(wěn)定性條件

在模型計算中，空間網(wǎng)格尺寸和時間間隔的選擇至關(guān)重要，既要保證計算精度，還得兼顧計算效率。對于空間十階差分格式，保證整個計算區(qū)域能穩(wěn)定計算的高階交錯網(wǎng)格差分顯式的穩(wěn)定性條件為

特別地，對于Δx＝Δz的情況，穩(wěn)定性條件為

2 典型地質(zhì)缺陷地球物理數(shù)值模擬及分析

為考查不同地質(zhì)缺陷的地震波場響應特征，檢驗不同地震方法的探測能力，設計了一個綜合地質(zhì)地球物理模型，如圖3所示。該模型包含凹陷、斷層、地裂縫、地下低速體、軟弱夾層及滑坡等典型地質(zhì)缺陷。

模型長400m，深80m，共分為7層。在地面坐標50m處有一個隱伏地裂縫，裂縫寬0.25m，從地下2m開裂到地下40m。在地面坐標182m、地下20m處有一個寬和高均為4m的規(guī)則地下低速體。另外，在地面坐標224m、地下45m處有一個寬7 m、高10m的不規(guī)則低速體。利用彈性波場有限差分模擬程序，計算了多分量地震記錄和面波地震記錄各150炮，之后分別進行了P－P、P－SV波成像和高階面波反演橫波速度成像。

模擬計算時，模型物性參數(shù)縱橫波速度之比相對于實際淺層地質(zhì)情況略微偏小，但不影響各種物探方法成像效果的驗證。

2.1 反射地震波場模擬及成像

波場模擬時采用高階交錯網(wǎng)格有限差分和PML吸收邊界技術(shù)，模型長×寬為400m×80m，取網(wǎng)格間距Δx＝Δz＝0.25m，時間步長Δt＝0.1 ms，取主頻為80Hz的Ricker子波作為震源函數(shù)。合成的單炮記錄為縱波源激發(fā)，單邊接收，每炮48道，采樣率為0.5ms，采樣長度為512ms，每道1 024個采樣點，道間距為2m，炮間距為2m，最小偏移距為4m。

圖4為數(shù)值模擬產(chǎn)生的典型單炮地震記錄（前兩炮），由圖4可知，來自各層的反射可以清晰識別。

圖5和圖6分別是采用常規(guī)反射疊加處理技術(shù)產(chǎn)生的P－P、P－SV波成像剖面。由圖5和圖6可知，由于模型復雜，在水平疊加剖面上，還存在繞射現(xiàn)象，經(jīng)疊后偏移處理，成像剖面的效果大大改善。再經(jīng)時深轉(zhuǎn)換，反射層位與真實模型層位之間的最大誤差不超過2m。在進行P－SV波偏移時，采用等效偏移距偏移（EOM）方法，其中P－SV波的偏移速度模型是通過轉(zhuǎn)換波高精度非雙曲時距曲線“雙曲化”處理后，再按常規(guī)速度分析方法進行速度掃描分析得到的，具體方法原理見文獻［19］。

圖3 典型地質(zhì)缺陷綜合地質(zhì)地球物理模型Fig.3 A comprehensive geological and geophysical model with typical geological defects

圖4 前兩炮典型單炮反射地震記錄Fig.4 Typical reflection seismic records of the first two shots

真實模型中的凹陷、斷層、淺層低速體、軟弱夾層及滑坡等主要地質(zhì)缺陷在反射地震成像剖面上均能清晰成像，且對應關(guān)系良好。真實模型中的地裂縫由于尺度較小且近垂直展布，不能形成有效反射，在成像剖面中反映不明顯。同時，較深層的地下低速體在成像剖面中反映也不明顯，其原因可能是該低速體沒有有效的水平反射界面、埋深較深、側(cè)面傾角較陡，因此，在地表無法接收來自低速體的有效反射信息所致。

圖5 模型地震記錄P－P波成像剖面Fig.5 Section from P－P wave imaging of the model

2.2 瑞利波場模擬及成像

彈性波模擬時，上、下、左、右4個邊界均采用吸收邊界，因此模擬記錄中不出現(xiàn)面波（被邊界吸收）。面波模擬時，仍采用高階交錯網(wǎng)格有限差分技術(shù)，取上邊界為自由邊界，其他邊界為PLM吸收邊界。模型長×寬為400m×80m，取網(wǎng)格間距Δx＝Δz＝0.25m，時間步長Δt＝0.05ms，取主頻為25Hz的Ricker子波作為震源函數(shù)。合成的單炮記錄為縱波源激發(fā)，單邊接收，每炮48道，采樣率為0.5ms，采樣長度為512ms，每道1 024個采樣點，道間距為2m，炮間距為2m，最小偏移距為4m。

圖7為數(shù)值模擬產(chǎn)生的典型單炮面波地震記錄：x分量中波場成分簡單，面波發(fā)育；z分量中波場成分相對復雜，但面波能量仍占主要地位。

圖8是z分量記錄的頻散曲線提取結(jié)果與理論頻散曲線對比圖。由圖8可知，由波場記錄提取得到的頻散曲線與理論頻散曲線基本上是一致的，說明了面波記錄的正確性。二者存在微小誤差的原因是波場模擬時所建模型在橫向方向界面是變化的，而理論頻散曲線的計算是基于水平層狀介質(zhì)基礎上得到的。另外，由圖8可知，由模擬記錄提取得到的頻散曲線最低頻率成分為4Hz左右，對應的瑞利波相速度約為400m/s；因此瑞利波的最大波長（波長等于相速度/頻率）為100m左右。根據(jù)半波長原理，可估算出本文所給模型瑞利波的最大探測深度為40～50m。

圖6 模型地震記錄P－SV波成像剖面Fig.6 Section from P－SV wave imaging of the model

圖7 第一炮典型單炮面波地震記錄Fig.7 Typical surface wave records of the first shot

圖8 z分量記錄頻散曲線的提取結(jié)果與理論頻散曲線對比圖Fig.8 Extracted dispersion from z component and its theoretical dispersion curves

圖9是由z分量面波記錄提取的頻散曲線反演得到的橫波速度斷面圖。由圖9可知，模型中的凹陷、地裂縫、斷層、淺層低速體、軟弱夾層及滑坡等主要地質(zhì)缺陷均有反映，尤其是反射剖面無法成像的地裂縫在面波反演速度斷面中有明顯的反映。與反射地震成像相比，面波成像在淺層結(jié)構(gòu)探測中具有更好的縱、橫向分辨能力，分辨率可達1m左右，但面波的反演深度有限。因此，對于不同深度、不同規(guī)模的地質(zhì)體，應采用不同的地球物理探測方法，如利用多道瞬態(tài)面波解決淺層問題，用反射地震探測更深的目標，用多波地震解決裂縫等精細構(gòu)造問題，以達到對地質(zhì)缺陷進行立體探測與定位的效果。

3 結(jié)論與討論

1）典型路基下伏地質(zhì)缺陷的數(shù)值模擬，是考查不同地質(zhì)缺陷地震波場響應特征、檢驗不同地震方法的探測能力、驗證其可行性和有效性的重要技術(shù)手段。因此，在正式布置勘探之前，對地質(zhì)缺陷進行波場數(shù)值模擬是非常必要的，可以最大限度地降低探測成本。

2）常規(guī)反射地震P波和SV波成像技術(shù)對真實模型中的凹陷、斷層、淺層低速體、軟弱夾層及滑坡等主要地質(zhì)缺陷均能清晰成像，且對應關(guān)系良好，最大誤差不超過2m。真實模型中的地裂縫由于尺度較小且近垂直展布，不能形成有效反射，在成像剖面中反映不明顯。

3）與反射地震成像相比，面波成像在淺層結(jié)構(gòu)探測中具有更好的縱、橫向分辨能力，但探測深度有限。因此，對于不同深度、不同規(guī)模的地質(zhì)體，應采用不同的地球物理探測方法，如：利用多道瞬態(tài)面波解決淺層問題，用反射地震探測更深的目標，用多波地震解決裂縫等精細構(gòu)造問題，以達到對地質(zhì)缺陷進行立體探測與定位的效果。

4）數(shù)值模擬結(jié)果驗證了地球物理方法探測路基下伏地質(zhì)缺陷的可行性，表明不同的地球物理探測方法對地質(zhì)缺陷的探測效果和適用范圍不同。試驗結(jié)果對實際施工中地球物理方法技術(shù)組合以及觀測系統(tǒng)的優(yōu)化設計具有重要的指導意義。

圖9 面波反演橫波速度斷面圖Fig.9 Shear wave velocity section from surface wave inversion

（References）：

［1］童立元，劉松玉，邱鈺，等.高速公路下伏采空區(qū)問題國內(nèi)外研究現(xiàn)狀及進展［J］.巖石力學與工程學報，2004，23（7）：1198－1202.Tong Liyuan，Liu Songyu，Qiu Yu，et al.Current Research State of Problems Associated with Mined－Out Regions Under Expressway and Future Development［J］.Chinese Journal of Rock Mechanics and Engineering，2004，23（7）：1198－1202.

［2］蔣正紅，蔣偉宏.鐵路路基病害地球物理層析成像方法簡介［J］.路基工程，1996（3）：25－26.Jiang Zhenghong，Jiang Weihong.Introduction to Geophysical Tomography Methods on Rail－Bed Disasters［J］.Subgrade Engineering，1996（3）：25－26.

［3］孫忠弟.高等級公路下伏空洞勘探、危害程度評價及處治研究報告集［M］.北京：科學出版社，1999.Sun Zhongdi.Reports of Cavity Exploration Under High－Grade Highway， Hazard Evaluation and Treatment［M］.Beijing：Science Press，1999.

［4］李慶春，劉永華，馮兵，等.濕陷性黃土暗穴的綜合地球物理探測試驗研究［J］.公路交通科技，2005，22（6）：138－142.Li Qingchun，Liu Yonghua，F(xiàn)eng Bing，et al.Integrated Trial Studies on Geophysical Exploration for Hidden Loess Cavities［J］.Journal of Highway and Transportation Research and Development，2005，22（6）：138－142.

［5］李志輝，張立，劉爭平.CT技術(shù)在鐵路路基病害勘探中的應用［J］.物探化探計算技術(shù)，2006，28（2）：137－141.Li Zhihui，Zhang Li， Liu Zhengping.The CT Technique Application to the Survey of Rail－Bed Disasters［J］.Computing Techniques for Geophysical and Geochemical Exploration，2006，28（2）：137－141.

［6］韓永琦.地震折射層析法在青藏鐵路路基開裂調(diào)查中的應用［J］.工程地球物理學報，2009，6（1）：1－5.Han Yongqi.The Application of Investigating the Roadbed Cracks of the Qinghai－Tibet Railway Using the Seismic Refraction Tomography Method［J］.Chinese Journal of Engineering Geophysics，2009，6（1）：1－5.

［7］韓金河.地震折射層析法在青藏鐵路石膏巖溶勘察中的應用［J］.鐵道勘察，2009（4）：18－20.Han Jinhe.Application of Seismic Refraction Tomography Method in Geologic Investigation of Anhydrite Karst Region on Qinghai－Tibet Railway［J］.Railway Investigation and Surveying，2009（4）：18－20.

［8］羅維斌，李慶春.電法在公路下伏地質(zhì)缺陷精細探測中的應用［J］.工程勘察，2010（6）：78－85.Luo Weibin，Li Qingchun.Electrical Methods Applied for Road Failures Fine Detection［J］.Journal of Geotechnical Investigation ＆Surveying，2010（6）：78－85.

［9］張華，陳小宏，劉松.地震勘探技術(shù)在公路采空區(qū)調(diào)查中的應用［J］.地下空間與工程學報，2011，7（2）：371－374.Zhang Hua，Chen Xiaohong，Liu Song.Application of Seismic Prospecting Technology in Investigation on the Goaf Under Highway［J］.Chinese Journal of Underground Space and Engineering，2011，7（2）：371－374.

［10］王秋生.公路路基下伏空洞的數(shù)值分析［J］.巖土力學，2005，26（增刊）：215－218.Wang Qiusheng.Numerical Aanalysis of an Excavation Under Roadbed［J］.Rock and Soil Mechanics，2005，26（Sup.）：215－218.

［11］鄧楊.隱伏溶洞山區(qū)公路路堤穩(wěn)定性數(shù)值模擬分析與處治［J］.北方交通，2013（6）：23－26.Deng Yang.Numerical Simulation Analysis and Treatment on the Embankment Stability of Mountain Road Hidden Karst Caves［J］.Northern Communications，2013（6）：23－26.

［12］周學明，李慶春，李志華，等.地震波交錯網(wǎng)格高階差分數(shù)值模擬研究［J］.鐵道工程學報，2011（8）：1－12.Zhou Xueming，Li Qingchun，Li Zhihua，et al.Research on Staggered－Grid High－Order Finite－Difference Numerical Simulation of Seismic Wave［J］.Journal of Railway Engineering Society，2011（8）：1－12.

［13］Bohlen T，Saenger E H.Accuracy of Heterogeneous Staggered－Grid Finite－Difference Modeling of Rayleigh Waves［J］.Geophysics，2006，71（4）：109－115.

［14］周竹生，劉喜亮，熊孝雨.彈性介質(zhì)中瑞雷面波有限差分法正演模擬［J］.地球物理學報，2007，50（2）：567－573.Zhou Zhusheng，Liu Xiliang，Xiong Xiaoyu.Finite－Difference Modelling of Rayleigh Surface Wave in Elastic Media［J］.Chinese Journal of Geophysics，2007，50（2）：567－573.

［15］周曉華，陳祖斌，曾曉獻，等.交錯網(wǎng)格有限差分法模擬微動信號［J］.吉林大學學報：地球科學版，2012，42（3）：852－857.Zhou Xiaohua，Chen Zubin，Zeng Xiaoxian，et al.Simulation of Microtremor Using Staggered－Grid Finite Difference Method［J］.Journal of Jilin University：Earth Science Edition，2012，42（3）：852－857.

［16］楊慶節(jié)，劉財，耿美霞，等.交錯網(wǎng)格任意階導數(shù)有限差分格式及差分系數(shù)推導［J］.吉林大學學報：地球科學版，2014，44（1）：375－385.Yang Qingjie，Liu Cai，Geng Meixia，et al.Staggered Grid Finite Difference Scheme and Coefficients Deduction of Any Number of Derivatives［J］.Journal of Jilin University：Earth Science Edition，2014，44（1）：375－385.

［17］Madariaga R.Dynamics of an Expanding Circular Fault［J］.Bulletin of the Seismological Society of America，1976，66（3）：639－666.

［18］Collino F，Tsogka C.Application of the Perfectly Matched Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media［J］.Geophysics，2001（1）：294－307.

［19］沈鴻雁，李慶春.轉(zhuǎn)換波高精度非雙曲時距曲線“雙曲化”處理［J］.石油地球物理勘探，2007，42（3）：251－255.Shen Hongyan， Li Qingchun.High－Precision Hyperbolic Simplified Processing for Non－Hyperbolic Time－Distance Curve of Converted Wave［J］.Oil Geophysical Prospecting，2007，42（3）：251－255.