分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程的最小二乘逆時(shí)偏移

陳漢明 周 輝* 田玉昆

(①油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102249； ②CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京 102249；③中國(guó)石油大學(xué)(北京)地球物理學(xué)院，北京 102249； ④中國(guó)地質(zhì)調(diào)查局油氣資源調(diào)查中心，北京 100083)

0 引言

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

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

此外，吸收衰減補(bǔ)償?shù)哪鏁r(shí)偏移方法近年來(lái)也是研究熱點(diǎn)之一。根據(jù)互相關(guān)成像條件[19]，為了補(bǔ)償?shù)卣饠?shù)據(jù)在傳播路徑上的振幅衰減和相位畸變，震源波場(chǎng)和伴隨波場(chǎng)都是隨時(shí)間呈指數(shù)增長(zhǎng)的，該過(guò)程可通過(guò)求解振幅衰減和相位畸變解耦的黏滯波動(dòng)方程實(shí)現(xiàn)[20-22]。然而，數(shù)值求解振幅增長(zhǎng)的黏滯波動(dòng)方程存在數(shù)值不穩(wěn)定問(wèn)題，通常采用低通濾波器濾掉高頻成分[23-25]； 同時(shí)又因截止頻率不易選取，會(huì)導(dǎo)致高頻丟失，進(jìn)而降低了分辨率。Wang等[26]和Zhao等[27]分別提出了不同的穩(wěn)定化策略，但都無(wú)法從根本上避免數(shù)值不穩(wěn)定問(wèn)題。

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

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

1 黏滯聲波最小二乘逆時(shí)偏移

1.1 分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程

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

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

(1)

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

(2)

(3)

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

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

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

(4)

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

1.2 Born正演模擬

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

(5)

可得到

=s(t)l(x-xs)

(6)

將式(6)減去式(1)，并經(jīng)簡(jiǎn)單整理得到

Aδu=mDu

(7)

其中

(8)

(9)

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

1.3 伴隨算子和梯度

利用拉格朗日乘子法定義擴(kuò)展的目標(biāo)函數(shù)[41]

?μ(Aδu-mDu0)dtdx

(10)

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

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

(11)

其中

(12)

注意A*與A(式(2))的區(qū)別僅在于最后一項(xiàng)的符號(hào)相反。正演模擬過(guò)程中，當(dāng)此項(xiàng)的符號(hào)為正時(shí)，對(duì)應(yīng)振幅衰減； 當(dāng)此項(xiàng)符號(hào)為負(fù)時(shí)，表示振幅增長(zhǎng)[23]。然而，由于伴隨方程式(11)中的殘差數(shù)據(jù)是逆時(shí)輸入的，式(12)中對(duì)時(shí)間一階偏導(dǎo)數(shù)的差分離散由向后差分變成了向前差分，因此伴隨波場(chǎng)μ沿時(shí)間外推的差分格式與正演模擬所用的格式一致，故伴隨波場(chǎng)也是振幅衰減的，不存在數(shù)值不穩(wěn)定問(wèn)題。

進(jìn)一步計(jì)算目標(biāo)函數(shù)對(duì)反射率模型的偏導(dǎo)數(shù)，即目標(biāo)函數(shù)的梯度可由背景波場(chǎng)和伴隨波場(chǎng)的零延遲互相關(guān)求取

(13)

式(11)和式(13)的具體推導(dǎo)過(guò)程見(jiàn)附錄A。

本文采用L-BFGS局部尋優(yōu)算法更新模型

mn+1=mn+αzn

(14)

2 數(shù)值例子

2.1 合成數(shù)據(jù)反演

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

圖2為模擬的第30炮共炮點(diǎn)記錄，炮點(diǎn)位于x=2475m。對(duì)比黏滯聲波記錄(圖2b)與聲波記錄(圖2a)，可見(jiàn)存在明顯振幅衰減?；谡鎸?shí)反射率模型計(jì)算的黏滯Born模擬記錄(圖2c)中主要包含一次反射波形，說(shuō)明Born模擬不同于全波模擬。觀察該炮第100道(x=2160m)記錄的波形(圖2d)，可見(jiàn)黏滯聲波數(shù)據(jù)相較于聲波數(shù)據(jù)存在明顯的振幅衰減和相位滯后，同時(shí)黏滯Born模擬數(shù)據(jù)與黏滯全波數(shù)據(jù)也有較明顯差異，這緣于線性化波動(dòng)方程誤差。

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

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

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

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

圖4進(jìn)一步比較了三種最小二乘逆時(shí)偏移方法所對(duì)應(yīng)的目標(biāo)函數(shù)的收斂速度。從該圖易知，三者中黏滯Born數(shù)據(jù)的黏滯聲波最小二乘逆時(shí)偏移的收斂速度最快，其次是黏滯全波數(shù)據(jù)的黏滯聲波最小二乘逆時(shí)偏移，而黏滯全波數(shù)據(jù)的傳統(tǒng)聲波最小二乘逆時(shí)偏移的收斂速度相對(duì)最慢。

圖4 三種最小二乘逆時(shí)偏移目標(biāo)函數(shù)的收斂速度對(duì)比

2.2 實(shí)際數(shù)據(jù)反演

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

圖7顯示了由最小二乘逆時(shí)偏移得到的反射率剖面，其中包括聲波最小二乘逆時(shí)偏移第1(圖7a)、第20(圖7c)、第40(圖7e)次和黏滯聲波最小二乘逆時(shí)偏移第1(圖7b)、第20(圖7d)和第40(圖7f)次的迭代結(jié)果，且第1次迭代相當(dāng)于逆時(shí)偏移。由第1次迭代的結(jié)果(圖7a和圖7b)可知，聲波和黏滯聲波逆時(shí)偏移均未能補(bǔ)償介質(zhì)的黏滯性，其偏移剖面深層的能量明顯弱于淺層，特別是黏滯聲波逆時(shí)偏移中因?yàn)檎鹪春桶殡S波場(chǎng)的能量都是衰減的，導(dǎo)致其偏移剖面(圖7b)的分辨率低于聲波偏移剖面(圖7a)。當(dāng)?shù)螖?shù)增至20時(shí)，黏滯聲波最小二乘逆時(shí)偏移逐漸補(bǔ)償了介質(zhì)的黏滯性，成像結(jié)果(圖7d)優(yōu)于傳統(tǒng)聲波最小二乘逆時(shí)偏移(圖7c)。對(duì)比圖7c與圖7d中解釋的3條斷層及若干同相軸可知，補(bǔ)償吸收衰減后的斷層面變得較突出，位于兩條斷層之間的小斷塊也能清晰可辨。

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

圖6 實(shí)際共炮點(diǎn)記錄

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

當(dāng)?shù)螖?shù)進(jìn)一步增至40時(shí)，黏滯聲波最小二乘逆時(shí)偏移結(jié)果進(jìn)一步凸顯了高頻細(xì)節(jié)，其成像剖面(圖7f)呈現(xiàn)更多具有一定連續(xù)性的細(xì)軸。可對(duì)比圖7d與圖7f中位于5000m

3 討論

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

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

4 結(jié)論

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

附錄A 黏滯聲波最小二乘逆時(shí)偏移伴隨算子和梯度表達(dá)式

二維情況下目標(biāo)函數(shù)(式(10))的具體形式為

(A-1)

(A-2)

對(duì)式(A-1)右端第二項(xiàng)先利用分部積分公式化簡(jiǎn)，再求其導(dǎo)數(shù)，例如

(A-3)

假設(shè)零值邊界條件δu(t=0,T)=0，μ(t=0,T)=0，則式(A-3)可進(jìn)一步化簡(jiǎn)為

(A-4)

同理，假設(shè)空間無(wú)窮遠(yuǎn)處滿足零值邊界條件，則有

(A-5)

由式(A-4)和(A-5)的推導(dǎo)可知，當(dāng)微分算子l 對(duì)稱且滿足零值邊界條件時(shí)，下式成立

(A-6)

從式(A-6)可得

(A-7)

(A-8)

(A-9)

式中下標(biāo)i和j為二維模型的網(wǎng)格索引。將其寫(xiě)成向量形式即為梯度表達(dá)式

(A-10)