栗學(xué)磊, 魏彥杰, 歐陽威

1 中國(guó)科學(xué)院深圳先進(jìn)技術(shù)研究院， 深圳 518055 2 中國(guó)科學(xué)院精密測(cè)量科學(xué)與技術(shù)創(chuàng)新研究院， 武漢 430077

0 引言

基于Born近似的線性反演方法是一類當(dāng)前流行的定量化地震成像和反演方法，可對(duì)物性反演和巖性解釋提供有效的支持，比如最小二乘偏移(LSM)和廣義Radon變換(GRT)反演.LSM在近20年一直是非?；钴S的研究領(lǐng)域(Stanton and Sacchi, 2017)，它在提高分辨率和壓制腳印等方面比常規(guī)偏移成像方法具有明顯優(yōu)勢(shì)(劉國(guó)章等，2020；劉玉敏等，2021).LSM通過迭代偏移和反偏移算子，擬合重建數(shù)據(jù)與原始數(shù)據(jù)，來求解基于Born近似的擾動(dòng)參數(shù).Nemeth等(1999)實(shí)現(xiàn)了最小二乘 Kirchhoff 偏移方法，壓制了由于數(shù)據(jù)不完整導(dǎo)致的成像腳印.當(dāng)前，LSM偏移算子已擴(kuò)展到雙程波動(dòng)方程(Dai et al., 2012)、多震源偏移(Dai, 2012)、彈性多參數(shù)偏移(Duan et al., 2016; Ren and Li, 2020)、成像域最小二乘偏移(IDLSM)(Fletcher et al., 2016; Guo and Wang, 2020).

然而，LSM仍然處于試驗(yàn)探索階段(楊勤勇和段心標(biāo), 2018)，面臨多種問題有待解決.當(dāng)前LSM存在兩個(gè)主要問題難以避免：(1) LSM需要迭代執(zhí)行偏移和反偏移算子，迭代計(jì)算成本過高；(2) LSM要求背景模型是準(zhǔn)確的，否則迭代收斂性較差.這兩個(gè)問題導(dǎo)致LSM在大規(guī)模應(yīng)用中難以執(zhí)行.

為了解決背景模型依賴問題，同時(shí)減少迭代次數(shù)，擴(kuò)展域Born偽逆算子在近期獲得了發(fā)展(Symes, 2008; ten Kroode, 2012; 劉玉金和李振春, 2015; Hou and Symes, 2015; Chauris and Cocher, 2017)，且當(dāng)前已推廣到2D變密度聲波介質(zhì)(Farshad and Chauris, 2020).該偽逆Born算子是在地下偏移距域逆時(shí)偏移(RTM)的基礎(chǔ)上，近似建立的一種地震數(shù)據(jù)與擴(kuò)展域模型之間的轉(zhuǎn)換算子.擴(kuò)展域模型與地震數(shù)據(jù)維度相同，可以保存大部分的原始地震數(shù)據(jù)信息.因此，該擴(kuò)展域方法能夠精確重建原始地震數(shù)據(jù)，即使背景模型是錯(cuò)誤的.同時(shí)，該擴(kuò)展域模型可通過Radon變換提取角度域信息.

擴(kuò)展域Born偽逆算子的主要弊端在于計(jì)算量明顯高于常規(guī)RTM.其計(jì)算成本與地下偏移距方向的離散采樣數(shù)量線性相關(guān).如果背景模型是已知且準(zhǔn)確的，該離散采樣數(shù)量可以減少，計(jì)算成本可以緩解(Hou and Symes, 2015).然而，對(duì)于3D模型情況，該擴(kuò)展域方法計(jì)算耗時(shí)將是巨大的，因?yàn)?D中的地下偏移距會(huì)在x和y方向同時(shí)離散采樣，采樣數(shù)量是2D模型的平方.因此，該方法也不適合未來3D大規(guī)模推廣.

GRT反演是另一種基于Born(Kirchhoff)近似的線性反演方法.與LSM相比，GRT反演是一種非迭代高效的線性反演方法，且已推廣到彈性多參數(shù)反演.Beylkin(1984, 1985)通過構(gòu)造 Beylkin 行列式建立了聲波單參數(shù)GRT逆變換算子理論.Miller等(1987)進(jìn)一步分析了GRT的反演特征與應(yīng)用.Bleistein(1987)基于Kirchhoff近似將GRT推廣到反射系數(shù)反演多參數(shù).之后，多參數(shù)GRT反演方法獲得發(fā)展.Beylkin和Burridge(1990)發(fā)展了聲波和彈性多參數(shù)GRT反演.Ikelle等(1992)進(jìn)一步分析了阻抗和速度擾動(dòng)參數(shù)反演分布情況.后來多參數(shù) GRT 反演推廣到彈性各向異性介質(zhì)(de Hoop et al., 1994, 1999; de Hoop and Bleistein, 1997; Burridge et al., 1998; Ursin, 2004).另外，Protasov和Tcheverda(2011)結(jié)合GRT和高斯束方法發(fā)展了真振幅成像.為解決大擾動(dòng)參數(shù)問題，基于GRT的雙散射反演方法獲得了發(fā)展(Ouyang et al., 2015, 2017; Ouyang and Mao, 2018).

然而，在過去20年GRT反演方法并未獲得充足的研究與推廣.現(xiàn)有GRT中的多參數(shù)提取存在明顯的不確定性和不穩(wěn)定性問題.多參數(shù)提取算子是一N階方陣的逆(N為提取參數(shù)數(shù)量)，與散射角積分直接相關(guān).而該散射角積分范圍隨方位角、傾角變化，在實(shí)際中是不確定且非唯一的.Beylkin和Burridge(1990)通過出現(xiàn)的最大散射角來確定散射角范圍，這種方法提取的多參數(shù)往往不準(zhǔn)確.同時(shí)，該多參數(shù)提取算子也經(jīng)常不穩(wěn)定.受參數(shù)數(shù)量影響，該提取算子常為病態(tài)，尤其是參數(shù)數(shù)量過多時(shí)(如三個(gè)以上).因此，利用現(xiàn)有多參數(shù) GRT 反演方法獲取精確和穩(wěn)定的反演結(jié)果非常困難.另外，與LSM類似，現(xiàn)有GRT反演同樣要求背景模型準(zhǔn)確.

為解決現(xiàn)有線性反演方法存在的主要問題，本文提出了角度域GRT(AD-GRT)多參數(shù)反演方法(2D聲波介質(zhì)).角度域信息是當(dāng)前地震成像與數(shù)據(jù)分析中的最有效信息(張宇, 2018)，可直接應(yīng)用于速度建模和AVA分析，且明顯降低成像與反演的非唯一性(Xu et al., 2001).我們利用局部走時(shí)和幅值近似搭建了角度域積分變換方程，實(shí)現(xiàn)了原始地震數(shù)據(jù)與角度域模型之間的精確變換.與擴(kuò)展域Born偽逆算子類似，該角度域模型也是一種擴(kuò)展域模型，可支持精確的地震數(shù)據(jù)重建，即使背景模型不準(zhǔn)確.角度域模型可用于提取準(zhǔn)確的擾動(dòng)多參數(shù).不同于擴(kuò)展域Born偽逆算子，AD-GRT計(jì)算耗時(shí)和常規(guī)GRT反演耗時(shí)相當(dāng).角度方向的采樣點(diǎn)數(shù)量對(duì)AD-GRT的計(jì)算量影響微弱.為提高角度域信息的精度，AD-GRT的角度方向采樣比常規(guī)角度域成像方法(Zhang et al., 2007; Tang et al., 2013)要明顯密集.同時(shí)，為解決傳統(tǒng)角度域離散方法(Xu et al., 2001)存在的幅值震蕩問題，我們?cè)O(shè)計(jì)了兩種有效的角度域離散分裂方法.實(shí)驗(yàn)驗(yàn)證了新提出的離散方法在提取平滑連續(xù)的角度域信息方面的優(yōu)勢(shì).另外，為提高AD-GRT反演的分辨率，我們將震源子波反濾波用于反演.

1 聲波廣義Radon正變換

該部分回顧在2D聲波介質(zhì)中基于Born近似的GRT正變換.2D非均勻聲波介質(zhì)頻率域無源波動(dòng)方程可表示為：

(1)

其中p(x,ω)為點(diǎn)x=(x1,x3)的聲波壓力波場(chǎng)，κ(x)=1/(ρc2)為壓縮率，σ(x)=1/ρ為密度倒數(shù)，c=c(x)和ρ=ρ(x)分別為波速和密度.在2D情況下，下標(biāo)j=1,3.本文公式中僅對(duì)全局坐標(biāo)系中的方向性下標(biāo)(即x1,x3方向)執(zhí)行愛因斯坦求和約定，其他上下標(biāo)符號(hào)不執(zhí)行求和操作.依照散射理論，如果將原始介質(zhì)屬性參數(shù)分解為背景參數(shù)和擾動(dòng)參數(shù)，則原始波場(chǎng)可分解為背景波場(chǎng)和散射波場(chǎng).所以設(shè)置參數(shù)：

σ(x)=σ0(x)+σ′(x),κ(x)=κ0(x)+κ′(x),

(2)

其中σ0和κ0為背景參數(shù)，一般是已知的，σ′和κ′為對(duì)應(yīng)的擾動(dòng)參數(shù).這種情況下，原始波場(chǎng)可分解為：

p(x,ω)=p0(x,ω)+p′(x,ω)，

(3)

其中p0和p′分別為背景波場(chǎng)和散射波場(chǎng).假設(shè)擾動(dòng)參數(shù)非常小，則散射波場(chǎng)以單散射為主，多散射波場(chǎng)可以忽略.基于聲波格林函數(shù)和Born近似，散射波場(chǎng)線性方程可表示為：

(4)

其中，G0=G0(x,r,ω)為背景模型中的聲波格林函數(shù)，且具有互易屬性G0(r,x,ω)=G0(x,r,ω)，r=(r1,r3)為散射波場(chǎng)接收點(diǎn).

基于式(4)，對(duì)應(yīng)的2D聲波GRT正變換可獲取(Beylkin和Burridge, 1990).式(4)中背景波場(chǎng)p0與震源位置和子波有關(guān)，在此用子波與格林函數(shù)F(ω)G0(x,s,ω)代替p0(x,ω)，其中F(ω)為震源子波頻率域公式，對(duì)應(yīng)的時(shí)間域子波用F(t)表示.如果子波為脈沖子波，則有F(ω)=1.本文中的時(shí)間和頻率域(或空間和波數(shù)域)函數(shù)使用相同函數(shù)符號(hào)表示，通過函數(shù)變量具體區(qū)分.將背景模型中的格林函數(shù)表示為：

(5)

其中，s=(s1,s3)為背景波場(chǎng)震源點(diǎn)，下標(biāo)s和r分別代表震源端和接收端，As和φs分別代表震源端振幅和走時(shí)，Ar和φr分別為接收端振幅和走時(shí).基于高頻近似，格林函數(shù)的空間偏導(dǎo)數(shù)可近似表示為：

(6)

將式(5)、(6)代入式(4)，整理可得：

×F(ω)f(x,θ)，

(7)

其中，f(x,θ)為一種角度域模型：

(8)

依據(jù)式(8)，可以利用兩個(gè)擾動(dòng)參數(shù)f1和f2來線性描述f(x,θ)，即：

(9)

θ=θ(s,x,r)為散射角，是s到x射線和r到x射線在x處的夾角，且滿足：

θ=αs-αr,

(10)

×F(t-φs-φr)，

(11)

其中，A(s,x,r)為振幅因子，其與As和Ar滿足|ω|AsArκ0=isgn(ω)A(s,x,r)，此處sgn(ω)為符號(hào)函數(shù).A(s,x,r)可用幾何擴(kuò)散物理量表示：

(12)

其中，q2為一個(gè)動(dòng)力學(xué)射線參數(shù).ρ0和c0分別為背景密度和波速.由于常規(guī)射線理論在焦散區(qū)誤差明顯，本文臨時(shí)忽略焦散區(qū)影響.

式(7)、(11)分別為頻率域和時(shí)間域2D聲波廣義Radon正變換.該變換相當(dāng)于角度域模型f(x,θ)到波場(chǎng)數(shù)據(jù)p′(r,s,t)的積分變換，此處的s和r僅分布在2D擾動(dòng)區(qū)域的邊界線上.理論上，f(x,θ)和p′(r,s,t)的均具有三個(gè)維度分布，因此，對(duì)應(yīng)的從p′(r,s,t)到f(x,θ)的逆變換也可搭建.

2 角度域廣義Radon逆變換

該部分在式(7)、(11)的基礎(chǔ)上，建立對(duì)應(yīng)的角度域廣義Radon逆變換(AD-GRT)，以支持從p′(r,s,t)到f(x,θ)的反向積分變換，并在此基礎(chǔ)上實(shí)現(xiàn)聲波多參數(shù)提取和數(shù)據(jù)重建擬合.

AD-GRT可以在Fourier變換的基礎(chǔ)上建立起來.f(x,θ)的2D Fourier變換可寫為：

(13a)

×exp[ikj(xj-yj)]，

(13b)

其中，k=(k1,k3)為2D波數(shù)矢量.x和y分別為擾動(dòng)點(diǎn)和反演點(diǎn).為了在式(13b)內(nèi)部組建與式(7)相同的表示形式，需要對(duì)比分析式(13b)、(7).式(7)、(13b)均包含積分項(xiàng)dx1dx3f(x,θ)，且分別包含復(fù)相位項(xiàng)exp[iω(φs+φr)]和exp[ikj(xj-yj)].剩余部分分別與(r,s,ω)和dk1dk3相關(guān).為了組建p′(r,s,ω)的逆變換，需要建立(r,s,ω)和dk1dk3之間的轉(zhuǎn)換關(guān)系.基于以上分析，下面通過以下三步完成建立AD-GRT.

(14)

k′和νj的取值范圍主要有兩種有效設(shè)置，第一種為0≤k′<∞，νj分布在整個(gè)單位元，第二種為-∞

(15a)

(15b)

其中psj(y)和prj(y)為慢度矢量ps(y)和pr(y)的分量.如果設(shè)置：

(16)

則可得出：

(17a)

(17b)

由式(17a)、(17b)可知，k′/ω是慢度矢量ps+pr的模，νj是ps+pr的單位方向矢量，也是走時(shí)φs(y)+φr(y)等值線傾角.

將式(15)、(17a)代入式(14)，整理可得：

×f(x,θ)exp[iω(φs(x)+φr(x))]

×exp[-iω(φs(y)+φr(y))]，

(18)

其中：

(19)

第三，建立積分項(xiàng)drdsdω相關(guān)的積分表達(dá)式.式(18)現(xiàn)有的積分項(xiàng)dωdν是二維空間積分，不能直接轉(zhuǎn)換到drdsdω.同時(shí)，傾角散射角坐標(biāo)系(ν,θ)可以映射到邊界線上的震源接收點(diǎn)坐標(biāo)系(s,r).所以式(18)對(duì)于組建p′(r,s,ω)到f(y,θ)逆變換尚缺少散射角積分項(xiàng)dθ.為此，對(duì)f(y,θ)添加一層積分dθ：

(20)

其中，θ0代表散射角θ的采樣點(diǎn).將式(18)代入式(20)，可得：

×f(x,θ)exp[iω(φs(x)+φr(x))]

×exp[-iω(φs(y)+φr(y))].

(21)

在2D單位圓上，積分項(xiàng)滿足轉(zhuǎn)換dνdθ=dαsdαr，其中，αs(y)和αr(y)為可映射震源s和接收點(diǎn)r的方向角，且滿足雅可比轉(zhuǎn)換：

dαs=Js(y)ds, dαr=Jr(y)dr，

(22)

dνdθ=Js(y)Jr(y)dsdr.

(23)

將關(guān)系式(23)代入式(21)，整理可得：

×Js(y)Jr(y)f(x,θ)exp[iω(φs(x)+φr(x))]

exp[-iω(φs(y)+φr(y))]，

(24)

至此，式(24)已包含了式(7)中相關(guān)的各項(xiàng)積分項(xiàng)和復(fù)相位項(xiàng).在此基礎(chǔ)上，補(bǔ)充其他因子項(xiàng)，即可建立逆變換.假設(shè)x在y附近周圍，則振幅滿足局部近似：

As(x,ω)≈As(y,ω),Ar(x,ω)≈Ar(y,ω).

(25)

對(duì)比式(24)、(7)，逆變換可建立：

(26)

類似于式(11)，頻率域式(26)可以轉(zhuǎn)換為對(duì)應(yīng)的時(shí)間域方程.設(shè)定：

(27)

(28)

式(28)即為2D聲波AD-GRT，與式(11)為對(duì)應(yīng)的p′(r,s,t)和f(y,θ0)之間的積分變換.在AD-GRT建立過程中，我們?cè)趦商幨褂昧司植拷?式(15)、(25))，只要輸入數(shù)據(jù)具有較好的連續(xù)性，這種近似導(dǎo)致的誤差是微弱的.基于反演的f(y,θ0)，聲波雙擾動(dòng)參數(shù)可以利用線性關(guān)系式(8)和最小二乘擬合方法獲得有效提取，原始波場(chǎng)數(shù)據(jù)可以利用正變換式(11)有效重建擬合.

3 角度域離散方法

雖然當(dāng)前已獲得AD-GRT最終表達(dá)式(28)，但是式中的δ(θ-θ0)并不方便直接參與計(jì)算，需要利用其他的近似插值或求和算法代替.并且，θ0離散采樣點(diǎn)與震源s和接收點(diǎn)r的離散分布并不能直接對(duì)應(yīng).這給AD-GRT的精確數(shù)值計(jì)算帶來了一定麻煩.為此，我們?cè)诔Ｒ?guī)角度域離散算法的基礎(chǔ)上，設(shè)計(jì)AD-GRT的準(zhǔn)確離散算法.角度積分表達(dá)式(20)可以近似表示為：

(29)

為了解決傳統(tǒng)角度域算法的離散問題，我們對(duì)θ相關(guān)的積分求和做進(jìn)一步分析，設(shè)計(jì)合理的數(shù)值求解方案.離散后的微元滿足關(guān)系式：

Δαs≈JsΔs,Δαr≈JrΔr,

(30)

其中Δs和Δr是炮點(diǎn)和接收點(diǎn)離散間隔，實(shí)際計(jì)算中是已知的.Δαs和Δαr分別為s和r傳播到y(tǒng)的射線方向間隔，可通過雅克比系數(shù)Js和Jr計(jì)算求得.因此，一個(gè)(s,r)的離散采樣點(diǎn)在y對(duì)應(yīng)的是一個(gè)離散單元，單元面積大小為Δαs×Δαr.而常規(guī)算法忽略該單元面積大小，僅依靠采樣點(diǎn)上的θ值來判斷累加位置是不合理的.當(dāng)采樣點(diǎn)θ非?？拷謪^(qū)邊界θ0±Δθ/2時(shí)，該單元有一部分應(yīng)屬于θ0附近分區(qū).這也是常規(guī)角度域計(jì)算方法出現(xiàn)幅值起伏震蕩的主要原因.

基于以上分析，我們將分區(qū)邊界θ0±Δθ/2附近的離散單元分裂為兩個(gè)離散單元，并求解分裂后的兩個(gè)離散單元的面積大小S1和S2.精確計(jì)算兩個(gè)分裂單元大小的方法可能有許多，在此我們?cè)O(shè)計(jì)兩種不同的離散單元分裂算法.

3.1 矩形分裂法

如圖1所示，每一個(gè)(αs,αr)采樣點(diǎn)的離散單元由一個(gè)矩形表示，矩形面積大小為Δαs×Δαr，采樣點(diǎn)位于長(zhǎng)和寬的中間位置.θbound=θ0±Δθ/2為θ0分區(qū)邊界.θmid為(αs,αr)采樣點(diǎn)對(duì)應(yīng)的θ值.當(dāng)θbound穿過一個(gè)離散單元時(shí)，將該單元大小拆分為S1和S2兩部分，并代替式(28)離散算法中的JsJrΔsΔr部分，分別累加到θ0及其附近的f(y,θ0)反演結(jié)果中.其中的S1和S2通過以下方式求解.

圖1 矩形離散單元與散射角等值線分布

如果ds≥dr，S1有表達(dá)式：

(31)

如果ds

(32)

基于S1，S2有表達(dá)式：

S2=dsdr-S1.

(33)

以上即為離散單元矩形分裂法.主要依據(jù)平面幾何方法分解矩形單元，形象且計(jì)算準(zhǔn)確.但該方法也有缺點(diǎn)：(1)多次分支對(duì)比影響計(jì)算耗時(shí)；(2)在較高維度(四維空間)推廣復(fù)雜.為此，我們?cè)O(shè)計(jì)了便于使用和推廣的第二種分裂法.

3.2 平行分裂法

圖2 平行四邊形離散單元與散射角等值線分布

(34a)

(34b)

此處同樣滿足S1+S2=dsdr.該方法也可以通過αs單獨(dú)控制離散單元，調(diào)換公式中的所有dr和ds即可實(shí)現(xiàn).該分裂方法計(jì)算簡(jiǎn)潔方便，且適用于未來高維離散計(jì)算.

4 數(shù)值試驗(yàn)

我們利用三種模型數(shù)據(jù)來驗(yàn)證我們提出的AD-GRT在角度域模型反演、多參數(shù)提取、數(shù)據(jù)重建擬合的有效性.同時(shí)，我們也測(cè)試提出的離散單元分裂法的計(jì)算效果，子波反濾波對(duì)分辨率的影響，以及錯(cuò)誤背景模型下的數(shù)據(jù)重建擬合.

三個(gè)模型均選擇Ricker子波作為震源子波(Wang, 2015)，子波主頻率fM均設(shè)為15 Hz.第一個(gè)模型地震記錄使用波動(dòng)方程有限差分法合成，另外兩個(gè)利用GRT正變換(式(11))合成，以便驗(yàn)證AD-GRT正反變換的準(zhǔn)確性.角度域模型在角度方向的采樣間隔均設(shè)置為Δθ=1.5°，在[0°, 180°]范圍內(nèi)共取120個(gè)采樣.負(fù)角度范圍在本文臨時(shí)不做分析.

4.1 水平單界面模型

水平單界面模型如圖3a所示.水平界面在深度z=1000 m，并且界面上下聲波波速分別為2000 m·s-1和2100 m·s-1，依照Gardner關(guān)系式ρ=0.31×c0.25(波速單位: m·s-1, 密度單位: g·cm-3)設(shè)定模型密度.模型在水平和垂直方向的網(wǎng)格間隔均為5 m.炮點(diǎn)和檢波點(diǎn)均分布于邊界線z=0 m上，且炮點(diǎn)間隔和檢波點(diǎn)間隔均為10 m，從炮點(diǎn)到檢波點(diǎn)的最大偏移距為2000 m.該模型地震記錄使用聲波波動(dòng)方程有限差分方法合成，以驗(yàn)證準(zhǔn)確波場(chǎng)數(shù)據(jù)反演結(jié)果的有效性.用于AD-GRT的背景波速和密度模型通過平滑原始模型來獲得.圖3b顯示了震源x=2000 m的單炮記錄，這里只顯示了主要反射信息，排除了直達(dá)波信息.

圖3 (a) 水平單界面模型； (b) 合成單炮地震記錄

圖4給出了利用AD-GRT反演方法(式(28))求解的角度域模型f(x,θ)反演結(jié)果(x=2000 m).為排除濾波導(dǎo)致的幅值震蕩，此處設(shè)置震源子波為F(ω)=1，忽略子波反濾波操作.反演的主要信息分布在水平線z=1000 m附近，且同相軸具有明顯的水平分布特征.同時(shí)，在反演結(jié)果邊緣位置存在明顯的傾斜噪聲(見區(qū)域A)，這與地震記錄最大偏移距邊界處的數(shù)據(jù)截?cái)嘤嘘P(guān).記錄邊界具有數(shù)據(jù)不連續(xù)特征，而不連續(xù)性不符合提供的AD-GRT的局部近似假設(shè)(式(15)、(25)).

圖4 角度域模型反演(忽略子波反濾波)

為了驗(yàn)證提供的離散單元分裂法在數(shù)值計(jì)算方面的有效性，我們分別利用傳統(tǒng)離散方法、矩形分裂法和平行分裂法數(shù)值計(jì)算反演結(jié)果(如圖4).三種反演結(jié)果在水平線z=1000 m上的幅值分布如圖5所示.圖5的橫軸方向使用角度的余弦值表示，以顯示角度域模型的線性傾斜分布特征(式(8)).其中的真實(shí)值是界面z=1000 m上下擾動(dòng)值的差值(跳躍值)，并且為方便對(duì)比，反演值乘了一個(gè)常數(shù).由此可見，傳統(tǒng)離散方法在角度域信息提取中易形成明顯的鋸齒狀震蕩，尤其在小角度范圍內(nèi).這種震蕩在地震記錄稀疏觀測(cè)情況下將會(huì)更為嚴(yán)重.相比之下，我們提出的兩種離散單元分裂法均能提供非常連續(xù)平滑的幅值分布效果，并且反演幅值分布與模型真實(shí)值分布在有效角度范圍內(nèi)均能高度擬合.由此可證明AD-GRT反演結(jié)果是準(zhǔn)確的，并且離散單元分裂法是有效的.

圖5 角度域反演幅值分布

聲波擾動(dòng)參數(shù)f1和f2(式(8))可以通過數(shù)據(jù)擬合從角度域模型(圖3)中提取.圖6顯示了雙參數(shù)提取結(jié)果.圖中f1提取結(jié)果非常靠近真實(shí)值(峰值相對(duì)誤差0.46%，絕對(duì)誤差0.00048).f2提取結(jié)果相對(duì)誤差比較明顯(峰值相對(duì)誤差7.3%，絕對(duì)誤差0.00078).這主要與兩點(diǎn)有關(guān)：(1)f2絕對(duì)值遠(yuǎn)小于f1，在雙參數(shù)共同提取過程中易形成相對(duì)誤差，而f1和f2的絕對(duì)誤差相差并不明顯.(2)f2提取的是幅值分布的斜率或?qū)?shù)，f1主要與幅值平均分布有關(guān)(見式(8))，而幅值分布的斜率(導(dǎo)數(shù))更容易受高頻震蕩噪聲的干擾.因此，圖6驗(yàn)證了基于角度域模型反演的多參數(shù)提取是合理的.

圖6 角度域反演中的多參數(shù)提取結(jié)果

如圖7所示，基于角度域模型反演結(jié)果(圖4)，地震數(shù)據(jù)可以利用GRT正變換(式(11))進(jìn)行重建.圖8顯示了不同偏移距的地震道對(duì)比.由此可知，重建數(shù)據(jù)中的主要擬合誤差是由數(shù)據(jù)截?cái)鄬?dǎo)致的，該誤差與圖4區(qū)域A處的傾斜噪聲相互對(duì)應(yīng).如果角度域中的傾斜噪聲可以消除，則在重建數(shù)據(jù)中的擬合誤差也可以壓制.除了該部分?jǐn)M合誤差外，重建數(shù)據(jù)與原始數(shù)據(jù)可以高精度擬合.由此可見，我們提供的AD-GRT方法可在非迭代情況下實(shí)現(xiàn)高精度數(shù)據(jù)重建.

圖7 利用角度域反演重建的地震數(shù)據(jù)

圖8 重建地震數(shù)據(jù)道

該測(cè)試驗(yàn)證了我們提出的AD-GRT在角度域模型反演、多參數(shù)提取、數(shù)據(jù)重建的準(zhǔn)確性，即使原始數(shù)據(jù)是由波動(dòng)方程有限差分法合成.同時(shí)，該測(cè)試也證明了提供的離散單元分裂法在去除傳統(tǒng)角度域離散方法存在的鋸齒狀震蕩的有效性.

4.2 水平層狀模型

圖9顯示了水平層狀模型的波速和密度模型.為驗(yàn)證AD-GRT正逆變換的準(zhǔn)確性，由此以后的模型測(cè)試均使用GRT正變換(式(11))合成地震記錄.圖9a、b顯示了波速和密度原始模型，背景模型通過平滑原始模型獲得.圖9c、d為擾動(dòng)參數(shù)模型.由此可知，界面z=1000 m和z=1300 m分別為f1和f2單參數(shù)界面.模型網(wǎng)格間隔和震源接收點(diǎn)排布情況均與之前的水平單界面模型相同.圖10顯示了合成地震記錄(震源x=2000 m).

圖9 水平層狀模型

圖10 合成單炮地震記錄

圖11顯示了角度域模型反演結(jié)果，其中(a)忽略了子波反濾波操作，而(b)執(zhí)行了該操作.圖11a中的紅點(diǎn)是界面上下真實(shí)值的差值.對(duì)比圖11a、b可知，子波反濾波明顯提高了反演結(jié)果的分辨率，但也形成了一定程度的震蕩，這對(duì)穩(wěn)定性有一定影響.圖11a、b中的反演值均能高精度的擬合真實(shí)值.另外，淺層反演結(jié)果形成了明顯的傾斜噪聲，這將影響之后的參數(shù)提取和數(shù)據(jù)重建.

圖12顯示了從角度域反演結(jié)果(圖11)中提取的多參數(shù)(式(8)).圖12a中的f2提取值誤差比f1明顯，圖12b中的f2提取值震蕩也較嚴(yán)重，這與之前的水平單界面模型情況相似.除了這些誤差和震蕩，圖12a、b的參數(shù)提取值和真實(shí)值都能較好的擬合.

地震數(shù)據(jù)可以利用角度域模型反演結(jié)果(圖11)進(jìn)行重建(式(11))，如圖13所示.圖13a、b分別為忽略和使用子波反濾波情況下的重建數(shù)據(jù)，兩種情況的重建擬合誤差基本相同，這說明震源子波是否已知對(duì)AD-GRT的數(shù)據(jù)重建沒有明顯影響，這與AD-GRT正逆變換中存在震源子波F(ω)補(bǔ)償有關(guān).重建數(shù)據(jù)的主要誤差為淺層強(qiáng)震動(dòng)噪聲，這與角度域反演結(jié)果(圖11)中的淺層傾斜噪聲直接相關(guān).這是由于震源和接收點(diǎn)附近的波場(chǎng)不符合局部近似假設(shè)(式(15)、(25))，易形成畸形震動(dòng).除了淺層和邊界的噪聲以外，兩種重建數(shù)據(jù)均能精確擬合原始地震數(shù)據(jù).圖14給出了忽略子波反濾波的不同偏移距的重建數(shù)據(jù)和原始數(shù)據(jù)對(duì)比，大部分位置重建和原始數(shù)據(jù)都能高精度擬合.

圖11 角度域模型反演

圖12 參數(shù)提取

圖13 重建地震數(shù)據(jù)

圖14 重建地震數(shù)據(jù)道

該模型測(cè)試驗(yàn)證了角度域模型與原始地震數(shù)據(jù)之間的對(duì)等關(guān)系，角度域模型保留了原始地震數(shù)據(jù)的絕大多數(shù)有效信息，并可支持穩(wěn)定的多參數(shù)提取和精確的數(shù)據(jù)重建.

4.3 Marmousi2模型

最后，我們提供復(fù)雜的Marmousi2模型測(cè)試，以驗(yàn)證AD-GRT在復(fù)雜模型中的有效性，討論焦散區(qū)存在的問題，測(cè)試錯(cuò)誤背景模型下的數(shù)據(jù)重建.圖15a、b顯示了波速和密度模型.地震記錄使用GRT正變換(式(11))合成.震源和接收點(diǎn)分布在地表z=0 m，間隔均為10 m.最大偏移距為4000 m.圖16顯示了震源x=4000 m和x=8000 m的地震記錄.可以看出，震源x=8000 m的地震記錄分布明顯復(fù)雜，這與震源附近區(qū)域模型結(jié)構(gòu)復(fù)雜有關(guān).

圖15 Marmousi2模型

圖16 合成地震記錄

圖17給出了位于x=6000 m和x=10000 m的角度域模型反演結(jié)果，該反演測(cè)試執(zhí)行了子波反濾波以提高反演分辨率.由圖可知，x=6000 m的角度域模型反演值能與真實(shí)值很好的擬合，而x=10000 m的反演值在許多位置不能匹配真實(shí)值.這與常規(guī)射線理論(ART)的應(yīng)用有關(guān).當(dāng)前的AD-GRT利用ART獲取波場(chǎng)傳播信息，易在復(fù)雜結(jié)構(gòu)區(qū)域形成焦散區(qū)，導(dǎo)致大量幅值和走時(shí)信息不連續(xù).該不連續(xù)性不符合AD-GRT的局部近似假設(shè)(式(15)、(25)).在x=10000 m附近介質(zhì)結(jié)構(gòu)復(fù)雜，因此這里的反演結(jié)果誤差明顯.

圖18顯示了從角度域模型反演值(圖17)中提取的多參數(shù)，圖19為x=6000 m和x=10000 m的多參數(shù)提取值與真實(shí)值的對(duì)比.由圖可知，提取的f1參數(shù)大部分都能非常靠近真實(shí)值，在復(fù)雜結(jié)構(gòu)附近誤差比較明顯，這是角度域反演值(圖17)誤差導(dǎo)致的.而f2參數(shù)多數(shù)不能靠近真實(shí)值，該誤差的形成原因與之前模型類似：(1) 一般f1參數(shù)值是f2的5倍以上，相對(duì)誤差明顯.(2)f2提取的是反演值對(duì)角度余弦的斜率(導(dǎo)數(shù))，易受局部震蕩噪聲的干擾.所以f2提取的穩(wěn)定性較差.有關(guān)參數(shù)提取的優(yōu)化技術(shù)有必要進(jìn)行進(jìn)一步研究.

圖17 角度域模型反演(使用子波反濾波)

圖18 參數(shù)提取值對(duì)比

圖19 參數(shù)提取值對(duì)比

圖20顯示了利用角度域模型反演值重建的震源x=4000 m和x=8000 m的地震數(shù)據(jù).由圖20a、c可知，震源x=4000 m數(shù)據(jù)的主要噪聲是由角度域反演值的淺層傾斜噪聲導(dǎo)致的.另外，地震記錄數(shù)據(jù)不連續(xù)也會(huì)造成重建誤差.除此之外，震源x=4000 m的數(shù)據(jù)大部都能擬合原始地震數(shù)據(jù).震源x=8000 m數(shù)據(jù)受復(fù)雜結(jié)構(gòu)焦散區(qū)影響，波場(chǎng)傳播幅值和走時(shí)不連續(xù)，導(dǎo)致大部分重建數(shù)據(jù)都有明顯誤差.

圖20 重建地震數(shù)據(jù)

為了測(cè)試錯(cuò)誤背景模型對(duì)角度域模型反演和數(shù)據(jù)重建的影響，同時(shí)分析波場(chǎng)傳播信息的不連續(xù)性，我們利用錯(cuò)誤背景模型進(jìn)行角度域模型反演和數(shù)據(jù)重建.如圖15c、d所示，錯(cuò)誤背景模型是簡(jiǎn)單線性梯度模型.圖21顯示了錯(cuò)誤背景模型下的角度域模型反演值.對(duì)比正確背景模型反演值(圖17)，x=6000 m的同相軸明顯彎曲，而x=10000 m的角度域信息是混亂的.即使如此，該反演值也能保存大部分原始地震數(shù)據(jù)信息.

圖22提供了利用錯(cuò)誤角度域反演值(圖21)重建的地震數(shù)據(jù).有圖22a、c可知，除了淺層的誤差以外，大部分重建數(shù)據(jù)能夠很好的擬合原始數(shù)據(jù).圖22b、d的重建數(shù)據(jù)誤差也不明顯，即使該震源數(shù)據(jù)是在復(fù)雜結(jié)構(gòu)區(qū)域合成的.對(duì)比正確背景模型下的重建數(shù)據(jù)(圖20)，圖22的重建數(shù)據(jù)反而能夠更好的擬合原始數(shù)據(jù).這是因?yàn)樵阱e(cuò)誤的背景模型中不存在焦散區(qū)和傳播信息不連續(xù)問題.圖22淺層噪聲的分布范圍比圖20偏大，這與背景模型近地表區(qū)域的波速較大有關(guān).圖23給出了使用正確和錯(cuò)誤背景模型重建的數(shù)據(jù)對(duì)比(震源x=8000 m，偏移距1200 m).由圖23可知，利用錯(cuò)誤但簡(jiǎn)單的背景模型重建的數(shù)據(jù)比利用正確但復(fù)雜的背景模型重建的數(shù)據(jù)能更好的擬合原始數(shù)據(jù).

圖21 角度域模型反演(使用錯(cuò)誤背景模型)

圖22 重建地震數(shù)據(jù)(利用錯(cuò)誤背景模型)

圖23 重建地震數(shù)據(jù)對(duì)比(震源x=8000 m, 偏移距1200 m)

Marmousi2模型測(cè)試說明了重建地震數(shù)據(jù)的精確度主要與波場(chǎng)傳播信息的連續(xù)性有關(guān)，而背景模型是否正確對(duì)數(shù)據(jù)重建并不重要.這也很好的驗(yàn)證了角度域模型能保留絕大部分原始地震信息，不管該角度域模型反演值是否正確.

5 討論

三個(gè)數(shù)值模型驗(yàn)證了新提出的AD-GRT反演方法的有效性.角度域模型的反演結(jié)果與真實(shí)值能夠高精度擬合，雙擾動(dòng)參數(shù)可獲得有效穩(wěn)定的提取，重建地震數(shù)據(jù)能與原始數(shù)據(jù)很好的擬合.

AD-GRT不要求背景模型是準(zhǔn)確的.當(dāng)背景模型是錯(cuò)誤的，角度域模型反演結(jié)果也是錯(cuò)誤的，多參數(shù)不能有效提取.但是，該角度域模型依舊很好的保存了原始地震數(shù)據(jù)的信息，能夠非常精確的重建地震數(shù)據(jù).并且，盡管該角度域模型反演結(jié)果是錯(cuò)誤的，其包含的角度域信息可直接用于進(jìn)一步改進(jìn)背景速度模型.

AD-GRT引入了子波反濾波(式(27)).如果地震數(shù)據(jù)的震源子波已知，子波反濾波操作能夠明顯提高角度域模型反演和多參數(shù)提取的分辨率和精度.但是子波反濾波對(duì)地震數(shù)據(jù)重建沒有明顯影響.如果子波未知，使用脈沖函數(shù)代替子波，地震數(shù)據(jù)仍然能夠精確重建.這與式(7)、(26)中子波的正反濾波有關(guān).

另外，AD-GRT當(dāng)前的主要弊端在于傳統(tǒng)射線理論(ART)的使用.ART的高頻近似和介質(zhì)平滑假設(shè)不符合復(fù)雜結(jié)構(gòu)模型特征.AD-GRT的許多反演噪聲與ART在焦散區(qū)生成的非連續(xù)波場(chǎng)傳播信息有關(guān)，該非連續(xù)性不符合AD-GRT中的走時(shí)和幅值局部近似(式(15)、(25)).因此，有必要討論頻率相關(guān)的波場(chǎng)傳播問題，為AD-GRT反演方法提供高精度的波場(chǎng)傳播信息.當(dāng)前，頻率相關(guān)的射線理論已獲得一定程度的發(fā)展(朱天飛, 2007; Protasov and Gadylshin, 2017).我們希望在此基礎(chǔ)上進(jìn)一步發(fā)展適用于AD-GRT的波場(chǎng)傳播方法.

6 結(jié)論

本文發(fā)展了一種新型的廣義Radon變換(GRT)反演方法(2D聲波介質(zhì))，即角度域GRT(AD-GRT).利用局部走時(shí)和幅值近似搭建的AD-GRT實(shí)現(xiàn)了原始地震數(shù)據(jù)與角度域模型之間的準(zhǔn)確轉(zhuǎn)換.角度域模型能夠記錄大部分地震數(shù)據(jù)信息，支持精確的地震數(shù)據(jù)重建，即使背景模型不準(zhǔn)確.并且，該角度域模型支持與角度相關(guān)的多擾動(dòng)參數(shù)的準(zhǔn)確穩(wěn)定提取.同時(shí)，本文提出的兩種角度域離散方法有效解決了傳統(tǒng)角度域離散方法存在的幅值震蕩問題.震源子波反濾波的引入明顯提高了AD-GRT的反演分辨率.另外，AD-GRT保持了常規(guī)GRT方法的非迭代高效運(yùn)行特征，角度域的添加對(duì)AD-GRT的計(jì)算耗時(shí)影響微弱，這與現(xiàn)有的地下偏移距擴(kuò)展域反演方法相比具有明顯優(yōu)勢(shì).因此，新提出的AD-GRT有效解決了現(xiàn)有線性反演存在的計(jì)算成本過高、背景模型依賴、多參數(shù)提取不確定性等問題.在此基礎(chǔ)上，我們下一步將擴(kuò)展AD-GRT反演方法到2D彈性各向同性介質(zhì)，進(jìn)一步靠近實(shí)際地下介質(zhì)模型.