劉煒，王彥春，畢臣臣，徐仲博

(1. 成都理工大學(xué) 地球物理學(xué)博士后科研流動站，四川 成都 610059；2.中國地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院，北京 100083)

0 引言

近年來，隨著油氣資源勘探開發(fā)程度的不斷加深，勘探對象也日益復(fù)雜化，因此薄互層、高陡構(gòu)造、微型構(gòu)造、隱蔽性油氣藏和巖性油氣藏等復(fù)雜地質(zhì)目標(biāo)的精確識別也愈發(fā)重要。相比于常規(guī)地面地震觀測方式，VSP觀測方式一般采用地表激發(fā)、井中接收的地震數(shù)據(jù)采集方式，因此VSP資料具有波場信息豐富、分辨率和信噪比高以及環(huán)境噪聲小等特點(diǎn)，有利于探明地下復(fù)雜地質(zhì)構(gòu)造以及提取儲層彈性參數(shù)等信息[1]。逆時偏移方法是一種基于雙程波波動方程的偏移成像方法，其不受地層傾角限制同時能夠適應(yīng)任意復(fù)雜速度區(qū)域，因此被認(rèn)為是目前成像精度最高的地震資料偏移成像方法[2]。VSP地震數(shù)據(jù)和逆時偏移方法的有機(jī)結(jié)合，有利于精確識別地下復(fù)雜構(gòu)造，因此研究VSP地震數(shù)據(jù)的逆時偏移方法具有一定的現(xiàn)實意義。

Whitemore于1983年首次提出了逆時偏移方法[3]，同年其他一些學(xué)者也對該方法進(jìn)行了詳細(xì)的探討[4-5]。經(jīng)過三十多年的發(fā)展，逆時偏移方法的理論體系逐漸完善，應(yīng)用對象和范圍也逐漸廣泛，出現(xiàn)了疊后逆時偏移和疊前逆時偏移、聲波逆時偏移[6]和彈性波逆時偏移[7]、各向同性介質(zhì)逆時偏移[8]、各向異性介質(zhì)逆時偏移[9]以及地面地震數(shù)據(jù)逆時偏移和VSP地震數(shù)據(jù)逆時偏移[10-12]等。逆時偏移的實現(xiàn)過程主要涉及到以下幾個關(guān)鍵方面：波場延拓、邊界條件、存儲策略、成像條件以及低頻噪聲壓制方法等。波場延拓是逆時偏移方法的核心，其實質(zhì)是離散求解波動方程，目前有限差分法最為常用[13]。在利用有限差分法求解波動方程時，提高其數(shù)值模擬精度是需要重點(diǎn)考慮的內(nèi)容，主要措施有增大差分階數(shù)和尋求優(yōu)化差分系數(shù)。根據(jù)頻散關(guān)系求解有限差分法的差分系數(shù)時，通常有泰勒級數(shù)展開法和最優(yōu)化方法兩種方法，一般基于后者所得到的差分系數(shù)的數(shù)值模擬精度更高。Dablain基于泰勒級數(shù)展開法求取了高階差分系數(shù)，有效地提高了有限差分法的數(shù)值模擬精度[14]。Liu和Sen利用泰勒級數(shù)法展開基于聲波方程的頻散關(guān)系分別獲得了傳統(tǒng)規(guī)則網(wǎng)格和交錯網(wǎng)格有限差分法的差分系數(shù)[15-16]。Liu采用最小二乘法優(yōu)化根據(jù)頻散關(guān)系建立的目標(biāo)函數(shù)分別求取了傳統(tǒng)規(guī)則網(wǎng)格和交錯網(wǎng)格有限差分法的優(yōu)化差分系數(shù)，進(jìn)一步提高了有限差分法的數(shù)值模擬精度[17-18]。在邊界條件和存儲策略方面，Clapp詳細(xì)探討了逆時偏移方法的存儲需求，提出了邊界存儲策略，有效地減少了震源波場存儲量[19]。在此基礎(chǔ)上，王保利等提出了有效邊界存儲策略，進(jìn)一步有效地減少了逆時偏移方法的存儲需求[20]。段沛然等提出了基于優(yōu)化算子邊界存儲策略的高效逆時偏移方法，其能夠兼顧逆時偏移方法的存儲量和計算時間[21]。在成像條件以及噪聲壓制方面，不少學(xué)者也進(jìn)行了大量的研究。常用的逆時偏移方法成像條件有：反褶積成像條件[22]、零延遲互相關(guān)成像條件[23]、震源歸一化零延遲互相關(guān)成像條件[24]等；常用的低頻噪聲壓制方法有：拉普拉斯濾波法[25]和高階拉普拉斯濾波法[26]等。

本文在二維變密度聲波波動方程的基礎(chǔ)上，基于優(yōu)化交錯網(wǎng)格有限差分法研究了VSP地震數(shù)據(jù)的高精度逆時偏移方法。優(yōu)化交錯網(wǎng)格有限差分法采用最小二乘法來確定其差分系數(shù)，因此相比于傳統(tǒng)的交錯網(wǎng)格有限差分法，當(dāng)差分算子長度相等時，它能夠更有效地壓制數(shù)值頻散，從而獲得更高精度的地震波場數(shù)值模擬結(jié)果。為了進(jìn)一步提高VSP數(shù)據(jù)逆時偏移方法的成像精度和實用性，采用PML吸收邊界條件和有效邊界存儲策略在有效地壓制地震波場邊界反射的同時極大地減少震源波場的存儲需求，采用震源歸一化零延遲互相關(guān)成像條件進(jìn)行成像計算，同時采用高階拉普拉斯濾波法壓制成像結(jié)果的低頻成像噪聲，進(jìn)而獲得高精度的成像結(jié)果。不同的模型測試結(jié)果表明：VSP地震數(shù)據(jù)逆時偏移方法可以實現(xiàn)對復(fù)雜地質(zhì)構(gòu)造的精確成像；由于VSP地震資料具有波場信息豐富、分辨率和信噪比高等特點(diǎn)，因此相比于常規(guī)地面地震數(shù)據(jù)的逆時偏移方法，VSP數(shù)據(jù)逆時偏移對井旁構(gòu)造、高陡構(gòu)造、微型構(gòu)造等復(fù)雜地質(zhì)構(gòu)造進(jìn)行成像更具有優(yōu)勢。

1 方法原理

逆時偏移方法的實現(xiàn)過程大體可以分為兩步：首先進(jìn)行震源波場的正向延拓和檢波點(diǎn)波場的逆向延拓，然后采用成像條件對震源波場和檢波點(diǎn)波場進(jìn)行成像計算，進(jìn)而得到最終的逆時偏移成像結(jié)果[27]。逆時偏移實現(xiàn)過程一般會涉及到以下幾個關(guān)鍵性問題：波場延拓、邊界條件、存儲策略、成像條件[28]以及低頻噪聲壓制方法等。

1.1 優(yōu)化交錯網(wǎng)格有限差分法

本文采用優(yōu)化交錯網(wǎng)格有限差分法離散求解變密度聲波波動方程實現(xiàn)逆時偏移方法的波場延拓過程。二維變密度聲波波動方程可以表示為[29]：

(1)

式中:ρ為介質(zhì)密度；p為波場值；K=ρv2為體積模量，v為介質(zhì)速度。

采用交錯網(wǎng)格有限差分法對式(1)左右兩邊進(jìn)行離散求解可得[30]：

(2)

式中：Δt為時間步長；h為空間網(wǎng)格步長，本文設(shè)定x和z方向的空間網(wǎng)格步長相等；M為差分算子長度，其值等于差分階數(shù)的一半；am為差分系數(shù)。

將式(2)代入式(1)，根據(jù)鏈?zhǔn)椒▌t并進(jìn)行簡化可得到二維變密度聲波波動方程時間2階、空間2M階的交錯網(wǎng)格差分格式[31]：

(3)

式中，r=vΔt/h；al和am為差分系數(shù)。

為了進(jìn)一步提高波場延拓的模擬精度，采用基于最小二乘法的優(yōu)化交錯網(wǎng)格有限差分法來計算每個時刻的波場值。該方法主要是采用最小二乘法優(yōu)化根據(jù)空間域頻散關(guān)系構(gòu)建得到的目標(biāo)函數(shù)從而得到優(yōu)化差分系數(shù)，其具體計算公式為[18]：

(n=2,3,…,M)

(5)

其中:

ψm(β)=2{sin[(m-0.5)β]-

2(m-0.5)sin(0.5β)}g(β)=β-2sin(0.5β)。(6)

式中：β=kh，k為波數(shù)；b為介于(0,π)之間的常數(shù)。

與傳統(tǒng)基于泰勒級數(shù)展開法的交錯網(wǎng)格有限差分法相比，基于最小二乘法的優(yōu)化交錯網(wǎng)格有限差分法可以進(jìn)一步壓制數(shù)值頻散，提高地震波場模擬精度，尤其是在大波數(shù)范圍內(nèi)。傳統(tǒng)交錯網(wǎng)格有限差分法的差分系數(shù)計算公式為[16,31]：

(7)

式中，m=1,2,…,M。為了比較分析本文優(yōu)化交錯網(wǎng)格有限差分法的優(yōu)越性，定義下列兩個參數(shù)[31]：

(8)

其中:

(9)

式中：δ(β,θ)為二維聲波波動方程數(shù)值模擬的相速度頻散；vFD為有限差分?jǐn)?shù)值求解的相速度；θ為平面波的傳播方向；ε(β,θ)為二維正演模擬時一個網(wǎng)格間距內(nèi)地震波傳播時間的相對誤差。由式(8)可知，當(dāng)δ越接近于1或ε越接近于0時，二維地震波場數(shù)值模擬的數(shù)值頻散就越小，即模擬精度越高。

為了比較分析傳統(tǒng)和優(yōu)化交錯網(wǎng)格有限差分法的數(shù)值頻散特征，設(shè)計一個均勻模型：介質(zhì)速度為 1 500 m/s，密度為2.0 g/cm3，模型大小為4 000 m×4 000 m，空間網(wǎng)格步長為20 m，時間步長為1 ms。圖1分別展示了不同差分算子長度時傳統(tǒng)和優(yōu)化交錯網(wǎng)格有限差分法的數(shù)值頻散誤差隨波數(shù)的變化規(guī)律。從圖中可以看出，隨著差分算子長度增加，基于優(yōu)化和傳統(tǒng)交錯網(wǎng)格有限差分法的地震波場數(shù)值模擬的頻散誤差均會減小，即模擬精度均會提高；在差分算子長度相等時，優(yōu)化交錯網(wǎng)格差分法的模擬精度優(yōu)于傳統(tǒng)交錯網(wǎng)格差分法。圖2分別展示了相同差分算子長度時傳統(tǒng)和優(yōu)化交錯網(wǎng)格有限差分法的頻散誤差隨地震波傳播方向的變化規(guī)律。從圖中可以看出，對于兩種交錯網(wǎng)格差分法，地震波場的模擬精度均與地震波的傳播方向有關(guān)；在地震波傳播方向相同時，優(yōu)化交錯網(wǎng)格差分法的模擬精度優(yōu)于傳統(tǒng)交錯網(wǎng)格差分法。設(shè)定震源子波是主頻為15 Hz的Ricker子波且位于模型正中心，地震記錄總時長為3 s，對該模型進(jìn)行正演模擬進(jìn)而對比分析兩種交錯網(wǎng)格差分法的模擬效果。圖3為不同差分算子長度時傳統(tǒng)和優(yōu)化交錯網(wǎng)格差分法在不同時刻的波場快照。從圖中可以看出，隨著差分算子長度增大，交錯網(wǎng)格差分法的數(shù)值頻散降低，模擬精度提高；在相同的差分算子長度條件下，優(yōu)化交錯網(wǎng)格差分法具有更高的數(shù)值模擬精度。

a—傳統(tǒng)交錯網(wǎng)格有限差分法;b—優(yōu)化交錯網(wǎng)格有限差分法a—conventional staggered-grid finite-difference;b—optimal staggered-grid finite-difference圖1 不同差分算子長度時交錯網(wǎng)格有限差分法的頻散誤差隨波數(shù)的變化規(guī)律Fig.1 Variations of dispersion error with wavenumber and different operator lengths by using different staggered-grid finite-difference methods

a—傳統(tǒng)交錯網(wǎng)格有限差分法;b—優(yōu)化交錯網(wǎng)格有限差分法a—conventional staggered-grid finite-difference;b—optimal staggered-grid finite-difference圖2 相同差分算子長度時(M=5)交錯網(wǎng)格有限差分法的頻散誤差隨傳播方向的變化規(guī)律Fig.2 Variations of dispersion error with different propagation angles by using different staggered-grid finite-difference methods (M=5)

a—傳統(tǒng)交錯網(wǎng)格有限差分法，M=5;b—優(yōu)化交錯網(wǎng)格有限差分法，M=5;c—傳統(tǒng)交錯網(wǎng)格有限差分法，M=10;d—優(yōu)化交錯網(wǎng)格有限差分法，M=10；從左至右依次為1 s和2.5 sa—conventional staggered-grid finite-difference for M=5;b—optimal staggered-grid finite-difference for M=5;c—conventional staggered-grid finite-difference for M=10;d—optimal staggered-grid finite-difference for M=10;from left to right time at 1s and 2.5s圖3 不同差分算子長度時交錯網(wǎng)格有限差分法在不同時刻的波場快照Fig.3 Snapshots by using different staggered-grid finite-difference methods and different operator lengths

算法穩(wěn)定性也是衡量有限差分?jǐn)?shù)值模擬方法優(yōu)劣的重要因素之一。交錯網(wǎng)格有限差分法的穩(wěn)定性條件可以表示為[16,31]：

(10)

式中s為穩(wěn)定性因子，其值越接近于1表示有限差分算法穩(wěn)定性越好。

利用上述均勻模型對本文交錯網(wǎng)格有限差分法進(jìn)行穩(wěn)定性分析。圖4給出了基于泰勒級數(shù)展開法的傳統(tǒng)交錯網(wǎng)格有限差分法以及基于最小二乘法的優(yōu)化交錯網(wǎng)格有限差分法的穩(wěn)定性因子曲線。由圖可知，優(yōu)化交錯網(wǎng)格有限差分法的穩(wěn)定性略差于傳統(tǒng)交錯網(wǎng)格有限差分法。通過以上分析可知，優(yōu)化交錯網(wǎng)格有限差分法在滿足穩(wěn)定性的條件下更加有利于壓制數(shù)值頻散，因此本文在后續(xù)的逆時偏移方法研究過程中主要采用時間2階、空間16階的優(yōu)化交錯網(wǎng)格有限差分法進(jìn)行波場延拓。

圖4 不同交錯網(wǎng)格有限差分法的穩(wěn)定性因子曲線Fig.4 Curves of stability factor of different staggered-grid finite-difference methods

1.2 PML邊界條件及有效邊界存儲策略

在地震波場的數(shù)值模擬過程中，由于計算機(jī)的計算區(qū)域有限，在模型四周存在的人工截斷邊界會造成強(qiáng)烈的邊界反射從而影響數(shù)值模擬結(jié)果，因此在地震波場的數(shù)值模擬過程中需采用有效的邊界條件來壓制邊界反射，從而提高地震波場數(shù)值模擬結(jié)果的準(zhǔn)確性。PML吸收邊界條件可以吸收任意方向、任意頻率的波，因此目前被廣泛應(yīng)用于地震波場的數(shù)值模擬過程[32-33]。

二維變密度聲波波動方程的PML吸收邊界條件的控制方程可以表示為：

(11)

式中，A為吸收衰減系數(shù)矩陣，其值從PML吸收邊界的內(nèi)邊界開始由內(nèi)到外依次增大，具體取值由吸收衰減因子決定。式(11)的交錯網(wǎng)格離散格式為：

(12)

PML吸收邊界條件的衰減因子分布情況如圖5所示。圖中區(qū)域E為實際模型區(qū)域，吸收衰減因子為零；區(qū)域A、B、C、D、F、G、H、I為吸收衰減區(qū)域，吸收衰減因子不為零。吸收衰減因子的種類繁多，本文采用余弦型衰減因子[34]，其具體表達(dá)式為：

(13)

式中：βx和βz分別表示沿x和z方向的余弦型衰減因子；B為衰減幅度因子；Lx和Lz分別表示x和z方

圖5 PML吸收邊界條件簡易示意Fig.5 Simple sketch of PML absorbing boundary condition

向完全匹配層的總網(wǎng)格數(shù)；lx和lz分別表示x和z方向距離完全匹配層外邊界的網(wǎng)格數(shù)。

PML吸收邊界條件能夠有效地壓制由人工截斷邊界造成的邊界反射，從而極大地削弱邊界反射對地震波場數(shù)值模擬結(jié)果中有效信息的影響，因此能夠有效地滿足逆時偏移成像的精度要求。然而，逆時偏移方法要求事先已知每個時刻的震源波場，存儲所有時刻的震源波場對計算機(jī)的存儲要求很高，一般計算機(jī)無法滿足。因此本文在PML吸收邊界條件的基礎(chǔ)上引進(jìn)有效邊界存儲策略，以求在保證波場延拓精度的條件下極大地減少逆時偏移震源波場的存儲需求。由圖5可知，在波場延拓的過程中，實際的正演模型包含內(nèi)部有效模型區(qū)域和外部吸收衰減邊界區(qū)域，即如圖6所示，區(qū)域A(深灰色)為內(nèi)部有效模型區(qū)域，區(qū)域B(白色)和區(qū)域C(淺灰色)為PML吸收邊界區(qū)域。有效邊界存儲策略只需存儲每個時刻區(qū)域C的震源波場值以及最后兩個時刻所有區(qū)域的震源波場值就可以精確地逆時重構(gòu)出所有時刻的震源波場，從而有效地降低逆時偏移的震源波場存儲需求。區(qū)域C的厚度與差分算子長度有關(guān)，在本文中波場延拓采用變密度聲波波動方程，因此區(qū)域C的網(wǎng)格點(diǎn)數(shù)等于差分算子長度的兩倍，即等于差分階數(shù)。

圖6 有效邊界存儲策略簡易示意Fig.6 Simple sketch of effective boundary storage strategy

為了驗證上述PML吸收邊界條件以及基于有效邊界存儲策略重構(gòu)震源波場的正確性和有效性，設(shè)計如圖7所示的多層層狀模型。該模型大小為4 000 m×4 000 m，空間網(wǎng)格步長為10 m，時間步長為1 ms，地震記錄總時長為3 s，PML吸收邊界的網(wǎng)格數(shù)為50，震源子波是主頻為15 Hz的Ricker子波且位于(2 000 m，1 800 m)處。圖8分別展示了該模型在不同時刻的正傳波場快照、應(yīng)用有效邊界存儲策略時所重構(gòu)的相對應(yīng)時刻的震源重構(gòu)波場快照以及它們二者之間的差異剖面。從圖中可以看出，在地震波場正向延拓過程中，當(dāng)?shù)卣鸩▊鞑ブ聊Ｐ瓦吔鐣r，地震波會繼續(xù)傳播到PML吸收邊界內(nèi)，然后PML吸收邊界對這些地震波進(jìn)行吸收衰減，因此由人工截斷邊界造成的邊界反射得到了極大地削弱，說明PML吸收邊界條件可以有效地壓制邊界反射，降低其對有效信號的影響。對比分析圖8a、圖8b和圖8c可以發(fā)現(xiàn)，在同一時刻震源的正傳波場和重構(gòu)波場的波形基本一致，二者之間的數(shù)值差異也很小，證明了采用有效邊界存儲策略在減小震源波場存儲需求的同時能夠根據(jù)所存儲的部分震源波場信息精確地重構(gòu)出震源波場。綜上可知，在逆時偏移方法中，綜合應(yīng)用PML吸收邊界條件以及有效邊界存儲策略，可以在有效壓制邊界反射的同時極大地降低震源波場的存儲需求。

圖7 多層層狀模型Fig.7 A multilayer model

a—正傳波場;b—重構(gòu)波場;c—正傳波場和重構(gòu)波場之間的差異;從左至右依次為0.4、0.8、1.2 s時刻a—forward wavefields;b—reconstructed wavefields;c—difference between the forward and reconstructed wavefields;from left to right time at 0.4,0.8,1.2 s圖8 多層層狀模型在不同時刻的波場快照Fig.8 Snapshots at different time for the multilayer model

1.3 成像條件

成像條件是影響逆時偏移成像精度的關(guān)鍵因素之一。目前，逆時偏移方法的成像條件種類繁多，如激發(fā)時間成像條件、最大振幅到達(dá)時成像條件、反褶積成像條件和零延遲互相關(guān)成像條件等。其中，零延遲互相關(guān)成像條件由于既能充分利用全部地震波場信息又能在一定程度上壓制成像噪聲，因此被廣泛應(yīng)用于逆時偏移方法中。其數(shù)學(xué)表達(dá)式為[23]：

(14)

式中：I(x,z)為成像值；S(x,z,t)為震源波場；R(x,z,t)為檢波點(diǎn)波場；Tmax為地震記錄總長度。

由式(14)可知，零延遲互相關(guān)成像條件通過對震源波場和檢波點(diǎn)波場進(jìn)行互相關(guān)計算從而獲得最終的逆時偏移成像結(jié)果，但是無法削弱震源對成像結(jié)果的影響，同時最終成像結(jié)果也無法直接反映地下地層的反射系數(shù)。在此基礎(chǔ)上，有學(xué)者提出了震源歸一化零延遲互相關(guān)成像條件，能夠有效地彌補(bǔ)零延遲互相關(guān)成像條件的缺陷，進(jìn)一步提高逆時偏移的成像精度。震源歸一化零延遲互相關(guān)成像條件的具體數(shù)學(xué)表達(dá)式為[24]：

(15)

對比式(14)和式(15)可知，震源歸一化零延遲互相關(guān)成像條件相當(dāng)于在零延遲互相關(guān)成像條件的基礎(chǔ)上除以震源波場的互相關(guān)結(jié)果，因此能夠有效地削弱震源對成像結(jié)果的影響，同時所得成像值可以直接對應(yīng)于地下地層的反射系數(shù)。

1.4 噪聲壓制

由于逆時偏移方法本身的原因，在其成像結(jié)果中不可避免地會產(chǎn)生低頻成像噪聲，為了提高成像精度，需采用有效方法對其進(jìn)行壓制。拉普拉斯濾波法是一種常用的噪聲壓制方法，相當(dāng)于將一個濾波器直接作用于逆時偏移成像數(shù)據(jù)，具有簡單易實施、噪聲壓制效果較好等特點(diǎn)。但是經(jīng)過一些學(xué)者的研究證明，常規(guī)拉普拉斯濾波法無法完全地壓制逆時偏移的成像噪聲，在去噪后的剖面上仍然殘留部分低頻噪聲。在前人的研究基礎(chǔ)上，本文采用高階拉普拉斯濾波法來壓制逆時偏移的低頻成像噪聲，其表達(dá)式為[26]：

式中：N=4,6,…,2n，即為大于2的偶數(shù)。該公式的簡易離散形式可以表示為：

式中：f(x,y)表示二維數(shù)據(jù)在(x,y)的數(shù)值。式(17)的上角標(biāo)并非表示指數(shù)，而是表示拉普拉斯算子的階數(shù)。

利用傅里葉變換將式(16)轉(zhuǎn)換到波數(shù)域可得：

(18)

式中：ω為角頻率，θ為地震波入射角。從式(18)可以看出，拉普拉斯類濾波法相當(dāng)于對逆時偏移成像數(shù)據(jù)進(jìn)行角度域衰減，其衰減程度與拉普拉斯算子的階數(shù)有關(guān)，不同階數(shù)的拉普拉斯濾波法僅衰減系數(shù)存在差異，分析可知高階拉普拉斯濾波法對大角度的逆時偏移成像噪聲有更好的壓制效果，但是階數(shù)過大也會加大損害逆時偏移結(jié)果中的有效信號，因此本文令N值等于4。

2 模型測試

為了驗證本文VSP地震數(shù)據(jù)逆時偏移方法的有效性，采用兩個模型進(jìn)行測試，分別為二維SEG/EAGE鹽丘模型和Marmousi模型。

2.1 二維SEG/EAGE鹽丘模型

二維SEG/EAGE鹽丘模型的速度和密度如圖9所示，其模型大小為3 380 m×2 100 m，空間網(wǎng)格步長為10 m，時間步長為1 ms，震源為15 Hz的Ricker子波，地震記錄總長度為4 s。對于該模型，采用常規(guī)地面地震觀測方式和VSP觀測方式進(jìn)行逆時偏移方法測試，分別對比分析了它們的逆時偏移成像結(jié)果。常規(guī)地面地震觀測系統(tǒng)的參數(shù)為：震源均勻地分布在地表，炮間距為30 m，共113炮；檢波器同樣均勻地分布在地表，道間距為10 m，共338道。VSP觀測系統(tǒng)參數(shù)為：震源分布方式與常規(guī)地面地震觀測系統(tǒng)一致，但是其檢波器則均勻地分布在井口坐標(biāo)為(0 m，0 m)和(3 380 m，0 m)的兩口井中，道間距為10 m，兩口井共420道。針對該SEG/EAGE鹽丘模型，主要分析討論了零延遲互相關(guān)成像條件與震源歸一化零延遲互相關(guān)成像條件的成像效果、高階拉普拉斯濾波法的低頻成像噪聲壓制效果(包括其對成像結(jié)果振幅和相位的影響)以及常規(guī)地面地震數(shù)據(jù)逆時偏移和VSP數(shù)據(jù)逆時偏移的差異。

圖9 二維SEG/EAGE鹽丘模型Fig.9 2D SEG/EAGE salt model

圖10顯示了分別應(yīng)用零延遲互相關(guān)成像條件和震源歸一化零延遲互相關(guān)成像條件時該二維SEG/EAGE鹽丘模型第57炮的VSP數(shù)據(jù)逆時偏移成像結(jié)果。從圖中可以看出，在基于零延遲互相關(guān)成像條件的逆時偏移結(jié)果中,震源對其成像結(jié)果的影響較大，而在基于震源歸一化零延遲互相關(guān)成像條件的逆時偏移結(jié)果中,震源對其成像結(jié)果的影響得到了極大的削弱，說明震源歸一化零延遲互相關(guān)成像條件能夠有效地削弱震源對逆時偏移成像結(jié)果的影響，從而提高其成像結(jié)果的精度。圖11給出了應(yīng)用震源歸一化零延遲互相關(guān)成像條件時采用高階拉普拉斯濾波法進(jìn)行噪聲壓制前后該二維SEG/EAGE鹽丘模型的地面地震數(shù)據(jù)和VSP數(shù)據(jù)逆時偏移結(jié)果。從圖中可以看出，在噪聲壓制前，常規(guī)地面地震數(shù)據(jù)和VSP地震數(shù)據(jù)的逆時偏移結(jié)果中均存在明顯的低頻成像噪聲，采用高階拉普拉濾波法進(jìn)行噪聲壓制后，二者的成像剖面均變得更為清晰，即成像精度提高，說明高階拉普拉濾波法能夠有效地壓制逆時偏移結(jié)果的低頻成像噪聲。同時，相比于常規(guī)的地面地震數(shù)據(jù)逆時偏移結(jié)果，VSP數(shù)據(jù)的逆時偏移結(jié)果能夠更加精確地反映地下介質(zhì)的構(gòu)造形態(tài)和構(gòu)造特征，尤其是高陡構(gòu)造部位以及井旁部位，說明VSP數(shù)據(jù)的逆時偏移方法更有利于精確地識別地下復(fù)雜地質(zhì)構(gòu)造。另外，由拉普拉斯算子的原理可知，拉普拉斯類濾波法實質(zhì)是對成像結(jié)果的數(shù)據(jù)進(jìn)行角度域濾波處理，其會對成像結(jié)果的振幅和相位產(chǎn)生一定的破壞。圖12給出了應(yīng)用高階拉普拉斯濾波法前后該二維SEG/EAGE鹽丘模型VSP逆時偏移結(jié)果的幅值譜和相位譜，其橫縱坐標(biāo)均表示空間頻率。從圖中可以看出，高階拉普拉斯濾波法可以壓制逆時偏移結(jié)果的成像噪聲，噪聲的振幅得到了有效的削弱，但是同時也可以發(fā)現(xiàn)濾波前后逆時偏移成像結(jié)果中有效信號的幅值譜和相位譜發(fā)生了改變，說明高階拉普拉濾波法會在一定程度上損害逆時偏移結(jié)果的振幅信息和相位信息。

a—零延遲互相關(guān)成像條件;b—震源歸一化零延遲互相關(guān)成像條件a—cross correlation imaging condition;b—normalized cross correlation imaging condition of sources圖10 二維SEG/EAGE鹽丘模型的第57炮VSP數(shù)據(jù)逆時偏移結(jié)果Fig.10 RTM results of the 57th VSP gather for the 2D SEG/EAGE salt model

a—噪聲壓制前幅度譜;b—噪聲壓制后幅度譜;c—噪聲壓制前相位譜;d—噪聲壓制后相位譜a—amplitude spectrum before noise suppression;b—amplitude spectrum after noise suppression;c—phase spectrum before noise suppression;d—phase spectrum after noise suppression圖12 二維SEG/EAGE鹽丘模型VSP逆時偏移結(jié)果的幅值譜和相位譜Fig.12 Amplitude spectrum and phase spectrum of VSP RTM results for the 2D SEG/EAGE salt model

2.2 Marmousi模型

為了進(jìn)一步驗證本文VSP逆時偏移方法對于復(fù)雜速度模型的有效性，采用Marmousi模型對其進(jìn)行測試。圖13給出了該Marmousi模型的速度和密度參數(shù)，其模型大小為7 000 m×7 000 m，空間網(wǎng)格步長為20 m，時間步長為1 ms，震源子波為15 Hz的Ricker子波，地震記錄總長度為6 s。對于該Marmousi模型，本文同樣對比分析了其常規(guī)地面地震數(shù)據(jù)和VSP地震數(shù)據(jù)的逆時偏移結(jié)果。兩種觀測系統(tǒng)的震源均以60 m為間距均勻分布在地表，共117炮，常規(guī)地面地震觀測系統(tǒng)的檢波器以20 m為間距均勻分布在地表，共350道，而VSP觀測系統(tǒng)的檢波器則以20 m為間距均勻地分布在井口坐標(biāo)為(0 m，0 m)和(7 000 m，0 m)的兩口井中，兩口井共700道。針對該Marmousi模型，首先分析討論了兩種不同的互相關(guān)成像條件的成像效果、高階拉普拉斯濾波法的低頻成像噪聲壓制效果、常規(guī)地面地震數(shù)據(jù)逆時偏移和VSP數(shù)據(jù)逆時偏移的差異，另外還進(jìn)一步討論了逆時偏移方法對偏移速度的敏感性。

圖14展示了分別應(yīng)用零延遲互相關(guān)成像條件和震源歸一化零延遲互相關(guān)成像條件時該Marmousi模型第59炮的VSP數(shù)據(jù)逆時偏移成像結(jié)果。從圖中可以看出，與二維SEG/EAGE鹽丘模型的單炮逆時偏移結(jié)果類似，震源歸一化零延遲互相關(guān)成像條件極大地削弱了震源對VSP數(shù)據(jù)逆時偏移成像結(jié)果的影響，從而提高了VSP數(shù)據(jù)單炮逆時偏移成像結(jié)果的精度。

圖13 Marmousi模型Fig.13 Marmousi model

圖15顯示了應(yīng)用震源歸一化零延遲互相關(guān)成像條件時采用高階拉普拉斯濾波法進(jìn)行噪聲壓制前后該Marmousi模型的地面地震數(shù)據(jù)和VSP數(shù)據(jù)逆時偏移結(jié)果。從圖中可以看出，與二維SEG/EAGE鹽丘模型的最終逆時偏移結(jié)果類似，無論是常規(guī)地面地震數(shù)據(jù)的逆時偏移結(jié)果還是VSP數(shù)據(jù)的逆時偏移結(jié)果，高階拉普拉濾波法均有效地壓制了逆時偏移結(jié)果中的低頻成像噪聲，提高了逆時偏移成像結(jié)果的精度。與此同時，VSP數(shù)據(jù)的逆時偏移結(jié)果比常規(guī)地面地震數(shù)據(jù)的逆時偏移結(jié)果更加精確，能夠更加有效地預(yù)測地下復(fù)雜地質(zhì)構(gòu)造。綜上可知，VSP數(shù)據(jù)的逆時偏移方法更加有利于識別和預(yù)測地下復(fù)雜地質(zhì)目標(biāo)。

a—零延遲互相關(guān)成像條件;b—震源歸一化零延遲互相關(guān)成像條件a—cross correlation imaging condition;b—normalized cross correlation imaging condition of sources圖14 Marmousi模型的第59炮VSP數(shù)據(jù)逆時偏移結(jié)果Fig.14 RTM results of the 59th VSP gather for the Marmousi model

a—地面數(shù)據(jù)低頻噪聲壓制前;b—VSP數(shù)據(jù)低頻噪聲壓制前;c—地面數(shù)據(jù)低頻噪聲壓制后;d—VSP數(shù)據(jù)低頻噪聲壓制后a—surface data before low frequency noise suppression;b—VSP data before low frequency noise suppression;c—surface data after low frequency noise suppression;d—VSP data after low frequency noise suppression圖15 Marmousi模型的常規(guī)地面地震數(shù)據(jù)和VSP數(shù)據(jù)的逆時偏移結(jié)果Fig.15 RTM results of conventional surface seismic data and VSP data for the Marmousi model

為了分析本文逆時偏移方法對速度的敏感性，將該Marmousi模型的速度模型和密度模型進(jìn)行平滑處理，其結(jié)果如圖16a所示，然后利用該平滑后的模型進(jìn)行逆時偏移，其結(jié)果如圖16b所示。從圖中可以看出，當(dāng)偏移速度模型不準(zhǔn)確時，常規(guī)地面地震數(shù)據(jù)逆時偏移結(jié)果以及VSP數(shù)據(jù)逆時偏移結(jié)果的精確度均會降低，尤其是在速度變化劇烈的部位，說明偏移速度是影響逆時偏移結(jié)果精度的重要因素之一。同時可以發(fā)現(xiàn)，此時VSP數(shù)據(jù)的逆時偏移結(jié)果仍然比地面地震數(shù)據(jù)的逆時偏移結(jié)果更為精確。

綜上說明，即使偏移速度出現(xiàn)細(xì)微誤差，在相同的條件下，VSP數(shù)據(jù)逆時偏移也能更為精確地識別地下地質(zhì)構(gòu)造。

a—速度;b—密度;c—地面地震數(shù)據(jù)逆時偏移結(jié)果;d—VSP數(shù)據(jù)逆時偏移結(jié)果a—velocity;b—density;c—RTM results of conventional surface seismic data;d—RTM results of VSP data圖16 平滑后的Marmousi模型及其逆時偏移結(jié)果Fig.16 The smoothed Marmousi model and its corresponding RTM results

3 結(jié)論與認(rèn)識

在前人研究的基礎(chǔ)上，研究了基于優(yōu)化交錯網(wǎng)格有限差分法的VSP數(shù)據(jù)逆時偏移方法。針對逆時偏移的不同方面，采用了不同的策略和措施。對于波場延拓，采用基于最小二乘法的優(yōu)化交錯網(wǎng)格有限差分法進(jìn)行波場的數(shù)值模擬，有效提高了數(shù)值模擬精度；采用PML吸收邊界條件和有效邊界存儲策略相結(jié)合的方式，在極大壓制邊界反射的同時有效地降低逆時偏移震源波場的存儲需求，轉(zhuǎn)而采用波場重構(gòu)的方式精確地重構(gòu)震源波場。在成像的過程中，應(yīng)用震源歸一化零延遲互相關(guān)成像條件進(jìn)行VSP數(shù)據(jù)的成像過程，有效地削弱震源對逆時偏移成像結(jié)果的影響。最后，采用高階拉普拉斯濾波法壓制逆時偏移成像結(jié)果的低頻噪聲，提高成像結(jié)果的精度。通過不同的模型測試，得到了以下幾點(diǎn)結(jié)論和認(rèn)識：

1)相比于傳統(tǒng)的基于泰勒級數(shù)展開法的交錯網(wǎng)格有限差分方法，基于最小二乘法的優(yōu)化交錯網(wǎng)格有限差分法能夠有效地提高地震波場的數(shù)值模擬精度，尤其是在大波數(shù)范圍內(nèi)；交錯網(wǎng)格有限差分法的數(shù)值模擬精度隨著差分算子長度的增大而提高，隨著地震波傳播方向的變化而發(fā)生改變。在相同的條件下，優(yōu)化交錯網(wǎng)格有限差分法的模擬效果優(yōu)于傳統(tǒng)交錯網(wǎng)格有限差分法。

2)PML吸收邊界條件能夠有效壓制由于計算空間有限造成的邊界反射；在此基礎(chǔ)上，應(yīng)用有效邊界存儲策略可以有效降低逆時偏移方法震源波場的存儲需求，其重構(gòu)震源波場與正傳震源波場基本一致，滿足高精度逆時偏移的精度要求。

3)震源歸一化零延遲互相關(guān)成像條件能夠有效地削弱震源對逆時偏移結(jié)果的影響，同時高階拉普拉斯濾波法能夠有效壓制低頻成像噪聲，從而提高逆時偏移成像結(jié)果的精度，但是高階拉普拉斯濾波法會在一定程度上破壞逆時偏移結(jié)果的振幅信息和相位信息。

4)相比于傳統(tǒng)的地面地震數(shù)據(jù)逆時偏移，VSP數(shù)據(jù)的逆時偏移能夠更加精確地識別地下復(fù)雜地質(zhì)構(gòu)造，尤其是井旁構(gòu)造、高陡構(gòu)造、微型構(gòu)造以及速度變化劇烈構(gòu)造等。