分數階拉普拉斯算子黏滯聲波方程的最小二乘逆時偏移

2020-06-03 08:13:52陳漢明田玉昆

石油地球物理勘探 2020年3期

陳漢明周輝* 田玉昆

(①油氣資源與探測國家重點實驗室，北京 102249； ②CNPC物探重點實驗室，北京 102249；③中國石油大學(北京)地球物理學院，北京 102249； ④中國地質調查局油氣資源調查中心，北京 100083)

0 引言

地震偏移利用波場傳播算子將地表觀測的地震數據轉換為地下空間分布的像，是整個地震數據處理流程的重要一環。在眾多偏移方法中，逆時偏移方法因無成像傾角限制、保幅性能好及易于向復雜介質擴展等優點，已成為重要的地震成像方法。然而，受制于幾何擴散、有限的觀測孔徑和子波帶寬等因素影響，逆時偏移的成像結果仍存在振幅分布不均衡、噪聲干擾及分辨率低等問題。最小二乘逆時偏移方法是逆時偏移技術的進一步發展，它將地震偏移視為反演問題，通過匹配反射波數據反演地下的反射系數模型，利用多次迭代壓制成像噪聲，同時能在一定程度上消除子波頻帶不足的負面影響，提高成像分辨率。

近年來，隨著計算機性能的大幅提升，特別是GPU并行計算技術的推廣，具有諸多優點的最小二乘逆時偏移方法已成為研究熱點之一。最小二乘逆時偏移方法最早由Dai等[1]基于聲波方程提出。針對該方法應用中存在的問題，現階段許多學者從不同的角度對其做了進一步的完善。如Yao等[2-3]和Dai等[4]對最小二乘逆時偏移的成像精度進行了研究。Yao等[5]將最小二乘逆時偏移表述為矩陣—向量形式，討論了非線性的最小二乘逆時偏移方法。根據反射系數具有稀疏特性，Wu等[6]研究了L1范數約束的最小二乘逆時偏移方法。方修政等[7]使用逆散射成像條件替代傳統的成像條件，壓制了最小二乘逆時偏移中的成像噪聲。針對復雜介質的影響，諸多學者研究了基于彈性波動方程的最小二乘逆時偏移方法[8-10]。Fang等[11]利用卷積型目標函數提出了不依賴于震源的彈性波最小二乘逆時偏移方法。李振春等[12]和郭旭等[13]進一步研究了基于VTI介質的最小二乘逆時偏移方法。最小二乘逆時偏移需對波動方程進行Born近似，為提高近似精度，陳生昌等[14]提出反射波動方程，并討論了基于反射波動方程的最小二乘逆時偏移方法。在提高計算效率方面，正則化約束下的多震源最小二乘逆時偏移方法在近年得到了較廣泛的關注[15-18]。

此外，吸收衰減補償的逆時偏移方法近年來也是研究熱點之一。根據互相關成像條件[19]，為了補償地震數據在傳播路徑上的振幅衰減和相位畸變，震源波場和伴隨波場都是隨時間呈指數增長的，該過程可通過求解振幅衰減和相位畸變解耦的黏滯波動方程實現[20-22]。然而，數值求解振幅增長的黏滯波動方程存在數值不穩定問題，通常采用低通濾波器濾掉高頻成分[23-25]；同時又因截止頻率不易選取，會導致高頻丟失，進而降低了分辨率。Wang等[26]和Zhao等[27]分別提出了不同的穩定化策略，但都無法從根本上避免數值不穩定問題。

最小二乘逆時偏移方法采用迭代思路反演反射率。若考慮介質的衰減性質，其正演和伴隨算子都是能量衰減型的，因此不存在數值不穩定問題。Dutta等[28]最早提出利用最小二乘逆時偏移方法補償介質黏滯性對成像結果的影響，避免了衰減補償逆時偏移的不穩定問題；李振春等[29]也討論了基于黏滯聲學介質的最小二乘逆時偏移方法； Guo等[30]進一步研究了基于廣義標準線性體黏彈性波動方程的最小二乘逆時偏移方法；曲英銘等[31]在最小二乘逆時偏移中同時考慮了介質的黏滯性和各向異性。

目前大多數衰減補償最小二乘逆時偏移方法均采用線性體黏滯波動方程作為正演方程[32-33]，該方程難以準確描述常Q(地震品質因子)模型[34]，由其模擬的地震波衰減和頻散與實際測量結果有一定差距[35]。基于此，本文利用一種新穎的分數階拉普拉斯算子黏滯聲波方程作為正演方程，推導其Born正演模擬算子，由此提出了一種新的衰減補償型最小二乘逆時偏移方法，用于補償介質的黏滯性對逆時偏移成像結果的影響。

1 黏滯聲波最小二乘逆時偏移

1.1 分數階拉普拉斯算子黏滯聲波方程

Born正演模擬算子是通過線性化全波方程得到的，本文采用Chen等[36]提出的常分數階拉普拉斯算子黏滯聲波方程描述地震波在黏滯聲學介質中衰減和頻散。該方程可表述為

Au0=s(t)l(x-xs)

(1)

式中：u0為波場變量；s(t)為震源子波； l為單位脈沖函數；x表示二維空間坐標；xs為震源位置；且有

(2)

(3)

式中：Q為地震品質因子；c0為定義在參考角頻率ω0上的地震波速度；ωd為震源的主頻。

Chen等[36]證實，在Q≥15的情況下，黏滯聲波方程式(1)的精度能很好地匹配Zhu等[20]所提出的變分數階拉普拉斯算子常Q黏滯聲波方程的精度。與后者相比，式(1)中拉普拉斯算子的階數不隨空間變化，因此數值求解該波動方程更容易。

目前求解分數階拉普拉斯算子黏滯聲波方程的數值方法包括偽譜法[37]、矩陣轉換法[38]、矩陣低秩分解法[39]等。為簡單起見，采用偽譜法對式(1)進行離散，分別使用二階中心差分和一階向后差分算子近似式(1)中的二階和一階時間偏導數，同時利用快速傅里葉變換(FFT)計算分數階拉普拉斯算子，例如

(4)

式中： F和F-1分別表示正、反FFT；k為波數。為避免邊界反射，本文使用分裂式完美匹配層(PML)[40]作為吸收邊界條件。

1.2 Born正演模擬

為了推導全波方程(式(1))對應的Born正演模擬算子，對背景速度模型c0施加微小擾動c=c0+δc，相應地，波場也產生微小擾動u0→u=u0+δu。將擾動后的速度模型c和波場u代入式(1)，并利用以下泰勒近似

(5)

可得到

=s(t)l(x-xs)

(6)

將式(6)減去式(1)，并經簡單整理得到

Aδu=mDu

(7)

其中

(8)

(9)

顯然，Born正演模擬包含兩次黏滯聲波方程正演模擬：一次是利用背景速度模型c0計算背景波場u0；另一次是將背景波場作為震源，計算反射波場δu，且反射波場與反射率m之間是線性關系。因此，Born近似下的最小二乘逆時偏移是線性反演問題。

1.3 伴隨算子和梯度

利用拉格朗日乘子法定義擴展的目標函數[41]

?μ(Aδu-mDu0)dtdx

(10)

式中：m={mi,j|0≤i

A*μ=R(δuobs-δucal)

(11)

其中

(12)

注意A*與A(式(2))的區別僅在于最后一項的符號相反。正演模擬過程中，當此項的符號為正時，對應振幅衰減；當此項符號為負時，表示振幅增長[23]。然而，由于伴隨方程式(11)中的殘差數據是逆時輸入的，式(12)中對時間一階偏導數的差分離散由向后差分變成了向前差分，因此伴隨波場μ沿時間外推的差分格式與正演模擬所用的格式一致，故伴隨波場也是振幅衰減的，不存在數值不穩定問題。

進一步計算目標函數對反射率模型的偏導數，即目標函數的梯度可由背景波場和伴隨波場的零延遲互相關求取

(13)

式(11)和式(13)的具體推導過程見附錄A。

本文采用L-BFGS局部尋優算法更新模型

mn+1=mn+αzn

(14)

2 數值例子

2.1 合成數據反演

首先使用Marmousi模型檢驗本文所提最小二乘逆時偏移方法效果。模型尺寸為340(x)×156(z)，網格單元為15m×10m，最大和最小速度分別為4700、1028m/s。基于真實速度模型(圖1a)得到光滑后速度模型(圖1b)；再利用經驗公式Q=12(c0/1000)2.2由光滑速度模型生成Q模型(圖1c)，最小Q值為28；進而得到真實的反射率模型(圖1d)。假設速度模型(圖1a和圖1b)定義在參考頻率200Hz上。基于真實的速度模型求解黏滯聲波方程式(1)，生成60炮觀測記錄，檢波點位于地表，炮點位于深度20m處，炮間距75m，第1炮距離模型左邊界水平距離為300m，采用中間放炮兩邊接收的排列，最大炮檢距為1800m，道間距15m。使用主頻為20Hz雷克子波作為震源，記錄時長為2.5s，時間采樣間隔為1ms。

圖2為模擬的第30炮共炮點記錄，炮點位于x=2475m。對比黏滯聲波記錄(圖2b)與聲波記錄(圖2a)，可見存在明顯振幅衰減。基于真實反射率模型計算的黏滯Born模擬記錄(圖2c)中主要包含一次反射波形，說明Born模擬不同于全波模擬。觀察該炮第100道(x=2160m)記錄的波形(圖2d)，可見黏滯聲波數據相較于聲波數據存在明顯的振幅衰減和相位滯后，同時黏滯Born模擬數據與黏滯全波數據也有較明顯差異，這緣于線性化波動方程誤差。

不同反演方法獲得的反射率模型如圖3所示，所使用的初始模型為零，其中黏滯聲波數據聲波最小二乘逆時偏移的第1次迭代結果(圖3a)，相當于傳統聲波逆時偏移的結果，可見雖然淺層結構較清晰，但深層能量弱，振幅沿深度分布不均衡，且與真實反射率模型(圖1d)有明顯差異。黏滯聲波數據黏滯聲波最小二乘逆時偏移的第1次迭代結果(圖3c)，整體上與圖3a相似，但位于500m

當迭代次數增至60時，聲波最小二乘逆時偏移反演黏滯聲波數據的結果分辨率較低，存在明顯的低頻干擾(圖3b)。相對而言，黏滯聲波最小二乘逆時偏移方法很好地補償了介質的黏滯性，其反演結果(圖3d)在z<800m范圍內與真實反射率模型吻合很好，在該深度以下的背斜和不整合面的成像也較清晰，但也存在一些噪聲，這是由于正演與反演算子欠匹配造成的。實際上，觀測數據是由全波方程式(1)生成的，反射波與反射率之間是非線性關系，但反演時利用Born算子進行線性反演，因此反演結果與真實結果存在差異。

作為對比，使用黏滯Born數據作為觀測數據進行最小二乘逆時偏移，同樣迭代60次。從偏移結果可知：當正演和反演算子同為線性時，反演結果有明顯改善；聲波最小二乘逆時偏移反演結果(圖3e)對構造成像較清晰，但振幅仍未恢復，尚存在低頻噪聲；而黏滯聲波最小二乘逆時偏移反演結果(圖3f)與真實反射率模型(圖1d)十分接近。

圖3 針對黏滯全波和黏滯Born數據采用兩種偏移方法反演得到的反射率模型a)、(b)對應黏滯全波數據聲波最小二乘逆時偏移第1和第60次迭代； (c)、(d)對應黏滯全波數據黏滯聲波最小二乘逆時偏移第1和第60次迭代； (e)黏滯Born數據聲波最小二乘逆時偏移第60次迭代； (f)黏滯Born數據黏滯聲波最小二乘逆時偏移第60次迭代

圖4進一步比較了三種最小二乘逆時偏移方法所對應的目標函數的收斂速度。從該圖易知，三者中黏滯Born數據的黏滯聲波最小二乘逆時偏移的收斂速度最快，其次是黏滯全波數據的黏滯聲波最小二乘逆時偏移，而黏滯全波數據的傳統聲波最小二乘逆時偏移的收斂速度相對最慢。

圖4 三種最小二乘逆時偏移目標函數的收斂速度對比

2.2 實際數據反演

進一步使用一條二維實際測線證實本文所提黏滯聲波最小二乘逆時偏移方法的可行性。該測線共80炮，炮間距為200m，炮點深度為20m，記錄時長為3s，時間采樣間隔為1ms，道間距為20m，從頻譜分析知該數據主頻約為20Hz。速度(圖5a)和Q(圖5b)模型的網格數為680(x)×320(z)，網格間距為20m，假設速度模型定義在參考頻率20Hz上。從第40炮共炮點記錄(圖6)可見若干反射波同相軸、未消除徹底的面波及其他干擾噪聲。最小二乘逆時偏移中使用主頻為20Hz雷克子波作為震源，并假設初始的反射率模型為零。

圖7顯示了由最小二乘逆時偏移得到的反射率剖面，其中包括聲波最小二乘逆時偏移第1(圖7a)、第20(圖7c)、第40(圖7e)次和黏滯聲波最小二乘逆時偏移第1(圖7b)、第20(圖7d)和第40(圖7f)次的迭代結果，且第1次迭代相當于逆時偏移。由第1次迭代的結果(圖7a和圖7b)可知，聲波和黏滯聲波逆時偏移均未能補償介質的黏滯性，其偏移剖面深層的能量明顯弱于淺層，特別是黏滯聲波逆時偏移中因為震源和伴隨波場的能量都是衰減的，導致其偏移剖面(圖7b)的分辨率低于聲波偏移剖面(圖7a)。當迭代次數增至20時，黏滯聲波最小二乘逆時偏移逐漸補償了介質的黏滯性，成像結果(圖7d)優于傳統聲波最小二乘逆時偏移(圖7c)。對比圖7c與圖7d中解釋的3條斷層及若干同相軸可知，補償吸收衰減后的斷層面變得較突出，位于兩條斷層之間的小斷塊也能清晰可辨。

圖5 實際地震資料的速度模型(a)和Q模型(b)

圖6 實際共炮點記錄

圖7 實際數據的最小二乘逆時偏移剖面(a)、(c)、(e)對應聲波最小二乘逆時偏移迭代1、20、40次； (b)、(d)、(f)對應黏滯聲波最小二乘逆時偏移迭代1、20、40次

當迭代次數進一步增至40時，黏滯聲波最小二乘逆時偏移結果進一步凸顯了高頻細節，其成像剖面(圖7f)呈現更多具有一定連續性的細軸。可對比圖7d與圖7f中位于5000m

3 討論

為了考察Q模型不準確對最小二乘逆時偏移結果的影響，這里重復圖3f中合成數據測試，即利用黏滯聲波最小二乘逆時偏移方法反演黏滯Born數據。所使用的Q模型與圖1c中真實Q模型分別存在-50%、-30%、-10%、50%、30%和10%的相對誤差，最大迭代次數為60次。其反演結果如圖8所示，其中標注的RMS值是與參考解(圖3f)之間的均方根誤差。從該反演結果可知：當Q偏小時，黏滯聲波最小二乘逆時偏移剖面的振幅存在正向異常，表明補償過度；當Q偏大時，偏移剖面的振幅偏弱，表明補償不夠。當Q模型的相對誤差為-10%和10%時，其反演結果(圖8e和圖8f)與參考剖面(圖3f)較為接近，其誤差剖面(圖8g和圖8h)的幅度較弱，說明黏滯聲波最小二乘逆時偏移方法對于10%的Q模型誤差具有一定的穩健性。

圖8 不同相對誤差Q模型的黏滯Born數據使用黏滯聲波最小二乘逆時偏移迭代60次的結果(a)-50%； (b)50%； (c)-30%； (d)30%； (e)-10%； (f)10%； (g)圖e與參考解(圖3f)的差異； (h)圖f與參考解的差異RMS表示偏移結果與參考解之間的均方根誤差

4 結論

本文基于一種新穎的分數階拉普拉斯算子常Q黏滯聲波方程，提出了一種具有補償介質吸收衰減功能的最小二乘逆時偏移方法。該方法從常分數階拉普拉斯算子黏滯聲波方程出發，定義了相應的Born正演模擬算子，在拉格朗日乘子法框架下推導出伴隨算子和梯度表達式，并利用L-BFGS局部尋優算法迭代更新反射率模型。因正演和伴隨算子都是能量衰減型的，故本文的黏滯聲波最小二乘逆時偏移方法能避免常規Q補償逆時偏移方法數值不穩定的弊端。合成和實際數據的反演算例都證實了本文所提方法能較好地補償介質的黏滯性，提高成像分辨率。