基于數(shù)據(jù)相似性的不依賴子波的頻率域全波形反演

辛天亮 黃建平* 解 飛② 周 濱 盧子卓

(①中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580； ②中國石油化工股份有限公司華東分公司,江蘇南京 210011； ③中海石油(中國)有限公司天津分公司 天津 300452)

0 引言

全波形反演[1-4](FWI)依據(jù)地震數(shù)據(jù)的動力學(xué)特征提取速度參數(shù)，與旅行時層析等基于地震數(shù)據(jù)的運(yùn)動學(xué)特征的方法相比，可以獲得更精確的地下速度結(jié)構(gòu)參數(shù)。

盡管FWI具有恢復(fù)地下精細(xì)速度信息的能力，但是仍然存在許多限制，如易陷入局部極值、缺失低頻信息、對初始模型要求高等，反演實際數(shù)據(jù)時缺少準(zhǔn)確的震源子波信息是主要制約因素之一[5]。早期FWI在模型試算時通常假設(shè)震源子波是已知的——已知震源特性反演法(KSS)。在實際資料處理中最直接的方法是從直達(dá)波中提取子波特征，但是由于地球介質(zhì)中傳播路徑的影響，直達(dá)波并不是真實的震源子波。Song等[6]、Pratt[7]先后在頻率域提出了包含子波迭代估算的反演方法(IES)，假設(shè)初始速度模型以及子波都接近真實情況，在反演過程中同時迭代、更新速度模型和震源子波。如果反演的模型嚴(yán)重偏離真實模型，那么在迭代中將會得到不正確的子波，兩者相互影響，導(dǎo)致反演無法收斂到真實解。

為了避免對子波的估計，Zhou等[8]提出了頻率域歸一化振幅反演法(ATN)，其基本原理是在每個頻率中用單炮記錄平均振幅對炮集的每個接收點振幅進(jìn)行歸一化，以消除震源特性信息。Lee等[9]提出了標(biāo)準(zhǔn)道歸一化的FWI(STN)。與ATN相比，STN使用單一參考道，并且使用頻率域波場數(shù)據(jù)。Choi等[10]提出了兩種與子波無關(guān)的頻率域振幅型目標(biāo)函數(shù)：一種具反褶積作用，與ATN的目標(biāo)函數(shù)形式基本一致，但是使用單一的參考道振幅數(shù)據(jù)，利用傳統(tǒng)的直接計算雅可比矩陣的方式計算梯度；另一種具褶積作用，將模擬數(shù)據(jù)和觀測數(shù)據(jù)分別與參考道的觀測數(shù)據(jù)和模擬數(shù)據(jù)相乘，然后使用新得到的數(shù)據(jù)構(gòu)建基于L2范數(shù)的目標(biāo)函數(shù)。由于數(shù)據(jù)同時包含理論子波和實際子波的影響，在相減過程中自然消去了子波的影響，從而達(dá)到不依賴于子波的效果。Greenhalgh等[11]基于聲波方程僅利用數(shù)據(jù)的振幅譜對淺層地震模型數(shù)據(jù)進(jìn)行KSS、IES和ATN，認(rèn)為KSS反演效果最好。Xu等[12]利用井間模型的頻率域彈性波場數(shù)據(jù)進(jìn)行KSS、IES、ATN和STN，認(rèn)為無論對于多分量還是單分量數(shù)據(jù)，KSS和IES能得到較好的反演效果，ATN反演效果次之，STN對參考道的依賴性很大，反演效果最差。Choi等[13]利用彈性波對數(shù)波場數(shù)據(jù)測試了子波估計型目標(biāo)函數(shù)和頻率域不依賴子波型目標(biāo)函數(shù)對噪聲的敏感性，對于含隨機(jī)噪聲數(shù)據(jù)，子波估計型目標(biāo)函數(shù)效果更好，對于含相干噪聲數(shù)據(jù)則相反。Choi等[14]將頻率域基于褶積的不依賴于子波的FWI用于時間域聲波方程，認(rèn)為使用平均參考道的效果好于使用單一參考道，且在已知子波情況下的反演效果最好。敖瑞德等[15]結(jié)合時間域不依賴子波的FWI與包絡(luò)反演，通過對比包絡(luò)、包絡(luò)對數(shù)和包絡(luò)平方三種目標(biāo)函數(shù)反演結(jié)果，發(fā)現(xiàn)包絡(luò)對數(shù)目標(biāo)函數(shù)的深層反演效果最好。Zhang等[16]針對時間域FWI的參考道褶積過程中由褶積和相關(guān)運(yùn)算產(chǎn)生的噪聲問題，提出對參考道加時窗的方式加快收斂以改善反演效果。楊濤等[17]在彈性波混合域?qū)崿F(xiàn)了不依賴子波的FWI。

本文利用頻率域單頻數(shù)據(jù)特點，在復(fù)數(shù)域設(shè)計了一種類似于實數(shù)域表征數(shù)據(jù)相似性的目標(biāo)函數(shù)(SOD)。SOD與相干性衡量目標(biāo)函數(shù)的結(jié)構(gòu)相似[18]，可在一定程度上降低反演的非線性。基于互相關(guān)的新方法可以弱化模擬數(shù)據(jù)與觀測數(shù)據(jù)間的振幅匹配，更適合于實際應(yīng)用。在構(gòu)造目標(biāo)函數(shù)時，自然消除了震源子波的影響，而且不需要使用參考道，從而避免了由于STN參考道選擇不當(dāng)導(dǎo)致收斂緩慢或無法收斂的缺陷。ATN可以使用頻率域波場數(shù)據(jù)或只使用振幅數(shù)據(jù)。若使用頻率域波場數(shù)據(jù)，由于求和所用的參考道數(shù)據(jù)為復(fù)數(shù)型，求和后不能校正數(shù)據(jù)能量，甚至求和后結(jié)果的振幅接近零時會導(dǎo)致反演不穩(wěn)定；若只使用振幅數(shù)據(jù)雖然可以避免上述缺陷，但是缺少相應(yīng)的相位信息，也難以達(dá)到頻率域波場數(shù)據(jù)的反演精度[12]。文中的SOD反演法(下文簡稱SOD)使用頻率域波場數(shù)據(jù)，并且不需要參考道求和，因此避免了ATN的缺陷。通過模型試算，對比了STN、ATN和SOD的反演結(jié)果，其中STN和ATN使用振幅型數(shù)據(jù)，SOD使用頻率域波場數(shù)據(jù)。首先對STN、ATN和SOD三種方法做無噪數(shù)據(jù)測試，用不同的子波生成模擬記錄，以驗證方法的有效性； 然后用一個梯度模型測試SOD的穩(wěn)定性[19-20]； 最后在單炮記錄中加入不同信噪比的隨機(jī)噪聲，測試SOD對含噪數(shù)據(jù)的適應(yīng)性。

1 方法原理

1.1 頻率域波動方程求解

在頻率域，聲波波動方程為

(1)

式中：K(x)為體積模量，x為笛卡爾坐標(biāo)下的介質(zhì)空間位置；ρ(x)為密度；u(x,ω)和f(x,ω)分別為頻率域波場值和震源值，ω為角頻率。本文將式(1)用有限差分格式離散[21-22]，并加入合適的邊界條件，得到矩陣表達(dá)形式

S(ω)U(ω)=F(ω)

(2)

式中：S(ω)為與模型參數(shù)、頻率、離散格式以及邊界條件有關(guān)的阻抗矩陣，設(shè)模型中共有n個節(jié)點，則S為n階的大型稀疏矩陣，其非零元素對稱分布在主對角線及兩邊，但是由于兩邊的值大小不同，所以它并非是對稱矩陣；U(ω)為所有模型節(jié)點的頻域波場值按照一定順序排列得到的n維列向量；F(ω)為頻率域的點源或組合形式的震源構(gòu)成的n維列向量。

針對式(2)的大型稀疏方程組的求解，通常可以使用分解法和迭代法。使用分解法直接求解(如LU分解)，只需要對S(ω)一次分解，就可以利用分解結(jié)果求解多個震源的波場，體現(xiàn)了頻率域正演的高效性。

1.2 目標(biāo)函數(shù)

頻率域中，觀測記錄和模擬記錄的共炮點道集用格林函數(shù)可分別表示為

(3)

(4)

對于第i炮，定義

(5)

式中：ui和di分別為第i炮的檢波點處的模擬值和觀測值組成的列向量； 上角標(biāo)*表示復(fù)數(shù)共軛。將波場值的格林函數(shù)代入式(5)，得到

(6)

由式(6)可見，E1中并不包含震源特征信息，具有不依賴子波的性能。因此，本文構(gòu)建了不依賴子波的目標(biāo)函數(shù)

i=1,2,…,Ns

(7)

1.3 目標(biāo)函數(shù)梯度求取

記E1的分子和分母對某一模型參數(shù)m1的偏導(dǎo)數(shù)分別為ξ和η，則

(8)

(9)

若取實部，則有

(10)

所以目標(biāo)函數(shù)對某一模型參數(shù)m1的導(dǎo)數(shù)為

(11)

為使公式簡約，記

(12)

則式(11)變?yōu)?/p>

(13)

將上式寫成矩陣形式

(14)

如果計算的是全部模型參數(shù)的偏導(dǎo)數(shù)，則可以得到目標(biāo)函數(shù)對整個模型參數(shù)m的梯度

(15)

其中

(16)

由式(16)可見，SOD對應(yīng)的梯度與常規(guī)L2型目標(biāo)函數(shù)的梯度形式相似，都非常簡潔，因此并不會增加額外的計算量。

2 模型試算

本文采用CSEG速度模型(圖1)驗證基于數(shù)據(jù)相似性的不依賴子波的頻率域FWI的有效性，并與STN、ATN的反演效果對比。本文使用的優(yōu)化方法為共軛梯度法，并取一個較小值作為固定步長與歸一化后的共軛梯度相乘更新模型。

圖1 CSEG速度模型

水平網(wǎng)格數(shù)為332,垂直網(wǎng)格數(shù)為200，水平與垂直網(wǎng)格尺寸均為15m，模型尺寸為4.98km×3.00km。震源位于z=15m處，起始炮位于x=60m處，炮間距為30m，共110炮。檢波器設(shè)置在地表進(jìn)行全接收，時間采樣間隔為1ms，采樣點數(shù)為2000，總記錄長度為2s

圖2 平滑模型

2.1 不含噪聲時的合成數(shù)據(jù)測試

為了驗證SOD的有效性，本文用常密度聲波波動方程測試各反演方法的效果。采用頻率域有限差分正演模擬，差分格式采用二階最優(yōu)9點格式[25]。初始速度模型為平滑模型(圖2)，用于生成觀測記錄的準(zhǔn)確子波為雷克子波(圖3a)，主頻為15Hz，最大振幅為1.0。為了測試本文方法對不同子波的處理能力，選取主頻均為15Hz、不同振幅的子波(圖3b～圖3f，其中圖3b、圖3c、圖3e、圖3f的最大振幅為1.0，圖3d的最大振幅為2.0)作為生成模擬記錄的震源子波，反演的頻率范圍為1～20Hz，頻率間隔為1 Hz，共20個頻率，每個頻率迭代20次，低頻的反演結(jié)果作為相鄰高頻的輸入。同時，為了提高計算效率，采用振幅編碼的方式將110炮編碼成40個超級炮[26-27]，為避免由于隨機(jī)編碼對反演結(jié)果的影響，測試中統(tǒng)一使用同一套編碼矩陣。

圖4為選取復(fù)雜子波(圖3f)的不同方法的反演結(jié)果。由圖可見，不同方法在復(fù)雜子波的情況下仍能較好地反演速度模型，特別是淺層反演效果較好，其中STN的深層更新效果較差(圖4a)，ATN有所改進(jìn)(圖4b)，SOD的反演效果最好(圖4c)。上述結(jié)果表明，本文方法在子波主頻相同的情況下反演結(jié)果不依賴子波。圖5為由圖4 得到的1.5km、3.0km處單道速度曲線。由圖可見，與ATN和SOD相比，STN的擬合較差，3.0km處單道速度曲線(圖5b)在500～1500m深度范圍較平滑，因此STN更新不足； ATN則有較大改善； SOD曲線最接近真實速度曲線。圖6為模型誤差收斂曲線。由圖可見，SOD收斂更好。綜上所述，本文提出的SOD的反演結(jié)果不依賴子波，而且與另兩種不依賴子波的反演方法(ATN、STN)相比，反演精度更高。

圖3 主頻為15 Hz、不同振幅的子波

圖4 選取復(fù)雜子波(圖3f)的不同方法的反演結(jié)果

值得注意的是，選用STN反演時，參考道的選擇是一個非常重要的問題。理論上，子波具有時變和空變特征，炮檢距越小，能量越大，畸變越小，因此小炮檢距參考道優(yōu)于大炮檢距參考道。但是炮記錄經(jīng)過編碼后，參考道選取原則要相應(yīng)改變。因此，為了減小由于參考道的選擇而導(dǎo)致各個超級炮產(chǎn)生的梯度間能量差異，本文設(shè)計了四種參考道選取方式。

(1)鄰近參考道。以(12×i)m所在道作為參考道，i=1,…,40，即在計算每一個超級炮的梯度時依次從左向右以120m為間隔選取參考道。

(2)遠(yuǎn)炮檢距參考道。當(dāng)i=1,…,20時，以(120×i+166)m所在道作為參考道，即第1炮選取286m所在道為參考道，依次從左向右以120m為間隔選取參考道； 當(dāng)i=21,…,40時，以(120×i-166)m所在道作為參考道，即第21炮選取2354m所在道為參考道，依次從左向右以120m為間隔選取參考道。

(3)最小振幅參考道。編碼后觀測記錄中最小振幅所在道。

(4)最大振幅參考道。編碼后觀測記錄中最大振幅所在道。

圖7為四種參考道選取方式對應(yīng)的反演結(jié)果。由圖可見： 方式(1)～方式(3)的反演結(jié)果只能得到一個大致的背景模型，無法反演精細(xì)構(gòu)造；方式(4)反演結(jié)果與真實模型較吻合。圖8為四種參考道選取方式反演結(jié)果模型誤差收斂曲線。由圖可見：①不同的參考道選取方式對反演結(jié)果影響很大，體現(xiàn)了STN的不穩(wěn)定性，雖然本文的參考道選取方式(方式(4))在一定程度上保證了STN的收斂性，但是基于振幅數(shù)據(jù)的ATN通過對每道求和后取平均的方式，可以自然地將歸一化后的數(shù)據(jù)能量校正到同一個水平，利于反演結(jié)果穩(wěn)定、高效地收斂；②SOD不需要選取參考道，并且同樣具備ATN的優(yōu)點，同時由于使用了包含更多信息的頻率域波場數(shù)據(jù)，收斂性更好；③方式(4)的誤差曲線收斂，方式(1)～方式(3)的誤差曲線不收斂。因此在后文STN測試中使用方式(4)的參考道選取方式。

圖5 由圖4得到的1.5km(a)、3.0km(b)處單道速度曲線

圖6 模型誤差收斂曲線

圖7 四種參考道選取方式的反演結(jié)果

圖8 四種參考道選取方式反演結(jié)果模型誤差收斂曲線

2.2 不同初始模型測試

為了測試不同反演方法對初始模型的依賴性，選取一個初始速度模型(圖9a)，使用的頻率與圖3相同。圖9為初始速度模型及不同方法迭代400次后的反演結(jié)果。由圖可見： STN反演效果很差，兩邊和中間都出現(xiàn)了虛假的低速構(gòu)造(圖9b)； 與STN相比，ATN在淺層速度較準(zhǔn)確，但在深部更新不明顯(圖9c)； SOD反演結(jié)果基本正確，而且明顯比ATN更新效果更好，在深部尤其如此(圖9d)。圖10為由圖9抽取的1.5km、3.0km處單道速度曲線。由圖可見： 在深度小于500m時，三種方法的反演曲線基本與真實速度曲線一致；STN和ATN的中部速度變化異常，深部速度更新不足； SOD速度曲線在中、深部均很接近真實速度曲線，因此SOD對初始模型的依賴性較小，可更穩(wěn)定、高效地反演模型參數(shù)。

圖9 初始速度模型及不同方法迭代400次后的反演結(jié)果

圖10 由圖9抽取的1.5km(a)、3.0km(b)處單道速度曲線

2.3 含隨機(jī)噪聲時的合成數(shù)據(jù)測試

為了使合成數(shù)據(jù)測試更接近真實情況，在炮記錄中加入一定比例的隨機(jī)噪聲。不同目標(biāo)函數(shù)對噪聲的抗噪性不同，L2型目標(biāo)函數(shù)對隨機(jī)噪聲具有較強(qiáng)的穩(wěn)定性。為了消除子波的影響而引入?yún)⒖嫉篮螅碌囊?guī)則化數(shù)據(jù)對噪聲的敏感性呈現(xiàn)出明顯差異。

本文在單炮記錄中加入了不同信噪比的隨機(jī)噪聲，用以測試噪聲對三種目標(biāo)函數(shù)反演結(jié)果的影響，反演所使用的頻率與圖3相同，初始模型與圖2相同。圖11為在單炮記錄中加入信噪比分別為10、20、30dB隨機(jī)噪聲的不同方法反演結(jié)果。由圖可見：STN受隨機(jī)噪聲影響嚴(yán)重，深部反演效果很差(圖11a)；ATN受噪聲的影響也較明顯(圖11b)；SOD受噪聲的影響較小(圖11c)。圖12和圖13分別為由圖11抽取的1.5km、3.0km處單道速度曲線。由圖可見，由于受噪聲影響，不同方法的反演結(jié)果在淺層都存在明顯波動，表現(xiàn)為： ①STN速度曲線非常不穩(wěn)定(圖12a、圖13a)；②ATN速度曲線在深層也存在明顯跳動(圖12b、圖13b)；③SOD速度曲線對噪聲的敏感程度最低，加入不同隨機(jī)噪聲的速度曲線與無噪數(shù)據(jù)基本一致，表明SOD對噪聲具有更好的穩(wěn)定性(圖12c、圖13c)。圖14為加入信噪比為30dB隨機(jī)噪聲的模型誤差收斂曲線。由圖可見，STN的收斂性明顯變差，ATN次之，SOD收斂性最好，與Choi等[14]的認(rèn)識一致。

圖11 在單炮記錄中加入信噪比分別為10dB(上)、20dB(中)、30dB(下)隨機(jī)噪聲的不同方法反演結(jié)果

圖12 由圖11抽取的1.5km處單道速度曲線

圖13 由圖11抽取的3.0km處單道速度曲線

圖14 加入信噪比為30dB隨機(jī)噪聲的模型誤差收斂曲線

Xu等[12]分析了STN和ATN的抗噪性，認(rèn)為STN目標(biāo)函數(shù)的分子和分母同時包含噪聲影響項，ATN目標(biāo)函數(shù)由于對參考道求和從而消除了分母的噪聲影響項，因此比STN的抗噪性好。本文提出的SOD目標(biāo)函數(shù)類似于ATN，通過相乘、求和運(yùn)算使噪聲在分子和分母中的比重相當(dāng)，因此目標(biāo)函數(shù)值較穩(wěn)定。穩(wěn)定的目標(biāo)函數(shù)值可以產(chǎn)生穩(wěn)定的梯度，這在反演過程中至關(guān)重要。若局部炮記錄產(chǎn)生的梯度能量發(fā)生異常，其他炮記錄產(chǎn)生的梯度作用相對減弱，反演的收斂速度將變慢，尤其是異常梯度指向局部極小值時，可能導(dǎo)致反演失敗。

3 結(jié)論

準(zhǔn)確獲取震源子波一直是地震資料處理的難題，本文提出的SOD在構(gòu)造目標(biāo)函數(shù)的過程中消除了震源子波對反演結(jié)果的影響。模型測試發(fā)現(xiàn)： ①SOD的反演結(jié)果不依賴子波，而且與另兩種不依賴子波的反演方法(ATN、STN)相比，反演精度更高，對實際資料的處理具有非常重要的意義； ②SOD的原理是基于數(shù)據(jù)間的相似性，對數(shù)據(jù)間能量差異具有更大的包容性，所以對初始速度模型的依賴性不高； ③即使存在隨機(jī)噪聲，通過相乘、求和運(yùn)算作用在分子和分母中的噪聲比重相當(dāng)，因此SOD目標(biāo)函數(shù)值較穩(wěn)定； ④SOD還可以與炮編碼策略結(jié)合，有效提高計算效率。