舒濤， 楊鍇, 王華忠

同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院， 上海 200092

0 引言

在現(xiàn)代地震勘探理論框架下，全(寬)波數(shù)帶的速度(或其他彈性參數(shù))的估算是處理技術(shù)的核心.全波形反演(Full waveform inversion, 簡(jiǎn)記為FWI.Tarantola，1984；Pratt，1999)試圖估算出全(寬)波數(shù)帶的速度場(chǎng).但由于數(shù)據(jù)不完備，正演算子不能很好地模擬實(shí)測(cè)波場(chǎng)，地震震源子波未知等原因，經(jīng)典的FWI得不到預(yù)期的結(jié)果.基于分波型或時(shí)空局部特征波場(chǎng)的層析成像與疊前深度偏移成像組合是全波形反演成像方法(FWI)實(shí)用化的途徑之一(王華忠等，2015).而上述方法技術(shù)的關(guān)鍵前提是如何獲得具有時(shí)空局部特征的地震數(shù)據(jù).

目前利用高斯束或高斯波包的特征波成像和反演的方法主要是基于數(shù)據(jù)域進(jìn)行地震數(shù)據(jù)的分解.然而基于數(shù)據(jù)域分解的方法需要對(duì)數(shù)據(jù)進(jìn)行遍歷，通常效率比較低.此外，基于數(shù)據(jù)域的分解在物理上不夠直觀，分解效果受地震數(shù)據(jù)的信噪比影響較大.因此，如何快速有效地分解地震數(shù)據(jù)仍然是高斯波包實(shí)用化面臨的主要挑戰(zhàn)之一.

不同于傳統(tǒng)的基于數(shù)據(jù)域的分解方法，本文提出了一種基于像域分解的特征高斯波包數(shù)據(jù)模擬方法.該方法首先對(duì)成像剖面作Gabor分解，利用具有局部空間特征的Gabor函數(shù)特征表達(dá)地下的擾動(dòng)界面，然后應(yīng)用高斯波包的散射波模擬實(shí)現(xiàn)特征高斯波包數(shù)據(jù)模擬.該方法的特點(diǎn)和優(yōu)勢(shì)在于，一方面，在成像剖面上作分解，物理上更加直觀，效率也更高，利用具有局部空間特征的Gabor函數(shù)可以特征表達(dá)地下的擾動(dòng)界面，從而能夠從數(shù)據(jù)中分離不同的反射同相軸；另一方面，基于高斯波包模擬的一次反射波數(shù)據(jù)具有時(shí)空局部特征，包含了地震數(shù)據(jù)的更多屬性，比如同相軸的到達(dá)時(shí)、空間位置和斜率，為特征波成像和特征波層析反演提供了強(qiáng)有利的數(shù)據(jù)支持.為了更好地用Gabor函數(shù)表達(dá)地質(zhì)界面，本文用公式建立了波數(shù)域表征的Gabor函數(shù)和空間域表征的Gabor函數(shù)之間的關(guān)系，并分析了Gabor函數(shù)中各個(gè)參數(shù)代表的物理意義.基于像域分解的Gabor函數(shù)表達(dá)的擾動(dòng)，經(jīng)過高斯波包的散射波模擬，能較好的表達(dá)地震數(shù)據(jù).通過數(shù)值實(shí)驗(yàn)對(duì)本文提出的方法進(jìn)行了驗(yàn)證，實(shí)驗(yàn)結(jié)果表明了該方法的有效性.

1 理論

1.1 空間域和波數(shù)域表征的Gabor函數(shù)

二維空間域表征的Gabor函數(shù)：

(1)

x′=-xsinθ+zcosθ,(2)

z′=-xcosθ-zsinθ,(3)

其中，x和z表示相對(duì)控制點(diǎn)的空間位置，A表示振幅，λ表示正弦函數(shù)的波長(zhǎng)，θ表示Gabor核函數(shù)的方向，ψ表示相位偏移(本文假設(shè)相位偏移為零)，σ表示高斯函數(shù)的標(biāo)準(zhǔn)差，γ表示空間方向的展布比例，又稱為橢圓率.

取公式(1)的實(shí)數(shù)部分，結(jié)合公式(2)和(3)，可得二維空間域表征的Gabor函數(shù)的詳細(xì)計(jì)算公式：

(4)

掌握不同參數(shù)對(duì)Gabor函數(shù)形態(tài)的影響是用Gabor函數(shù)特征表達(dá)地下的擾動(dòng)模型的前提和關(guān)鍵，有助于提高像域分解的有效性和精確性.為了便于直觀地認(rèn)識(shí)和分析不同的Gabor函數(shù)參數(shù)對(duì)于Gabor函數(shù)形態(tài)的影響，圖1—4分別展示了λ、θ、σ、γ不同取值下的Gabor函數(shù)形態(tài).如圖1所示，正弦函數(shù)的波長(zhǎng)λ影響的是垂直Gabor核函數(shù)方向θ的頻率，波長(zhǎng)越大，頻率越小，震蕩衰減越小.圖2展示了不同Gabor核函數(shù)方向θ下Gabor函數(shù)形態(tài)，本文Gabor核函數(shù)方向θ取水平向右為正，隨順時(shí)針方向逐漸增大.圖3表示的是高斯函數(shù)的標(biāo)準(zhǔn)差σ對(duì)Gabor函數(shù)形態(tài)的影響，可以看到標(biāo)準(zhǔn)差σ越大，Gabor函數(shù)展布的空間范圍也越大，這跟理論上的分析是一致的.圖4展示了空間方向上的比例γ對(duì)Gabor函數(shù)形態(tài)的影響，空間方向上的比例在這里定義為Gabor函數(shù)在平行Gabor核函數(shù)方向θ與垂直Gabor核函數(shù)方向θ的空間展布比例.

二維波數(shù)域表征的Gabor函數(shù)：

g(x;A,k,K)=Aexp(ikTx-0.5xTKx),(5)

其中，x表示相對(duì)控制點(diǎn)的空間位置，A表示振幅，k表示Gabor函數(shù)的波數(shù)向量，K表示Gabor函數(shù)的衰減系數(shù)矩陣.

取公式(5)的實(shí)部并將其展開，可得二維波數(shù)域表征的Gabor函數(shù)的詳細(xì)計(jì)算公式：

+2Kxzxz)]cos(kxx+kzz).

(6)

通過對(duì)比二維空間域表征的Gabor函數(shù)公式(4)和二維波數(shù)域表征的Gabor函數(shù)公式(6)，可以建立不同域表征的Gabor函數(shù)之間參數(shù)的對(duì)應(yīng)關(guān)系：

(7)

(8)

(9)

(10)

(11)

圖1 正弦函數(shù)的波長(zhǎng)λ對(duì)Gabor函數(shù)形態(tài)的影響(a) λ=33.564; (b) λ=67.128; (c) λ=134.256.其他參數(shù)分別為θ=0°，σ=80，γ=1.Fig.1 The influence of the wavelength of sine function on the shape of Gabor function(a) λ=33.564; (b) λ=67.128; (c) λ=134.256 with θ=0°，σ=80，γ=1.

圖2 Gabor核函數(shù)的方向θ對(duì)Gabor函數(shù)形態(tài)的影響(a) θ=0°; (b) θ=45°; (c) θ=90°.其他參數(shù)分別為λ=67.128，σ=20，γ=1/6.325.Fig.2 The influence of the direction of Gabor kernel function on the shape of Gabor function(a) θ=0°; (b) θ=45°; (c) θ=90° with λ=67.128，σ=20，γ=1/6.325.

圖3 高斯函數(shù)的標(biāo)準(zhǔn)差σ對(duì)Gabor函數(shù)形態(tài)的影響(a) σ=20; (b) σ=40; (c) σ=80. 其他參數(shù)分別為λ=67.128，θ=0°，γ=1.Fig.3 The influence of standard deviation of Gaussian function on the shape of Gabor function(a) σ=20; (b) σ=40; (c) σ=80 with λ=67.128，θ=0°，γ=1.

圖4 空間方向上的比例γ對(duì)Gabor函數(shù)形態(tài)的影響(a) γ=1/1; (b) γ=1/2; (c) γ=1/4.其他參數(shù)分別為λ=67.128，θ=0°，σ=80.Fig.4 The influence of the proportion of spatial scale on the shape of Gabor function(a) γ=1/1; (b) γ=1/2; (c) γ=1/4 with λ=67.128，θ=0°，σ=80.

圖5 基于像域分解的特征高斯波包數(shù)據(jù)模擬的流程圖Fig.5 Flow chart of characteristic Gaussian wave packet data modeling based on image domain decomposition

通過公式(7—11)，參數(shù)向量k和參數(shù)矩陣K的物理意義更加明確.需要指出的是，本文采用基于高斯束與高斯波包的Gabor框架散射波模擬方法(李輝和王華忠，2015)來模擬擾動(dòng)波場(chǎng).因此，上述關(guān)系的構(gòu)建有助于我們將空間域表征的Gabor函數(shù)變換到波數(shù)域表征的Gabor函數(shù)，從而用上述方法實(shí)現(xiàn)擾動(dòng)波場(chǎng)的模擬.

1.2 基于像域分解的思路

本文基于像域分解的思路如圖5所示.首先，通過地震偏移獲得成像剖面，成像剖面的質(zhì)量會(huì)直接影響特征界面的拾取.然后，在成像剖面上拾取特征界面并計(jì)算對(duì)應(yīng)的地層傾角，用于重新表達(dá)地下擾動(dòng)模型時(shí)確定Gabor函數(shù)的位置和核函數(shù)的方向.如何準(zhǔn)確且自動(dòng)化地拾取地質(zhì)界面是地球物理領(lǐng)域的研究熱點(diǎn)和難點(diǎn)，近些年發(fā)展了不少基于自動(dòng)化的層位拾取方法(Wu and Hale，2015；Wu and Fomel，2018; Yan and Wu，2021).對(duì)于特定震相的成像和走時(shí)反演而言，拾取成像剖面中所有的地質(zhì)特征是不必要的.本文側(cè)重于拾取出特征較為明顯且連續(xù)性較好的反射界面，而并不尋求準(zhǔn)確充分地拾取成像剖面中的所有地質(zhì)界面，這在一定程度上降低了拾取的難度和復(fù)雜度.本文采用Wang等(2015)提出的基于結(jié)構(gòu)張量的自動(dòng)拾取算法來拾取特征界面并計(jì)算對(duì)應(yīng)的地層傾角.

接著，由公式(4)計(jì)算空間域表征的Gabor函數(shù)以特征表達(dá)成像剖面.此環(huán)節(jié)中特征表達(dá)成像剖面的Gabor函數(shù)的參數(shù)和個(gè)數(shù)的選取至關(guān)重要，前者主要影響特征表達(dá)成像剖面的精度，后者主要影響高斯波包傳播的精度和效率.基于像域分解的原則是利用Gabor函數(shù)稀疏表達(dá)特征成像剖面，即在較為準(zhǔn)確的特征表達(dá)成像剖面的基礎(chǔ)上，盡可能減少Gabor函數(shù)的個(gè)數(shù).對(duì)于Gabor函數(shù)參數(shù)的選取，本文在1.1節(jié)的測(cè)試基礎(chǔ)上考慮到平行層理方向的地質(zhì)界面特征更為明顯，將Gabor函數(shù)中對(duì)地質(zhì)界面水平方向影響較小的參數(shù)，如正弦函數(shù)的波長(zhǎng)λ、相位偏移ψ、高斯函數(shù)的標(biāo)準(zhǔn)差σ設(shè)定為不變量，根據(jù)經(jīng)驗(yàn)分別取λ=67.128，ψ=0，σ=20.而Gabor函數(shù)中余下的參數(shù)如相對(duì)控制點(diǎn)的空間位置x和z，振幅A，Gabor核函數(shù)的方向θ，空間方向的展布比例γ，分別根據(jù)選取的層位位置，對(duì)應(yīng)位置成像的能量大小，對(duì)應(yīng)層位的地層傾角以及線性度等來確定.利用盡可能少的Gabor函數(shù)稀疏表達(dá)特征成像剖面則可以看作是高斯波包傳播的精度和效率之間的一種折中，這可以當(dāng)作一個(gè)稀疏約束下的反問題進(jìn)行求解，找到滿足精度和效率的最優(yōu)解.

最后，基于Gabor函數(shù)特征表達(dá)的成像剖面進(jìn)行高斯波包擾動(dòng)波場(chǎng)的模擬，獲得基于高斯波包模擬的地震數(shù)據(jù).不同于波動(dòng)方程模擬，基于高斯波包模擬的地震數(shù)據(jù)其計(jì)算效率更高，內(nèi)存占用量更少，此外基于高斯波包重新模擬的地震數(shù)據(jù)包含的信息更豐富，比如反射波的到達(dá)時(shí)、空間位置和波場(chǎng)的傳播方向，可以用于后續(xù)的特征反射波成像和反演.

2 數(shù)值實(shí)驗(yàn)

在本節(jié)中，我們將通過斷層模型和Marmousi模型來驗(yàn)證基于像域分解的特征高斯波包數(shù)據(jù)模擬的可行性和有效性，最后將本文方法用于部分Sigsbee 2A模型的高斯波包走時(shí)層析反演實(shí)驗(yàn)以說明其應(yīng)用效果.

2.1 斷層模型實(shí)驗(yàn)

首先，利用一個(gè)斷層模型測(cè)試基于像域分解的特征高斯波包數(shù)據(jù)模擬方法的有效性.斷層模型橫向和縱向的網(wǎng)格大小為501×301，橫向和縱向的網(wǎng)格間距均為10 m，如圖6a所示.在地表均勻布置41炮，炮間距為100 m，炮點(diǎn)橫向范圍500～4500 m.每炮501道地表全接收，道間距為10 m，接收記錄長(zhǎng)度3.0 s，時(shí)間采樣間隔1 ms.震源信號(hào)是主頻為20 Hz的雷克子波.基于上述觀測(cè)系統(tǒng)進(jìn)行波動(dòng)方程有限差分正演，獲得的地震記錄作為觀測(cè)數(shù)據(jù)并用于逆時(shí)偏移成像(RTM)，得到的成像剖面如圖6b所示.基于結(jié)構(gòu)張量算法進(jìn)行層位拾取(見圖6d)并計(jì)算對(duì)應(yīng)的地層傾角(見圖6c).利用拾取的層位和地層傾角信息可以基于公式(4)計(jì)算空間域表征的Gabor函數(shù)以特征表達(dá)成像剖面，如圖7所示.對(duì)比圖6b，可發(fā)現(xiàn)基于Gabor函數(shù)特征表達(dá)的成像剖面很好地恢復(fù)了原成像剖面中的主要反射同相軸，僅在近地表附近略有不同.RTM成像剖面在近地表附近存在較多的偏移假象，這是由于成像疊加干涉不干凈導(dǎo)致的.而基于Gabor函數(shù)來表達(dá)的特征成像剖面僅僅保留原成像剖面中連續(xù)性較好和能量較強(qiáng)的反射同相軸，并不追求完整地表達(dá)原成像剖面中的所有信息，如無地質(zhì)意義的成像噪聲.

圖6 斷層模型(a) 速度模型； (b) 基于圖6(a)得到的RTM成像剖面； (c) 基于結(jié)構(gòu)張量算法計(jì)算的地層傾角場(chǎng)； (d) 基于結(jié)構(gòu)張量算法拾取的層位.Fig.6 Fault model(a) Velocity model; (b) RTM imaging profile based on Fig.6 (a); (c) The formation inclination field calculated based on the structure tensor algorithm; (d) The layer picked up by the structure tensor algorithm.

圖7 基于Gabor函數(shù)重新表達(dá)的擾動(dòng)界面Fig.7 Perturbation interface reconstruction based on Gabor function

最后，基于圖7所示的擾動(dòng)界面進(jìn)行高斯波包擾動(dòng)波場(chǎng)模擬，獲得地表記錄的擾動(dòng)波場(chǎng).如圖8b所展示的單炮記錄，炮點(diǎn)橫向位置為2500 m處.同波動(dòng)方程有限差分法模擬結(jié)果(圖8a)相比，基于像域分解的高斯波包模擬的地震數(shù)據(jù)完整地恢復(fù)了主要的反射同相軸，且反射同相軸特征明顯，連續(xù)性強(qiáng)，噪聲也相對(duì)較少.二者的抽道對(duì)比(見圖9)可以看到高斯波包模擬的擾動(dòng)波場(chǎng)同波動(dòng)方程有限差分法模擬得到的擾動(dòng)波場(chǎng)十分接近，僅在振幅上略有差異.總的來說，二者具有較好的一致性.測(cè)試結(jié)果證明了本文提出的基于像域分解的特征高斯波包數(shù)據(jù)模擬方法的可行性.

2.2 Marmousi模型實(shí)驗(yàn)

為驗(yàn)證本文方法對(duì)于復(fù)雜模型的適應(yīng)性，在此對(duì)Marmousi模型進(jìn)行測(cè)試.圖10展示了Marmousi模型的速度場(chǎng)，模型橫向和縱向的網(wǎng)格點(diǎn)數(shù)為501×351，橫向和縱向的網(wǎng)格間距均為10 m.在地表均勻布置51炮，炮間距為100 m，炮點(diǎn)橫向范圍0～5000 m.每炮501道接收，道間距為10 m.震源信號(hào)是主頻為20 Hz的雷克子波.基于上述觀測(cè)系統(tǒng)進(jìn)行波動(dòng)方程有限差分正演，獲得的地震記錄作為觀測(cè)數(shù)據(jù)并用于逆時(shí)偏移成像(RTM)，得到的成像剖面如圖11所示.利用拾取的層位和地層傾角信息可以基于公式(4)計(jì)算空間域表征的Gabor函數(shù)以特征表達(dá)成像剖面，如圖12所示.重新表達(dá)后的成像剖面雖然沒有完全匹配RTM成像剖面，但是其質(zhì)量更高，噪聲更少，主要的反射界面特征更明顯，便于在數(shù)據(jù)域中獲得與主要特征反射界面對(duì)應(yīng)的反射同相軸.

圖8 斷層模型合成數(shù)據(jù)的分解(a) 波動(dòng)方程有限差分法模擬的單炮數(shù)據(jù)，作為原始數(shù)據(jù)(切除了直達(dá)波)；(b) 基于像域分解的特征高斯波包模擬的地震數(shù)據(jù). 炮點(diǎn)橫向位置為2500 m處.Fig.8 Decomposition of synthetic data of fault model(a) Single shot data simulated by wave equation finite difference method, as original data (direct wave is removed); (b) Single shot data simulated by characteristic Gaussian wave packet based on image domain decomposition. The horizontal position of shot is 2500 m.

圖9 地震記錄抽道對(duì)比(a) 地表1.0 km處； (b) 地表2.5 km處； (c) 地表4.0 km處. 紅線代表波動(dòng)方程有限差分法；藍(lán)線代表基于像域分解的特征高斯波包.Fig.9 A comparison of three normalized traces of shot gathers in Fig.8 at (a) 1.0 km; (b) 2.5 km; (c) 4.0 km. Red curve: finite difference method; blue curve: characteristic Gaussian wave packet based on image domain decomposition

圖10 Marmousi速度模型Fig.10 The velocity field of Marmousi model

圖11 RTM成像剖面Fig.11 RTM imaging profile

圖12 基于Gabor函數(shù)重新表達(dá)的擾動(dòng)界面Fig.12 Perturbation interface reconstruction based on Gabor function

最后，基于圖12所示的擾動(dòng)界面做高斯波包擾動(dòng)波場(chǎng)模擬，獲得地表記錄的擾動(dòng)波場(chǎng).圖13展示了波動(dòng)方程有限差分法模擬和基于像域分解的高斯波包模擬的地震數(shù)據(jù)單炮記錄(炮點(diǎn)橫向位置為3000 m處)，二者在相位上具有較好的一致性，僅在振幅上略有差異.圖14展示了二者的抽道對(duì)比，高斯波包模擬的擾動(dòng)波場(chǎng)同波動(dòng)方程有限差分法模擬得到的擾動(dòng)波場(chǎng)在相位上較為一致，基本反映了數(shù)據(jù)的變化趨勢(shì)，僅在振幅上略有差異，表現(xiàn)出隨偏移距增大振幅差異也變大的規(guī)律.在計(jì)算上，同波動(dòng)方程有限差分法模擬相比，高斯波包模擬計(jì)算效率更高，內(nèi)存占用量更小，如表1.此外，基于高斯波包重新模擬的地震數(shù)據(jù)包含的信息更豐富，比如反射波的到達(dá)時(shí)、空間位置和波場(chǎng)的傳播方向，這在層析反演中用于數(shù)據(jù)測(cè)量和偏移成像等處理時(shí)更具優(yōu)勢(shì).

圖13 Marmousi模型合成單炮記錄的對(duì)比(a) 波動(dòng)方程有限差分法模擬的單炮數(shù)據(jù)； (b) 基于像域分解的特征高斯波包模擬的地震數(shù)據(jù). 炮點(diǎn)橫向位置為3000 m處.Fig.13 Comparison of single shot records simulated with Marmousi model(a) Single shot data simulated by wave equation finite difference method; (b) Seismic data represented by characteristic Gaussian wave packet based on image domain decomposition. The horizontal position of shot is 3000 m.

圖14 地震記錄抽道對(duì)比(a) 地表1.0 km處； (b) 地表2.0 km處； (c) 地表3.0 km處. 紅線代表波動(dòng)方程有限差分法；藍(lán)線代表基于像域分解的特征高斯波包.Fig.14 A comparison of three normalized traces of shot gathers in Fig.8 at (a) 1.0 km; (b) 3.0 km; (c) 5.0 km. Red curve: finite difference method; blue curve: characteristic Gaussian wave packet based on image domain decomposition.

表1 Marmousi模型不同正演方法計(jì)算效率與內(nèi)存占用量對(duì)比Table 1 Comparison of computational efficiency and memory consumption of different forward modeling methods for Marmousi model

2.3 反演實(shí)驗(yàn)

圖15 部分Sigsbee 2A速度模型Fig.15 The partial Sigsbee 2A velocity model

圖16 初始速度模型Fig.16 The initial velocity model

為進(jìn)一步驗(yàn)證本文方法的有效性，這里對(duì)不包含巖體構(gòu)造的Sigsbee 2A模型進(jìn)行測(cè)試并用于高斯波包走時(shí)層析反演.圖15展示了部分Sigsbee 2A模型的速度場(chǎng)，模型橫向和縱向的網(wǎng)格點(diǎn)數(shù)為601×301，橫向網(wǎng)格間距為10 m，縱向網(wǎng)格間距為10 m.在地表均勻布置51炮，炮間距為100 m，炮點(diǎn)橫向范圍500～5500 m.每炮601道接收，道間距為10 m.震源信號(hào)是主頻為20 Hz的雷克子波.初始速度為常梯度模型(見圖16).基于上述觀測(cè)系統(tǒng)進(jìn)行波動(dòng)方程有限差分正演，獲得的地震記錄作為觀測(cè)數(shù)據(jù)并用于逆時(shí)偏移成像(RTM)，得到的成像剖面如圖17所示.利用拾取的層位和地層傾角信息可以基于公式(4)計(jì)算空間域表征的Gabor函數(shù)以特征表達(dá)成像剖面，如圖18所示.最后，基于圖16所示的初始速度和圖18所示的擾動(dòng)界面做高斯波包擾動(dòng)波場(chǎng)模擬，獲得地表記錄的擾動(dòng)波場(chǎng)(圖19a).

圖17 基于初始模型的RTM成像剖面Fig.17 RTM imaging profile based on initial model

圖18 基于Gabor函數(shù)重新表達(dá)的擾動(dòng)界面Fig.18 Perturbation interface reconstruction based on Gabor function

由于地下真實(shí)速度未知，在速度不正確的情況下，高斯波包擾動(dòng)波場(chǎng)模擬的擾動(dòng)波場(chǎng)通常與地表接收到的地震記錄存在走時(shí)差，類似于反射波走時(shí)層析計(jì)算走時(shí)差，通過時(shí)差校正可以較好地消除相位差異.本文采用動(dòng)態(tài)規(guī)整算法(Hale，2013)計(jì)算走時(shí)差并對(duì)高斯波包擾動(dòng)波場(chǎng)記錄進(jìn)行了時(shí)差校正(圖19b).同基于真實(shí)速度的波動(dòng)方程有限差分法模擬結(jié)果(圖19c)相比，校正后的基于像域分解的高斯波包模擬的地震數(shù)據(jù)同相軸特征明顯，連續(xù)性強(qiáng).圖20展示了抽道對(duì)比情況，基于初始速度進(jìn)行高斯波包模擬的擾動(dòng)波場(chǎng)經(jīng)過時(shí)差校正后，同波動(dòng)方程有限差分法模擬得到的擾動(dòng)波場(chǎng)具有較好的一致性，基本反映了數(shù)據(jù)的變化趨勢(shì)，僅在振幅上略有差異，這與高斯波包模擬擾動(dòng)波場(chǎng)的精度有關(guān)，在大角度反射時(shí)誤差相對(duì)較大.但是，同波動(dòng)方程有限差分法模擬相比，高斯波包模擬計(jì)算效率更高，包含的信息更豐富，用于反射波走時(shí)層析反演更為直接和高效.

圖19 單炮記錄的對(duì)比(a) 基于初始速度(圖16)和擾動(dòng)界面(圖18)進(jìn)行高斯波包模擬獲得的單炮記錄； (b) 基于(a)進(jìn)行時(shí)差校正后的單炮記錄; (c) 基于真實(shí)模型進(jìn)行聲波方程正演獲得的單炮記錄. 炮點(diǎn)橫向位置為3000 m處.Fig.19 Single shot records obtained by Gaussian wave packet simulation based on initial velocity and disturbance interface; (b) Single shot record after time difference correction based on (a); (c) Single shot record obtained by forward modeling of acoustic equation based on real model. The horizontal position of shot is 3000 m.

將基于像域分解的特征高斯波包模擬的51炮地震數(shù)據(jù)用于特征高斯波包走時(shí)層析，初始速度為常梯度模型(見圖16).圖21展示了特征高斯波包走時(shí)層析迭代20輪的反演結(jié)果，其與真實(shí)速度模型的光滑背景部分吻合程度較好.進(jìn)一步對(duì)比反演結(jié)果，抽取不同地表位置(X=1000 m，3000 m, 5000 m)的縱向速度，如圖22所示，可以看到，反演結(jié)果與真實(shí)速度模型的變化趨勢(shì)較為吻合，較好地恢復(fù)了中-大尺度的背景速度結(jié)構(gòu).該實(shí)驗(yàn)表明了利用高斯波包重新模擬的地震數(shù)據(jù)在背景速度反演中的實(shí)用性.以上測(cè)試結(jié)果進(jìn)一步證明了本文提出的基于像域分解的特征高斯波包數(shù)據(jù)模擬方法的有效性.

圖21 反演結(jié)果Fig.21 The inversion result

圖22 不同地表位置(X=1000 m, 3000 m, 5000 m)處的速度抽線對(duì)比藍(lán)色曲線，綠色曲線和紅色曲線分別為初始速度模型，反演結(jié)果和真實(shí)速度模型.Fig.22 A comparison of the initial velocity model (blue curve)， inversion result (green curve) and true velocity model (red curve) at X=1000 m, 3000 m, 5000 m

3 結(jié)論

本文提出了一種基于像域分解的特征高斯波包數(shù)據(jù)模擬方法.該方法首先利用具有空間局部特征的Gabor函數(shù)對(duì)成像剖面上進(jìn)行特征分解，然后通過高斯波包擾動(dòng)波模擬獲得高斯波包模擬的能量較強(qiáng)、連續(xù)性較好的一次反射波數(shù)據(jù)，最后用于特定震相的反演與成像.本文的方法不僅可以利用成像剖面地質(zhì)意義更為直觀的特點(diǎn)來拾取特征反射界面，還能有效地克服數(shù)據(jù)域分解精度和效率不高的問題，保證了分解后的反射數(shù)據(jù)具有更加明確的物理意義.在實(shí)際問題中即便速度不正確的情況下，通過時(shí)差校正本文方法也同樣有效.不同于波動(dòng)方程有限差分模擬，除了模擬效率的提升和內(nèi)存占用量的減少，基于時(shí)空局部特征的高斯波包模擬的地震數(shù)據(jù)還包含更多的數(shù)據(jù)屬性，例如反射波的到達(dá)時(shí)、空間位置和波場(chǎng)的傳播方向等，為特征波成像和特征波層析反演提供了有利的數(shù)據(jù)支持，用于反射波走時(shí)反演更為直接和高效.數(shù)值實(shí)驗(yàn)測(cè)試表明了本文提出的分解方法的有效性，將像域分解的特征高斯波包數(shù)據(jù)用于反射波層析反演能夠獲得較好的背景速度模型.

致謝感謝審稿專家對(duì)本文提出的寶貴意見和建議.感謝華為技術(shù)有限公司的李輝博士，同濟(jì)大學(xué)的劉玉柱教授、馮波研究員和吳成梁博士在論文撰寫過程中的指導(dǎo)與建議.