霍科宇,邵廣周*,王國(guó)順,白帥,吳華
(1.長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院,陜西西安 710054; 2.長(zhǎng)安大學(xué)理學(xué)院,陜西西安 710054)
傳統(tǒng)的瑞雷面波勘探技術(shù)常以人工激發(fā)方式產(chǎn)生的振動(dòng)作為震源,通過(guò)對(duì)地面記錄的信號(hào)進(jìn)行分析,提取面波頻散曲線并反演探區(qū)橫波速度,以實(shí)現(xiàn)對(duì)探區(qū)的淺地表結(jié)構(gòu)成像。但因該技術(shù)探測(cè)深度較小,故在城市環(huán)境等外界噪聲干擾大的地區(qū)施工受到諸多條件限制。而以人類生產(chǎn)活動(dòng)引起的振動(dòng)與自然界海嘯、地震、火山噴發(fā)等天然振動(dòng)作為震源的被動(dòng)源面波勘探技術(shù),能夠探測(cè)到更深的地層,且具有野外施工方式靈活、數(shù)據(jù)采集成本低等優(yōu)點(diǎn),所受關(guān)注越來(lái)越多。
被動(dòng)源成像技術(shù)的研究和應(yīng)用主要集中在數(shù)據(jù)采集方式、頻散曲線提取方法和頻散曲線反演成像等方面[1-3]。隨著研究的不斷深入,被動(dòng)源波場(chǎng)模擬研究逐漸得到充分重視。劉軍[4]利用射線追蹤方法進(jìn)行了二維微地震多波場(chǎng)正演模擬研究,但只考慮了震源隨位置的變化,而未考慮震源隨時(shí)間的變化;容嬌君等[5]采用延時(shí)爆炸反射點(diǎn)代替微地震震源,采用單程波動(dòng)方程相移結(jié)合插值的波場(chǎng)延拓方法對(duì)設(shè)計(jì)的地質(zhì)模型進(jìn)行了數(shù)值模擬,但只考慮了二維情況下單分量接收和震源隨時(shí)間的變化; Zhao等[6]對(duì)流體注入儲(chǔ)層引發(fā)的微地震進(jìn)行了數(shù)值模擬,得到了與實(shí)際觀測(cè)數(shù)據(jù)一致的結(jié)果;朱恒等[7]基于地震波有限差分?jǐn)?shù)值模擬得到虛炮集記錄的疊加剖面,同時(shí)驗(yàn)證了虛炮集記錄的準(zhǔn)確性,證明了地震干涉技術(shù)的可行性;張喚蘭[8]根據(jù)水力裂縫特征建立微地震震源分布的數(shù)學(xué)模型,利用三維高階交錯(cuò)網(wǎng)格有限差分法對(duì)微地震波場(chǎng)進(jìn)行了數(shù)值模擬,實(shí)現(xiàn)了微地震的定位,但模擬過(guò)程中壓力震源需滿足正態(tài)分布;唐杰等[9]通過(guò)有限差分正演模擬探討了微地震波場(chǎng)的傳播特性,結(jié)果表明微地震震源特性對(duì)波形有顯著影響;徐宗博[10]采用面波理論公式模擬背景噪聲,實(shí)現(xiàn)了背景噪聲波場(chǎng)模擬,并將該方法的模擬結(jié)果與有限差分模擬結(jié)果對(duì)比,發(fā)現(xiàn)有限差分模擬結(jié)果細(xì)節(jié)信息更豐富,更符合實(shí)際情況;Wang等[11]對(duì)地下儲(chǔ)氣庫(kù)進(jìn)行微震檢測(cè)過(guò)程中進(jìn)行了微地震波傳播的數(shù)值模擬,并根據(jù)模擬的結(jié)果對(duì)微震事件進(jìn)行了分析;雷朝陽(yáng)等[12]在被動(dòng)源成像研究中,分別將震源均勻分布在地下某一深度和隨機(jī)分布在400~500 m 范圍內(nèi),采用聲波有限差分進(jìn)行波場(chǎng)數(shù)值模擬,同時(shí)利用隨機(jī)時(shí)間序列與Ricker子波褶積構(gòu)造噪聲源合成噪聲記錄,證實(shí)了從透射響應(yīng)中提取反射響應(yīng)的有效性;張晟瑞[13]采用兩種不同的震源加載方式,同時(shí)依據(jù)水力壓裂誘發(fā)微地震的震源激發(fā)方式不同,引入三種不同震源機(jī)制進(jìn)行微地震數(shù)值模擬,對(duì)正演地震記錄和不同時(shí)刻波場(chǎng)快照進(jìn)行對(duì)比分析,得到了不同震源機(jī)制情況下的波場(chǎng)特征。
總體而言,無(wú)論是基于已知模型對(duì)背景噪聲成像方法以及頻散曲線提取方法的有效性進(jìn)行驗(yàn)證,還是對(duì)背景噪聲數(shù)據(jù)采集方式進(jìn)行參數(shù)優(yōu)化,實(shí)現(xiàn)背景噪聲波場(chǎng)模擬均顯得尤為重要。本文通過(guò)將隨機(jī)震源加載于不同地質(zhì)模型上,模擬計(jì)算檢波器的疊加響應(yīng),從合成記錄中提取頻散曲線,并與Schwab-Knopoff快速計(jì)算方法[14]得出的理論頻散曲線進(jìn)行對(duì)比、研究,驗(yàn)證了本文三維背景噪聲模擬方法的可行性和有效性。
三維背景噪聲波場(chǎng)模擬過(guò)程中,隨機(jī)震源子波類型的產(chǎn)生和震源的加載都是必不可少的。常見(jiàn)子波類型有Ricker 子波、Fuchs-Muller 子波、正弦衰減子波和Gaussian導(dǎo)數(shù)子波等,四種子波表達(dá)式分別為
式中:f0為子波中心頻率;td為延遲時(shí)間;t為時(shí)間。延遲時(shí)間可調(diào)整子波起跳時(shí)間,從而控制隨機(jī)震源起震時(shí)刻[15]。首先,根據(jù)模擬區(qū)域的大小確定隨機(jī)震源的產(chǎn)生范圍,根據(jù)背景噪聲的特點(diǎn)確定隨機(jī)震源個(gè)數(shù)、空間坐標(biāo)、延遲時(shí)間、主頻、能量大小、子波類型(1 代表Ricker 子波,2 代表Fuchs-Muller 子波,3 代表正弦子波,4代表Gaussian導(dǎo)數(shù)子波)和震源加載方向(2代表X軸方向加載集中力源,3 代表Y軸方向加載集中力源,4 代表Z軸方向加載集中力源,5 代表自定義集中力源加載方向(自定義力源加載可通過(guò)設(shè)置力源角度α和β實(shí)現(xiàn)任意方向的源加載,其中α表示X軸偏向Z軸的偏向角,β表示Y軸偏向Z 軸的偏向角))的范圍,并設(shè)置背景噪聲波場(chǎng)模擬的相關(guān)參數(shù)、模型參數(shù)及檢波器的排列等,其中Y軸為垂直方向。然后,由子程序產(chǎn)生隨機(jī)震源的個(gè)數(shù)、空間坐標(biāo)、延遲時(shí)間、主頻、能量大小、子波類型和震源加載方向等參數(shù)。表1為隨機(jī)生成的9個(gè)震源的震源參數(shù)。
表1 隨機(jī)震源參數(shù)
隨機(jī)震源的產(chǎn)生是波場(chǎng)模擬的前提,同時(shí),隨機(jī)震源的加載在整個(gè)背景噪聲模擬過(guò)程中亦是至關(guān)重要的,只有在模擬的過(guò)程中有效加載了震源才能得到地震波在介質(zhì)中的傳播情況。本文采用的三維均勻介質(zhì)彈性波應(yīng)力—速度方程為
式中:ρ為介質(zhì)密度;Vi為速度分量;τij為應(yīng)力分量;fi為外力分量;λ、μ為介質(zhì)的拉梅常數(shù); 下標(biāo)x、y、z為坐標(biāo)方向。
對(duì)三維均勻介質(zhì)彈性波速度—應(yīng)力方程采用交錯(cuò)網(wǎng)格有限差分[16-18]進(jìn)行模擬。首先將密度、速度和應(yīng)力分量等參數(shù)交錯(cuò)布置與模型網(wǎng)格點(diǎn)或半網(wǎng)格點(diǎn)上;其次,將式(4)所有一階偏導(dǎo)數(shù)進(jìn)行差分近似,根據(jù)計(jì)算需要得到相應(yīng)階數(shù)的差分方程;最后,將時(shí)間變量離散,隨著時(shí)間的步進(jìn),根據(jù)差分方程循環(huán)交替地計(jì)算各網(wǎng)格點(diǎn)的應(yīng)力或速度,最終得到各網(wǎng)格點(diǎn)不同時(shí)刻的變量值。本文所采用的交錯(cuò)網(wǎng)格有限差分的優(yōu)點(diǎn)之一就是分開(kāi)計(jì)算速度和應(yīng)力,從而可根據(jù)不同數(shù)據(jù)集以不同方式便利地加載震源。
在地震勘探中,常見(jiàn)震源有集中、爆炸和剪切震源[19]。本文采用的差分算法可將震源加載到應(yīng)力和速度分量中,通過(guò)這兩種方法的組合能夠構(gòu)建各種震源。將震源加載到應(yīng)力τxx、τxy網(wǎng)格點(diǎn)的時(shí)間二階、空間四階差分形式如下
式中:空間步長(zhǎng)為Δx=Δy=Δz=Δh;時(shí)間步長(zhǎng)為Δt;st+1為t+1時(shí)刻的震源;為差分系數(shù);上標(biāo)t+1 為時(shí)間前進(jìn)一個(gè)步長(zhǎng);下標(biāo)分別為空間上x(chóng)、y、z方向前進(jìn)半個(gè)步長(zhǎng)。
將震源代入不同速度網(wǎng)格點(diǎn)可實(shí)現(xiàn)不同方向集中力源的加載過(guò)程,將震源加載到速度分量Vx網(wǎng)格點(diǎn)的時(shí)間二階、空間四階差分形式如下
同理,可得到震源加載到正應(yīng)力(τyy,τzz)、剪切應(yīng)力(τxz,τyz,)和速度分量(Vy,Vz)網(wǎng)格點(diǎn)上的差分形式。此外,為避免出現(xiàn)數(shù)值頻散以及提高數(shù)值模擬精度,網(wǎng)格間距需滿足下式
式中:λmin為最小波長(zhǎng);Vmin為模型最小波速;fmax為源信號(hào)最大頻率;n為每個(gè)最小波長(zhǎng)所含網(wǎng)格點(diǎn)數(shù),其值與差分算子階數(shù)有關(guān)。
同時(shí),為避免數(shù)值模擬過(guò)程中發(fā)生不穩(wěn)定,需進(jìn)行穩(wěn)定性判別,判別條件[20]采用
式中:Vmax為模型最大波速;βo為對(duì)應(yīng)差分算子階數(shù)的差分系數(shù),下標(biāo)o為差分算子階數(shù);q為穩(wěn)定性條件因子,其值與差分算子階數(shù)有關(guān)。
本文背景噪聲波場(chǎng)模擬在已知模型介質(zhì)結(jié)構(gòu)和參數(shù)的前提下實(shí)現(xiàn)。首先將產(chǎn)生的隨機(jī)震源加載于模型的交錯(cuò)網(wǎng)格系統(tǒng)之上,利用交錯(cuò)差分計(jì)算并研究地震波在地下介質(zhì)中的傳播規(guī)律; 然后提取各檢波器的疊加響應(yīng)得到頻散曲線。模擬得到的信號(hào)為無(wú)道頭微機(jī)格式,能夠作為二進(jìn)制文件直接進(jìn)行讀取,跳過(guò)文件卷頭和道頭,即為信號(hào)存儲(chǔ)的數(shù)據(jù)文件。通常在提取頻散曲線前需對(duì)模擬數(shù)據(jù)進(jìn)行如濾波、臺(tái)站道間距計(jì)算、時(shí)間域數(shù)據(jù)標(biāo)準(zhǔn)化和頻譜白化等預(yù)處理。整個(gè)隨機(jī)震源數(shù)值模擬及頻散曲線提取過(guò)程如圖1所示。
圖1 隨機(jī)震源數(shù)值模擬及頻散曲線提取流程
本文通過(guò)三組理論模型的被動(dòng)源數(shù)值模擬數(shù)據(jù)提取頻散曲線與理論頻散曲線或頻散能量圖進(jìn)行對(duì)比,以驗(yàn)證本文隨機(jī)震源的產(chǎn)生與加載在被動(dòng)源噪聲波場(chǎng)模擬中的可行性和有效性。
高速雙層模型尺寸為3120 m×3120 m×3120 m,具體參數(shù)如表2所示。
模型網(wǎng)格間距Δx=Δy=Δz=12 m; 震源個(gè)數(shù)為550; 震源強(qiáng)度在0.001~1.000 之間隨機(jī)生成; 子波延遲時(shí)間在0~4.0 s之間隨機(jī)生成; 子波類型通過(guò)隨機(jī)生成; 震源主頻范圍為0.5~20 Hz; 震源加載方向隨機(jī)生成。圖2和圖3分別為高速雙層模型隨機(jī)震源空間分布圖和隨機(jī)震源的相關(guān)參數(shù)柱狀圖。
圖2 高速雙層模型隨機(jī)震源空間分布圖
圖3 高速雙層模型隨機(jī)震源的相關(guān)參數(shù)柱狀圖
沿地表x方向布設(shè)221 個(gè)線性排列的檢波器,間距為12 m。采樣率為0.2 ms,采樣時(shí)間為4.0 s。模型四周采用吸收層厚度為20 個(gè)網(wǎng)格點(diǎn)的CPML 吸收邊界,地表采用鏡像法處理的自由邊界[19],正演得到有效的合成地震記錄(圖4)。
圖4 高速雙層模型噪聲合成記錄(每10 道顯示1 道)
模型尺寸為200 m×100 m×100 m,具體參數(shù)如表3所示。
表3 三層含低速層模型參數(shù)
模型網(wǎng)格間距Δx=Δy=Δz=0.5 m;震源個(gè)數(shù)為729; 震源強(qiáng)度在0.001~1.000 之間隨機(jī)生成;子波延遲時(shí)間在0~4.0 s 之間隨機(jī)生成;子波類型隨機(jī)生成;震源主頻范圍為0.5~30 Hz;震源加載方向隨機(jī)生成。在此參數(shù)下,既滿足背景噪聲的實(shí)際情況,也盡可能地避免了數(shù)值頻散。圖5 和圖6 分別為三層含低速層模型隨機(jī)震源空間分布圖和隨機(jī)震源的相關(guān)參數(shù)柱狀圖。
圖5 三層含低速層模型隨機(jī)震源空間分布圖
圖6 三層含低速層模型隨機(jī)震源的相關(guān)參數(shù)柱狀圖
沿地表x方向布設(shè)86個(gè)線性排列的檢波器,間距為2 m,采樣率為0.2 ms,采樣時(shí)間4.0 s。模型四周采用吸收層厚度為20 個(gè)網(wǎng)格點(diǎn)的CPML 吸收邊界,地表采用鏡像法處理的自由邊界。正演得到有效的合成地震記錄見(jiàn)圖7。
圖7 三層含低速層模型噪聲合成記錄(每5 道顯示1 道)
模型尺寸為200 m×100 m×100 m,具體參數(shù)如表4所示。
表4 垂直斷層模型參數(shù)
垂直斷層斷面位于x方向80 m 處,模型結(jié)構(gòu)如圖8所示。模型網(wǎng)格間距Δx=Δy=Δz=0.5 m;震源個(gè)數(shù)為500;震源強(qiáng)度在0.001~1.000 之間隨機(jī)生成;子波延遲時(shí)間在-0.2~12.0 s 之間隨機(jī)生成;子波類型隨機(jī)生成;震源主頻范圍為0.5~30.0 Hz?;谝陨蠗l件進(jìn)行模擬,可以得到符合實(shí)際背景噪聲的地震記錄,且有效防止了數(shù)值頻散的出現(xiàn)。
圖8 垂直斷層模型
圖9 和圖10 分別為垂直斷層模型隨機(jī)震源空間分布圖和隨機(jī)震源的相關(guān)參數(shù)柱狀圖。
圖9 斷層模型隨機(jī)震源空間分布圖
圖10 垂直斷層模型隨機(jī)震源相關(guān)參數(shù)柱狀圖
沿地表x方向線性排列的檢波器有86個(gè),道間距為2 m,采樣率為0.2 ms,采樣時(shí)間為1.0 s。模型四周采用吸收層厚度為20 個(gè)網(wǎng)格點(diǎn)的CPML 吸收邊界,地表采用鏡像法處理的自由邊界[21]。正演得到有效的合成地震記錄見(jiàn)圖11。
為驗(yàn)證三維背景噪聲模擬的可行性和有效性,本文通過(guò)理論合成噪聲記錄提取面波頻散曲線,并與頻散方程得到的理論頻散曲線進(jìn)行對(duì)比、研究。
在提取面波頻散曲線之前,需要對(duì)數(shù)值模擬生成的被動(dòng)源地震記錄進(jìn)行預(yù)處理。首先,計(jì)算各道之間的距離,為隨后的互相關(guān)計(jì)算提供參數(shù)。然后,依據(jù)不同設(shè)計(jì)地質(zhì)模型采用不同通帶、阻帶頻段進(jìn)行帶通濾波,并對(duì)濾波后的地震記錄采用滑動(dòng)絕對(duì)平均法進(jìn)行時(shí)間域標(biāo)準(zhǔn)化,以消除采集臺(tái)站周圍異常信號(hào)的影響。經(jīng)多次試驗(yàn),發(fā)現(xiàn)時(shí)間窗口長(zhǎng)度為40 s時(shí)得到的數(shù)據(jù)結(jié)果最佳。再次,經(jīng)過(guò)傅里葉變換將信號(hào)從時(shí)間域轉(zhuǎn)換到頻率域,對(duì)時(shí)間域標(biāo)準(zhǔn)化后的數(shù)據(jù)進(jìn)行頻譜白化。頻譜白化的目的是使頻譜曲線更均勻,避免之后計(jì)算出現(xiàn)“譜孔現(xiàn)象”。最后,利用傅里葉逆變換將頻譜白化后的數(shù)據(jù)再次轉(zhuǎn)換到時(shí)間域。
本文提取頻散曲線采用基于Aki 公式的兩道法[22],該方法基于互相關(guān)理論,利用互相關(guān)譜和第一類零階貝塞爾函數(shù)進(jìn)行擬合。而互相關(guān)處理是對(duì)預(yù)處理的頻率域的實(shí)部進(jìn)行負(fù)半軸的反序和正半軸的疊加,再進(jìn)行對(duì)稱化處理的過(guò)程。
在高速雙層模型的互相關(guān)處理過(guò)程中,選取第1道與第10 道的隨機(jī)震源合成噪聲記錄進(jìn)行相關(guān)處理和傅里葉變換,得到互相關(guān)結(jié)果和頻率域互相關(guān)譜。同時(shí),為了驗(yàn)證三維背景噪聲波場(chǎng)模擬在頻散曲線提取效果上的優(yōu)勢(shì),在相關(guān)參數(shù)保持不變的前提下,本文基于二維和三維波動(dòng)方程模擬背景噪聲,并選取第1 道與第55 道的二維和三維隨機(jī)震源合成噪聲記錄進(jìn)行相關(guān)處理。為驗(yàn)證三維背景噪聲波場(chǎng)模擬方法在復(fù)雜模型的適用性,使用垂直斷層模型并選取其第1 道與第78 道的三維隨機(jī)震源合成噪聲記錄進(jìn)行相關(guān)處理。
兩道法[22]提取頻散曲線的基本原理是對(duì)數(shù)值模擬數(shù)據(jù)經(jīng)預(yù)處理和互相關(guān)處理后,得到的互相關(guān)譜實(shí)部,進(jìn)行三次多項(xiàng)式插值獲得相應(yīng)零點(diǎn)fi,聯(lián)合第一類零階貝塞爾函數(shù)的零點(diǎn)zi,代入c(wi)=2πfir/zi計(jì)算出對(duì)應(yīng)的相速度c,其中wi為互相關(guān)譜的第i個(gè)零點(diǎn),r為臺(tái)站間距。由于計(jì)算得到的互相關(guān)譜可能存在環(huán)境干擾,使互相關(guān)譜實(shí)部在某一頻段上存在缺少或增加零點(diǎn),因此用式c(wi)=2πfir/zi+2m替代上式,m為缺失或多余的零點(diǎn)個(gè)數(shù),其取值為0,1,-1,2,-2,…??衫没ハ嚓P(guān)函數(shù)計(jì)算不同m值擬合出不同頻散曲線,然后再依據(jù)設(shè)定模型橫波速度范圍選擇一條合適的頻散曲線[23-25]。
3.3.1 高速雙層模型
在高速雙層模型頻散曲線提取時(shí),發(fā)現(xiàn)m=0 時(shí)頻散曲線提取效果最佳,最低頻率可以提取到1.9 Hz左右。圖12即為m=0時(shí),模擬數(shù)據(jù)利用兩道法提取的頻散曲線與理論頻散曲線對(duì)比圖。由圖可以看出,利用兩道法計(jì)算的頻散曲線(離散點(diǎn))與理論頻散曲線擬合較好,只有零星的離散點(diǎn)偏離理論頻散曲線。
圖12 理論頻散曲線與合成噪聲記錄提取頻散曲線對(duì)比
3.3.2 三層含低速層模型
三層含低速層模型頻散曲線提取過(guò)程中,發(fā)現(xiàn)m=0時(shí)頻散曲線提取效果同樣最為理想。為了驗(yàn)證三維背景噪聲波場(chǎng)模擬的優(yōu)勢(shì),本文對(duì)三維模擬數(shù)據(jù)和二維模擬數(shù)據(jù)分別利用兩道法提取頻散曲線并進(jìn)行比較。圖13 即為m=0 時(shí)三維和二維模擬數(shù)據(jù)對(duì)應(yīng)的頻散曲線與理論頻散曲線對(duì)比圖。由圖可知,利用兩道法計(jì)算的頻散曲線(離散點(diǎn))能較好地?cái)M合理論頻散曲線,三維模擬數(shù)據(jù)提取的頻散曲線在低頻段提取精度高于二維模擬數(shù)據(jù),且高頻段對(duì)應(yīng)的頻散曲線更光滑。
圖13 理論頻散曲線與二維(a)、三維(b)波動(dòng)方程合成噪聲記錄提取的頻散曲線對(duì)比
3.3.3 垂直斷層模型
由于理論頻散曲線的計(jì)算模型局限于層狀介質(zhì),無(wú)法計(jì)算斷層模型的理論頻散曲線,故設(shè)置一個(gè)相同斷層模型使用主動(dòng)源方法作為對(duì)照,主動(dòng)源方法炮檢距為0,震源子波采用主頻為20 Hz 的Ricker 子波。
將主動(dòng)源提取的頻散能量圖和背景噪聲提取的頻散曲線進(jìn)行對(duì)比。在提取垂直斷層模型頻散曲線時(shí),發(fā)現(xiàn)m=0 時(shí)頻散曲線提取效果最佳,圖14 即為m=0 時(shí),對(duì)背景噪聲模擬數(shù)據(jù)利用兩道法提取的頻散曲線與主動(dòng)源方法頻散能量圖。由圖可見(jiàn),利用兩道法計(jì)算的頻散曲線(離散點(diǎn))與主動(dòng)源方法模擬頻散能量圖擬合較好,只有少量的離散點(diǎn)偏離理論頻散曲線。
圖14 背景噪聲頻散曲線與主動(dòng)源方法頻散能量圖
本文以三維波動(dòng)方程和二維波動(dòng)方程為理論基礎(chǔ),用隨機(jī)震源產(chǎn)生的振動(dòng)疊加模擬背景噪聲,實(shí)現(xiàn)了不同地質(zhì)模型的背景噪聲波場(chǎng)模擬。對(duì)比、分析模擬數(shù)據(jù)提取的頻散曲線與理論頻散曲線表明,當(dāng)震源參數(shù)處于合理范圍時(shí),模擬數(shù)據(jù)提取的頻散曲線與理論頻散曲線擬合程度較為理想。在無(wú)法得到理論頻散曲線的模型中,將背景噪聲提取的頻散曲線與相同模型使用主動(dòng)源方法得到的頻散能量圖進(jìn)行對(duì)比,也可得到較理想結(jié)果,從而驗(yàn)證了背景噪聲模擬對(duì)較復(fù)雜模型的適用性。此外,三維模擬數(shù)據(jù)頻散曲線提取精度在低頻段比二維模擬數(shù)據(jù)更高,且高頻段的頻散曲線更光滑。