偽深度域彈性波有限差分數值模擬及逆時偏移

鄒 強 黃建平* 李慶洋 雍 鵬 辛天亮

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

0 引言

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

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

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

Alkhalifah等[23]首先提出將深度坐標轉換成垂向旅行時坐標的思想，并應用于VIT介質聲波方程正演模擬，驗證了方法的有效性，之后研究了偽深度域坐標下的速度分析，證明在偽深度域進行速度分析要比在深度域更穩定[16]； Li等[24]和Sun等[25]完成了偽深度域下的聲波最小二乘逆時偏移，該方法能有效降低最小二乘逆時偏移算法對速度的依賴性； Plessix[26]研究了在偽深度域中的聲波全波形反演，認為用垂向時間代替深度進行速度分析，能克服反演的零空間問題，從而提供了一個更穩健的算法。“偽深度域”策略是一種能有效改善局部過采樣現象、節約內存占用、提高計算效率的方法。

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

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

1 方法原理

1.1 不同坐標系間的映射關系

假設坐標系A(x,z)到坐標系B(ξ,η)的映射關系為

(1)

定義Fa、Fab為任意變量F的一階導數和二階導數，有

(2)

式中a、b為坐標變量。根據坐標變換的鏈式法則有

(3)

求解式(3)，可以得到兩種坐標系下一階導數的轉換關系為

(4)

定義g=xξzη-xηzξ，則任意變量F在兩個坐標系下的一階導數映射關系為

(5)

二階導數映射關系為

(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)的右端第二項和式(9)的右端第三項相等且符號相反，則式(7)可化簡為

(10)

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

(11)

(12)

式(10)～式(12)即為從坐標系A到坐標系B的二階導數映射關系。

1.2 PDD-EWE公式推導

在笛卡爾坐標系外推地震波場時，每個子波波長的網格數可表示為

(13)

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

(14)

偽深度域垂向網格點數為

(15)

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

(16)

式中vsm為平滑速度場。通常，偽深度域垂向網格點數要少于深度域，能節約計算內存占用。同理從偽深度域坐標到笛卡爾坐標的映射關系為

(17)

則深度域到偽深度域的一階、二階導數映射關系為

(18)

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

(19)

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

(20)

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

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

(21)

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

(22)

二階導數的離散格式為

(23)

對于PDD-EWE方程，雖然偏導數項的增加提高了計算的復雜度，但由于該方程沒有中間變量，與一階方程相比能節約計算內存。

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

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

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

(24)

式中h為網格橫向間距。

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

(25)

(26)

當給定f、h后，誤差函數只與速度v和差分算子半長度M有關。當給定誤差函數ε，則不同的速度會對應不同的差分算子長度，速度越大，需要的差分算子長度越短。根據速度和限定誤差，調整差分算子長度，能有效提高計算效率。

2 數值模擬

2.1 層狀介質

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

(27)

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

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

圖1 層狀模型

2.2 Marmousi 模型

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

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

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

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

圖4 Marmousi模型

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

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

由此可見，偽深度域變長度差分算子彈性波模擬方法適用于復雜介質模型，且能節省約28%內存占用。

3 偏移成像

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

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

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

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

4 結論

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

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

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

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