馮浩宇， 俞海龍， 王 典， 劉 財

(吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，長春 130026)

0 引言

隨著地震勘探程度的不斷深入以及地震勘探方法技術(shù)的發(fā)展，人們對地震數(shù)據(jù)采集的質(zhì)量、數(shù)據(jù)處理方法和數(shù)據(jù)解釋的合理性具有越來越高的要求，其中野外地震數(shù)據(jù)的質(zhì)量尤為重要，要想采集到準確的、高質(zhì)量的野外地震數(shù)據(jù)，就應(yīng)該針對不同的地質(zhì)目標體，設(shè)置恰當?shù)挠^測方式，采用合理的采集參數(shù)[1]，才能更多的獲取到所要研究的目標體的地質(zhì)信息，對于那些地質(zhì)構(gòu)造復(fù)雜的地區(qū)，應(yīng)用地震定向照明分析技術(shù)，可以有效的解決在地震數(shù)據(jù)采集中所存在的問題。

積分方程法、射線追蹤法和波動方程法是目前地震數(shù)值模擬三種常用的方法，射線追蹤法由于缺少地震波動力學(xué)信息以及高頻近似和奇點問題，當?shù)刭|(zhì)構(gòu)造變得復(fù)雜時會產(chǎn)生很大的計算誤差，積分方程法由于受幾何形態(tài)地影響較大，使其難以在實際中得到應(yīng)用，而同時用到了地震波運動學(xué)和動力學(xué)信息的波動方程的地震數(shù)值模擬方法，能夠很好地適應(yīng)橫向速度強烈變化介質(zhì)，與實際中所觀測到的地震波場非常接近[2]。因此，我們進行聲波的波動方程數(shù)值模擬[3-4]，進而實現(xiàn)基于聲波波動方程的定向照明分析。

地震照明分析是通過地震波場數(shù)值模擬的方法來研究地震波在地下介質(zhì)的傳播規(guī)律和能量分布的定量分析方法[5]，通常采用對地震波場的角度域分解或震源直接激發(fā)產(chǎn)生定向地震波場兩種方法來實現(xiàn)定向照明分析，對于波場的角度域分解方法，Wu等[6]利用Radon 變換方法分解平面波并獲取了波場的角度域信息，Xie等[7]利用小波變換分解波場的方法來獲得角域信息，并將地震波在地下介質(zhì)中的能量分布作為照明度，用于地震定向照明分析，裴正林[8]根據(jù)地震波能量理論，對地震波場進行方向性分解，并給出了地震波定向照明度的計算方法，計算出了基于聲波方程的震源定向照明度和源-檢定向照明度，Mao等[9]提出了一種基于局部指數(shù)標架的局部角度域地震定向照明分析方法，通過局部余弦和正弦變換的線性組合將不同發(fā)射接收系統(tǒng)的空間域波場分解為局部角度域，該方法相比于傳統(tǒng)的方法具有更高的效率，耿瑜等[10]在局部指數(shù)標架角度域方向照明的基礎(chǔ)上，將方向照明分析應(yīng)用在鹽丘下部目標區(qū)域上，研究采集系統(tǒng)對成像質(zhì)量的影響，從而達到提高成像質(zhì)量的目的，馮偉[11-12]、陳生昌[13]、謝小碧[15]、Tao等[14]也采用不同的變換方法對地震波場進行角度域分解，實現(xiàn)了對地下介質(zhì)的定向照明分析。對于直接產(chǎn)生定向地震波場方法，王忠仁等[16]利用有限差分法對相控震源chirp信號掃描的地震響應(yīng)進行了數(shù)值模擬，得到了相比于常規(guī)組合激發(fā)方法更好的信號，為了提高原始地震數(shù)據(jù)的信噪比，姜弢等[17-18]提出了一種基于可控震源陣列的定向照明控制方法，Tomas等[19]對組合震源激發(fā)原理進行了系統(tǒng)的總結(jié)，并給出了其在提高地震分辨率的應(yīng)用，蔡紀琰等[20]利用4種常用的震源進行了波動方程有限差分正演模擬，從理論上研究了不同激發(fā)方式在震源藥量相同的情況下的波場傳播特征，以及所對應(yīng)的能量照明關(guān)系，鞏向博等[21]給出了在地表起伏條件下，組合震源所產(chǎn)生的定向地震波場的計算方法，秦龍[22]提出通過增加炮點的數(shù)目來提高地震照明的均勻性[22]。

對于角度域分解法，由于需要對地震波場逐點進行角度域分解以便獲得角度信息，因此計算量大、計算效率低，對于直接產(chǎn)生定向地震波場方法，大多采用一階壓力-速度聲波方程，相比于二階聲波方程，涉及公式多、編程實現(xiàn)難度大、計算效率低，鑒于此筆者直接基于二階聲波波動方程，利用規(guī)則網(wǎng)格有限差分技術(shù)對聲波波場進行正演模擬，通過改變組合震源的數(shù)目來控制地震波主波束的集中程度、通過改變組合震源之間的激發(fā)時間，來調(diào)整地震波主波束的方向，然后計算地震波場能量在地下的分布情況，進而實現(xiàn)定向照明分析。

1 方法原理

基于二維聲波方程采用時間2階空間10階的規(guī)則網(wǎng)格有限差分進行正演模擬，在二維空間域內(nèi)，聲波方程為式(1)。

(1)

式中：u為聲波波場；v為聲波速度；t為時間；x、z為空間坐標。

由于實際正演模擬的區(qū)域是有限的，當?shù)卣鸩▊鞑サ竭吔鐣r會發(fā)生強烈的邊界反射，嚴重影響著計算結(jié)果的準確度，完全匹配層(Perfectly Matched Layer，簡稱PML)吸收邊界條件能能夠很好的消除邊界反射的影響，其基本思想是在所研究區(qū)域的邊界上加上吸收衰減層，使地震波傳播到邊界時不會發(fā)生反射，并且在吸收區(qū)域地震波會快速衰減，以此來減弱邊界反射所帶來的影響。

我們采用的是王維紅等[3]提出的完全匹配層吸收邊界條件，其給出的二維聲波方程完全匹配層的控制方程為式(2)。

(2)

式中，A(x,y)為吸收衰減因子，在實際正演模擬區(qū)域中，衰減因子A(x,y)為零，式(2)和式(1)完全相同，此時波場傳播滿足聲波方程，在PML區(qū)域中，衰減因子A(x,y)不為零，地震波無反射的進入該區(qū)域后會快速衰減。

圖1為完全匹配層示意圖，內(nèi)部區(qū)域是實際地震波場正演模擬的區(qū)域，在該區(qū)域中x和z方向的衰減因子均為零，在PML區(qū)域1、2、3中衰減因子均大于零，陳可洋等[23]提出的余弦衰減因子，其表達式為式(3)。

lx=0,1,2,…PML

lz=0,1,2,…PML

(3)

式中:B為衰減幅度因子;αz表示沿z方向即PML區(qū)域1中的衰減因子;αx表示沿x方向即PML區(qū)域2中的衰減因子;PML區(qū)域3中的衰減因子為αx+αz;PML表示吸收層厚度;lx和lz分別沿x和z方向網(wǎng)格點到邊界的距離，從式(3)中可以看出，衰減因子隨著點到邊界距離的增加而逐漸增大，對應(yīng)圖1中的顏色逐漸加深。

由此可得到基于PML邊界條件的時間2階空間高階(2K階)有限差分方程：

(4)

圖1 完全匹配層示意圖Fig.1 Schematic diagram of complete matching layer

圖2 未加邊界條件的波場快照Fig.2 Snapshot of wave field without boundary conditions

圖3 加邊界條件的波場快照Fig.3 Snapshot of wave field with boundary conditions

圖4 組合震源的方向性Fig.4 Directivity of the combined source

圖2、圖3分別是無邊界條件和加入PML吸收邊界條件的波場快照，從圖2、圖3中可以看出，當無邊界條件時，波場傳到邊界會發(fā)生強烈的邊界反射，對計算結(jié)果產(chǎn)生很大的影響，當加入PML吸收邊界條件時，波場傳到邊界時無邊界反射，具有很好地吸收效果。

如圖4所示，地表有n個震源，分別為O1、O2、…、On，相鄰兩個震源之間的間距為d且保持恒定的相位差β，S是距離震源較遠的位置處的一點，因此各震源到S點的射線可以看成是平行的，θ是射線方向與地表的夾角，由此可以得到S點來自相鄰兩個震源的地震波的相位差為ψ=kdcosθ+β，假設(shè)震源O1在S點產(chǎn)生的場強為E1，則震源O1、O2、…、On在S點產(chǎn)生的總場強為E=E1+E2+E3+…+En=E1(1+ejψ+ej2ψ+…+ej(n-1)ψ)=E1(sin(nψ/2)/sin(ψ/2))ejξ，式中ξ=((n-1)/2)ψ是總場強的相位角，與總場強的方向無關(guān)，因此可以取fa(θ)=sin(nψ/2)/sin(ψ/2)，我們將fa(θ)稱為組合震源地震波場的方向因子，當ψ=0時，fa(θ)取最大值，即

(5)

于是可以得到歸一化的組合震源地震波場的方向因子F為：

(6)

其中：n為震源的數(shù)目;k為波數(shù);d為相鄰震源的間距;θ為傳播角度;β為相鄰震源間的恒定相位差。

當相鄰兩震源間的相位差(kdcosθ+β)發(fā)生變化時，地震波主波束方向會隨之變化，主波束方向?qū)?yīng)最大能量指向，所以最大能量的指向可以通過改變恒定相位差β來控制

θmax=arccos(-β/dk)=arccos(-βv/πfd)

(7)

圖5 組合震源不同時間延遲的地震波場理論方向圖Fig.5 Seismic wave field theory with different time delays in combination of earthquake source(a)n=5,τ=0.001 5 s；(b)n=5,τ=0.003 5 s ；(c)n=5,τ=-0.002 5 s(d)n=10,τ=0.001 5 s；(e)n=10,τ=0.003 5 s；(f)n=10,τ=-0.002 5 s

當相鄰兩震源間依次延時τ，則任意時刻引起的相位差為β=-2πfτ，由此引起的地震波主波束方向為：

(8)

由式(8可以得到選定地震波方向后，對應(yīng)的延時時間為：

(9)

圖5是震源個數(shù)分別為5和10，震源間距為10 m，地下介質(zhì)速度為2 000 m/s，震源從左到右的延遲時間分別為0.001 5 s、0.003 5 s、-0.002 5 s所對應(yīng)的地震波方向圖。從圖5中可以看出，不同的延遲時間，對應(yīng)著不同的地震波主波束指向，且隨著震源數(shù)目的增多，主波束指向更加集中。

假設(shè)位于(xn,zn)處的單個震源激發(fā)產(chǎn)生的地震波場在地下任意一點(x,z)處的分布為Us(xn,zn;x,z,t)，那么多個震源組合激發(fā)產(chǎn)生的地震波場Us(θ;x,z,t)為式(10)。

圖6 水平層狀介質(zhì)模型Fig.6 Horizontal layered medium model

(10)

根據(jù)地震照明理論，組合震源在二維模型空間點(x,z)處的照明度為[17]式(11)。

(11)

2 數(shù)值算例

2.1 水平層狀介質(zhì)模型

圖6是四層水平層狀介質(zhì)模型，速度從上到分別為1 000 m/s、1 500 m/s、2 000 m/s、2 500 m/s，模型大小為2 000 m×2 000 m，圖7(a)、圖7(b)、圖7(c)是炮點個數(shù)為9，分別位于960 m到1 040 m，炮點間隔為10 m，激發(fā)延時為0 ms、3 ms、5 ms的照明結(jié)果。圖7(d)、圖7(e)、圖7(f)為對應(yīng)的波場快照，從圖中可以看出，當炮點無延時激發(fā)時，波場垂直地表向下傳播，能量延垂直地表向下分布，且隨著與炮點距離的增加，能量分布范圍逐漸變大；當炮點激發(fā)有延遲時，隨著延遲時間的增加，能量傳播方向與垂向角度逐漸增大，且地層速度越高，對應(yīng)的角度越大。

圖7 不同激發(fā)延時能量照明及波場快照Fig.7 Energy illumination and wave field snapshots with different excitation time delay(a)無時間延遲；(b)從左到右延遲3 ms；(c)從左到右延遲5 ms(d)無時間延遲波場；(e)從左到右延遲3 ms波場；(f)從左到右延遲5 ms波場

圖8 逆沖推覆構(gòu)造模型及照明結(jié)果Fig.8 Thrust nappe structure model and lighting results(a)逆沖推覆構(gòu)造模型；(b)從左到右延遲7 ms；(c)從右到左延遲7 ms；(d)從右到左延遲8.7 ms

圖9 Marmousi速度模型Fig.9 The marmousi model

圖10 不同激發(fā)延時照明結(jié)果Fig.10 Illumination results of different excitation delays(a)從左到右延遲5 ms； (b)從左到右延遲6 ms ；(c)從左到右延遲7.5 ms； (d)從左到右延遲9 ms

2.2 逆沖推覆構(gòu)造模型

圖8 (a)是逆沖推覆構(gòu)造模型，從上到下各層的速度分別為1 000 m/s、1 500 m/s、1 800 m/s、2 200 m/s、2 500 m/s，模型大小為2 070 m×1 270 m，本文所要研究的是斷層所在的區(qū)域，圖(b)是炮點個數(shù)為9，分別位于460 m到540 m，炮點間隔為10 m，從左到右相鄰炮點激發(fā)延時為7 ms的能量分布圖，波場延與地表垂向向右夾角大約為45。方向傳播，能量主要集中在此方向上，使其與斷層方向大致平行，圖8(c)是炮點個數(shù)為9，分別位于960 m到1 040 m，炮點間隔為10 m,從右到左相鄰炮點激發(fā)延時為7 ms的能量分布圖，波場延與地表垂向向左夾角大約為45。方向傳播，能量主要集中在此方向上，使能量盡可能多的分布在斷層所在的區(qū)域；圖8(d)是炮點個數(shù)為9，分別位于1 460 m到1 540 m，炮點間隔為10 m,從右到左相鄰炮點激發(fā)延時為8.7 ms的能量分布圖，波場延與地表垂向向左夾角大約為60.5。方向傳播，能量主要集中在此方向上，使能量盡可能多的分布在斷層所在的區(qū)域，從中可以看出，采用上述方法設(shè)計炮點位置，能量主要分布在斷層區(qū)域，可以獲得更多的該區(qū)域的反射信息，進而能夠很好的對該區(qū)域進行深入的研究。

圖11 目標層不同激發(fā)延時能量曲線Fig.11 Different excitation delay energy curves of the target layer

2.3 marmousi模型

圖9是二維復(fù)雜的marmousi速度模型，模型大小為1 220 m×3 840 m，模型上部有三條大的斷裂，中部有一鹽丘，深部靠近模型的左右兩側(cè)有兩個高速異常體，在深度1000 m ～1 100 m、水平2 500 m ～3 000 m位置處存在一低速目標層，網(wǎng)格點數(shù)為122×384，網(wǎng)格間距為10 m，時間步長為0.5 ms，采樣時間為2.0 s。圖10是采用雷克子波(主頻為30 Hz)，炮點個數(shù)為9，分別位于1 380 m到1 540 m，炮點間隔為20 m，從左到右相鄰炮點激發(fā)延時間分別為5 ms、6 ms、7.5 ms、9 ms的能量分布圖，所要研究的目標區(qū)域為圖中黑色虛線所在的低速區(qū)域，從圖9中可以看出，在相同震源數(shù)目的情況下，組合震源不同的延遲時間，目標區(qū)域的能量分布是不同的，當延遲時間為6 ms時，目標區(qū)域的能量分布明顯高于延遲時間為5 ms、7.5 ms和9 ms所對應(yīng)的能量分布，因此，在實際野外勘探時，為了得到目標區(qū)域而更多的反射信息，獲取高信噪比的數(shù)據(jù)，延遲時間應(yīng)設(shè)置為6 ms。

圖11是在深度為1 040 m、水平方向為2 500 m～3 000 m處的組合震源不同激發(fā)延時的能量曲線，從圖11中可以看出，延遲時間為6 ms時目標層具有最高的能量分布，其余能量分布由大到小的延遲時間依次為5 ms、7.5 ms和9 ms，這也進一步說明延遲時間為6 ms時，可以獲取更多的與目標層相關(guān)的反射信息，能夠得到最佳的勘探效果。

3 結(jié)論

1)基于二階標量聲波波動方程，利用規(guī)則網(wǎng)格有限差分技術(shù)，對聲波波場進行正演模擬，相比于一階壓力-速度方程容易編程實現(xiàn)、計算效率高。

2)采用組合震源的激發(fā)方式，通過改變相鄰兩個震源之間的激發(fā)延遲時間來產(chǎn)生定向地震波場，進而實現(xiàn)組合震源定向照明分析，從層狀介質(zhì)模型試驗中可以看出，延遲時間越大，能量傳播方向與垂向夾角越大，地層速度越高對應(yīng)的角度也越大。

3)將定向照明方法應(yīng)用于復(fù)雜的逆沖推覆構(gòu)造模型和marmousi速度模型中，對于不同的炮點位置，通過改變激發(fā)時間來實現(xiàn)對目標區(qū)域的照明分析，當炮點位置相同時，可以通過調(diào)整激發(fā)時間來獲取最佳的照明效果。

4)筆者是在二維水平地表的情況進行研究的，當?shù)乇砥鸱鼤r該方法不再適用，此后會將該方法推廣到起伏地表和三維模型中，以便更好地滿足實際勘探的要求。