基于非均勻介質(zhì)的譜比曲線正演技術(shù)及應(yīng)用

劉鐵華,劉鐵,張邦,卞友艷,張占榮,化希瑞

(中鐵第四勘察設(shè)計(jì)研究院有限公司，湖北 武漢 430063)

0 引言

隨著城鎮(zhèn)化和城市圈的發(fā)展，城市人口、交通等壓力越來越大，許多城市為建設(shè)集約型城市，亟須開發(fā)地下空間，實(shí)現(xiàn)城市開發(fā)的可持續(xù)發(fā)展，而開展針對(duì)城市地下空間的勘探是地下空間開發(fā)的前提。在復(fù)雜城市環(huán)境中，常用地球物理勘探技術(shù)受建筑交通阻礙、干擾強(qiáng)的因素限制明顯，微動(dòng)技術(shù)因其具有狹小空間適應(yīng)性和強(qiáng)抗機(jī)械干擾能力[1]，在復(fù)雜城市環(huán)境中具有突出應(yīng)用優(yōu)勢(shì)。

微動(dòng)是指地層表面發(fā)育的天然微弱振動(dòng)，其信號(hào)為包含體波(P波和S波)和面波的綜合波場(chǎng)，在城市環(huán)境中，面波占微動(dòng)信號(hào)總能量的70%以上[2]，通常利用微動(dòng)信號(hào)中的面波信息來研究地下橫波速度結(jié)構(gòu)工程。應(yīng)用中，常利用垂直分量微動(dòng)信號(hào)中的瑞利波成分，開展基于頻散曲線的面波勘探(微動(dòng)勘探)，按照裝置和數(shù)據(jù)采集方式不同，面波勘探存在主動(dòng)源面波和被動(dòng)源面波兩種技術(shù)路線[3]，其均是基于采集的數(shù)據(jù)提取頻散曲線后實(shí)現(xiàn)地層結(jié)構(gòu)反演的。此外，還可利用三分量微動(dòng)信號(hào)的水平分量與垂直分量的比值關(guān)系(譜比曲線)，根據(jù)譜比曲線特征開展微動(dòng)譜比法勘探[4-5]。微動(dòng)譜比法僅需單點(diǎn)三分量微動(dòng)數(shù)據(jù)即可，相較于面波法具有更強(qiáng)的場(chǎng)地適應(yīng)性，兩類方法在數(shù)據(jù)采集和處理流程上具有一定的相似性，實(shí)際應(yīng)用中結(jié)合微動(dòng)譜比法和面波法的多源頻率域地震勘探技術(shù)既保留了面波勘探技術(shù)的速度穩(wěn)定性，又繼承了微動(dòng)譜比法的強(qiáng)場(chǎng)地適應(yīng)性和高橫向分辨率特點(diǎn)[6]。

在開展微動(dòng)勘探工作時(shí)，核心技術(shù)體現(xiàn)在微動(dòng)信號(hào)的特征曲線(頻散曲線和譜比曲線)反演過程中，最終獲得勘探范圍內(nèi)的速度結(jié)構(gòu)?，F(xiàn)有的特征曲線正演方法均是建立在水平層狀介質(zhì)模型基礎(chǔ)上，其基本假設(shè)為：①多層層狀介質(zhì)中，瑞利波的相速度存在頻散特性，微動(dòng)水平分量和垂直分量的比值關(guān)系可依據(jù)面波和體波在不同方向的貢獻(xiàn)表示；②面波的能量和探測(cè)深度主要集中于一個(gè)波長(zhǎng)范圍內(nèi)，且不同波長(zhǎng)具有不同的穿透深度；③瑞利波的質(zhì)點(diǎn)振動(dòng)在傳播面上作橢圓運(yùn)動(dòng)，同時(shí)具有水平分量和垂直分量的位移。目前，微動(dòng)特征曲線正演主要基于水平層狀介質(zhì)的假設(shè)前提，但實(shí)際地層，尤其是淺表地層應(yīng)該看作為非均勻介質(zhì)，在存在橫向變化較大情況下導(dǎo)致理論譜比曲線與實(shí)測(cè)譜比曲線擬合度較低，最終導(dǎo)致勘探精度降低。例如，基巖面起伏顯著的地層中，觀測(cè)到的譜比曲線為寬單峰特征，與理論窄單峰譜比曲線存在明顯差異。所以，有必要基于非均勻介質(zhì)假設(shè)條件進(jìn)行微動(dòng)信號(hào)的特征曲線的正演。

1 微動(dòng)勘探技術(shù)

1.1 層狀介質(zhì)中的面波傳播

地震波有兩種基本的類型：一種是在地球介質(zhì)內(nèi)傳播的地震波(體波)，另一種是沿介質(zhì)表面?zhèn)鞑サ牡卣鸩?面波)。體波包括縱波(P波)和橫波(S波)兩種，瑞利波、勒夫波、導(dǎo)波和斯通利波是4種不同類型的面波，工程應(yīng)用中主要基于瑞利波和勒夫波的頻散特性進(jìn)行地質(zhì)勘探。

瑞利波最早是由英國(guó)物理學(xué)家Rayleigh于1885年在理論上推導(dǎo)出，后在天然地震記錄上得到證實(shí)，分布于自由表面上或者表面為疏松的覆蓋層內(nèi)[7]。Haskell發(fā)現(xiàn)層狀介質(zhì)中瑞利波存在頻散特性[8]，即不同頻率成分的瑞利波具有不同的相速度，并推導(dǎo)了瑞利波的頻散方程。瑞利波是通過P和SV波在自由表面上的疊加形成的，能量主要集中在一個(gè)波長(zhǎng)的范圍內(nèi)且不同波長(zhǎng)的瑞利波穿透深度也不同，與體波相比，瑞利波具有能量衰減慢、振幅大、質(zhì)點(diǎn)振動(dòng)存在橢圓極化等特征[9]。此后，人們開始利用其頻散特性來研究地殼、上地幔等深部介質(zhì)中橫波速度結(jié)構(gòu)[10-16]，近年來瑞利波探測(cè)的應(yīng)用領(lǐng)域逐漸擴(kuò)展到淺部工程地球物理勘探。

勒夫面波產(chǎn)生于介質(zhì)表面的低速覆蓋層以及該層與下面介質(zhì)的分界面上，是一種SH型波，具有頻散現(xiàn)象，是面波中傳播速度最快的一種波，由波與自由表面相互作用的結(jié)果。勒夫波使介質(zhì)在水平方向發(fā)生形變，且只有在剪切剛度隨深度增加而增大的情況下才會(huì)出現(xiàn)。

1.2 層狀介質(zhì)中的頻散特征曲線正演

在均勻半空間中，自由界面上的邊界條件是法向與切向的應(yīng)力必須為零，界面一側(cè)為自由空間，無介質(zhì)約束質(zhì)點(diǎn)的振動(dòng)，當(dāng)SV波大于臨界角入射到自由表面時(shí)，與反射P波沿著自由表面前進(jìn),相互干涉便形成了瑞利波。對(duì)于不同的介質(zhì)，瑞利波水平和垂直方向的主要能量均大部分集中在深度小于一個(gè)波長(zhǎng)內(nèi)，即認(rèn)為瑞利波的穿透深度約為一個(gè)波長(zhǎng)。在自由表面上，瑞利波質(zhì)點(diǎn)的振動(dòng)軌跡為一逆進(jìn)的橢圓。大約在0.2倍波長(zhǎng)的深度時(shí)，水平位移極性反轉(zhuǎn)，這就是說從此處看是其動(dòng)軌跡為一順時(shí)針轉(zhuǎn)動(dòng)的橢圓。

目前計(jì)算多層層狀介質(zhì)中瑞利波頻散曲線的方法主要包括[17-20]：Haskell算法、Knopoff 算法、Abo-Zena 算法、廣義反射—透射系數(shù)算法和快速矢量傳遞算法?？紤]到實(shí)際地層的分層性，可以利用圖1所示的水平層狀介質(zhì)模型來近似模擬實(shí)際地下介質(zhì)，在層狀介質(zhì)模型中，沿水平方向介質(zhì)是均勻的而沿垂直方向介質(zhì)是非均勻的。瑞利波在n層水平層狀介質(zhì)模型中傳播，x軸與瑞利波的傳播方向一致，z軸垂直向下并從自由表面起算，界面序號(hào)及層序號(hào)如圖1所示。

圖1 層狀介質(zhì)模型示意Fig.1 Schematic diagram of layered medium model

對(duì)于層狀介質(zhì)模型中，第m層介質(zhì)，其參數(shù)分別為：ρm，第m層的密度；dm，第m層的厚度；λm、μm，第m層的拉梅系數(shù)；c，相速度；ω，角頻率；αm、βm，第m層的縱波和橫波速度；γm=2(βm/c)2?；诮橘|(zhì)中應(yīng)力應(yīng)變連續(xù)原則，使用迭代矩陣求解，由地表已知邊界條件開始，不斷由地層參數(shù)遞推出新的矩陣，即根據(jù)m層的矩陣遞推得到m+1層的矩陣。考慮到自由表面處應(yīng)力為0的條件，滿足以下公式：

(1)

其中，

(2)

(3)

V(0)=0，

(4)

W(0)=(γ1-1)2，

(5)

(6)

S(0)=0。

(7)

式(1)為相速度與頻率的函數(shù)，求解時(shí)根據(jù)式(4)～(7)條件使用層狀模型遞推公式建立針對(duì)當(dāng)前模型的矩陣方程組，再逐一求解一定范圍內(nèi)的各頻率對(duì)應(yīng)的相速度值，通常情況下可以求解出至少1個(gè)解。對(duì)于某一個(gè)固定頻率ω，其最低的相速度值c被稱為基階波相速度，稱具有較高相速度值c的瑞利波為高階模式瑞利波，根據(jù)相速度的大小依次被命名為一階波、二階波等，即在多層層狀介質(zhì)傳播時(shí)瑞利波存在多種模式，如圖2所示。

圖2 層狀介質(zhì)模型各階理論頻散曲線示意Fig.2 Schematic diagram of theoretical dispersion curves of layered media model

1.3 層狀介質(zhì)中的譜比特征曲線正演

微動(dòng)信號(hào)的譜比特征曲線表示水平和垂直分量的能量在各個(gè)頻率上的比值，大多數(shù)研究者傾向于微動(dòng)譜比特征的垂向信號(hào)主要成分為瑞利面波，水平成分為瑞利面波和勒夫波面波及它們的高階分量，微動(dòng)譜比可以用瑞利面波基階橢圓率解釋，用體波解釋微動(dòng)譜比的觀點(diǎn)將被放棄。此外，Arai和Tokimatsu基于擴(kuò)散場(chǎng)理論，假設(shè)環(huán)境噪聲為擴(kuò)散場(chǎng)，在彈性介質(zhì)中，點(diǎn)P處的譜比曲線計(jì)算公式為[21]：

(8)

E(P,ω)=-2πμEsk-1Im[Gmm(P,P)]，

(9)

其中，橫波波數(shù)k=ω/β，橫波平均能量密度Es=ρω2S2，S2為橫波平均譜密度。

當(dāng)瑞利波沿自由表面?zhèn)鞑r(shí)，介質(zhì)質(zhì)點(diǎn)運(yùn)動(dòng)軌跡是橢圓形的，水平分量和垂直分量的振幅比稱為橢圓率，在半無限均勻彈性介質(zhì)中傳播時(shí)，自由表面上瑞利波的質(zhì)點(diǎn)運(yùn)動(dòng)軌跡呈現(xiàn)出逆向橢圓，垂直分量幅度約為水平分量振幅的1.5倍。然而，當(dāng)瑞利波在多層層狀介質(zhì)中傳播時(shí)，橢圓率不是固定值，不同頻率的瑞利波所對(duì)應(yīng)的橢圓率是不同的，例如當(dāng)水平層狀介質(zhì)中存在一個(gè)波阻抗比(相鄰兩層介質(zhì)中橫波速度比值)較大(比值大于3)的界面時(shí)，橢圓率曲線會(huì)在特定的頻率值上呈現(xiàn)奇異點(diǎn)，如峰值和谷值，橢圓率曲線與地下介質(zhì)中橫波速度結(jié)構(gòu)密切相關(guān)。張立對(duì)多層層狀介質(zhì)中瑞利波在水平分量和垂直分量的位移公式進(jìn)行推導(dǎo)[22]，根據(jù)已給定的地層模型參數(shù)，利用瑞利面波速度頻散方程計(jì)算相應(yīng)頻率下的相速度，即獲得頻散數(shù)據(jù)，然后將不同頻率所對(duì)應(yīng)的相速度求得所研究頻帶范圍每個(gè)單一頻率成分所對(duì)應(yīng)的每層瑞利面波質(zhì)點(diǎn)位移的水平分量和垂直分量，進(jìn)而求得質(zhì)點(diǎn)橢圓運(yùn)動(dòng)軌跡(由水平分量和垂直分量合成)隨深度的變化特征。z=0 m處的地表，其橢圓極化特征曲線相當(dāng)于野外采集數(shù)據(jù)水平分量和垂直分量的譜比特征曲線。

2 基于非均勻介質(zhì)的特征曲線正演技術(shù)

2.1 層狀介質(zhì)特征曲線正演缺陷及對(duì)策

上文介紹了微動(dòng)信號(hào)特征曲線的正演方法，現(xiàn)有技術(shù)均是建立在水平層狀介質(zhì)的基礎(chǔ)上，是對(duì)現(xiàn)實(shí)地質(zhì)結(jié)構(gòu)的近似。實(shí)際地層并非為水平地層結(jié)構(gòu)時(shí)，實(shí)測(cè)特征曲線往往與一維水平模型計(jì)算的理論特征曲線存在差異。如圖3所示，為某地面水平下覆基巖為單斜構(gòu)造情況下的實(shí)測(cè)譜比曲線，其表現(xiàn)為較寬的單峰形態(tài)，并非理論上的狹窄單峰特征。

圖3 實(shí)測(cè)譜比曲線示意Fig.3 Schematic diagram of measured spectral ratio curve

出現(xiàn)實(shí)測(cè)特征曲線與理論特征曲線差異的情況，本質(zhì)上是理論特征曲線計(jì)算時(shí)采用的一維水平層狀模型假設(shè)導(dǎo)致的，需要構(gòu)建一套基于非均勻介質(zhì)模型的特征曲線正演方法，本文在研究面波傳播特征和曲線形態(tài)特征的基礎(chǔ)上，提出了動(dòng)態(tài)優(yōu)選技術(shù)的正演方法，可適用于非均勻介質(zhì)模型下的譜比曲線正演，頻散曲線正演的適用性還有待進(jìn)一步研究。

2.2 動(dòng)態(tài)優(yōu)選法正演原理

由面波的波場(chǎng)傳播理論可知，面波的能量和探測(cè)深度主要集中于一個(gè)波長(zhǎng)范圍內(nèi)，且不同波長(zhǎng)具有不同的穿透深度。即可認(rèn)為，非均勻介質(zhì)中地表某點(diǎn)的特征曲線是由來自不同深度范圍的介質(zhì)在不同頻率按椎體形態(tài)由一定比例組成的。

如圖4所示，在求取C7點(diǎn)處特征曲線時(shí)，對(duì)其存在貢獻(xiàn)信號(hào)的空間范圍如圖中綠色陰影區(qū)域所示，該點(diǎn)的特征曲線可表示為：

圖4 動(dòng)態(tài)優(yōu)選法原始示意Fig.4 Schematic diagram of equivalent dynamic weighting method

(10)

其中，xb=x0-0.25tan (θ)V(x,z)/f，xe=x0+0.25tan (θ)V(x,z)/f。A(x,f)為特征曲線，是關(guān)于頻率f和空間坐標(biāo)x的函數(shù)；V(x,z)是模型給出的橫波速度；z為當(dāng)前測(cè)點(diǎn)不同深度的垂向坐標(biāo)值；x0為當(dāng)前測(cè)點(diǎn)的橫向坐標(biāo)值；θ為孔徑系數(shù)，為小于90的角度，當(dāng)θ取45時(shí)表示影響區(qū)域?yàn)橐詼y(cè)點(diǎn)為圓心的向下45°角的圓錐體；λ為優(yōu)選系數(shù)，小于1且趨于1，通常位于0.9～1.0之間。

2.3 非均勻介質(zhì)模型的特征曲線計(jì)算流程

微動(dòng)信號(hào)特征曲線均是關(guān)于頻率的函數(shù)，為測(cè)點(diǎn)處地下空間一定范圍內(nèi)波場(chǎng)的綜合反應(yīng)，基于非均勻介質(zhì)的動(dòng)態(tài)優(yōu)選法有效實(shí)現(xiàn)了該條件下的特征曲線正演，因微動(dòng)譜比曲線具有更高的橫向分辨率，受不均勻性影響更為明顯。本文以非均勻條件下微動(dòng)譜比曲線的正演為例，說明動(dòng)態(tài)優(yōu)選法實(shí)施的具體思路：首先，計(jì)算測(cè)點(diǎn)附近一定影響范圍內(nèi)不同位置的等效水平層狀模型的譜比曲線；其次，計(jì)算隨頻率變化的孔徑范圍和加權(quán)系數(shù)，通過動(dòng)態(tài)優(yōu)選法的方式計(jì)算測(cè)點(diǎn)處的綜合譜比曲線，計(jì)算流程如下：

1)建立等效水平層狀模型。將非均勻介質(zhì)按一定距離dx分割成多個(gè)獨(dú)立模型，再將每個(gè)獨(dú)立模型等效為水平層狀模型。如圖4，可建立兩層模型，上層覆蓋層速度設(shè)為v0，下層基巖速度設(shè)為v1。

2)計(jì)算各等效水平層狀模型的譜比曲線。計(jì)算步驟1)中按固定間隔dx分割的等效水平層狀模型的理論譜比曲線；圖4中 C5、C6、C7、C8、C9五個(gè)等效水平層狀模型點(diǎn)對(duì)應(yīng)的理論譜比曲線如圖5所示。

圖5 不同測(cè)點(diǎn)譜比曲線Fig.5 Spectral ratio curves of different measuring points

3)計(jì)算模擬點(diǎn)頻率f的影響孔徑。以當(dāng)前模擬點(diǎn)為中心，計(jì)算不同頻率f的影響孔徑，計(jì)算公式為R=0.25×k×v0/f，其中R為孔徑的半徑；v0為覆蓋層速度；k為空間轉(zhuǎn)換系數(shù)，為一常數(shù)，通?？扇?。對(duì)于C7測(cè)點(diǎn)的影響孔徑，通過計(jì)算可能為C5、C6、C7、C8、C9五處等效水平層狀模型點(diǎn)。

4)搜索頻率f對(duì)應(yīng)孔徑內(nèi)模擬點(diǎn)數(shù)。對(duì)孔徑內(nèi)各模擬點(diǎn)的對(duì)應(yīng)頻點(diǎn)f，以步驟3)中計(jì)算得到的空間R為半徑，搜索該范圍內(nèi)符合條件的模擬點(diǎn)的模擬譜比曲線，取頻率f對(duì)應(yīng)的譜比值Qi，1<=i<=n，n為孔徑內(nèi)的模擬點(diǎn)個(gè)數(shù)。

5)計(jì)算f對(duì)應(yīng)孔徑內(nèi)各點(diǎn)權(quán)系數(shù)。對(duì)步驟4)中各點(diǎn)，按各點(diǎn)到所計(jì)算模擬點(diǎn)的距離計(jì)算權(quán)值Wi=1-λr/R，其中r為各點(diǎn)到模擬點(diǎn)的距離，i為孔徑內(nèi)點(diǎn)編號(hào)，i<=n，λ為優(yōu)選系數(shù)。

7)綜合特征曲線計(jì)算。重復(fù)步驟3)～步驟6)的過程，計(jì)算所有頻率f對(duì)應(yīng)的綜合特征值，得到綜合特征曲線，如圖6所示。

圖6 C7測(cè)點(diǎn)綜合譜比曲線Fig.6 C7 comprehensive spectral ratio curve

以上步驟為雙層非均勻介質(zhì)條件下的特征曲線計(jì)算簡(jiǎn)化流程，在實(shí)際正演計(jì)算中通常更為復(fù)雜，不能簡(jiǎn)單地按照頻率進(jìn)行優(yōu)選，需先根據(jù)模型的層速度信息將頻率域的特征曲線轉(zhuǎn)化到深度域，再在深度域進(jìn)行優(yōu)選，最終將優(yōu)選的特征曲線轉(zhuǎn)換到頻率域。

3 應(yīng)用案例分析

3.1 項(xiàng)目概況

應(yīng)用項(xiàng)目為我國(guó)南方某濱海灘涂地貌和海沖積平原地貌，因城市建設(shè)沿線場(chǎng)地經(jīng)填、挖、整平等人工改造，形成現(xiàn)在的較平坦的地形地貌，地層為第四系殘積土、震旦系變粒巖等。地表多為居民民房和工業(yè)廠房，周邊交通繁忙，電磁波干擾和機(jī)械干擾嚴(yán)重。

勘探工作處于工程前期定測(cè)階段，根據(jù)前期勘探成果已經(jīng)確認(rèn)該段基巖深度出現(xiàn)有規(guī)律的變化，采用被動(dòng)源面波與微動(dòng)譜比法聯(lián)合的多源頻率域勘探技術(shù)進(jìn)行勘探，旨在獲取勘探范圍內(nèi)的基巖起伏情況?？碧焦ぷ鞑捎?道1排列線性多源頻率地震數(shù)據(jù)采集裝置進(jìn)行數(shù)據(jù)采集，點(diǎn)距約10 m(根據(jù)場(chǎng)地條件作適當(dāng)調(diào)整)。數(shù)據(jù)采集過程中，改善耦合條件，克服現(xiàn)場(chǎng)隨機(jī)干擾，水平角度控制在2°范圍內(nèi)，授時(shí)精度控制在1×10-6s內(nèi)，數(shù)據(jù)質(zhì)量良好。

3.2 勘探成果應(yīng)用分析

本文選取該區(qū)域內(nèi)某典型單斜—水平結(jié)構(gòu)段進(jìn)行分析(圖7)，測(cè)線長(zhǎng)度200 m，由左往右基巖逐漸變深，在基巖埋深漸變段譜比曲線(圖7b)為較窄的單峰結(jié)構(gòu)，而末端近似水平段的譜比曲線(圖7a)表現(xiàn)為較寬單峰形態(tài)，這一特征印證了前文分析結(jié)果。

a—水平段；b—單斜段a—the horizontal section;b—the single section圖7 典型實(shí)測(cè)數(shù)據(jù)譜比曲線Fig.7 Spectral ratio curves of typical measured data

應(yīng)用中，首先提取了3個(gè)排列的頻散曲線和譜比曲線(每個(gè)排列由7道三分量長(zhǎng)周期微動(dòng)數(shù)據(jù)組成)，兩種特征曲線采用聯(lián)合反演的方式獲得地層橫波速度剖面，綜合解釋橫波速度475 m/s為土石界面(W2與W3風(fēng)化界線)，如圖8所示。

完成地面微動(dòng)勘探后，在里程2 565 m和2 760 m處布置了兩個(gè)鉆孔，兩處揭示的土石界面標(biāo)高如圖中“十”字標(biāo)所示。對(duì)比圖8a、b，可明顯發(fā)現(xiàn)非均勻介質(zhì)模型的反演橫波速度剖面對(duì)橫向的不均勻性靈敏度更高，對(duì)基巖探測(cè)具有更高的探測(cè)精度。

a—基于非均勻介質(zhì)反演的橫波速度剖面；b—基于水平水層層狀介質(zhì)反演的橫波速度剖面a—shear wave velocity profile based on inhomogeneous medium inversion;b—shear wave velocity profile based on horizontal water layered medium inversion圖8 反演后的橫波速度剖面成果Fig.8 S-wave velocity profile after inversion

4 結(jié)論

本文針對(duì)微動(dòng)信號(hào)處理中的正演問題，提出了一種新的正演方法，取得以下結(jié)論：

1)水平層狀介質(zhì)假設(shè)的地質(zhì)模型無法滿足微動(dòng)特征曲線的高精度正演要求，需基于非均勻模型進(jìn)行正演；

2)本文提出的動(dòng)態(tài)優(yōu)選法的微動(dòng)特征曲線正演技術(shù)，在等效水平層狀介質(zhì)理論曲線的基礎(chǔ)上，計(jì)算隨頻率變化的孔徑范圍和優(yōu)選系數(shù)，可實(shí)現(xiàn)非均勻介質(zhì)的高精度譜比曲線正演；

3)通過實(shí)際應(yīng)用表明，采用動(dòng)態(tài)優(yōu)選法微動(dòng)特征曲線正演技術(shù)后，能夠更好地反演非均勻介質(zhì)的橫波分布，勘探成果具有更高的橫向分辨率和勘探精度。