三維聲波全波形反演的實現(xiàn)與驗證

胡光輝，賈春梅，夏洪瑞，賀劍波，宋 林，沈忠秋

（1.中國石油化工股份有限公司石油物探技術研究院，江蘇南京211103；2.中國石油大學（華東），山東青島266580）

隨著石油工業(yè)的發(fā)展和勘探開發(fā)的不斷深入，石油勘探的難度越來越大。這就需要不斷尋找新的技術方法以適應當今石油勘探的需求。全波形反演（FWI）就是順應這種要求而迅速發(fā)展起來的一項新技術，并已成為當今地球物理界的研究熱點。20世紀80年代，Lailly［1］和 Tarantola［2］分別提出了時間域的全波形反演理論與方法，Gauthier等［3］和 Mora［4］實現(xiàn)了二維的實際資料處理。20世紀90年代初，Pratt等［5］將全波形反演推廣到頻率域，使得全波形反演在三維實際資料處理中具有了可行性。受計算機技術的限制，直到近年全波形反演技術才應用于三維實際資料處理中［6－9］，但也只局限于海上資料的處理，陸上實際資料的應用還存在很大的挑戰(zhàn)。

全波形反演是一個基于地震全波場模擬的數(shù)據(jù)擬合過程，幾乎使用了地震記錄中所有有效信息［1］，而不像其它傳統(tǒng)的方法（如旅行時層析成像等技術）僅使用地震記錄中的部分信息。全波形反演是一個局部尋優(yōu)的過程，通過更新迭代初始模型進而減小計算數(shù)據(jù)和觀測數(shù)據(jù)之間的誤差，逐步逼近真實模型。全波形反演是一種建立高精度速度模型的有效手段［10］，尤其是在地下地質(zhì)構造復雜、地震成像難的區(qū)域，全波形反演能夠充分利用地震剖面上所記錄的全部波形及相位信息。實際上，在全波形反演的局部尋優(yōu)過程中，關鍵技術在于高效、精確地獲取梯度信息［11－12］，它可以提供精度為半波長的高精度物理參數(shù)模型［13－14］。但是全波形反演是一個“病態(tài)”問題，因此要求初始模型足夠接近真實地質(zhì)模型，以避免在迭代過程中陷入局部最小。此外，全波形反演對計算資源的要求也是制約其發(fā)展的瓶頸之一，這也是三維彈性波反演至今受到限制的原因。近年來，隨著計算機技術的快速發(fā)展和一系列優(yōu)化算法的出現(xiàn)，三維聲波全波形反演的應用開始陸續(xù)出現(xiàn)［6－9］。

我們采用頻率域反演聯(lián)合時間域正演的混合優(yōu)化算法實現(xiàn)三維聲波的全波形反演。在全波形反演方法原理與算法介紹的基礎上，具體論述了有效計算梯度的伴隨狀態(tài)法和適應計算資源的并行算法，并使用SEG／EAGE目標區(qū)塊模型對這一算法進行驗證。

1 全波形反演方法

全波形反演方法的原理是：給定一個初始模型，通過正演模擬得到其傳播波場，將結果與實際采集資料進行對比，若兩者的誤差沒有滿足給定的誤差范圍，則修正模型，重復上述操作，直到誤差達到設定要求。全波形反演是在實際資料約束下逐步尋優(yōu)的過程。由于時間域的數(shù)據(jù)量巨大，需要巨大的存儲空間，而當前的計算機集群往往不能滿足要求，尤其是在三維情況下，這一矛盾顯得尤為突出，因此頻率域反演算法應時而生。

頻率域反演使用幾個離散頻率的傅里葉級數(shù)取代龐大的時間采樣點信息，大大節(jié)省了計算機存儲資源。其次，頻率域迭代算法從低頻開始反演到高頻，有效利用低頻信息降低全波形反演的非線性，較時間域反演而言，頻率域迭代算法更利于這一算法的實現(xiàn)。尤其是對衰減系數(shù)等物理參數(shù)的反演，頻率域反演有時間域反演不可比擬的優(yōu)越性。因此頻率域全波形反演日益受到全球地球物理學家的重視。頻率域全波形算法提高了全波形反演方法處理實際問題的能力，尤其在處理三維大數(shù)據(jù)量時更顯示了其不可取代的地位。因此，頻率域全波形反演得到了迅猛的發(fā)展，在正演波場模擬、反演頻率選擇策略、目標函數(shù)設置、梯度計算等方面都相繼取得了重大突破。在上述算法中，正演波場模擬是反演的基礎，正演的精度和效率在很大程度上決定了反演的精度和效率。

1.1 正演方法的選擇

正演是反演的重要組成部分，正演方法的選擇決定著反演的精度和效率。因此，能否正確選擇正演方法決定了全波形反演的可行性，尤其是三維多炮點勘探情況下。在二維情況下，頻率域正演有不可替代的優(yōu)越性。

對于某一頻率，波場u和震源s是線性關系：

式中：A（x，ω）為阻抗矩陣；u（x，ω）為某一頻率波場；s（x，ω）為某一頻率震源函數(shù)。

矩陣A是高度稀疏、非對稱、非正定的，直接求逆非常困難。因此，常采用稀疏矩陣因式分解直接求解法。

利用直接求解法，得到

式中：LU是阻抗矩陣A的分解。從（2）式可以看出，LU完全獨立于震源s，即一旦LU建立，對不同的震源是一個線性求解過程，LU阻抗矩陣分解結果可以對不同震源重復利用，從而大大降低了計算量。但每次正演只能計算一個頻率。在三維情況下則完全不同，Operto等［15］和 Virieux等［16］提出在三維情況下LU的建立對于內(nèi)存存儲空間需求極高，甚至難以實現(xiàn)。Operto等［15］指出，在使用SEG／EAGE三維推覆體模型進行波場模擬（空間步長為50m，正演頻率為10Hz）時，需要320GB內(nèi)存存儲單精度的LU分解結果。一般的計算機集群不能達到這一要求。為了避免這一問題，我們在處理三維情況下的全波形反演時，采取時間域正演方法有

式中：u（x，t）為波場值；B（x，t）為阻抗矩陣；s（x，t）為震源。

同時靈活運用離散傅里葉變換，將時間域波場轉換到頻率域來實現(xiàn)頻率域的全波形反演。這樣我們既得到了頻率域的波場值，又可以在時間域對數(shù)據(jù)進行所需要的預處理。一次正演可以得到多個頻率的波場值。表1給出了時間域與頻率域正演的優(yōu)越性對比結果。

表1 時間域正演與頻率域正演優(yōu)越性對比

從表1中不難看出，時間域正演比頻率域正演節(jié)省內(nèi)存，但如果需要存儲整個波場值，時間域正演的內(nèi)存要求為O（N5）。而如果計算資源允許，時間域正演使用震源并行算法，計算時間可降低到O（N4），頻率域正演在多節(jié)點并行計算方面則提速有限。此外，時間域正演可適用于任何采集系統(tǒng)，而直接求解法則受限，比如拖攬采集，一旦檢波器位置改變則需要花費大量時間重新構建LU矩陣。

1.2 反演方法選擇

與正演方法相同，反演也可以在時間域和頻率域實現(xiàn)。時間域反演，可以對資料進行所需要的預處理，或者針對某種波進行處理，加快收斂速度。但全波形反演是一個高度非線性問題，時間域反演同時反演所有頻率，不利于迭代算法的實現(xiàn)，而且時間域反演數(shù)據(jù)量龐大，尤其對三維大工區(qū)來說數(shù)據(jù)存儲量和計算量都尤為巨大甚至不可實現(xiàn)。另外，時間域資料包含所有頻率信息，不利于多尺度反演策略的實現(xiàn)。

如果使用所有頻率，頻率域反演完全等價于時間域反演，但是對于數(shù)據(jù)的預處理存在局限性。頻率域反演的一個優(yōu)勢在于：該方法利用一個復數(shù)（傅里葉級數(shù)）代替整個時間域各采樣點，同時使用少數(shù)的幾個離散頻率即可得到高分辨率速度模型，從而大大減小了存儲空間和計算量。此外，頻率域反演可以很容易地從低頻向高頻反演實現(xiàn)多尺度反演策略。從盡可能低的頻率開始反演，利用低頻數(shù)據(jù)對周期調(diào)頻不敏感的反演特性使問題更加線性化。我們使用頻率域反演聯(lián)合時間域正演的混合算法，實現(xiàn)三維全波形反演。

1.3 梯度的有效計算

梯度計算是全波形反演技術最核心的部分。如何高效率地計算梯度，決定著反演的計算量及程序開發(fā)的復雜程度。為了避免支付Frechet矩陣導數(shù)的計算所需的巨大代價，我們使用伴隨狀態(tài)法來計算梯度［11－12］。這種方法使用炮點正向傳播波場與接收點逆時傳播殘差波場的互相關生成梯度，避免了Frechet導數(shù)的直接計算。此外，偽保守形式的一階速度－壓力場聲波方程的使用，使梯度計算更加簡單有效。由于其算子的對稱性，使得炮點正向傳播波場和接收點逆時傳播殘差波場可以使用同一算子，大大提高梯度計算的效率，簡化了程序開發(fā)的復雜度。梯度的散射內(nèi)核獨立于正演算子且為對角矩陣，這一性質(zhì)使得梯度計算適用于不同的正演算子且大大降低了計算的復雜程度。

定義誤差函數(shù)為

式中：u表示模擬波場值；R表示算子；dobs為觀測記錄值。

正演算子的偽保守形式一階速度－壓力場聲波方程可寫為

（5）式可以簡寫成

其中，

需要指出的是，這里B′是對稱的，因為所有的物理參數(shù)都在方程的左側。這一改進將大大簡化梯度的計算。

使用伴隨狀態(tài)法求取梯度的公式為

式中：mi和λ＊分別表示模型參數(shù)和回傳波場。顯而易見，在計算炮點正向傳播波場與接收點逆時傳播殘差波場時，這種方法可以使用同一個正演算子，由于正演算子B′的對稱性，在求取對角散射內(nèi)核?Λ／?mi時完全獨立于正演算子，且Λ是對角矩陣，這一性質(zhì)使我們在計算梯度時可以使用不同的正演算子，且正演波場模擬和殘差回傳波場模擬時可以使用相同的正演模擬手段，使得梯度的計算高效精確。比較經(jīng)典的做法是在同一個反演算法里，靈活應用有限差分和有限元差分正演算子［17］，以適應不同要求。

從數(shù)學角度可以分析伴隨狀態(tài)法計算梯度的有效性。下面將驗證這種方法的準確性和精確度。使用有限差分法計算梯度的公式為

我們采用41m×41m×41m均質(zhì)立方體模型，模型中心加一半徑為50m的繞射體。x，y和z方向具有相同的步長（50m），震源采用5Hz主頻的雷克子波。圖1為采用伴隨狀態(tài)法計算得到的梯度水平方向切片（z＝800m）。該梯度使用頻率域反演聯(lián)合時間域正演得到。正演使用空間4階、時間2階有限差分交錯網(wǎng)格法。梯度計算僅使用5Hz單頻。圖2為采用有限差分法計算得到的梯度水平方向切片（z＝800m）。圖3為二者的一維對數(shù)的直接對比結果。

圖1 采用伴隨狀態(tài)法計算的均質(zhì)立方體梯度

在圖3中，實線和虛線分別為伴隨狀態(tài)法和有限差分法得到的梯度。從圖3不難看出，二者形狀大體一致，存在可以接受的誤差。這是因為有限差分對步長敏感，當采用不同步長時，該誤差隨步長而改變，因此說這個誤差來自有限差分的離散化。

2 SEG／EAGE目標區(qū)塊模型驗證

2.1 模型簡介

SEG／EAGE推覆體模型大小為20.00km×20.00km×4.65km，速度 范 圍 為 2 200～6 000m／s（圖4）。該模型模擬野外復雜地質(zhì)構造，主要由受構造運動影響拉張分裂的基巖和推覆于其上的沉積巖序列兩部分構成。其頂界面為一不整合面。在此不整合面之上，發(fā)育有12層連續(xù)沉積層序，每層厚度設計尺度與野外實際觀測相同，沉積序列被其底部發(fā)育的巖體分開。推覆體模型的巖相及層序地層學特征模擬了區(qū)域張應力作用下的海相與構造混合沉積分布于某些沉積層位中的河道及裂縫。在此區(qū)域中，所有封閉的圈閉均含有油氣。

圖4 SEG／EAGE推覆體模型

推覆體模型構造復雜，縱向層位多變并具有強烈的橫向速度變化，可以用來很好地檢驗全波形反演的適應性和其對復雜構造的成像能力。斷層的平面分布等構造現(xiàn)象可以考察全波形反演對構造細節(jié)刻畫的能力。尤為重要的是，中深層較小的含油氣圈閉，可以有效地體現(xiàn)全波形反演的高分辨能力以及儲層刻畫能力。

2.2 全波形反演參數(shù)設置

全波形反演要求大偏移距、寬方位角的觀測系統(tǒng)，因此我們設計了如圖5所示的固定檢波器采集系統(tǒng)，圖中黑線點和紅線分別為檢波器和炮線。

圖5 推覆體模型全波形反演－固定排列地表采集觀測系統(tǒng)

三維全波形反演發(fā)展的難點就在于計算量和單核內(nèi)存需求的問題。為了節(jié)省計算資源，測試全波形反演的效果，本次實驗采用最小75m網(wǎng)格，使用盡可能少的震源，最高反演頻率為7Hz。反演參數(shù)如表2所示，其中，Ngroup為頻率組標號；Nx，Ny，Nz分別為模型東西、南北及深度方向采樣點數(shù)；Δh為正演模型空間步長；Δt為正演時間步長；f為頻率組內(nèi)頻率范圍；Nsour和Nrec分別為震源和檢波器個數(shù)；Nite為迭代次數(shù)。

表2 SEG／EAGE推覆體模型全波形反演參數(shù)設置

本次實驗共需至少259 200次正演過程，因此，如果不解決計算量問題，該實驗就會無法完成。得益于雙并行算法的有效利用［9］，使得三維全波形反演的計算效率大幅提高。受計算資源限制，本次實驗使用了720個CPU，共耗時84.5h。

2.3 全波形反演結果

觀測數(shù)據(jù)由時間域正演獲得。在SEG／EAGE模型反演過程中，需要注意：①對2.6km深度處古河道重建以及其附近圈閉的準確刻畫（圖6a）；②對斷層面附近復雜構造地帶，常規(guī)成像手段難以成像區(qū)域的準確描述；③對模型各沉積層的準確標定（圖6b）。初始模型建立采用高斯平滑手段（圖7），目的是平滑掉所有沉積層和河道附近的圈閉。由圖7可見，在初始模型z＝2.6km處，古河道及其附近所有圈閉已全部平滑掉。由于速度對比度太大，且面積較大，其上部的高速異常體沒有被完全平滑掉，仍然存在（圖7a）。圖7b是該模型在y＝5.3km處的縱切面，從圖7b中可以發(fā)現(xiàn)，斷層附近復雜構造帶已被摧毀，較薄的沉積層也已全部消失，整個速度場的速度值都發(fā)生了變化，可以近似認為為一個梯度場。

圖8為第1組頻率（3.2～4.0Hz）的反演結果，不難看出，在4Hz已經(jīng)初步恢復了古河道的大致形態(tài)，但是細節(jié)部分還沒有任何顯示。縱向上比較厚的沉積層位以及斷層的大致走向得到了重建，但是較薄層位還沒有準確找到（圖9）。這是因為全波形反演的分辨率跟使用的頻率有極大的關系。全波形反演能提供最高1／2波長的分辨率，這里使用的最大頻率為4Hz，按照平均速度3 500m／s計算，只能找到大于437.5m的構造。這也體現(xiàn)了多尺度迭代法反演的優(yōu)越性，有效地降低了全波形反演的非線性。如圖9所示，在2.0km及2.5km處的沉積層，由于儲層較薄，使用第1組頻率沒有得到任何有效的改善。另外，此次計算進行了10次迭代，結果不能令人滿意，因為可以看到淺層的速度得到了有效恢復，而中深層僅僅只找到了正確方向，并沒有準確標定速度信息。但是由于本次實驗的目的是準確標定目標層，而準確重建速度需要更多次的迭代。綜合以上因素，此次的計算止步于此。

圖8 第1組頻率反演結果

圖9 第1組頻率反演模型與初始模型在x＝2.3km，y＝5.3km處的一維對數(shù)對比

使用第1組頻率反演的模型，作為第2組頻率（4.0～5.0Hz）反演的初始模型，得到了如圖10所示的結果。運用同樣的方法可以得到第3組頻率（5.0～7.0Hz）的反演模型（圖11）。這兩組頻率的初始速度模型與對應的最終反演所得模型的一維對數(shù)對比分別如圖12a和圖12b所示。

圖10 第2組頻率反演結果

第2組頻率的最大頻率為5Hz，從圖10a可以看出，古河道的形態(tài)得到有效恢復。河道內(nèi)部高速異常帶完全準確重建；從圖10b中可以看出，在第1組反演模型中沒有任何改善的2.0km及2.5km處的較薄沉積層開始向正確方向收斂，但是速度信息還不準確。第3組頻率已經(jīng)基本找出所有沉積層位（圖11），刻畫出比較精確的河道信息。河道附近的圈閉信息也有明顯顯示。中淺層速度信息得到了有效準確的恢復，但是深層部分只是收斂到了正確的方向。

如果反演繼續(xù)進行，更多的迭代次數(shù)或者更高的頻率將使深度信息得到有效改善，速度信息更加準確，但也意味著將占用更加巨大的計算資源。為了說明問題，驗證全波形反演所能達到的精度，我們縮小了模型，使用更少的炮點和檢波點，反演到更高的頻率，結果如圖13至圖15所示。不難看出，全波形反演可以精準描述儲層信息，精細刻畫復雜構造，可以基本恢復真實速度信息，達到常規(guī)速度建模無法實現(xiàn)的精度，與真實速度場幾乎完全吻合。圖16為實際模型、初始模型和反演模型在x＝5.00km，y＝3.37km 處的一維對數(shù)對比結果。從圖16可以看出，除了采集腳印帶來的影響，全波形反演基本恢復了速度值的真實信息，尤其是中淺層部分反演模型和真實模型幾乎完全重合。

圖16 實際模型、初始模型和反演模型在x＝5.00km，y＝3.37km處的一維對數(shù)對比

3 結束語

通過三維聲波全波形反演方法的理論論述及模型測試，可以得到以下認識。

1）全波形反演可以滿足高精度速度建模的要求，對復雜構造、較小儲層等具有精細刻畫的能力。

2）頻率域全波形反演聯(lián)合時間域正演的優(yōu)化算法的提出，使得三維全波形反演逐步走向實用化。

3）并行計算的應用大幅提高了全波形反演的計算效率，突破了全波形反演巨大計算量的瓶頸，使得三維全波形反演技術的實用化成為可能。

4）本次研究采用SEG／EAGE推覆體模型，模擬野外復雜地質(zhì)構造，對儲層及薄層沉積層的重建幾乎和實際模型完全吻合，達到了對模型高精度重建的目的。

全波形反演作為新興的技術，對資料的要求比較高（如信噪比，低頻信息等），這些原始資料的品質(zhì)問題需在野外采集過程中以予克服，或通過反演前的預處理過程使之提高，才能滿足全波形反演的應用要求。

［1］Lailly P.The seismic inverse problem as a sequence of before stack migrations［C］∥Conference on inverse scattering：theory and applications.Philadelphia：SIAM，1984：206－220

［2］Tarantola A.Inversion of seismic reflection data in the acoustic approximation［J］.Geophysics，1984，49（8）：1259－1266

［3］Gauthier O，Virieux J，Tarantola A.Two－dimensional nonlinear inversion of seismic waveforms：numerical results［J］.Geophysics，1986，51（7）：1387－1403

［4］Mora P R.Nonlinear two－dimensional elastic inversion of multi－offset seismic data［J］.Geophysics，1987，52（9）：1211－1228

［5］Pratt R G，Worthington M H.Inverse theory applied to multi－source cross－h(huán)ole tomography，part I：acoustic wave－equation method［J］.Geophysical Prospecting，1990，38（3）：287－310

［6］Sirgue L，Barkved O I，Dellinger J，et al.Full waveform inversion：the next leap forward in imaging at Valhall［J］.First Break，2010，28（1）：65－70

［7］Plessix R E，Perkins C.Full waveform inversion of a deep water ocean bottom seismometer dataset［J］.First Break，2010，28（1）：71－78

［8］Ben Hadj Ali H，Operto S，Virieux J.An efficient frequency－domain full waveform inversion method using simultaneous encoded sources［J］.Geophysics，2011，76（4）：R109

［9］Hu G，Etieen V，Castellanos C，et al.Assessment of 3d acoustic isotropic full waveform inversion of wide－azimuth obc data from valhall［J］.Expanded Abstracts of 82ndAnnual Internat SEG Mtg，2012，1－6

［10］Pratt R G.Velocity models from frequency－domain waveform tomography：past，present and future［J］.Expanded Abstracts of 66thEAGE Annual Conference，2004，Z99

［11］Plessix R E.A review of the adjoint－state method for computing the gradient of a functional with geophysical applications［J］.Geophysical Journal International，2006，167（2）：495－503

［12］Castellanos C，Etienne V，Hu G，et al.Algorithmic and methodological developments towards full waveform inversion in 3delastic media［J］.Expanded Abstracts of 81stAnnual Internat SEG Mtg，2011，2793－2798

［13］Sirgue L，Pratt R G.Efficient waveform inversion and imaging：a strategy for selecting temporal frequencies［J］.Geophysics，2004，69（1）：231－248

［14］Ravaut C，Operto S，Improta L，el al.Multi－scale imaging of complex structures from multi－fold wide－aperture seismic data by frequency－domain full－wavefield inversions：application to a thrust belt［J］.Geophysical Journal International，2004，159（3）：1032－1056

［15］Operto S，Virieux J，Amestoy P，et al.3Dfinitedifference frequency－domain modeling of visco－acoustic wave propagation using a massively parallel direct solver：a feasibility study［J］.Geophysics，2007，72（5）：SM195－SM211

［16］Virieux J，Operto S.An overview of full－waveform inversion in exploration geophysics［J］.Geophysics，2009，74（6）：WCC1－WCC26

［17］Hu G.Three－dimensional acoustic full waveform inversion：method，algorithms and application to the Valhall petroleum field［D］.Grenoble：Universite de Universite Joseph Fourier，2012