定向延遲組合激發(fā)正演照明方法及應(yīng)用

馮浩宇， 俞海龍， 王 典， 劉 財(cái)

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

0 引言

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

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

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

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

1 方法原理

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

(1)

式中：u為聲波波場(chǎng)；v為聲波速度；t為時(shí)間；x、z為空間坐標(biāo)。

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

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

(2)

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

圖1為完全匹配層示意圖，內(nèi)部區(qū)域是實(shí)際地震波場(chǎng)正演模擬的區(qū)域，在該區(qū)域中x和z方向的衰減因子均為零，在PML區(qū)域1、2、3中衰減因子均大于零，陳可洋等[23]提出的余弦衰減因子，其表達(dá)式為式(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)格點(diǎn)到邊界的距離，從式(3)中可以看出，衰減因子隨著點(diǎn)到邊界距離的增加而逐漸增大，對(duì)應(yīng)圖1中的顏色逐漸加深。

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

(4)

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

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

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

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

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

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

(5)

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

(6)

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

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

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

(7)

圖5 組合震源不同時(shí)間延遲的地震波場(chǎng)理論方向圖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

當(dāng)相鄰兩震源間依次延時(shí)τ，則任意時(shí)刻引起的相位差為β=-2πfτ，由此引起的地震波主波束方向?yàn)椋?/p>

(8)

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

(9)

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

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

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

(10)

根據(jù)地震照明理論，組合震源在二維模型空間點(diǎn)(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)是炮點(diǎn)個(gè)數(shù)為9，分別位于960 m到1 040 m，炮點(diǎn)間隔為10 m，激發(fā)延時(shí)為0 ms、3 ms、5 ms的照明結(jié)果。圖7(d)、圖7(e)、圖7(f)為對(duì)應(yīng)的波場(chǎng)快照，從圖中可以看出，當(dāng)炮點(diǎn)無延時(shí)激發(fā)時(shí)，波場(chǎng)垂直地表向下傳播，能量延垂直地表向下分布，且隨著與炮點(diǎn)距離的增加，能量分布范圍逐漸變大；當(dāng)炮點(diǎn)激發(fā)有延遲時(shí)，隨著延遲時(shí)間的增加，能量傳播方向與垂向角度逐漸增大，且地層速度越高，對(duì)應(yīng)的角度越大。

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

圖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ā)延時(shí)照明結(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)是炮點(diǎn)個(gè)數(shù)為9，分別位于460 m到540 m，炮點(diǎn)間隔為10 m，從左到右相鄰炮點(diǎn)激發(fā)延時(shí)為7 ms的能量分布圖，波場(chǎng)延與地表垂向向右夾角大約為45。方向傳播，能量主要集中在此方向上，使其與斷層方向大致平行，圖8(c)是炮點(diǎn)個(gè)數(shù)為9，分別位于960 m到1 040 m，炮點(diǎn)間隔為10 m,從右到左相鄰炮點(diǎn)激發(fā)延時(shí)為7 ms的能量分布圖，波場(chǎng)延與地表垂向向左夾角大約為45。方向傳播，能量主要集中在此方向上，使能量盡可能多的分布在斷層所在的區(qū)域；圖8(d)是炮點(diǎn)個(gè)數(shù)為9，分別位于1 460 m到1 540 m，炮點(diǎn)間隔為10 m,從右到左相鄰炮點(diǎn)激發(fā)延時(shí)為8.7 ms的能量分布圖，波場(chǎng)延與地表垂向向左夾角大約為60.5。方向傳播，能量主要集中在此方向上，使能量盡可能多的分布在斷層所在的區(qū)域，從中可以看出，采用上述方法設(shè)計(jì)炮點(diǎn)位置，能量主要分布在斷層區(qū)域，可以獲得更多的該區(qū)域的反射信息，進(jìn)而能夠很好的對(duì)該區(qū)域進(jìn)行深入的研究。

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

2.3 marmousi模型

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

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

3 結(jié)論

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

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

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

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