基于優(yōu)化交錯(cuò)網(wǎng)格有限差分法的VSP逆時(shí)偏移

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

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

0 引言

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

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

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

1 方法原理

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

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

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

(1)

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

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

(2)

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

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

(3)

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

為了進(jìn)一步提高波場(chǎng)延拓的模擬精度，采用基于最小二乘法的優(yōu)化交錯(cuò)網(wǎng)格有限差分法來計(jì)算每個(gè)時(shí)刻的波場(chǎng)值。該方法主要是采用最小二乘法優(yōu)化根據(jù)空間域頻散關(guān)系構(gòu)建得到的目標(biāo)函數(shù)從而得到優(yōu)化差分系數(shù)，其具體計(jì)算公式為[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)基于泰勒級(jí)數(shù)展開法的交錯(cuò)網(wǎng)格有限差分法相比，基于最小二乘法的優(yōu)化交錯(cuò)網(wǎng)格有限差分法可以進(jìn)一步壓制數(shù)值頻散，提高地震波場(chǎng)模擬精度，尤其是在大波數(shù)范圍內(nèi)。傳統(tǒng)交錯(cuò)網(wǎng)格有限差分法的差分系數(shù)計(jì)算公式為[16,31]：

(7)

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

(8)

其中:

(9)

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

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

a—傳統(tǒng)交錯(cuò)網(wǎng)格有限差分法;b—優(yōu)化交錯(cuò)網(wǎng)格有限差分法a—conventional staggered-grid finite-difference;b—optimal staggered-grid finite-difference圖1 不同差分算子長(zhǎng)度時(shí)交錯(cuò)網(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)交錯(cuò)網(wǎng)格有限差分法;b—優(yōu)化交錯(cuò)網(wǎng)格有限差分法a—conventional staggered-grid finite-difference;b—optimal staggered-grid finite-difference圖2 相同差分算子長(zhǎng)度時(shí)(M=5)交錯(cuò)網(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)交錯(cuò)網(wǎng)格有限差分法，M=5;b—優(yōu)化交錯(cuò)網(wǎng)格有限差分法，M=5;c—傳統(tǒng)交錯(cuò)網(wǎng)格有限差分法，M=10;d—優(yōu)化交錯(cuò)網(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 不同差分算子長(zhǎng)度時(shí)交錯(cuò)網(wǎng)格有限差分法在不同時(shí)刻的波場(chǎng)快照Fig.3 Snapshots by using different staggered-grid finite-difference methods and different operator lengths

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

(10)

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

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

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

1.2 PML邊界條件及有效邊界存儲(chǔ)策略

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

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

(11)

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

(12)

PML吸收邊界條件的衰減因子分布情況如圖5所示。圖中區(qū)域E為實(shí)際模型區(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吸收邊界條件簡(jiǎn)易示意Fig.5 Simple sketch of PML absorbing boundary condition

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

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

圖6 有效邊界存儲(chǔ)策略簡(jiǎn)易示意Fig.6 Simple sketch of effective boundary storage strategy

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

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

a—正傳波場(chǎng);b—重構(gòu)波場(chǎng);c—正傳波場(chǎng)和重構(gòu)波場(chǎng)之間的差異;從左至右依次為0.4、0.8、1.2 s時(shí)刻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 多層層狀模型在不同時(shí)刻的波場(chǎng)快照Fig.8 Snapshots at different time for the multilayer model

1.3 成像條件

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

(14)

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

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

(15)

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

1.4 噪聲壓制

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

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

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

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

(18)

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

2 模型測(cè)試

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

2.1 二維SEG/EAGE鹽丘模型

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

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

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

2.2 Marmousi模型

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

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

圖13 Marmousi模型Fig.13 Marmousi model

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

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

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

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

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

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

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

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

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

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