李海山 楊午陽 雍學善

(中國石油勘探開發(fā)研究院西北分院，甘肅蘭州 730020)

1 引言

隨著地震勘探的深入，逆時偏移等深度域成像技術得到了快速發(fā)展，但這些方法的成像質量很大程度上取決于速度模型的精度[1，2]。目前存在大量的速度建模方法，如疊加速度分析、偏移速度分析、旅行時層析反演和全波形反演方法等[3]。全波形反演(Full waveform inversion，FWI)充分利用了地震波的運動學和動力學特征，通過最優(yōu)化算法不斷地更新介質參數模型，使模擬地震記錄與實際觀測地震記錄之間的殘差達到最小，從而得到與實際介質參數最接近的模型[4]。相對于旅行時層析等方法，FWI精度更高，不僅能夠滿足復雜構造成像方法對速度參數的精度要求，而且具備巖性反演與儲層預測的潛力[3，5]。

自從Tarantola[6]和Mora[7]提出FWI方法以來，眾多學者致力于該類方法的理論及應用研究，已被廣泛應用于海上資料[4，8]，并逐漸應用于陸上資料[9，10]。FWI包括三個關鍵環(huán)節(jié)，即炮點正向傳播波場和檢波點逆時傳播殘差波場的計算、模型參數梯度(或Frechét導數)的計算和模型參數的更新，其中梯度計算是最關鍵的環(huán)節(jié)，梯度的計算方式及精度決定了反演的計算量及收斂性。由于直接求取Frechét導數往往較困難，需要巨大的計算量，因此許多學者采用間接方法避免Frechét導數的直接求取，如Tarantola[6]利用格林函數法推導出目標函數關于地下介質參數的梯度。伴隨狀態(tài)法是求解目標函數關于模型參數梯度的有效方法之一，已廣泛用于地球物理反演問題的求解[11，12]。胡光輝等[13]利用伴隨狀態(tài)法導出了頻率域梯度計算公式，實現了時間域正演與頻率域反演相結合的三維聲波FWI。

FWI可在不同域中實現，如時間域[6，7]、頻率域[14，15]和拉普拉斯域[16]等。盡管在頻率域可實現從低頻到高頻的多尺度FWI，但頻率域正演內存需求巨大，特別是三維情況下矩陣求逆運算對于內存需求極高，甚至難以實現[17]；與頻率域反演方法相比，拉普拉斯域反演方法不易陷入局部極值，對初始模型的依賴程度更低，但只能恢復模型中的長波長構造信息，對模型細節(jié)的刻畫能力不足[18]。由于在不同域中實現FWI具有不同的優(yōu)缺點，因此出現了將各自優(yōu)點相結合的混合域FWI方法，如時間—頻率域[14，19]、拉普拉斯—頻率域[18]、時間—拉普拉斯—頻率域[20]等FWI方法。

盡管時間域反演方法由于同時反演所有頻率成分，增大了問題的非線性程度，對初始模型的依賴性較強，但實際資料反演表明反射層析初始建模與時間域多尺度反演策略相結合可以顯著降低這種影響[21，22]； 此外，時間域波場正向傳播及逆時外推時較直接、快速，可以更好地估計震源的延遲時間，可以靈活地對數據進行必要的預處理； 同時，時間域正演可以方便地使用震源并行算法，能夠較好地適應各種采集系統。

綜合前述分析，本文基于三維一階速度—應力聲波方程，采用基于攝動理論的伴隨狀態(tài)法，推導出時間域三維一階速度—應力聲波方程的伴隨方程及相應的縱波速度的梯度計算公式； 在此基礎上實現了時間域三維聲波FWI，并通過模型數據檢驗了方法的可行性和有效性。

2 全波形反演目標函數

非均勻各向同性介質中，三維聲波方程的一階速度—應力方程為[23]

(1)

式中：p為壓力場；vi(i=x，y，z)為波場速度分量；ρ為介質密度；fp為震源項；K為體積模量，與縱波速度V的關系為

K=ρV2

(2)

全波形反演是利用全波場信息，通過優(yōu)化方法估計地下介質參數模型，使正演地震記錄與觀測地震記錄達到最佳擬合[5]。設目標函數為波場殘差能量，即

δ(x-xr)dΩdt

(3)

式中：u(m；x，t)為正演地震記錄；T為地震記錄長度；uobs(x，t)為觀測地震記錄；m為模型參數向量；u=(vx，vy，vz，p)T為地震波場向量；xr為接收點位置；Ω為模型的邊界。

利用FWI方法求取地下介質模型參數時，如果采用局部最優(yōu)化方法，例如共軛梯度法，需要計算目標函數對模型參數的梯度，且反演結果的有效性依賴于梯度的有效性和準確性。自Plessix[11]將源于控制論的伴隨狀態(tài)法引入反演問題，許多學者將該方法引入地球物理反演問題的求解，用于計算依賴于一組狀態(tài)變量的目標函數的梯度[24]。依據伴隨狀態(tài)法原理，對式(3)施加式(1)的約束，可以得到如下的最優(yōu)化問題

(4)

式中：λ(x，t)=(λx，λy，λz，λp)T為伴隨波場向量；E(m)=0為式(1)的簡寫形式。

3 伴隨方程及模型參數梯度

將式(4)展開，有

(5)

令式(5)對地震波場向量的偏導數為零，即

(6)

結合邊界條件

(7)

及終止條件

(8)

可以得到下面的伴隨方程(詳細推導見附錄A)

(9)

該方程即為波場殘差從檢波點位置開始的逆時反傳方程。

式(5)對體積模量K求偏導，可得

(10)

定義分辨率核函數為

(11)

式(11)即為體積模量的梯度計算公式。由式(2)可得關于縱波速度的梯度計算公式為

(12)

考慮到地震成像中速度模型為首要需求，本文暫不考慮密度變化的影響，即只進行三維常密度聲波方程的FWI。

4 全波形反演策略

目前多采用局部優(yōu)化反演方法使目標函數達到最小[25]。在局部最優(yōu)化方法中，雖然共軛梯度法只利用一階導數信息，卻克服了最速下降法收斂慢的缺點，避免了牛頓法需要存儲和計算Hessian矩陣并求逆的缺點，因此得到廣泛使用[3，26]。共軛梯度法的基本思想是利用目標函數上某一點的梯度構建出一組共軛方向，沿著共軛方向進行搜索，尋求目標函數的極小值點。共軛梯度法縱波速度更新的表達式離散化后寫成向量形式為

Vn+1=Vn+αφn

(13)

式中：n為迭代次數；α為共軛方向上的搜索步長；φn為共軛方向，即

(14)

式中βn為標量參數。共軛梯度法根據βn的不同求取形式分類，本文采用Fletcher-Reeves方法計算[26]

(15)

共軛方向上的搜索步長α對反演收斂有重要的影響，步長過小會增加迭代次數、降低反演效率；步長過大易導致反演陷于局部極值。通常采用線性搜索方法和拋物線擬合法確定模型參數更新步長，拋物線擬合法需要額外的二次正演模擬，因此本文采用非精確線性搜索方法求取步長[27]。

在每次迭代過程中，炮點正向傳播波場用式(1)計算，檢波點逆時傳播殘差波場用式(9)計算，用式(12)計算縱波速度梯度，接著用式(13)更新縱波速度模型。由于交錯網格有限差分法比中心網格有限差分法具有更好的數值頻散抑制能力[24]，因此波場傳播采用交錯網格有限差分法模擬；此外，在時間域波場傳播較直接、快速，結合GPU并行計算的優(yōu)點[28]，可以顯著提高計算速度。

5 模型測試

3D SEG/EAGE逆掩推覆速度模型縱向層位多且具有強烈的橫向速度變化，可以用于驗證FWI的適應性，檢驗FWI對斷層、河道等細節(jié)的刻畫能力?？紤]到計算量和單個GPU節(jié)點的內存限制，在測試時對3D SEG/EAGE逆掩推覆速度模型相關參數進行了修改(圖1)，其中網格數為160×160×120，空間步長為15m。

采用正交觀測系統，利用速度模型(圖1)生成三維聲波觀測記錄，震源子波是主頻為15Hz的雷克子波，炮點與接收點均置于地表，共激發(fā)256炮進行全方位接收，每炮160×160個檢波點，炮點距和炮線距均為150m，檢波點距和檢波線距均為15m，記錄長度為1.5s，時間采樣間隔為1ms，最大炮檢距為3400m。反演初始速度模型通過三維高斯平滑得到，共迭代60次得到最終反演結果。

圖1 3D SEG/EAGE速度模型

圖2為225m深度反演速度切片，可見反演后構造邊界， 特別是淺層河道邊界及內部速度結構得到了準確刻畫。圖3為375m深度反演速度切片，盡管初始模型嚴重偏離準確模型，但反演后斷層及局部細節(jié)都得到了準確刻畫。需要注意的是，在圖2和圖3中，邊界區(qū)域的速度反演精度遠低于中心區(qū)域，主要是受覆蓋次數低及炮檢距較小的影響。

圖4為第80條聯絡測線準確速度、初始速度及反演速度剖面對比，可見淺層河道、構造等得到了較好刻畫，而隨著深度的增加，反演準確性逐漸降低，部分原因是深層反射能量不足，另一原因是需要采用更優(yōu)的梯度預條件算子及梯度更新優(yōu)化方法提高反演精度[29]。同時可見，逆掩推覆體右側陡傾構造的反演精度高于左側， 這是因為逆掩推覆體右側的觀測炮檢距更大，為陡傾構造的反演提供了更好的照明。因此，FWI要求盡可能提供寬方位大炮檢距資料。

圖2 真實速度(a)、初始速度(b)和反演結果(c)225m深度切片對比

圖4 真實速度(a)、初始速度(b)和反演結果(c)第80條聯絡測線對比

圖5為(0.9km，0.9km)、(1.2km，1.2km)及(1.5km，1.5km)三個不同位置準確速度、初始速度及反演速度單道曲線對比，由于(0.9km，0.9km)處于逆掩推覆復雜構造上，反演速度曲線與準確速度曲線偏差較大(圖5a)；(1.2km,1.2km)、(1.5km,1.5km)位置處于逆掩推覆復雜構造翼部，反演速度的準確性明顯提高(圖5b、圖5c)。下一步，可以通過實現基于模型分割的GPU/CPU異構協同并行反演方法突破計算和內存的限制，增大模型和炮檢距提高復雜高陡構造區(qū)的速度反演精度；同時，為提高反演曲線(圖5)低頻趨勢的準確性，可采用時間域多尺度反演策略。

圖6為FWI時歸一化波場殘差隨迭代次數變化曲線，可見本文反演方法穩(wěn)定性和收斂性均較好，迭代60次后觀測記錄與擬合記錄之間的殘差較小。需要注意的是，本文反演中采用準確子波，反演方法的噪聲抑制能力也有待提高。為進一步提高反演方法的性能，可采用褶積型目標函數[30]，并引入模型參數的正則化項[31]。

圖5 不同位置真實速度(藍線)、初始速度(黑線)和反演結果(紅線)單道曲線對比

圖6 歸一化波場殘差隨迭代次數的變化曲線

本文方法基于GPU在時間域實現，相比單節(jié)點CPU，波場正向傳播及逆時外推時較直接、快速。表1給出了模型網格數不同時單個CPU與GPU的正演計算時間對比，表2給出了單節(jié)點上GPU卡的數量與正演計算時間對比(GPU卡為Tesla K10)，可見，與CPU相比，GPU具有較高的計算效率，尤其是對于較大模型的三維FWI。

表1 單個CPU與GPU的正演計算時間(s)對比

表2 不同GPU卡數量的正演計算時間(s)對比

6 結論

本文采用基于攝動理論的伴隨狀態(tài)法，實現了基于一階速度—應力聲波方程的三維全波形反演方法。具體計算時，炮點正向傳播波場用時間域三維一階速度—應力聲波方程計算，檢波點逆時傳播殘差波場用推導出的伴隨方程計算，波場傳播采用交錯網格有限差分法模擬。模型測試結果證明了方法的可行性和有效性。

本文方法由于在時間域實現，波場正向傳播及逆時外推時較直接、快速，結合GPU并行計算的優(yōu)點，可以顯著提高計算速度。本文方法只是實現了在單個節(jié)點上多個GPU卡的并行計算，在計算效率、反演準確性的提高、對噪聲的抑制能力、對實際資料的適應性等方面有待進一步改進。

附錄A 伴隨方程的推導

式(5)對波場速度分量vx求導，可得

(A-1)

對上式分部積分，有

(A-2)

式(A-2)滿足終止條件及邊界條件：λx|t=T=0，λp|x=x1=λp|x=x2=0，其中x1和x2分別為地下介質參數模型在x方向的兩個邊界。結合終止條件及邊界條件，式(A-2)可以化簡為

(A-3)

同理，式(5)對速度分量vy和vz求導，結合終止條件及邊界條件，可得

(A-4)

(A-5)

式(5)對壓力場p求偏導，可得

(A-6)

對式(A-6)分部積分， 結合終止條件及邊界條件：λp|t=T=0，λx|x=x1=λx|x=x2=0，λy|y=y1=λy|y=y2=0,λz|z=0=λz|z=Z=0, 其中y1、y2和Z分別為地下介質參數模型在y方向的兩個邊界和在z方向的下邊界，得

(A-7)

令式(A-3)、式(A-4)、式(A-5)和式(A-7)等于0，可得伴隨方程(式(9))。