李志遠，梁光河，谷丙洛

1中國科學(xué)院礦產(chǎn)資源研究重點實驗室中國科學(xué)院地質(zhì)與地球物理研究所，北京 100029

2中國科學(xué)院大學(xué)，北京 100049

1 引 言

多波多分量地震勘探技術(shù)是具有極大的科學(xué)價值和發(fā)展前途的勘探地震學(xué)前緣學(xué)科之一.它在提高地震資料的分辨率，提高成像精度，精細預(yù)測巖性，研究地下介質(zhì)方位各向異性等方面有著常規(guī)縱波勘探無法比擬的優(yōu)勢［1］.多波多分量地震勘探技術(shù)將應(yīng)用于所有的陸地或是海底地震勘探［2］.

隨著多波多分量地震勘探技術(shù)的發(fā)展，多波多分量地震資料處理技術(shù)也得到了快速發(fā)展.基于波場分離的多波多分量處理技術(shù)便是其中很重要的一類處理方法.這類方法以標量波場理論為基礎(chǔ)，首先將多波多分量地震數(shù)據(jù)分解為縱波數(shù)據(jù)和轉(zhuǎn)換橫波數(shù)據(jù)，然后對這兩種數(shù)據(jù)分別進行處理［3］.因此，縱、橫波的分離是多波多分量地震資料處理中很重要的一步.

縱、橫波的分離方法有許多種，例如：偏振濾波法［4－5］，Radon變換法［6］，波動方程法［7］等.本文依據(jù)Helmholtz分解理論分離縱、橫波利用的是這兩種波場偏振特性的差異，屬于偏振濾波法.

由Helmholtz分解理論可知，一個矢量場可以被分解成一個無旋的標量場和一個無散的矢量場.而在各向同性介質(zhì)中，縱波是無旋場，橫波是無散場，因此根據(jù)Helmholtz分解理論，對地震波場分別求取散度和旋度，可以提取出純縱波和純橫波［8］.但是在各向異性介質(zhì)中，縱波的偏振方向并不完全平行于波的傳播方向，橫波的偏振方向也不完全垂直于波的傳播方向.因此對于各向異性介質(zhì)中的波場，直接求散度和旋度無法得到純的縱、橫波.為了在各向異性介質(zhì)中也能夠分離出縱、橫波，Dellinger和Etgen［9］以Helmholtz分解理論為基礎(chǔ)，在頻率－波數(shù)域求解Christoffel方程，得到縱、橫波質(zhì)點的振動方向，繼而得到所謂的波型分離算子（類似于散度和旋度算子），然后再反變換回時間－空間域，從而實現(xiàn)了二維VTI介質(zhì)中的縱、橫波分離.Dellinger［10］又將此方法推廣到三維各向異性介質(zhì)中.但是此分離算子對于均勻介質(zhì)效果很好，而對于非均勻介質(zhì)，其適用性不佳.

Yan和Sava［11－12］在Dellinger理 論 的 基 礎(chǔ) 上，利用Christoffel方程在每個網(wǎng)格節(jié)點上計算分離算子，從而將其適用性從均勻VTI介質(zhì)推廣到非均勻的TTI介質(zhì)和VTI介質(zhì)中.Zhang和McMechan［13］對比研究了前幾種方法，提出了適用性更廣的波場分離方法.該方法以Helmholtz分解理論和Christoffel方程為基礎(chǔ)，推導(dǎo)出可以直接獲得縱、橫波各分量的公式.其物理意義明確，而且可以很好的適用于3D各向異性介質(zhì)中.Yan和Sava［14］針對該方法計算成本高、效率低的問題，提出了更為高效的計算方法.郭鵬等［15］也研究了非均勻VTI介質(zhì)中的縱、橫波分離.

上述的縱、橫波的分離方法，只適用于正演模擬.對于給定的地震記錄，堯德中等［16］從Helmholtz分解理論出發(fā)，在頻率－波數(shù)域?qū)SP地震記錄求取散度和旋度，從而解決了地震記錄中有一個方向上的偏導(dǎo)數(shù)無法求取的問題.然后再反變回時間－空間域，得到純的縱、橫波地震記錄.但是該方法要求各地層的縱、橫波速度為已知.Sun［17］將Helmholtz分解理論與波場延拓相結(jié)合，先將地面地震記錄在均勻空間中延拓到某一位置，然后利用求散度和旋度的方法得到純的縱、橫波，最后再將分離后的波場分別逆向延拓回地面，從而實現(xiàn)地面地震記錄的縱、橫波的分離.Sun等［18］將此方法應(yīng)用到3D介質(zhì)中.但是此方法對波場進行了散度和旋度的運算，使分離后的波場產(chǎn)生了相移，并且使得縱、橫波的振幅比也發(fā)生了改變.因此需要對產(chǎn)生的相移進行校正［19］，并恢復(fù)真實的縱、橫波振幅比［20］.

本文依據(jù)Helmholtz分解理論，在頻率－波數(shù)域求解散度和旋度來進行縱、橫波分離.針對地表速度不準確造成縱、橫波分離不徹底的問題，提出了準確估算地表速度值的方法.針對進行散度和旋度運算后，縱、橫波的相位和振幅比會發(fā)生改變的問題，提出了相應(yīng)的相位校正以及縱、橫波振幅比校正的方法.

2 利用散度和旋度進行縱橫波分離的原理

堯德中等［16］指出在頻率－波數(shù)域進行散度和旋度運算可以解決地震記錄中某個方向的偏導(dǎo)數(shù)無法求取的問題.

根據(jù)Helmholtz分解理論，位移場U可表示為一個標量勢的梯度與一個矢量勢的散度之和［8］.用公式表示為

根據(jù)（2）、（3）兩式對（1）式兩邊分別求旋度和散度，可得：

由于S只包含了橫波分量，P只包含了縱波分量，因此可以用S和P來分別表示橫波波場和縱波波場.

在二維情況下，令U＝［u，w］，其中u表示位移場的水平分量，w表示位移場的垂直分量.將（4）、（5）兩式展開，可得（二維情況下S是一標量）：

在時間－空間域中，地面地震記錄z方向的偏導(dǎo)數(shù)是無法求得的.將（6），（7）兩式變換到頻率－波數(shù)域求解［16，21］，有：

由于在地面接收到的是上行波，所以（10）式取負號.對于（8）式，v取縱波波速，對于（9）式，v取橫波波速［16］.求得和后，再反變換回時間－空間域，便可得到分離的縱、橫波.

圖1為合成的地面地震記錄.模型為一均勻介質(zhì)，速度是vp＝3000m／s，vs＝1700m／s，所用震源為集中力震源，置于模型中心.

圖1 合成地震記錄（a）水平分量；（b）垂直分量.Fig.1 Synthetic seismic data（a）horizontal component（b）vetical component

從圖1中抽取一道地震數(shù)據(jù)顯示在圖2中，圖2a是水平分量（用u表示），圖2b是垂直分量（用w表示）.圖2c中的實線是u對z的偏導(dǎo)數(shù)，虛線是w對x的偏導(dǎo)數(shù)；圖2e是圖2c中兩條線相減的結(jié)果.從圖2c和圖2e中可以看出，（6）式中之所以能將縱波成分消除掉，是因為?u／?z和?w／?x中的縱波成分相等，橫波成分不相等，兩者相減便可將縱波成分消除，只保留橫波成分.圖2d中的實線是u對x的偏導(dǎo)數(shù)，虛線是w對z的偏導(dǎo)數(shù)；圖2f是圖2d中兩條線相加的結(jié)果.從圖2d和圖2f中可以看出，（7）式中之所以能將橫波成分消除，是因為?u／?x和?w／?z中的橫波成分振幅大小相等，相位相差π，而縱波成分振幅不等，相位相差很小，兩者相加便可將橫波成分消除，只保留縱波成分.

在上述計算過程中，給出的速度為真實速度.如果給出的速度不準確，結(jié)果如何？我們知道，（8）式和（9）式中的kx是與速度無關(guān)的，由（10）式可知，kz與速度有關(guān).也就是說，速度影響的只是u和w對z的偏導(dǎo)數(shù)（?u／?z和?w／?z）.

圖3a為縱波速度分別為2500、3000、3500m／s時，?u／?z的波形圖.從圖中可以看出速度對?u／?z的振幅影響很明顯，而且速度越快振幅越小，速度越慢振幅越大.速度對相位也有影響，但是影響比較小，在圖中很難看出.這表明，如果速度給的不準確，?u／?z－?w／?x是不能夠?qū)⒖v波成分消除掉的，如圖3c所示，縱波有殘留.圖3b為橫波速度分別為1200、1700、2200m／s時，?w／?z的波形圖.從圖中可以看出速度對?w／?z的影響與對?u／?z的影響一樣有著相同的規(guī)律.在速度不準確的情況下，?u／?x＋?w／?z也無法將橫波成分消除掉，如圖3d所示，橫波有殘留.從上面分析可以知，只有在速度準確的情況下，縱橫波才能分離開.因此，需要估算地表的縱、橫波速度.

圖2 圖1中的一道地震數(shù)據(jù)（a）水平分量u；（b）垂直分量w；（c）?u／?z和?w／?x的波形圖；（d）?u／?x和?w／?z的波形圖；（e）分離出的橫波；（f）分離出的縱波.Fig.2 A trace data of Fig.1（a）Horizontal component u；（b）Vetical component w；（c）Waveform of?u／?zand?w／?x；（d）Waveform of?u／?xand?w／?z；（e）The separated P－waves；（f）The separated S－waves.

3 地表縱、橫波速度的估算

令uP（ω）和wP（ω）分別表示單頻ω上，平面縱波的水平分量和垂直分量，則（x，z）處的uP（ω）和wP（ω）可表示為［20］

圖3 （a）縱波速度取不同值時，?u／?z的波形圖；（b）橫波速度取不同值時，?w／?z的波形圖；（c）縱波速度為3500m／s時，分離出的橫波；（d）橫波速度為1200m／s時，分離出的縱波Fig.3 （a）Waveforms of?u／?zon different P velocities；（b）Waveforms of?w／?zon different S velocities；（c）The separated S－waves using P velocity 3500m／s；（d）The separated P－waves using S velocity 1200m／s

其中，φu（ω）和φw（ω）分別表示單頻ω上，縱波水平分量以及垂直分量的振幅；kxP和kzP分別為水平方向和垂直方向上的波數(shù).

令UP和WP分別為平面縱波的水平分量和垂直分量

通過前面的分析可知，若要消除掉縱波，（15）、（16）兩式應(yīng)相等，于是有：

由（24）、（25）式可知，（?UP／?z）＊和?WP／?x的振幅比近似等于真實縱波速度與所給縱波速度之比.（?WS／?z）＊和?US／?x的振幅比近似等于真實橫波速度與所給橫波速度之比的負數(shù).因此利用（24）、（25）式，可以估算出地表的真實速度值.但這兩個式子是一個粗略的近似，可利用多次迭代提高速度估算的精度.具體的做法是（以估算縱波速度為例）：

（2）計算出（?UP／?z）＊和?WP／?x，求得兩者振幅之比；

（3）判斷振幅比是否滿足要求（接近于1），若滿足，停止計算，若不滿足由（24）式，計算出，重復(fù)步驟2；

同理，可估算出地表的橫波速度值.

在實際計算時，我們選取一個子波長度的時間窗口，該窗口內(nèi)的縱波子波或橫波子波與其它子波有很小的干涉.在該窗口內(nèi)計算（?UP／?z）＊和?WP／?x最大振幅之比或（?WS／?z）＊和?Us／?x最大振幅之比.

例如，對于上述地震地面記錄，給定初始的縱波速度值為3500m／s，經(jīng)兩次迭代后，求得（?UP／?z）＊和?WP／?x最大振幅之比為0.9628，利用（24）式求得縱波速度為3012.7m／s.再以此速度計算（?UP／?z）＊和?WP／?x最大振幅之比為0.9956，已經(jīng)非常接近于1，此時估算的縱波速度為2999.4m／s，與真實速度只差0.6m／s.同樣，當給定初始的橫波速度值為2200m／s時，經(jīng)過兩次迭代，求得（?WS／?z）＊和?US／?x最大振幅之比－0.9955，估算的橫波速度為1701.1m／s，與真實速度相差1.1m／s.圖4（a，b）分別為縱、橫波的初始速度為3500m／s以及1200m／s時分離的橫波和縱波，可以看到圖3c中殘留的縱波以及圖3d中殘留的橫波已經(jīng)消除掉.

本文也驗證了所給速度與真實速度相差非常大的情況.如當縱波初始速度給出8000m／s，經(jīng)過5次迭代，估算出的縱波速度為3021.7m／s；橫波初始速度給出6700m／s時，經(jīng)過4次迭代，估算出的橫波為1704.6m／s.分離的縱、橫波如圖4（c，d）所示.又如，當縱波初始速度為20m／s時，經(jīng)過5次迭代，估算出的縱波速度為3022m／s；橫波初始速度為10m／s時，經(jīng)過4次迭代，估算出的橫波速度為1704.6m／s.分離的縱、橫波如圖4（e，f）所示.從圖4（c－f）可以看出在所給出的縱、橫波速度與真實速度相差很大的情況下，也能夠很好地分離縱、橫波.

4 相位和縱橫波振幅比的校正

4.1 相位的校正

Sun等［19］指出在對波場進行散度和旋度計算后，縱、橫波的相位將產(chǎn)生π／2的移動.從（8）式和（9）式可以看出，對波場求取散度和旋度后，等式的右邊被乘上了一個i因子（i＝exp（iπ／2）），從而使分離后的P波和S波產(chǎn)生了π／2的相移.由于Sun是在空間域求取的散度和旋度，因此無法直接在計算的過程中對相位進行校正，需要在求取散度和旋度之后，將分離的縱、橫波分別通過對時間t做Hilbert變換［22］來進行校正.而本文是在頻率波數(shù)域求取的散度和旋度，因此可以在計算過程中將（8）式和（9）式兩邊同乘以－i（－i＝exp（－iπ／2）），直接將產(chǎn)生的π／2的相移校正回去，即：

4.2 縱、橫波振幅比的校正

Sun等［20］指出在對波場進行散度和旋度計算后，使得縱、橫波的振幅比產(chǎn)生了vP／vS的變化，需要在分離后的橫波波場乘以vS／vP，從而將縱、橫波振幅比校正到真實的縱、橫波振幅比上.而對于本文來說，由于在計算（8）式時，為了保證中不含有縱波，（10）式中的速度v取vP，在計算（9）式時，為了保證中不含有橫波，（10）式中的速度v取所以利用（8）式和（9）式計算出的S～和P～，再反變換回時間－空間域時，縱橫波的振幅比實際上是產(chǎn)生了vS／vP的變化（見附錄A）.因此，我們在分離后的橫波波場上乘以vP／vS來校正縱、橫波的振幅比.

圖5是做了相位與縱橫、波振幅比校正的縱、橫波波場圖.在成圖時將分離后的縱、橫波振幅進行了規(guī)格化［20］.通過圖5a與圖2a中的橫波（橫波能量主要集中在水平分量上），圖5b與圖2b中的縱波（縱波能量主要集中在垂直分量上）對比可以看到相位和振幅比得到了恢復(fù).

5 模型驗證

圖6為Sigsbee鹽丘模型.模型表層的縱波速度為1517m／s，橫波速度為876m／s.正演模擬時，震源放在地下50m，采用爆炸震源，檢波器置于地表，中間放炮，兩邊接收，模型四周使用了混合吸收邊界［23］.

圖7是利用圖6所示模型正演出的地面地震記錄.從圖中可以看出記錄中波場的主要能量是來自鹽丘表面的反射，淺部地層的反射能量很微弱.波場很復(fù)雜，縱、橫波混疊在一起，尤其是在1.5s以后，幾乎無法將兩者區(qū)分出來.

在進行縱、橫波分離時，假定地表速度未知，給出的初始縱波速度為3000m／s，初始橫波速度為1800m／s.利用上述估算速度值的方法，分別經(jīng)過10次和9次迭代，估算出的地表縱波速度為1517m／s，橫波速度為877m／s，與表層的縱、橫波速度非常接近.圖8（a，b）為最終分離的縱波和橫波，在頻率波數(shù)域計算散度和旋度時，使用的是（26）式和（27）式，并進行了縱、橫波振幅比的校正，以及規(guī)格化［19］.從圖8中可以看出，原本混疊在一起無法區(qū)分的縱、橫波得到了很好的分離.

圖6 Sigsbee鹽丘模型（a）縱波速度，從1517～4564m／s；（b）橫波速度，876～2635m／s.Fig.6 The sigsbee salt model（a）The P－wave velocity ranging from 1517to 4564m／s；（b）The S－wave velocity ranging from 876to 2635m／s.

6 結(jié) 論

本文依據(jù)Helmholtz分解定理，利用對地面地震記錄求散度和旋度的方法，分離出了縱、橫波.分析了地表縱、橫波的速度在波場分離中的重要性，針對速度值不準確會使縱、橫波分離不徹底的問題，提出了估算地表速度值的方法.利用此方法，可以在地表速度值未知時，將其估算出來.雖然（24）式和（25）式只是近似式，但通過鹽丘模型驗證可以看出，經(jīng)過幾次迭代后依然可以比較準確地估算出地表速度值，并且利用估算出的速度值很好地分離出了縱、橫波，說明了本文方法的有效性.估算速度值時，子波的選取很關(guān)鍵，選取比較純凈的縱波子波或橫波子波，估算的速度值才能更加接近真實值.根據(jù)本文波場分離的方法提出了相應(yīng)的相位校正以及縱、橫波振幅比校正的公式.在實際計算時，直接利用（26）和（27）兩式進行計算，不需要再額外進行相位校正，節(jié)省了計算成本.本文只對位移形式表示的波場進行了分離，對于速度形式表示的波場，本文的方法依然是適用的.

附錄A

這部分證明主要參考了文獻［20］.與其不同的是，在對波場求取完散度和旋度后，文獻［20］中的縱、橫波振幅比變?yōu)閏vP／vS，而本文中變?yōu)閏vS／vP.主要原因是，本文需要人為選取波數(shù)，而文獻［20］則不需要.

圖7 （a）模擬地震記錄的水平分量；（b）模擬地震記錄的垂直分量；（c）切除直達波后模擬地震記錄的水平分量；（d）切除直達波后模擬地震記錄的垂直分量.Fig.7 （a）Horizontal component of synthetic seismogram；（b）Vetical component of synthetic seismogram；（c）Horizontal component of synthetic seismogram after removing direct arrival；（d）Vertical component of synthetic seismogram after removing direct arrival.

令UP和US分別表示在坐標（x，z）處的縱波和橫波的位移矢量，則有

其中，φ0（ω）和θ0（ω）分別是頻率為ω時縱波和橫波的振幅；m和n是單位矢量，分別平行于縱波和橫波的偏振方向；kx，kz分別為縱波x方向和z方向上的波數(shù)；lx，lz分別為橫波x方向和z方向上的波數(shù).

令φ0（ω）和θ0（ω）滿足如下關(guān)系式：

c是常量.

彈性波波場U可以表示為縱波波場和橫波場的線性疊加：

圖8 分離出的縱波和橫波（a）分離出的縱波；（b）分離出的橫波.Fig.8 The separated P－and S－wavefield（a）The separated P－waves；（b）The separated S－waves.

同樣，在求旋度時，為保證得到純橫波，求取波數(shù)lz時，β應(yīng)該取縱波速度，因此l＝ω／vP，則

對比（A6）式和（A10）式可以看出，經(jīng)過了散度和旋度運算后，橫波與縱波的振幅比由c變成了cvS／vP.

（

）

［1］ 趙邦六.多分量地震勘探在巖性氣藏勘探開發(fā)中的應(yīng)用.石油勘探與開發(fā)，2008，35（4）：397－423.

Zhao B L.Application of multi－component seismic exploration in the exploration and production of lithologic gas reservoirs.PetroleumExplorationandDevelopment（in Chinese），2008，35（4）：397－423.

［2］ Davis T L.Multicomponent seismology：The next wave.Geophysics，2001，66（1）：49.

［3］ 孫歧峰，杜啟振.多分量地震數(shù)據(jù)處理技術(shù)研究現(xiàn)狀.石油勘探與開發(fā)，2011，38（1）：67－72.

Sun Q F，Du Q Z.A review of the multi－component seismic data processing.PetroleumExplorationandDevelopment（in Chinese），2011，38（1）：67－72.

［4］ 胡天躍，張廣娟，趙偉等.多分量地震波波場分解研究.地球物理學(xué)報，2004，47（3）：504－508.

Hu T Y，Zhang G J，Zhao W，et al.Decomposition of multicomponent seismic wavefileds.ChineseJ.Geophys.（in Chinese），2004，47（3）：504－508.

［5］ 李英康，崔作舟.分離縱波和橫波的偏振旋轉(zhuǎn)法.地球物理學(xué)報，1994，37（增刊1）：372－382.

Li Y K，Cui Z Z.P and S－waves separated by polarization revolving method.ActaGeophysicaSinica（ChineseJ.Geophys.）（in Chinese），1994，37（Supp.1）：372－382.

［6］ 馮晅，張先武，劉財?shù)?帶有多道相關(guān)的拋物線Radon變換法分離P－P、P－SV波.地球物理學(xué)報，2011，54（2）：304－309.Feng X，Zhang X W，Liu C，et al.Separating P－P and P－SV wave by parabolic Radon transform with multiple coherence.ChineseJ.Geophys.（in Chinese），2011，54（2）：304－309.

［7］ Devaney A J，Oristagliot M L.A plane－wave decomposition for elastic wave fields applied to the separation of P－waves and S－waves in vector seismic data.Geophysics，1986，51（2）：419－423.

［8］ Aki K，Richards P G.Quantitative Seismology（2nd ed）.California：University Science Books，2002.

［9］ Dellinger J，Etgen J.Wave－field separation in two－dimensional anisotropic media.Geophysics，1990，55（7）：914－919.

［10］ Dellinger J.Anisotropic Seismic Wave Propagation［Doctor′s thesis］.California：Stanford University，1991：65－69.

［11］ Yan J，Sava P.3Delastic wave mode separation for TTI media.79th Annual International Meeting，SEG，Expanded Abstracts，2009：4294－4298.

［12］ Yan J，Sava P.Elastic wave－mode separation for VTI media.Geophysics，2009，74（5）：WB19－WB32.

［13］ Zhang Q S，McMechan G A.2Dand 3Delastic wavefield vector decomposition in the wavenumber domain for VTI media.Geophysics，2002，75（3）：D13－D26.

［14］ Yan J，Sava P.Improving the efficiency of elastic wave－mode separation for heterogeneous tilted transverse isotropic media.Geophysics，2011，76（4）：T65－T78.

［15］ 郭鵬，何兵壽，張一凡.VTI介質(zhì)中的波場分離.工程地球物理學(xué)報，2011，8（3）：261－268.

Guo P，He B S，Zhang Y F.Wave－field Separation in VTI Media.ChineseJournalofEngineeringGeophysics（in Chinese），2011，8（3）：261－268.

［16］ 堯德中，周熙襄，鐘本善.VSP記錄的縱、橫波分離方法與應(yīng)用.石油地球物理勘探，1993，28（5）：623－628.

Yao D Z，Zhou X X，Zhong B S.Method for separating out P－wave or S－wave in VSP data，and its application.OGP（in Chinese），1993，28（5）：623－628.

［17］ Sun R.Separating P－and S－waves in a prestack 2－dimensional elastic seismogram.61th Ann.Mtg.，Eur.Assn.Geosci.Eng.，Extended Abstracts，1999：6－23.

［18］ Sun R，McMechan G A，Hsiao H H，et al.Separating Pand S－waves in prestack 3Delastic seismograms using divergence and curl.Geophysics，2004，69（1）：286－297.

［19］ Sun R，Chow J D，Chen K J.Phase correction in separating P－and S－waves in elastic data.Geophysics，2001，66（5）：1515－1518.

［20］ Sun R，McMechan G A，Chuang H H.Amplitude balancing in separating P－and S－waves in 2Dand 3Delastic seismic data.Geophysics，2011，76（3）：S103－S113.

［21］ 堯德中.單程彈性波逆時偏移和相移偏移方法.石油地球物理勘探，1994，29（4）：449－455，461.

Yao D Z.Reverse－time and phase－shift migrations of singleway elastic waves.OGP（in Chinese），1994，29（4）：449－455，461.

［22］ Claerbout J F.Fundamentals of Geophysical Data Processing：With Applications to Petroleum Prospecting.California：Blackwell Scientific Publications，1985：20－24.

［23］ Liu Y，Sen M K.A hybrid absorbing boundary condition for elastic wave modeling with staggered－grid finite difference.80th SEG meeting，Denver，Colorado，USA，Expanded Abstracts，2010：2945－2949.