高 啟 朱 偉 郝 龍 鄧小虎 許 可 楊俊杰

（1.長江大學地球物理與石油資源學院，湖北 武漢 430100；2.武漢拓盟能源科技有限公司，湖北武漢 430074；3.武漢地震工程研究院有限公司，湖北 武漢 430071）

0 引言

瑞雷面波勘探是工程物探中的常用方法。在城市、鐵路和公路附近，地表振動強烈，主動源面波勘探可能難以應用，而被動源面波勘探利用場地背景噪聲中的面波，數(shù)據(jù)采集的環(huán)境要求低，是當前工程物探領(lǐng)域重點研究的方法［1-4］。

20 世紀50 年代，國外學者開始研究被動源面波記錄探測地層結(jié)構(gòu)［5］。在國內(nèi)，被動源勘探研究始于20 世紀60 年代，2000 年以后應用于探測地下橫波速度結(jié)構(gòu)［6-7］。李翀等［8］利用被動源面波勘探方法劃分南昌某地區(qū)地層結(jié)構(gòu)，確定了古河道的大致位置。劉國峰［9］利用被動源面波勘探對內(nèi)蒙古淺覆蓋區(qū)的礦區(qū)厚度進行劃分。邵廣周［10］采用主、被動源面波聯(lián)合勘探，準確探測黃土覆蓋區(qū)的分層結(jié)構(gòu)?；谀Ｐ偷恼⒎囱菅芯繉Ρ粍釉磾?shù)據(jù)的采集、處理、反演和解釋具有重要的指導作用，趙東［11］證明空間自相關(guān)方法提取頻散曲線與實際頻散曲線擬合度更高。白帥［12］對比分析不同速度模型的模擬結(jié)果，為被動源地震成像提供了有力的理論依據(jù)。

花崗巖在我國東南地區(qū)廣泛出露,與原巖相比，風化花崗巖的力學性質(zhì)變化較大，增加了市政和道路工程施工的難度。花崗巖的不均勻風化表現(xiàn)出不同的形式,孤石是一種常見的不均勻風化產(chǎn)物，對地鐵盾構(gòu)施工的影響較大［13-15］。對山區(qū)道路橋梁、隧道和邊坡施工而言，探測不同風化程度的花崗巖的分布具有重要意義。

在道路改擴建工程中，場地情況復雜，背景噪聲大，花崗巖風化程度的探測與評價是工程物探的難題。目前，被動源勘探在這方面的應用相對較少。本研究根據(jù)道路工程勘探和建設(shè)施工中的實際情況，設(shè)計花崗巖不均勻風化模型，模擬生成被動源地震記錄，分析不同模型的地震記錄的差異。

1 波動方程正演原理

在具有線狀分布的隨機強震動情況下，被動源勘探有可能采用與主動源勘探相似的線性觀測系統(tǒng)，可用二維模型模擬被動源地震記錄。假設(shè)二維模型的上表面為地面，設(shè)置自由表面邊界條件［16］；在其他邊界應用完全匹配層邊界條件［17］。波動方程采用一階速度—應力格式的彈性波方法，數(shù)值模擬方法為交錯網(wǎng)格有限差分法。數(shù)量眾多的隨機震源設(shè)置在地表附近。

1.1 彈性波方程

一階速度—應力格式的彈性波方程如公式（1）［18-19］。

式中：ρ為密度；vx、vz表示速度；σxx、σxz、σzz表示應力；f x、f z表示體力量；λ、μ為拉梅系數(shù)。

1.2 邊界條件

模型邊界條件的設(shè)置如圖1 所示。假設(shè)模型的上表面為地表，采用自由表面邊界條件。其他三個邊界設(shè)置完全匹配層（PML）邊界條件。

圖1 邊界條件示意

1.2.1 吸收邊界條件。PML 邊界條件是當前彈性波數(shù)值模擬中應用最多、效果最好的吸收邊界條件［20-21］。通過設(shè)置合理的參數(shù)，它可以衰減絕大部分入射波的能量。應用PML 邊界條件時，對波場變量沿坐標軸方向分解，在邊界層內(nèi)利用阻尼衰減相應方向的波場變量的分量。阻尼系數(shù)的表達式為式（2）［22］。

式中：R為理論反射系數(shù)；δ為吸收層厚度，取值越大吸收效果越好；x為吸收層網(wǎng)格點到邊界的距離。

1.2.2 自由表面邊界條件。自由表面邊界條件的方法有多種，目前模擬精度較好的是聲學彈性邊界近似法。其基本思想是利用聲學—彈性邊界近似代替自由表面，令正應力σzz在自由表面處直接為零，同時還考慮自由表面上下橫向應力保持連續(xù)的條件。在二維數(shù)值模擬時的處理方式為式（3）［21］。

式中：σzz、ρ分別為自由表面上的垂向正應力、密度；λ、μ為拉梅常數(shù)；ρ1和μ1分別為自由表面下介質(zhì)的密度和拉梅常數(shù)。

1.3 交錯網(wǎng)格有限差分法

交錯網(wǎng)格有限差分法由Vireux 在進行SH 和PSV 波正演時建立，能穩(wěn)定地應用在非均勻模型中，對流固邊界的適應性也較好［23-24］。交錯網(wǎng)格有限差分法可以與PML 邊界條件和自由表面邊界條件相結(jié)合，實現(xiàn)在有限大小的模型中模擬瑞雷面波的傳播。交錯網(wǎng)格有限差分示意如圖2 所示，其主要特點是應力和位移分布在網(wǎng)格的不同位置，且能滿足波場變量差分計算的要求。

圖2 交錯網(wǎng)格有限差分示意

1.4 震源

被動源面波勘探常采用傅里葉變換法求取地震記錄的頻波譜，提取頻散曲線，反演橫波速度。這說明時間域的波形不是重要影響因素。隨機噪聲在地震記錄中的波形可能沒有特定的形態(tài)。因此，震源子波可以用雷克子波表示為式（4）。

式中：A為振幅；f m為主頻；td為延遲時間。

模擬被動源記錄，震源設(shè)置較為關(guān)鍵。震源是一組隨機震源的組合。各個震源的加載位置、振幅、主頻和延遲時間都應具有一定的隨機性。

1.5 均勻模型測試

本研究建立了一個大小為1 000 m×600 m 的均勻模型測試模擬方法。地層的縱波速度為1 000 m/s，橫波速度為400 m/s，密度為2 000 kg/m3。網(wǎng)格大小0.6 m，時間步長0.08 ms。單一震源設(shè)置在模型中部地表附近。50 ms時刻垂直速度分量的波場快照如圖3所示，可以清楚地觀測到瑞雷面波、橫波和縱波。

圖3 50 ms時刻垂直速度分量的波場快照

2 不均勻風化花崗巖地層模型的被動源面波正演與分析

2.1 地層模型

本研究構(gòu)建了三個不均勻風化花崗巖地層模型，如圖4 所示。三個模型的大小均為3 000 m×600 m，網(wǎng)格大小均為0.6 m?；◢弾r的風化程度分為全風化、強風化和中風化三類。模型一［圖4（a）］的第一層為全風化層，厚度為50 m；第二層為強風化層，厚度100 m；第三層為中風化層，厚度為450 m。模型二［圖4（b）］的第一層為全風化層，厚度為20 m；第二層為強風化層，厚度為80 m；第三層為中風化層，厚度為500 m。模型三［圖4（c）］的第一層為全風化層，厚度為20 m；第二層為強風化層，厚度為480 m；在強風化層上部殘留一個球形的中風化花崗巖體，垂向厚度約為200 m，頂部與全風化層相連。各層花崗巖的縱波速度、橫波速度和密度見表1。在工程地震勘探中，檢波器排列長度一般較短。本研究模型的長度遠遠大于排列長度，主要是為了布置數(shù)量眾多的震源。

表1 不同風化程度花崗巖的縱波速度、橫波速度和密度

圖4 三個不均勻風化花崗巖模型

2.2 隨機震源

被動源地震模擬需要大量震源。這些震源激發(fā)的空間位置，頻率、出現(xiàn)時間和振幅須滿足一定的隨機性。在有限的空間和時間內(nèi)模擬地震波傳播并記錄，必須對這些參數(shù)的分布范圍進行控制。5 000 個隨機震源的參數(shù)統(tǒng)計結(jié)果如圖5所示。由圖5可知，震源的激發(fā)位置在水平方向分布較均勻，在垂直方向分布在地表以下50 m 范圍內(nèi)；震源主頻相對較低；延遲時間分布范圍約為250～750 ms。這5 000 個隨機震源在模擬之前確定，因此在三個模型的模擬中，震源參數(shù)是一致的。

圖5 震源參數(shù)直方圖

2.3 波場快照

波場快照能夠展示地震波在模型中的傳播特征，對觀測系統(tǒng)的布置具有重要的指導意義。三個花崗巖模型在50 ms 時刻的垂直速度分量的波場快照如圖6所示。由圖6可知，在不同風化程度的花崗巖地層和巖體中地震波場具有明顯不同。這說明地層速度的變化，對波傳播具有重要影響。

圖6 三個模型在50 ms時刻的垂直速度分量的波場快照

2.4 時間記錄

在地表附近記錄速度垂直分量形成時間記錄。在數(shù)值模擬中，地表每一個網(wǎng)格點都可以是接收點，從而形成高密度的地震記錄。三個模型垂直速度分量的時間記錄如圖7所示，由于模型長度達3 000 m，記錄道數(shù)多，圖中按一定的間隔顯示了模型中部約2 000 m 范圍內(nèi)的記錄。由圖7 可知，不同的花崗巖地層模型，時間記錄存在差異，但在波形圖中并不明顯。三個模型記錄兩兩相減結(jié)果的波形如圖8 所示，表明三個時間記錄存在明顯差異。

圖7 三個模型垂直速度分量的時間記錄

圖8 三個模型時間記錄的兩兩差值

3 結(jié)語

在我國東南部山地丘陵地區(qū)，淺層風化花崗巖地質(zhì)體的勘探是道路工程建設(shè)中的重要工作內(nèi)容。在道路改擴建工程等背景噪聲強烈的地區(qū)，被動源面波勘探可能是一種有效的方法。本研究構(gòu)建了三個不均勻風化的花崗巖地層二維模型，設(shè)計隨機震源，采用彈性波方程正演，獲得了被動源地震記錄和波場快照。當模型結(jié)構(gòu)發(fā)生變化時，波場快照和地震記錄均存在明顯的變化，表明被動源記錄對模型的響應特征不同，能夠反映不同模型的特征。