一種基于改進共軛梯度法的彈性波全波形反演速度分層建模方法

王立德，王小衛，周 輝，吳 杰，張志強，王建樂，王德英，馮 剛

（1.中國石油勘探開發研究院西北分院，蘭州 730020；2.中國石油大學（北京），北京昌平 102249；3.中國石油大慶油田有限責任公司采氣分公司，黑龍江大慶 163000）

0 引言

隨著油氣資源的開發利用，勘探開發成本低、類型相對簡單的構造型油氣藏大多已被發現和開發，而復雜構造區往往是一些成藏質量較差的區域，地震成像復雜，勘探難度大，因此高精度地震成像需求越來越迫切。彈性波全波形反演（EFWI）是利用多分量地震記錄的多種波場信息進行巖性參數反演，相比其他的建模方法具有更高的精度［1］。以全波形反演為代表的地震反演成像技術在地震地質建模綜合研究中也扮演著越來越重要的角色［2］，但是EFWI 存在很強的非線性，初始模型的建立、低頻數據的缺失等問題均可能導致反演陷入局部極小值，從而制約了EFWI 在實際資料（尤其是陸上實際資料）中的推廣應用。

全波形反演（FWI）需要直接求取目標函數對模型參數的導數從而進行速度更新，這將帶來巨大的計算量和存儲量。為克服大計算量問題，Lailly［3］提出了基于震源波場和殘差波場的梯度求取方法，其具體做法是對2 種波場進行一次互相關。該方法對于每個單獨激發的震源只要正演2 次就可進行1次更新，能夠減少全波形反演的計算量，由于該方法是基于線性理論的，對模型的大尺度背景恢復起來很困難。Tarantola［4］基于時間域縱波方程將非線性原理應用到最小二乘反演中，并推導出了利用伴隨方程法求取介質參數更新方向的公式。Tarantola［5］基于彈性介質中的波動方程形式提出并實現了利用彈性波進行全波形反演的思路，并將FWI 理論推廣到了EFWI。任浩然等［6］利用多種地震信息進行建模，但多重因素的疊加加大了反演的非線性并降低了實用性。為降低EFWI 的非線性程度，許多學者做了針對性的研究，最常用的是自地震數據低頻成分到高頻成分逐漸匹配的多尺度反演策略，該策略可以在時間域［7］、頻率域［8］和拉普拉斯域［9］實現，但后2種域的計算過程復雜，對計算資源的硬件性能要求更高。在時間域可以通過對地震數據的濾波實現多尺度反演，由于地震數據本身就是時空域的，對地震數據進行時間和空間域相關處理更方便。從小炮檢距數據到大炮檢距數據的反演，通過限制參與反演數據的炮檢距來實現。大炮檢距可以有效重建低波數信息，小炮檢距則對高波數成分具有很好的恢復作用［10-12］。另外，比較常用的降低EFWI非線性程度的方法是利用初至波進行近地表速度反演，因為EFWI 的非線性大部分由反射波引起的，當初始速度精度不夠時，反射點很難歸位［13］，因此，避免反射波參與反演可有效地降低反演的非線性。盡管很多學者針對EFWI的非線性問題進行了大量研究，并在理論模型上得到了驗證，但對于實際資料（尤其是陸上資料）而言，EFWI 的速度建模方法依然難以見到效果，這是因為陸上資料更容易受震源能量、激發響應、噪音污染、近地表傳播方式等的影響［14］，從而加劇EFWI 的非線性程度。王華忠等［15］指出EFWI的非線性取決于介質模型的復雜性和描述地震波物理傳播過程的正算子的復雜性，并在分析了陸上和海上地震數據特點的基礎上，提出了全波形反演走向實用化的策略，即利用特征波場反演而不是全波場，同時增加相位信息，盡量充分考慮初始模型的先驗信息。

EFWI 的非線性很大部分來自于深層反射波，這將直接影響縱、橫波速度場的反演精度。為提高速度反演的精度，提出一種分層反演策略：①利用時窗選取包含淺層速度信息的較早時間到達的地震數據，用小偏移距早至波數據，自低頻到高頻多尺度重建淺層速度場，以減少中深層反射波及噪聲帶來的強非線性影響，提高淺層速度的反演質量。早至波與初至波一樣，并不是特指某一類型的波，而是包含了很多波的信息。小偏移距早至波主要包括淺層直達波、反射波及小折射波等。篩選的具體標準是根據反演淺層界面的深度以及相應深度的初始速度場，估算出地震記錄時間，這樣選擇的偏移距范圍內的地震記錄時間能夠約束地震波穿透深度，從而獲得包含淺層速度場的波場記錄。在EFWI 中，淺層速度的反演精度將直接決定深層反演的好壞。在這個過程中用小尺度交錯網格有限差分模擬淺層波場。②固定好淺層速度，用全波場、全偏移距數據對中深層速度進行從低頻到高頻的多尺度反演，用大尺度交錯網格有限差分模擬中深層波場，在減少計算量的同時保證波場模擬的精度。③最后使用伴隨波場法求取縱、橫波速度梯度，引入深度學習優化器RMSprop［16］的梯度處理方式，提高反演的穩定性和收斂速率，進一步提高深層速度場的反演精度。將該分層反演策略應用到理論模型和陸上實際資料，以期有效提高淺層和中深層速度場的反演精度。

1 方法原理

1.1 時間域EFWI 基本原理

全波形反演的基本思想是用模擬數據不斷地逼近觀測獲得的實際數據，使它們的誤差最小化，此時擬合結果對應的模型參數最接近真實的模型參數。基于模擬和觀測記錄，全波形反演目標泛函往往可以被定義成多種不同的誤差形式，不同形式的目標泛函采用不同的參考標準匹配模擬數據和觀測數據。理論上如果目標泛函能夠達到一定范圍內的極小值，此時的模型參數就會更加接近實際情況。此次研究采用的EFWI 目標函數表達式如下：

式中：uobs(x,t)為觀測記錄；u(m;x,t)為正演模擬記錄，在分層反演方法中，模擬記錄按照分層時窗大小進行調整；dΩ 表示對檢波點接收區域各道的積分求和；dt表示對記錄時間的積分求和；x為地震道；m為模型參數，本文中指淺層及中深層的縱、橫波速度；T表示地震記錄的時間，s；Ω 為檢波點接收區域。

直接求取目標函數關于速度的梯度非常困難，比較常用的方法是通過伴隨波場求取的伴隨方法。Zhang 等［17］從Lagrange 乘子法原理出發，通過定義褶積型目標函數，得出了基于彈性波波動方程的伴隨方程，并給出了相應模型參數的梯度公式，本文采用此種方法求取模型參數的梯度。二維彈性波波動方程通過伴隨方程法得到的縱、橫波速度場的梯度表達式如下：

式中：ρ為密度，g/cm3，此次研究不涉及密度反演，因而為常數；vP，vS分別為縱、橫波速度，m/s；λ為伴隨波場值；v為震源波場值。

1.2 基于RMSProp 的共軛梯度優化方法

1.2.1 基本理論

EFWI 主要利用非線性最優化方法最小化觀測記錄和理論記錄之間的誤差，以重建地下介質的物性參數。為了最小化目標函數，通常使用梯度類優化算法來解決非線性問題。牛頓法具有二階收斂的優勢，相較于一階收斂的最速下降法和介于一階與二階收斂的共軛梯度法，它的收斂更快，但是需要直接求取目標函數的海森矩陣，難度較大，且需要巨大的計算量和巨大的內存資源成本，因此該方法在EFWI 中很難推廣應用。共軛梯度法（CG）的收斂速率介于一階和二階方法之間，可求解無約束二次規劃問題，由于其表達式簡單，實現起來容易，在大規模數據的無約束優化問題中具有獨特的優勢［18］。其主要邏輯是找到一組共軛正交的方向向量，對每個方向上的模型參數依次進行迭代更新，當每個方向上的誤差達到或接近0 時，迭代結束。該方法的收斂速度大于最速下降法等一階的優化方法。它的模型參數更新方式可以用下式表示：

式中：gn為公式（2）和（3）所計算的縱波和橫波速度的梯度場；n為迭代次數；dn為通過共軛梯度法計算的更新方向，將其代入下式進行模型參數的更新

式中：α為更新步長。

RMSProp 來源于深度學習，本文將這種優化理論應用到共軛梯度法中，相比于傳統共軛梯度方法，新方法可使EFWI 收斂更快、反演精度更高。優化后需對縱、橫波速度場每個位置的梯度gn做如下處理：

式中：gn,i為縱、橫波速度梯度gn的第i個元素；為對gn,i移動加權平均的結果；η為常數，一般取值為0.99；為處理后的縱、橫波速度梯度元素，所有形成了新的梯度場。將代入到公式（4）和（5）便可進行縱、橫波速度的迭代更新。相對于傳統的共軛梯度法，優化后的方法具有如下優勢：①可使當前參與更新的梯度包含先前若干次迭代時的梯度信息，可保證整個反演過程的穩定性并提高收斂效率。②公式（7）的作用相當于實現了搜索步長的自適應調整，當反演梯度gn,i非常大時，計算得到的也很大，將變小，代入公式（4）計算的dn,i相應變小，按照公式（5）的步長α進行速度更新時更新量亦相應減小，可避免發生因更新過量而越過目標泛函極小點且在極小點附近擺動的現象，從而提高了反演精度。反之，梯度gn,i非常小時，通過增大梯度可以加速目標函數的下降，使收斂更快。③在EFWI 中，深層的梯度能量很弱，導致深層區域縱、橫波速度場的重建效果不理想。本文提出的公式（7）的處理方式隱含淺層和中深層速度場梯度能量的均衡化處理，使深層的梯度能量相對提升，提高了深層速度的反演質量。

1.2.2 數值算例

為驗證基于RMSProp 改進的CG 方法對EFWI效果的提升，選擇了二維Overthrust 模型進行反演測試。在反演參數相同的情況下，分別用傳統的CG 法和改進后的CG 法重建該模型（圖1）。

圖1 Overthrust 模型的真實速度和初始速度Fig.1 True and initial velocity of Overthrust model

該二維Overthrust 模型共計800×170 個網格點，交錯網格有限差分空間步長為8 m；使用雷克子波生成震源，模型設計50 炮，每炮800 個檢波器；為降低非線性影響，選擇多尺度反演，共設置4 個頻帶，各頻帶的主頻依次是4 Hz，7 Hz，11 Hz 和15 Hz，每個頻帶迭代20 次。對比2 種方法對應的反演結果可以看出，無論是縱波速度還是橫波速度，改進后的CG 方法重構的速度場精度更高，尤其在深層，構造更加清晰，速度值更接近真實值，說明本文方法能夠有效平衡淺層和深層的梯度能量，從而改善深層反演的結果（圖2）。

圖2 不同方法反演的Overthrust 模型速度場對比Fig.2 Comparison of velocity field of Overthrust model inversed by different methods

為了更加直觀地對比2 種方法的反演精度，分別求取了2 種方法在整個反演過程中的速度誤差值，該誤差是基于真實速度模型的相對誤差。在4個頻帶共80 次的迭代過程中，改進CG 反演的速度誤差曲線均位于傳統方法誤差曲線下方，誤差曲線下降更快，說明優化后反演的速度場精度更高，收斂更快，通過更少的迭代次數即可達到傳統CG 反演的精度（圖3）。

圖3 不同頻帶反演速度的相對誤差Fig.3 Relative errors of inversion velocity in different frequency bands

1.3 分層反演策略

1.3.1 基本理論

EFWI 強非線性是由資料采集等客觀因素引起的地震數據與反演參數之間的強不確定性，因此通過篩選不同數據子集，分步驟、多尺度的反演是降低EFWI 非線性程度的有效策略（圖4）?；谶@一規律的時間域EFWI 分層反演策略，可降低反演的非線性程度，提高反演質量。具體反演過程如下：

圖4 速度場分層反演流程Fig.4 Flowsheet of layered inversion of velocity field

（1）利用反演淺層界面的深度以及相應深度的初始速度，估算出相應地震記錄時間，用該時間窗口選取較早到達的地震數據，同時選定偏移距范圍。用包含豐富淺層速度信息的小偏移距直達波、小折射波［19］、淺層反射波等早至波反演淺層速度，規避主要由中深層反射波帶來的強非線性干擾。

（2）通過從低頻到高頻逐漸匹配的多尺度反演手段，重建淺層速度。反演時采用小網格剖分的交錯網格有限差分方法進行波場模擬，以保證淺層速度波場模擬的精度，提高淺層反演質量。

（3）固定反演重建的淺層速度場，約束中深層反演，降低反演的不確定性，利用全采樣時間的全偏移距數據反演中深層速度。

（4）由低頻到高頻多尺度反演中深層速度，由于中深層速度較大，采用大網格的交錯網格有限差分模擬波場，利用含有中深層速度場信息的地震記錄段反演中深層速度場。該方法在減少計算量的同時可獲得高精度的中深層反演結果。

分層反演時速度場梯度需要按照以下公式進行處理

式中：d1表示淺層速度場的深度，m。梯度gn經過公式（9）定義的分段處理函數p(z) 處理后，淺層梯度值為0，不再更新，而深層梯度的能量得到提升。

1.3.2 數值算例

為驗證分層反演策略的有效性，利用二維Marmousi 模型進行測試（圖5）。速度更新時，對二維Marmousi 模型的初始速度場分別進行常規EFWI和分層反演策略下的EFWI。整個速度場的地震記錄采樣時長為4 s，采樣間隔為1 ms，共40 炮激發，用0～250 ms 的時間窗篩選地震記錄以確定偏移距范圍，進行多尺度反演的雷克子波主頻依次為4 Hz，8 Hz，12 Hz，18 Hz 和25 Hz。淺層深 度 為300 m，真實縱、橫波速度場深度為0～300 m 時為低勻速層，有限差分網格為5 m×5 m。淺層反演結束后固定淺層速度，用整個4 s 的全偏移距地震記錄進行多頻帶反演中深層速度場，有限差分網格為10 m×10 m。對比基于2 種反演策略得到的縱、橫波速度場，可見淺層低勻速層反演結果在反演參數完全相同的情況下，無論是縱波速度還是橫波速度，分層反演的低速層均更干凈，更接近真實速度場，而常規反演在低勻速層存在干擾；同時，在深度為0～1.5 km，距離為0～6 km 的區域內，分層反演得到的縱、橫波速度場均更加清晰，局部構造的分辨率更高，弱反射層得到了較好的恢復（圖6）。

圖5 二維Marmousi 模型的真實速度和初始速度Fig.5 True and initial velocity of 2D Marmousi model

2 實際資料應用

為了進一步測試新方法對實際資料的有效性和適應性，對四川某三維區塊的二維測線進行EFWI速度分層建模測試。該三維區塊為正交觀測系統采集，最小炮檢距為28.28 m，覆蓋次數為11，道距為40 m，子波主頻為25 Hz，三分量地震記錄的采樣時長為5 s，采樣間隔為1 ms。X分量地震記錄信噪比低，很難看到有效信息，Y分量信噪比更低，Z分量信噪比相對較高（圖7）。通過層析反演得到該區塊的初始速度，待反演的二維速度場橫向距離為5 km，縱向深度為7 km。

分層反演時，給定縱、橫波初始速度場（圖8），采用0～1 s 時窗反演深度為0～2 500 m 的速度場。反演該深度時有限差分的空間采樣間隔為15 m，反演更深的速度場時空間采樣間隔為25 m。共選取300 炮近偏移距數據參與淺層速度的反演，反演設置了5個頻帶，采用主頻依次為6 Hz，10 Hz，15 Hz，20 Hz 和25 Hz 的雷克子波進行波場模擬。從淺層速度場的反演結果（圖9a，9b）來看，淺層速度場得到了有效更新。由于淺層速度已經被重建，與待反演的深層速度在交界位置產生了不耦合層，這將在下一步的反演中產生虛假反射影響深層速度建模，因此在反演深層速度場時，需要對該不耦合層進行光滑處理以降低虛假反射的影響。在平滑不耦合層后，固定淺層速度場，迭代更新深層速度，由于深度2 000 m 附近的速度場被重新平滑，因此反演深層速度時需要從這一深度開始反演。深層速度場反演時采用了300 炮的全偏移距全記錄時長的數據，除空間采樣間隔和所用地震數據范圍外，所有的反演參數與淺層一致。從全速度場最終的反演結果（圖9c，9d）可以看出，淺、中、深層速度場均被更新，可見水平形態的構造趨勢。

圖8 四川某三維區塊的初始速度場Fig.8 Initial velocity field of a 3D block in Sichuan

圖9 四川某三維區塊分層反演后的速度場Fig.9 Velocity field after layered inversion of a 3 D block in Sichuan

高斯束彈性波偏移共成像點（CIP）道集的拉平程度可用于驗證速度場的精度，分別利用初始速度場和反演重建速度場進行偏移并抽取相同位置的共成像點道集進行對比，同相軸拉平程度高的對應速度場越精確。對比四川某工區實際資料的縱、橫波共成像點道集可以看出，相比初始速度場，使用分層反演后的縱、橫波速度偏移得到的共成像點道集同相軸拉平程度均更高、同相軸更清晰，能量亦更為集中，說明優化后速度場的精度更高（圖10）。

圖10 四川某三維區塊高斯束彈性波偏移共成像點（CIP）道集Fig.10 Common imaging point（CIP）gathers of Gaussian beam elastic wave migration of a 3D block in Sichuan

3 結論

（1）相較于常規的基于走時的反演方法，EFWI通過近似模擬地震波在彈性介質中的波動方程，充分利用地震信號的運動學和動力學特征進行反演，速度建模的精度更高。

（2）利用分層反演的策略，降低制約EFWI 的非線性因素影響，提高淺層和中深層縱、橫波速度反演的精度；改進的共軛梯度法可提高反演的穩定性，提高收斂速率。同時，改進的共軛梯度法隱含淺層和中深層速度場梯度能量的歸一化，使深層的梯度相對提升，提高了深層速度的反演精度。

（3）分層反演策略，由淺到深，由近到遠地選取地震資料，可減少反演過程中的非線性影響，提高反演準確性；還可兼顧淺層低速和深層高速采用大小不同的有限差分網格，在減少計算量的同時保證波場模擬的精度，從而提高反演速度場的精度。