鄒 強(qiáng) 黃建平* 李慶洋 雍 鵬 辛天亮

(①中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院地球物理系，山東青島 266580；②中國石化中原油田分公司物探研究院，河南濮陽 457001)

0 引言

隨著地震勘探技術(shù)和計算機(jī)性能的發(fā)展，彈性波數(shù)值模擬與成像技術(shù)在勘探領(lǐng)域得到了越來越多的關(guān)注[1-3]。利用彈性波偏移算法對多分量數(shù)據(jù)進(jìn)行建模和成像，可以提供更多的地下地質(zhì)構(gòu)造信息[4-5]，如驗(yàn)證亮點(diǎn)反射、評估介質(zhì)參數(shù)和檢測裂縫等[6-7]。

地震波場外推通常采用有限差分法[8-12]，利用恒定間距的笛卡爾坐標(biāo)網(wǎng)格近似波動方程的時間和空間偏導(dǎo)數(shù)，該方法具有編程簡單、運(yùn)算速度快、能夠分析處理各種復(fù)雜地質(zhì)構(gòu)造中的波動問題等優(yōu)點(diǎn)[13]。但傳統(tǒng)有限差分法忽略了波場的變化，在實(shí)際應(yīng)用中會出現(xiàn)幾個方面的問題： 第一，地下介質(zhì)速度隨深度變化，地震子波波長也隨深度變化，如果網(wǎng)格間距不變，導(dǎo)致波場空間采樣不均勻，在低速層會出現(xiàn)欠采樣，在高速層會出現(xiàn)過采樣，在橫向變速劇烈的區(qū)域則更嚴(yán)重；第二，地震勘探所面對的區(qū)域通常涉及到起伏地表[14]、低速異常體、孔隙裂縫等特殊地質(zhì)構(gòu)造，采用固定網(wǎng)格間距的常規(guī)深度域采樣可能會導(dǎo)致算法的不穩(wěn)定和精度的降低；第三，當(dāng)存在地層倒轉(zhuǎn)時，在傳統(tǒng)笛卡爾坐標(biāo)系下，波場延拓方向與波場能量傳播的方向不一致[15]。

物理上，地震波波長隨介質(zhì)速度的擾動而變化，因此，網(wǎng)格的設(shè)計應(yīng)該基于真實(shí)地下介質(zhì)的速度擾動[16]。與傳統(tǒng)規(guī)則的笛卡爾網(wǎng)格相比，根據(jù)地下介質(zhì)速度設(shè)計的不規(guī)則網(wǎng)格能更精確地描繪復(fù)雜地質(zhì)構(gòu)造[17-19]。為此，黃超等[20]提出了變網(wǎng)格間距的方法，根據(jù)地下介質(zhì)的速度擾動設(shè)計不同的空間網(wǎng)格步長，提高了計算效率；張慧等[21]和曲英銘[22]提出了時空雙變的方法，對速度場進(jìn)行局部加密，同時對時間層也進(jìn)行了局部變換，提高了對低降速帶等復(fù)雜地區(qū)的計算效率和適用性。這些方法雖然能有效改善過采樣現(xiàn)象，提高計算效率，但在過渡區(qū)插值時會引入相應(yīng)的誤差，且實(shí)施過程較為復(fù)雜。

Alkhalifah等[23]首先提出將深度坐標(biāo)轉(zhuǎn)換成垂向旅行時坐標(biāo)的思想，并應(yīng)用于VIT介質(zhì)聲波方程正演模擬，驗(yàn)證了方法的有效性，之后研究了偽深度域坐標(biāo)下的速度分析，證明在偽深度域進(jìn)行速度分析要比在深度域更穩(wěn)定[16]； Li等[24]和Sun等[25]完成了偽深度域下的聲波最小二乘逆時偏移，該方法能有效降低最小二乘逆時偏移算法對速度的依賴性； Plessix[26]研究了在偽深度域中的聲波全波形反演，認(rèn)為用垂向時間代替深度進(jìn)行速度分析，能克服反演的零空間問題，從而提供了一個更穩(wěn)健的算法?！皞紊疃扔颉辈呗允且环N能有效改善局部過采樣現(xiàn)象、節(jié)約內(nèi)存占用、提高計算效率的方法。

從偏移算法中所采用的速度函數(shù)看，地震偏移可分為時間偏移和深度偏移[27]。時間偏移算法具有較高的計算效率，對偏移速度的依賴性小，然而該方法基于速度橫向均勻性的假設(shè)，當(dāng)遇到嚴(yán)重的速度橫向變化或復(fù)雜的速度界面時，成像效果差； 深度偏移算法能較好地適應(yīng)強(qiáng)橫向變速的復(fù)雜介質(zhì)，然而該算法對速度的依賴性很強(qiáng)，且計算量大于時間偏移。偽深度域偏移算法對速度做了加權(quán)處理，因此能降低對速度的依賴性； 而且偽深度域算法本質(zhì)上也是基于雙程波動方程理論，因此其成像精度接近疊前逆時深度偏移。

基于前人的研究成果，本文首先利用坐標(biāo)變換的鏈?zhǔn)椒▌t，推導(dǎo)了波場的一階和二階導(dǎo)數(shù)從深度域到偽深度域的轉(zhuǎn)換形式，然后推導(dǎo)了偽深度域二階彈性波方程(PDD-EWE)；為了提高算法的計算效率，引入了變長度差分算子[28]的等效交錯網(wǎng)格[29-31]有限差分法，實(shí)現(xiàn)了偽深度域彈性波數(shù)值模擬和逆時偏移。結(jié)果表明，與深度域方法相比，偽深度域變算子長度差分算法既能保證計算精度，又能顯著提高計算效率、降低內(nèi)存占用。

1 方法原理

1.1 不同坐標(biāo)系間的映射關(guān)系

假設(shè)坐標(biāo)系A(chǔ)(x,z)到坐標(biāo)系B(ξ,η)的映射關(guān)系為

(1)

定義Fa、Fab為任意變量F的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)，有

(2)

式中a、b為坐標(biāo)變量。根據(jù)坐標(biāo)變換的鏈?zhǔn)椒▌t有

(3)

求解式(3)，可以得到兩種坐標(biāo)系下一階導(dǎo)數(shù)的轉(zhuǎn)換關(guān)系為

(4)

定義g=xξzη-xηzξ，則任意變量F在兩個坐標(biāo)系下的一階導(dǎo)數(shù)映射關(guān)系為

(5)

二階導(dǎo)數(shù)映射關(guān)系為

(6)

將Fxx展開，可化為

(Fξzη-Fηzξ)ηzξ]

(7)

定義Fxx=F1,xx+F2,xx，則

xξzξ ηzη-xηzξ ξzη)

(8)

xξzξ ηzη-xηzξ ξzη)

(9)

可以看出，式(8)的右端第二項(xiàng)和式(9)的右端第三項(xiàng)相等且符號相反，則式(7)可化簡為

(10)

同理，可得Fxz、Fzz的展開形式為

(11)

(12)

式(10)～式(12)即為從坐標(biāo)系A(chǔ)到坐標(biāo)系B的二階導(dǎo)數(shù)映射關(guān)系。

1.2 PDD-EWE公式推導(dǎo)

在笛卡爾坐標(biāo)系外推地震波場時，每個子波波長的網(wǎng)格數(shù)可表示為

(13)

式中：λ為波長；v為波的傳播速度；h為網(wǎng)格間距；T為子波周期。h和T通常是固定的，因此，當(dāng)速度增大時，每個子波波長的網(wǎng)格點(diǎn)數(shù)會增加，導(dǎo)致過采樣，加重內(nèi)存負(fù)擔(dān)。偽深度域的思想是將垂向深度轉(zhuǎn)換成垂直旅行時坐標(biāo)，然后進(jìn)行均勻重采樣，從而避免過采樣。值得注意的是，對垂直旅行時進(jìn)行均勻重采樣時，為避免空間假頻，垂向采樣間隔 應(yīng)滿足

(14)

偽深度域垂向網(wǎng)格點(diǎn)數(shù)為

(15)

式中：vmin為最小速度；fmax為地震波最大頻率；tmax為最大單程旅行時。定義從笛卡爾坐標(biāo)到偽深度域坐標(biāo)的映射關(guān)系為

(16)

式中vsm為平滑速度場。通常，偽深度域垂向網(wǎng)格點(diǎn)數(shù)要少于深度域，能節(jié)約計算內(nèi)存占用。同理從偽深度域坐標(biāo)到笛卡爾坐標(biāo)的映射關(guān)系為

(17)

則深度域到偽深度域的一階、二階導(dǎo)數(shù)映射關(guān)系為

(18)

將式(18)代入式(10)～式(12)，可以得到變量F在偽深度域中的二階空間導(dǎo)數(shù)

(19)

在均勻各向同性介質(zhì)中，以位移表示的二階常速度彈性波方程在深度域的表達(dá)式為

(20)

式中：u、w分別為x、z方向的位移分量；vP、vS分別為縱波速度和橫波速度。

將式(19)代入式(20)，可得PDD-EWE方程

(21)

值得注意的是，本文是利用縱波速度實(shí)現(xiàn)從深度域到偽深度域的轉(zhuǎn)換。測試結(jié)果表明，利用偽深度域的方法轉(zhuǎn)換得到的橫波速度場與利用縱橫波比得到的橫波速度場基本一致。以分量u為例說明式(21)中的空間偏導(dǎo)數(shù)等效交錯網(wǎng)格有限差分法格式。其一階導(dǎo)數(shù)的離散格式為

(22)

二階導(dǎo)數(shù)的離散格式為

(23)

對于PDD-EWE方程，雖然偏導(dǎo)數(shù)項(xiàng)的增加提高了計算的復(fù)雜度，但由于該方程沒有中間變量，與一階方程相比能節(jié)約計算內(nèi)存。

1.3 變長度差分算子基本原理

為了進(jìn)一步提升算法的計算效率，引入變長度差分算子[23]的思想。在一定的精度要求下，根據(jù)速度的不同設(shè)計不同長度的差分算子，可有效改善速度的非均勻性導(dǎo)致的計算消耗。

本文的波場延拓方法是基于等效交錯網(wǎng)格有限差分法，與傳統(tǒng)一階的交錯網(wǎng)格具有相同的精度，因此u分量的一階空間偏導(dǎo)數(shù)的交錯網(wǎng)格離散形式為

(24)

式中h為網(wǎng)格橫向間距。

對式(24)進(jìn)行傅里葉變換，并代入平面波解u=u0exp(ikx)，有

(25)

(26)

當(dāng)給定f、h后，誤差函數(shù)只與速度v和差分算子半長度M有關(guān)。當(dāng)給定誤差函數(shù)ε，則不同的速度會對應(yīng)不同的差分算子長度，速度越大，需要的差分算子長度越短。根據(jù)速度和限定誤差，調(diào)整差分算子長度，能有效提高計算效率。

2 數(shù)值模擬

2.1 層狀介質(zhì)

為了說明本文變長度差分算子偽深度域彈性波數(shù)值模擬的有效性和優(yōu)勢，設(shè)計一個簡單的層狀模型，其主要模擬參數(shù)為：網(wǎng)格點(diǎn)數(shù)為301×401，橫、縱向網(wǎng)格間距均為5m，時間采樣間隔為0.3ms，總接收時間為2.4s，震源采用主頻為30Hz的雷克子波，加載在模型表面中間，振幅增益為104倍，縱橫波速比為1.73，深度域縱波速度模型如圖1a所示。首先將模型從深度域轉(zhuǎn)換到偽深度域，在滿足采樣定理的前提下，取Δη=2.5ms，偽深度域速度模型如圖1b所示，網(wǎng)格點(diǎn)數(shù)為301×300。對比圖1b和圖1a可以看出，偽深度域的速度場實(shí)際上在淺層低速時被拉伸，深層高速時被壓縮，使得垂向波場采樣變得均勻，從而避免了過采樣。針對模擬參數(shù)，根據(jù)上文的變長度差分算子計算方法，當(dāng)f=30Hz、h=5m、將誤差控制在10-9時，計算的不同速度對應(yīng)的不同差分算子長度如圖1c所示。根據(jù)圖1c，可以得到層狀模型對應(yīng)的差分算子長度，傳統(tǒng)深度域差分算子長度2M=12。另外，在數(shù)值模擬的過程中，為防止邊界反射帶來的干擾，采用海綿吸收邊界條件[31]，其衰減因子為

(27)

式中：np為吸收邊界層厚度，吸收層加在模型外層；i為計算點(diǎn)到外邊界的距離；a為衰減系數(shù)。在本文的數(shù)值模擬和逆時偏移中，吸收層厚度均設(shè)為50個網(wǎng)格。

圖2是利用深度域方法、偽深度域方法及引入變長度差分算子后的偽深度域方法得到的層狀模型單炮地震記錄，從圖中可以看出，三種方法模擬的單炮地震記錄基本一致。進(jìn)一步對三種方法模擬的單炮地震記錄做殘差分析。圖3a為偽深度域與深度域模擬結(jié)果的殘差，水平和垂直分量的相對誤差分別約為0.25%和0.12%。產(chǎn)生這種誤差的原因是：將圖1a轉(zhuǎn)換成圖1b時，需要對速度重采樣，本文選擇了三次樣條插值方法。圖3b為偽深度域變長度差分算子與深度域模擬結(jié)果的殘差，與圖3a基本相同；圖3c為偽深度域變差分算子長度與固定差分算子長度模擬結(jié)果的殘差，水平和垂直分量的相對誤差分別約為1.5×10-6和7.0×10-7，幾乎可以忽略不計。原始深度域速度模型網(wǎng)格數(shù)為301×401，轉(zhuǎn)換到偽深度域中的速度模型網(wǎng)格數(shù)為301×300，相當(dāng)于節(jié)約了25.2%的內(nèi)存。應(yīng)用CPU型號為Intel(R) Xeon(R) CPU E5-2630 v4、主頻為2.20GHz的計算機(jī)，三種方法運(yùn)行時間分別為3705、3466和2946s，偽深度域變長度差分算子模擬方法與深度域方法相比，效率提高了20.5%，與偽深度域固定長度差分算子模擬方法相比，提高了13.8%。

圖1 層狀模型

2.2 Marmousi 模型

選取Marmousi模型驗(yàn)證本文方法對復(fù)雜介質(zhì)的適用性。將原始Marmousi模型橫向按10∶1、縱向按3∶1抽稀，得到的深度域縱波速度模型，網(wǎng)格數(shù)為1360×700，如圖4a所示。將深度域速度場轉(zhuǎn)換到偽深度域，縱向采樣間隔Δη=2.5ms，網(wǎng)格數(shù)為1360×507，如圖4b所示。可以看出，偽深度域模型淺層被拉伸，深層被壓縮，對應(yīng)的有限差分算子長度如圖4c所示，而傳統(tǒng)深度域方法差分算子長度2M=12。模擬參數(shù)為：時間采樣間隔為0.3ms，空間采樣間隔為5m，總接收時間為2.1s，震源采用中心頻率為30Hz的雷克子波，振幅增益104倍，縱橫波速度比設(shè)為1.73。

圖5為傳統(tǒng)深度域方法和本文偽深度域變長度差分算子方法彈性波數(shù)值模擬單炮記錄對比，圖6為兩種方法模擬記錄在x=2.0km處的單道記錄，可以看出本文方法和深度域方法的單道地震記錄基本保持一致，沒有走時誤差，只存在振幅差異，水平和垂直分量最大相對誤差分別為0.3%、0.23%。

圖2 層狀模型不同方法模擬的水平分量(上)和垂直分量(下)單炮記錄

圖3 層狀模型不同方法模擬的水平分量(上)和垂直分量(下)單炮記錄殘差

圖4 Marmousi模型

圖5 Marmousi模型兩種方法模擬的水平分量(上)和垂直分量(下)單炮記錄及及殘差

圖6 Marmousi模型兩種方法模擬結(jié)果在x=2.0km處單道(左)及局部放大(右)

由此可見，偽深度域變長度差分算子彈性波模擬方法適用于復(fù)雜介質(zhì)模型，且能節(jié)省約28%內(nèi)存占用。

3 偏移成像

將本文偽深度域變長度差分算子彈性波模擬方法應(yīng)用于Marmousi模型的彈性波逆時偏移，驗(yàn)證本文方法的實(shí)用價值。模擬參數(shù)為：網(wǎng)格數(shù)為1360×700，縱橫向網(wǎng)格間距為5m，時間采樣間隔為0.3ms，總接收時間為3.0s，震源采用主頻為30Hz的雷克子波，地面布置61炮，炮點(diǎn)距為100m，采用全接收方式，道間距為5m，速度模型和差分算子長度如圖4所示，吸收邊界厚度為50個網(wǎng)格間距。

本文逆時偏移采用互相關(guān)的成像條件。圖7為常規(guī)深度域方法和本文方法Marmousi模型彈性波逆時偏移的PP波和PS波成像結(jié)果，可以看出，與傳統(tǒng)深度域方法一樣，本文方法偏移結(jié)果能較好地刻畫Marmousi模型的構(gòu)造層位，可見本文方法能適用于復(fù)雜介質(zhì)的彈性波逆時偏移成像。將圖7c和圖7d轉(zhuǎn)換到深度域，抽取x=2.0和6.0km處單道偏移結(jié)果與圖7a和圖7b對比，如圖8所示。可以看出，兩種方法的單道曲線之間基本沒有深度誤差，只存在較小的振幅誤差?；谙嗤挠嬎銠C(jī)資源(CPU型號為Intel(R) Xeon(R) CPU E5-2630 v4、主頻為2.20GHz，18個節(jié)點(diǎn))，傳統(tǒng)深度域方法Marmousi模型彈性波逆時偏移用時24.5h，本文的偽深度域變長度差分算子方法用時18.5h，效率提升了24.4%。

圖7 Marmousi模型兩種方法彈性波逆時偏移結(jié)果

圖8 Marmousi模型兩種方法彈性波逆時偏移結(jié)果不同位置處單道對比

4 結(jié)論

本文發(fā)展了一種基于偽深度域的變長度差分算子有限差分彈性波數(shù)值模擬方法，與傳統(tǒng)深度域方法相比，有效改善了由于速度差異過大導(dǎo)致的波場采樣不均勻現(xiàn)象。理論分析和模型測試結(jié)果表明：

(1)偽深度域彈性波數(shù)值模擬方法主要優(yōu)勢在于能節(jié)約內(nèi)存占用，變長度差分算子方法能有效提升計算效率。在保證精度的前提下，與常規(guī)深度域方法相比，偽深度域變長度差分算子方法能節(jié)省約25%～30%內(nèi)存，效率提升約20%～25%；

(2)對速度場在偽深度域重采樣處理會導(dǎo)致模擬結(jié)果的振幅誤差，與傳統(tǒng)深度域方法相比，本文方法的簡單模型和Marmousi模型模擬結(jié)果的振幅相對誤差不超過1%；

(3)將本文方法應(yīng)用于Marmousi模型彈性波逆時偏移，與常規(guī)深度域方法偏移結(jié)果的單道波形一致性較好，只存在較小的振幅誤差，驗(yàn)證了本文方法對復(fù)雜介質(zhì)模型的適用性。