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

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

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

0 引言

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

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

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

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

1 方法原理

1.1 頻率域波動方程求解

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

(1)

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

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

(2)

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

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

1.2 目標函數

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

(3)

(4)

對于第i炮，定義

(5)

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

(6)

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

i=1,2,…,Ns

(7)

1.3 目標函數梯度求取

記E1的分子和分母對某一模型參數m1的偏導數分別為ξ和η，則

(8)

(9)

若取實部，則有

(10)

所以目標函數對某一模型參數m1的導數為

(11)

為使公式簡約，記

(12)

則式(11)變為

(13)

將上式寫成矩陣形式

(14)

如果計算的是全部模型參數的偏導數，則可以得到目標函數對整個模型參數m的梯度

(15)

其中

(16)

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

2 模型試算

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

圖1 CSEG速度模型

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

圖2 平滑模型

2.1 不含噪聲時的合成數據測試

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

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

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

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

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

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

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

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

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

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

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

圖6 模型誤差收斂曲線

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

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

2.2 不同初始模型測試

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

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

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

2.3 含隨機噪聲時的合成數據測試

為了使合成數據測試更接近真實情況，在炮記錄中加入一定比例的隨機噪聲。不同目標函數對噪聲的抗噪性不同，L2型目標函數對隨機噪聲具有較強的穩定性。為了消除子波的影響而引入參考道后，新的規則化數據對噪聲的敏感性呈現出明顯差異。

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

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

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

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

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

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

3 結論

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