高 啟 朱 偉 郝 龍 鄧小虎 許 可 楊俊杰
(1.長江大學地球物理與石油資源學院,湖北 武漢 430100;2.武漢拓盟能源科技有限公司,湖北武漢 430074;3.武漢地震工程研究院有限公司,湖北 武漢 430071)
瑞雷面波勘探是工程物探中的常用方法。在城市、鐵路和公路附近,地表振動強烈,主動源面波勘探可能難以應(yīng)用,而被動源面波勘探利用場地背景噪聲中的面波,數(shù)據(jù)采集的環(huán)境要求低,是當前工程物探領(lǐng)域重點研究的方法[1-4]。
20 世紀50 年代,國外學者開始研究被動源面波記錄探測地層結(jié)構(gòu)[5]。在國內(nèi),被動源勘探研究始于20 世紀60 年代,2000 年以后應(yīng)用于探測地下橫波速度結(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ū)道路橋梁、隧道和邊坡施工而言,探測不同風化程度的花崗巖的分布具有重要意義。
在道路改擴建工程中,場地情況復雜,背景噪聲大,花崗巖風化程度的探測與評價是工程物探的難題。目前,被動源勘探在這方面的應(yīng)用相對較少。本研究根據(jù)道路工程勘探和建設(shè)施工中的實際情況,設(shè)計花崗巖不均勻風化模型,模擬生成被動源地震記錄,分析不同模型的地震記錄的差異。
在具有線狀分布的隨機強震動情況下,被動源勘探有可能采用與主動源勘探相似的線性觀測系統(tǒng),可用二維模型模擬被動源地震記錄。假設(shè)二維模型的上表面為地面,設(shè)置自由表面邊界條件[16];在其他邊界應(yīng)用完全匹配層邊界條件[17]。波動方程采用一階速度—應(yīng)力格式的彈性波方法,數(shù)值模擬方法為交錯網(wǎng)格有限差分法。數(shù)量眾多的隨機震源設(shè)置在地表附近。
一階速度—應(yīng)力格式的彈性波方程如公式(1)[18-19]。
式中:ρ為密度;vx、vz表示速度;σxx、σxz、σzz表示應(yīng)力;f x、f z表示體力量;λ、μ為拉梅系數(shù)。
模型邊界條件的設(shè)置如圖1 所示。假設(shè)模型的上表面為地表,采用自由表面邊界條件。其他三個邊界設(shè)置完全匹配層(PML)邊界條件。
圖1 邊界條件示意
1.2.1 吸收邊界條件。PML 邊界條件是當前彈性波數(shù)值模擬中應(yīng)用最多、效果最好的吸收邊界條件[20-21]。通過設(shè)置合理的參數(shù),它可以衰減絕大部分入射波的能量。應(yīng)用PML 邊界條件時,對波場變量沿坐標軸方向分解,在邊界層內(nèi)利用阻尼衰減相應(yīng)方向的波場變量的分量。阻尼系數(shù)的表達式為式(2)[22]。
式中:R為理論反射系數(shù);δ為吸收層厚度,取值越大吸收效果越好;x為吸收層網(wǎng)格點到邊界的距離。
1.2.2 自由表面邊界條件。自由表面邊界條件的方法有多種,目前模擬精度較好的是聲學彈性邊界近似法。其基本思想是利用聲學—彈性邊界近似代替自由表面,令正應(yīng)力σzz在自由表面處直接為零,同時還考慮自由表面上下橫向應(yīng)力保持連續(xù)的條件。在二維數(shù)值模擬時的處理方式為式(3)[21]。
式中:σzz、ρ分別為自由表面上的垂向正應(yīng)力、密度;λ、μ為拉梅常數(shù);ρ1和μ1分別為自由表面下介質(zhì)的密度和拉梅常數(shù)。
交錯網(wǎng)格有限差分法由Vireux 在進行SH 和PSV 波正演時建立,能穩(wěn)定地應(yīng)用在非均勻模型中,對流固邊界的適應(yīng)性也較好[23-24]。交錯網(wǎng)格有限差分法可以與PML 邊界條件和自由表面邊界條件相結(jié)合,實現(xiàn)在有限大小的模型中模擬瑞雷面波的傳播。交錯網(wǎng)格有限差分示意如圖2 所示,其主要特點是應(yīng)力和位移分布在網(wǎng)格的不同位置,且能滿足波場變量差分計算的要求。
圖2 交錯網(wǎng)格有限差分示意
被動源面波勘探常采用傅里葉變換法求取地震記錄的頻波譜,提取頻散曲線,反演橫波速度。這說明時間域的波形不是重要影響因素。隨機噪聲在地震記錄中的波形可能沒有特定的形態(tài)。因此,震源子波可以用雷克子波表示為式(4)。
式中:A為振幅;f m為主頻;td為延遲時間。
模擬被動源記錄,震源設(shè)置較為關(guān)鍵。震源是一組隨機震源的組合。各個震源的加載位置、振幅、主頻和延遲時間都應(yīng)具有一定的隨機性。
本研究建立了一個大小為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時刻垂直速度分量的波場快照
本研究構(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 三個不均勻風化花崗巖模型
被動源地震模擬需要大量震源。這些震源激發(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ù)直方圖
波場快照能夠展示地震波在模型中的傳播特征,對觀測系統(tǒng)的布置具有重要的指導意義。三個花崗巖模型在50 ms 時刻的垂直速度分量的波場快照如圖6所示。由圖6可知,在不同風化程度的花崗巖地層和巖體中地震波場具有明顯不同。這說明地層速度的變化,對波傳播具有重要影響。
圖6 三個模型在50 ms時刻的垂直速度分量的波場快照
在地表附近記錄速度垂直分量形成時間記錄。在數(shù)值模擬中,地表每一個網(wǎng)格點都可以是接收點,從而形成高密度的地震記錄。三個模型垂直速度分量的時間記錄如圖7所示,由于模型長度達3 000 m,記錄道數(shù)多,圖中按一定的間隔顯示了模型中部約2 000 m 范圍內(nèi)的記錄。由圖7 可知,不同的花崗巖地層模型,時間記錄存在差異,但在波形圖中并不明顯。三個模型記錄兩兩相減結(jié)果的波形如圖8 所示,表明三個時間記錄存在明顯差異。
圖7 三個模型垂直速度分量的時間記錄
圖8 三個模型時間記錄的兩兩差值
在我國東南部山地丘陵地區(qū),淺層風化花崗巖地質(zhì)體的勘探是道路工程建設(shè)中的重要工作內(nèi)容。在道路改擴建工程等背景噪聲強烈的地區(qū),被動源面波勘探可能是一種有效的方法。本研究構(gòu)建了三個不均勻風化的花崗巖地層二維模型,設(shè)計隨機震源,采用彈性波方程正演,獲得了被動源地震記錄和波場快照。當模型結(jié)構(gòu)發(fā)生變化時,波場快照和地震記錄均存在明顯的變化,表明被動源記錄對模型的響應(yīng)特征不同,能夠反映不同模型的特征。