趙雪然，劉四新，張健橋，侯國正，范運(yùn)嶺

1.吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長春 130026；2.山東省物化探勘查院，濟(jì)南 250013；3.廣東省地質(zhì)調(diào)查院，廣州 510000

0 引言

微動(dòng)是地球表層時(shí)刻存在著的非地震引起的微弱振動(dòng)。微動(dòng)在地震原始記錄中表現(xiàn)為連續(xù)擾動(dòng)，過去對(duì)于研究地震波來說是一種難以排除的干擾噪聲。然而，微動(dòng)中其實(shí)攜帶著與地殼淺部介質(zhì)密切相關(guān)的信息[1]。

微動(dòng)是一種由體波(P波和S波)和面波(瑞利波和拉夫波)組成的復(fù)雜振動(dòng)，面波的能量占總能量的70%以上[2]，并且它與體波的主要區(qū)別是在不均勻的介質(zhì)中傳播時(shí)面波發(fā)生頻散，而體波無頻散現(xiàn)象[3]。微動(dòng)在時(shí)間和空間上存在高度變化、無規(guī)律性及無重復(fù)性的特點(diǎn)，研究表明世界各地的微動(dòng)頻譜形態(tài)大體相似[4]。所以，微動(dòng)的頻譜特性反映了微動(dòng)在時(shí)間和空間上的變化，即微動(dòng)信號(hào)在一定時(shí)空范圍內(nèi)具有統(tǒng)計(jì)穩(wěn)定性，可以用時(shí)間和空間上的平穩(wěn)隨機(jī)過程描述。微動(dòng)探測(cè)就是以此理論為基礎(chǔ)從信號(hào)中提取面波的頻散曲線，然后反演出地下介質(zhì)的橫波速度結(jié)構(gòu)[5]。

微動(dòng)探測(cè)方法在國外，特別是在日本已廣泛應(yīng)用于盆地結(jié)構(gòu)調(diào)查、工程地質(zhì)調(diào)查和場(chǎng)地穩(wěn)定性評(píng)價(jià)等[6--7]。中國是從20世紀(jì)80年代末到90年代初對(duì)長周期微動(dòng)展開了研究。近幾年來利用微動(dòng)臺(tái)陣觀測(cè)技術(shù)界定地層速度結(jié)構(gòu)已經(jīng)應(yīng)用到城市地質(zhì)調(diào)查、斷層探測(cè)、地?zé)嵴{(diào)查和基巖埋深探測(cè)等多個(gè)領(lǐng)域[8--9]，探測(cè)深度從近地表到幾千米。

微動(dòng)勘探方法根據(jù)臺(tái)陣布設(shè)形式和頻散曲線提取方法主要分為頻率--波數(shù)法(frequency-wavenumber)、空間自相關(guān)法(spatial auto correlation)和背景噪聲互相關(guān)法(ambient noise correlation function)等。頻率--波數(shù)法需要相對(duì)大量的地震臺(tái)站，臺(tái)站布置越多，頻譜精度越高[10]；空間自相關(guān)法在現(xiàn)場(chǎng)數(shù)據(jù)采集和數(shù)據(jù)分析方面比頻率--波數(shù)法簡單，但需要布置圓形臺(tái)陣，在實(shí)際中受地理?xiàng)l件影響，要盡量減少臺(tái)站布置[11]；背景噪聲互相關(guān)法適用于大尺度地球探測(cè)，在小尺度面波淺勘領(lǐng)域應(yīng)用效果不及頻率--波數(shù)法和空間自相關(guān)法。

傳統(tǒng)的空間自相關(guān)法(以下簡稱SPAC法)通常采用圓形臺(tái)陣(圖1)，對(duì)臺(tái)陣圓周上的檢波器采集的微動(dòng)信號(hào)分別與圓心處采集的微動(dòng)信號(hào)做相關(guān)計(jì)算，進(jìn)而求取頻散曲線，再對(duì)曲線求平均。筆者針對(duì)常見的雙重圓形臺(tái)陣，研究其不同檢波距組合方式的頻散曲線有效頻帶范圍與測(cè)點(diǎn)之間相關(guān)性的好壞、是否共中心點(diǎn)等因素之間的關(guān)系，針對(duì)不同半徑的頻散曲線，用分段平均的方式取代原先直接平均的方式，得到了精度更高的頻散曲線。

圖1 微動(dòng)觀測(cè)臺(tái)陣示意圖Fig.1 Sketch of microtremor observation array

1 微動(dòng)勘探SPAC法的原理及方法

空間自相關(guān)法(以下簡稱SPAC法)最早由Aki提出其理論基礎(chǔ)和實(shí)驗(yàn)方法，其理論認(rèn)為：①隨機(jī)振動(dòng)在時(shí)間和空間上具有穩(wěn)定的特性；②該微動(dòng)波場(chǎng)主要由傳播在地表面的頻散性波動(dòng)組成的。在此條件下，計(jì)算的不同臺(tái)站對(duì)的方位平均的自相關(guān)系數(shù)近似為第一類零階貝塞爾(Bessel)函數(shù)[12]。

其基本原理是在自由表面上利用地震檢波器圓形排列臺(tái)陣觀測(cè)微動(dòng)，對(duì)每個(gè)臺(tái)站的波形記錄在相同時(shí)間段內(nèi)進(jìn)行快速傅里葉變換求功率譜和互功率譜：

(1)

(2)

(3)

式中：R1,2=∑u0(t)ur(t+τ)，其中u0為圓心點(diǎn)選取的波形記錄，ur為半徑為r的圓上選取的波形記錄；R1,1=∑u0(t)u0(t+τ)；R2,2=∑ur(t)ur(t+τ)。

S(r,θ,ω)為圓心點(diǎn)與半徑為r的圓上一點(diǎn)的互功率譜；S(0,ω)為圓心點(diǎn)的自功率譜；S(r,ω)為半徑r的圓上一點(diǎn)的自功率譜。得到功率譜之后，計(jì)算每個(gè)臺(tái)站對(duì)的空間自相關(guān)函數(shù)，最后將半徑相同的臺(tái)站對(duì)的空間自相關(guān)函數(shù)進(jìn)行方位平均，得到該半徑對(duì)應(yīng)的空間自相關(guān)系數(shù)：

(4)

根據(jù)Aki的理論假設(shè)，推導(dǎo)出空間自相關(guān)系數(shù)可以用零階貝塞爾函數(shù)表示，即：

(5)

式中：r為兩臺(tái)站間的距離；f為頻率；c(f)為瑞雷波相速度；ρ為方位平均后的自相關(guān)系數(shù)；J0為第一類零階貝塞爾函數(shù)。該式適用于以基階面波能量為主，利用垂向分量觀測(cè)數(shù)據(jù)進(jìn)行空間自相關(guān)，提取瑞雷面波頻散的情形。

2 微動(dòng)數(shù)據(jù)處理

2.1 臺(tái)站的設(shè)計(jì)和數(shù)據(jù)采集

SPAC方法在野外工作中通常采用雙重圓形臺(tái)陣(圖2)，即一個(gè)檢波器(X7)分布在圓心O，檢波器X2、X4、X6均勻分布在半徑為r的圓周上，檢波器X1、X3、X5均勻分布在半徑為2r的圓周上，三角形臺(tái)陣布設(shè)可有效改善波源的方位不均勻。

圖2 微動(dòng)雙重觀測(cè)臺(tái)陣示意圖Fig.2 Sketch of microtremor double observation array

2.2 數(shù)據(jù)處理的基本步驟

利用微動(dòng)數(shù)據(jù)反演S波速度結(jié)構(gòu)時(shí)，首先對(duì)原始數(shù)據(jù)預(yù)處理，然后進(jìn)行傅里葉變換，計(jì)算每個(gè)點(diǎn)的自功率譜(圖3)以及每兩個(gè)點(diǎn)之間的互功率譜；再根據(jù)公式(4)求出方位平均的自相關(guān)系數(shù)，最后利用SPAC法提取頻散曲線。

圖3 自功率譜圖Fig.3 Auto-power spectrum

2.3 利用SPAC法提取頻散曲線

通過上述步驟得到自功率譜和互功率譜之后，再利用公式(5)，將自相關(guān)系數(shù)與零階一類貝塞爾函數(shù)擬合，計(jì)算出相速度，提取頻散曲線?？臻g自相關(guān)系數(shù)ρ(r,f)有兩個(gè)變量：頻率f和半徑r，通常固定其中一個(gè)變量再進(jìn)行擬合，具體過程為：SPAC方法是固定半徑為r0，用標(biāo)準(zhǔn)零階貝塞爾函數(shù)擬合ρ(r0,f)-f曲線，求出貝塞爾函數(shù)的宗量x；進(jìn)而由x=2πfr0/c(f)得到瑞雷波相速度c。

圖4 多重觀測(cè)臺(tái)陣示意圖Fig.4 Sketch of microtremor multiple observation array

SPAC法中，通常采用多重三角臺(tái)陣(圖4)，對(duì)于雙重臺(tái)陣，有5種半徑組合方式(圖5)，每一種半徑組合方式可以得到一條頻散曲線，將5條曲線平均即可得到最終的頻散曲線[14]。實(shí)際中不同的半徑組合方式得到的自相關(guān)系數(shù)的可靠頻帶范圍并不同，因此直接將不同半徑對(duì)應(yīng)的頻散曲線平均會(huì)產(chǎn)生較大的誤差。本文針對(duì)這一問題提出了分段平均的方法。

圖5 雙重臺(tái)陣中不同半徑的組合方式Fig.5 Combination mode of different radii in a double array

圖6 不同半徑的自相關(guān)系數(shù)曲線Fig.6 Autocorrelation coefficient curves with different radii

在對(duì)實(shí)際數(shù)據(jù)進(jìn)行處理時(shí)，恰當(dāng)?shù)倪x擇擬合頻段范圍對(duì)于估算相速度的準(zhǔn)確性十分重要。根據(jù)擬合貝塞爾函數(shù)時(shí)誤差的大小，可以將貝塞爾函數(shù)劃分成單調(diào)遞減和振蕩衰減兩部分(圖7)。單調(diào)遞減部分即宗量x從零到第一個(gè)波谷點(diǎn)的部分，振蕩衰減部分是宗量大于第一個(gè)波谷點(diǎn)的部分，函數(shù)值在零附近振蕩且逐漸衰減。

根據(jù)Cho等人的研究，當(dāng)空間自相關(guān)系數(shù)在貝塞爾函數(shù)遞減部分，其擬合誤差較小，能得到準(zhǔn)確的相速度；當(dāng)空間自相關(guān)系數(shù)在貝塞爾函數(shù)振蕩部分，其擬合誤差較大，估算的相速度值誤差也大[15]。因此通常選擇自相關(guān)系數(shù)曲線的單調(diào)遞減部分與貝塞爾函數(shù)進(jìn)行擬合。

圖7 標(biāo)準(zhǔn)一類零階貝塞爾函數(shù)Fig.7 Standard class zero order Bessel function

通過對(duì)比不同半徑的空間自相關(guān)系數(shù)與標(biāo)準(zhǔn)零階貝塞爾函數(shù)不難發(fā)現(xiàn)，半徑為2r的臺(tái)站對(duì)的自相關(guān)系數(shù)的有效頻段大于半徑為r的臺(tái)站對(duì)自相關(guān)系數(shù)的有效頻段，而其他臺(tái)站對(duì)組合方式，由于臺(tái)站對(duì)并不共中心點(diǎn)，相關(guān)性較差，對(duì)平穩(wěn)隨機(jī)微動(dòng)源的要求更高，而實(shí)際工作中隨機(jī)源的方位并不均勻，且可能存在優(yōu)勢(shì)源，所以可利用的頻率范圍更小[16]。即在圖5中，雖然d、e組合方式的半徑更大，但有效頻帶范圍并沒有明顯大于前3種。第c種的臺(tái)站組合方式在共中心點(diǎn)的前提下，半徑也相對(duì)較大，因此有效頻帶范圍最寬，第a、b種在相對(duì)高頻部分有可靠地頻散曲線，而第e種作為半徑最大的組合方式，在低頻部分有可靠的頻散曲線。

綜上所述，針對(duì)雙重圓形觀測(cè)臺(tái)陣，在擬合頻散曲線時(shí)，首先剔除異常點(diǎn)，然后采取在低頻部分選取半徑為2r、3r、3.46r的曲線的平均，在較高頻部分選取半徑為r、1.73r、2r的曲線的平均，通過這種分段平均的方式不僅擴(kuò)寬了頻帶寬度，而且相比于通常把不同半徑的頻散曲線直接簡單做平均的處理，這種分段平均的方式精確度更高。

圖8 分段平均擬合出的頻散曲線Fig.8 Dispersion curve obtained by piecewise average fitting

最后將利用多種半徑組合、分段擬合方法得到的頻散曲線進(jìn)行反演，獲得地下介質(zhì)的S波速度結(jié)構(gòu)(圖9)。

圖9 反演的S波速度結(jié)構(gòu)Fig.9 Velocity structure of inversed S-wave

3 場(chǎng)地試驗(yàn)分析

本試驗(yàn)測(cè)線南北向布設(shè)，臺(tái)站點(diǎn)距10 m，臺(tái)站采用雙重觀測(cè)系統(tǒng)，由兩個(gè)同心圓的內(nèi)接正三角形和圓心組成，小三角形邊長為10 m，大三角形邊長為20 m。使用1 Hz垂直檢波器，德國Summit地震儀。采樣頻率1 ms，單次采樣時(shí)間16 s，重復(fù)采樣若干次。

本次試驗(yàn)由于臺(tái)站間距較小，所以采用固定距離(即半徑)，擬合ρ-f曲線。根據(jù)前人經(jīng)驗(yàn)分析[17--18]，空間自相關(guān)法最大勘探深度為臺(tái)陣半徑的2至7倍之間，除了臺(tái)站距離，具體勘探深度還受多種因素影響，如檢波器周期、數(shù)據(jù)質(zhì)量等。

分別采用傳統(tǒng)的SPAC法和改進(jìn)后的SPAC法對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行處理，對(duì)比最后獲得的二維橫波速度--深度圖(圖10)，可以看出，通過多種半徑組合方式，分段平均的方法探測(cè)深度可以達(dá)到80 m，而傳統(tǒng)的SPAC法只能探測(cè)到30 m，說明本文提出的方法對(duì)提高探測(cè)深度和精度有一定效果。

從二維橫波速度--深度圖(圖10)可以看到，從地表到35 m±，橫波速度在100～300 m/s，為覆蓋土層，35 m以下速度明顯增大，判斷為花崗巖基巖。通常用同一巖體中風(fēng)化與未風(fēng)化巖體的縱波速度之比可以將花崗巖的風(fēng)化程度進(jìn)行劃分(表1)[19]。類似的，可以利用同一巖體中風(fēng)化與未風(fēng)化巖體的橫波速度之比對(duì)巖體風(fēng)化程度進(jìn)行劃分，且橫波波速比與縱波波速比的值近似[20]，因此通過巖體中橫波波速可以得到波速比，進(jìn)而推斷巖體風(fēng)化程度。

表1 花崗巖風(fēng)化程度分級(jí)

結(jié)合表1分析，可以將本次試驗(yàn)獲得的橫波速度剖面進(jìn)行劃分：

(1)覆蓋土：從地表到40 m±，vs速度約在150～300 m/s，在10～20 m范圍內(nèi)出現(xiàn)部分低速異常，vs速度約100 m/s±，判斷為礫質(zhì)覆蓋土。

(2)強(qiáng)--中風(fēng)化花崗巖層：vs速度約在300～400 m/s，速度呈層狀變化，較均勻。頂界面埋深約在35～40 m，有一定起伏，厚度大約為20 m。

a.傳統(tǒng)SPAC法得到的速度--深度圖;b.改進(jìn)后得到的速度--深度圖。圖10 二維橫波速度--深度對(duì)比圖Fig.10 2D S-wave speed-depth comparison map

(3)微風(fēng)化花崗巖層：vs速度約在400～450 m/s，頂界面埋深在45～60 m，界面起伏較大，最厚處達(dá)到20 m，最薄僅5 m。

(4)未風(fēng)化花崗巖層：vs速度>500 m/s，頂界面埋深在70 m±，界面相對(duì)比較平緩，速度穩(wěn)定，巖性均勻。

根據(jù)測(cè)線起點(diǎn)附近的ZK0002鉆孔巖芯揭露，0～4 m為素填土；4～37.4 m為粉砂質(zhì)淤泥質(zhì)土，粉砂組主要為細(xì)粒石英砂，其中在17 m±有少量腐植質(zhì)及貝殼碎屑；37.4～44.8 m為石英砂及少量碎石土；44.8 m以下均為早奧陶紀(jì)黑云母二長細(xì)粒花崗巖，44.8～57 m花崗巖呈散體結(jié)構(gòu)，風(fēng)化嚴(yán)重；57～75 m的花崗巖中部分長石微風(fēng)化，少量中風(fēng)化；75 m以下的花崗巖及結(jié)構(gòu)清晰，幾乎未風(fēng)化。

從圖11的二維速度--深度圖與鉆孔巖芯對(duì)比可以看出，通過SPAC法得到的巖性劃分與鉆孔巖芯資料基本吻合，驗(yàn)證了本方法的正確性。

圖11 二維橫波速度--深度圖Fig.11 2D S-wave speed-depth map

4 結(jié)論

(1)通過分析雙重圓形觀測(cè)臺(tái)陣中不同半徑組合方式下的自相關(guān)系數(shù)曲線，在共中心點(diǎn)的情況下，半徑越大的臺(tái)站對(duì)求出的自相關(guān)系數(shù)曲線形態(tài)越好，擬合得到的頻散曲線可靠頻帶范圍也越寬，驗(yàn)證了共中心點(diǎn)觀測(cè)臺(tái)陣對(duì)面波入射方向均勻性的要求較低，效果較好。

(2)針對(duì)微動(dòng)信號(hào)處理，對(duì)傳統(tǒng)SPAC法進(jìn)行改進(jìn)，即根據(jù)其臺(tái)站分布方式和半徑大小，選擇各自可靠部分的曲線，最后將不同半徑的頻散曲線進(jìn)行綜合，主要以共中心點(diǎn)且距離最大的組合方式所得到的頻散曲線為主，大于此距離的組合方式對(duì)應(yīng)的頻散曲線取低頻部分，小于此距離的組合方式對(duì)應(yīng)的頻散曲線取高頻部分，最后平均得到最終的頻散曲線。通過此方法，提高了SPAC法求取頻散曲線的精度，進(jìn)一步提高了探測(cè)精度。

(3)通過對(duì)實(shí)測(cè)數(shù)據(jù)的處理，與鉆孔資料進(jìn)行對(duì)比，驗(yàn)證了本文提出的方法具有計(jì)算精度更高，探測(cè)范圍更廣的優(yōu)勢(shì)。