基于波動方程解析解的塊約束廣義聲阻抗反演

李遠強 霍志周 李景葉*

陳小宏①②③ 張 健①②③ 耿偉恒①②③

（①中國石油大學（北京）地球物理學院，北京102249；②油氣資源與探測國家重點實驗室，北京102249；③海洋石油勘探國家工程實驗室，北京102249；④中國石化石油勘探開發研究院，北京100083

0 引言

地震反演技術可以利用地震數據的旅行時、振幅和波形等信息估計地下巖石的彈性參數［1-2］，在儲層預測與評價、油藏特征描述等方面得到了廣泛的應用［3-4］。作為儲層預測與表征的核心技術，疊后阻抗反演是應用最廣、最成熟的反演方法［5-9］。該技術經歷了從最初直接遞推反演［5］到基于模型的迭代反演［7-10］的發展過程。

隨著油氣勘探開發的不斷深入，要求阻抗反演結果達到一定的分辨率，便于更好地識別薄層［11］。然而，下列兩個原因限制阻抗剖面的分辨率：①基于單一界面情形計算反射系數，正演方法不完備。假設地震數據僅包含一次波信息，忽略波的傳播效應以及地層厚度的影響（調諧效應）。傳播效應如透射損失以及調諧效應將造成地震記錄的振幅異常，進而對反演結果的數值準確度以及分辨率造成影響。而層間多次波將會當作有效反射，造成錯誤的反演解釋結果，特別是存在薄層且彈性參數差異較大的情況下［12］。同時，當前地震處理技術是無法完全消除透射損失、調諧效應、層間多次波的影響。②阻抗反演算法大多以最小二乘為基礎，并加入合適的先驗約束獲取穩定的反演結果［13］。然而這些約束往往會使得反演結果過度光滑，無法清晰刻畫地層的邊界［14-16］。因此只有突破上述基本假設，才能得到分辨率更高的反演結果。

相對于常規阻抗褶積正演，將波動方程的解作為正演工具獲取全波場響應，可以更為準確預測地下介質的彈性參數［17］。目前，波形阻抗反演方法有兩種思路。一種是使用數值求解方法（有限差分、有限元等）近似求解波動方程，從而進行反問題求解獲取阻抗結果，可稱之為全波形反演（FWI）［18-19］。由于FWI對計算機儲存能力以及計算效率要求極高，實際實施仍然存在困難。另一種是對層狀假設下的波動方程解析求解并執行一維介質波形反演［20-21］。Ursin等［22］證明常規阻抗正演方法是解析解的零階近似，因此利用解析解波形反演阻抗，結果將更為真實合理。

其次，可通過正則化約束獲取高分辨且穩定的阻抗結果［17，23］。稀疏和平滑約束在地震阻抗反演中廣泛應用。平滑約束如高斯約束、L2范數約束等往往會忽略地下介質的分層特性，使得解過度光滑，無法清晰刻畫地層的邊界。稀疏約束如柯西約束、L1范數約束雖然可以獲取高分辨率的解，但無法避免局部極小，從而引起某些地質假像［24］。借鑒圖像處理中的圖形銳化算法［25］，Theune等［24］引入微分拉普拉斯塊狀約束，通過多道反演改善結果的垂直分辨率，但由于多道反演需要求解大型稀疏矩陣，將該約束簡化到一維情況是更好的選擇。

另外，疊前道集的優化過程中，特別是最為重要的角度道集提取環節［26］，多數情況下基于聲學介質假設，勢必會導致道集趨向于聲學AVO 特征。而疊后成像、反演過程均忽略炮檢距對反射系數的影響，部分疊加剖面未得到相應的重視。根據聲波在界面的能量分配［27-28］可知，聲波反射系數隨入射角／炮檢距變化，且隨著入射角的增大而增大。借助于不同入射角度的反射信息，可更為準確反演彈性參數，尤其是速度信息。由于直接反演密度不穩定，需借鑒彈性阻抗反演技術［2］，對部分疊加資料進行廣義阻抗反演，進一步計算速度和密度信息，以避免密度反演的不穩定性。

借助上述三項關鍵技術，本文提出基于塊狀約束的一維聲波波形廣義聲阻抗反演技術，用于提取更為真實、可靠且穩定的速度和密度參數。首先推導了兩界面聲波反射系統傳播函數的遞歸表達式，最終獲取不同入射角的反射響應和Fréchet導數。基于貝葉斯反演框架，在原始目標函數中加入塊狀約束，在壓制病態數值噪聲（即觀測數據的小擾動對應的模型數據將產生巨大擾動的數值噪聲）的同時保持地層的邊界特征。針對目標函數是一個非線性問題，可利用高斯牛頓法求解，從而獲得一個穩定的高分辨率波阻抗結果。利用廣義聲阻抗與入射角、速度和密度關系求解穩定的速度和密度。最終，通過模型以及實際數據驗證新方法的有效性和實用性。

1 原理

1.1 聲波方程反射響應解析求解

聲波方程的解析求解往往基于矩陣遞歸框架，但該算法框架需額外計算除反射響應以外的其他響應（如透射響應等），計算效率較低。另外，該算法的輸入為速度和密度，無法直接輸入聲阻抗參數用于疊后反演。本文修改遞歸框架，將聲阻抗作為算法輸入，推導僅包含反射響應的遞歸計算公式，將矩陣運算簡化為單元素運算，由于沒有額外的計算量，可以迅速獲取總的反射響應。

圖1為兩個界面聲波反射系統示意圖，Di（ω）表示第i層頂界面的下行波，均勻介質中簡諧波傳播可表示為exp（－jωt），則第i層底界面的下行波可表示為

式中hi和vi分別為第i 層厚度、速度。同樣第i層底界面的上行波可表示為

式中Ui（ω）表示第i層頂界面的上行波。

圖1 兩個界面的上、下行波關系

根據Treitel等［29］關于聲波界面能量分配理論，可得

式中：Ri、Ti分別為下行垂直入射時第i 個反射界面的反射系數和透射系數；R′i、T′i分別為上行垂直入射時第i個反射界面的反射系數和透射系數。反射系數和透射系數滿足

聯立式（1）～式（4），可得

用式（6）除以式（5），有

根據傳播響應函數的定義ri（ω）＝Ui（ω）／Di（ω），表示反射界面i以下的總反射響應。因此可得反射響應的遞歸計算公式為

垂直入射時反射系數為

根據聲波在非垂直入射界面時能量分配規則［11］，有

式中入射角θ與透射角θ′滿足Snell定律p＝sinθ／vi＝sinθ′／vi＋1，其中p 為射線參數。

最深第M 層不存在反射，因此rM＝0。至此可由式（8）得到總反射響應函數r0。用r0乘以頻率域子波s（ω），并進行傅里葉反變換，可獲對應地震記錄

此時整個算法的輸入是速度、密度。將輸入參數從深度域轉換到時間域，并利用復雙程旅行時間修改式（8），可得

式中τi＝Δt，Δt為時間采樣間隔。

定義廣義聲阻抗為

界面的反射系數為

則可將廣義聲阻抗作為算法的輸入用于計算不同角度入射的地震響應。

將輸入參數從深度域轉換到時間域，將算法與常規疊后反演預處理流程（建立初始模型）有機地結合在一起，便于對比。此外，還可以提高反演問題的穩定性，同時降低該問題非線性程度。本文借助殘差函數圖（RFM）作為轉換前后反演問題的適定性評價標準［30］。通過對相鄰層的波阻抗擾動，根據深度域和時間域擾動所建立的RFM 如圖2所示。十字虛線的交點表示對應的真實值。等高線越密集，殘差的擾動變化率越大，反演問題的穩定性越差；對應橢圓的離心率越大，該問題非線性程度越高。由此可見，當輸入參數從深度域轉換到時間域，反演問題的適定性會有所改善。

1.2 Fréchet導數解析求解

上述反演問題的損失函數項（模型估計值與觀測值之間的差異度量）是非線性的，對于非線性最優問題一般可通過全局優化算法如模擬退火、遺傳算法等以及梯度類下降算法求解。全局優化算法可免于計算梯度但計算代價極大，本文使用梯度下降算法。該算法的關鍵在于求解Fréchet導數，即地震記錄相對于模型參數的導數

圖2 深度域（a）與時間域（b）RFM 對比

由于子波與模型參數無關，僅需計算反射系數的Fréchet導數。可通過鏈式法則獲取對應導數元素表達式

其中

根據式（14），有

式中

圖3 為常規線性正演矩陣與反射系數的Fréchet導數矩陣對比。線性矩陣（圖3a）呈現條帶狀，帶寬為2，隨著深度變化其數值并未發生變化，除條帶外其他元素均為0。這與線性假設一致，某點反射系數僅與相鄰層的彈性參數有關。恰恰相反，導數矩陣（圖3b）呈現條帶狀，帶寬為6～8個元素；受傳播效應影響，隨著深度變化矩陣數值逐漸變小；由于計算的全波場響應，條帶周圍元素不為0，說明某點反射系數與整個層狀介質的彈性參數都有關，只是該點周圍的幾個層占主要因素。

圖3 線性矩陣（a）與導數矩陣（b）的對比

1.3 貝葉斯反演框架構建

正演過程可表示為

式中：d 為觀測數據向量；m 為模型參數向量；G為m 映射到d 的非線性算子，本文中G（m）為解析法正演的結果；n 表示二者之間的誤差。模型參數的先驗概率分布為

式中：S（m）為模型參數的正則化項；P0m為模型參數的歸一化常數。假設噪聲相互獨立并且服從高斯分布，利用貝葉斯公式可得模型參數的后驗概率分布為

當模型的先驗信息服從高斯分布時，有

高斯分布約束是一個過度平滑的約束。地下介質的成層性要求反演結果的垂向梯度在層內為零，而層與層之間有個較大的值。這意味著概率分布函數集中在零點附近，而正負無窮處較小，微分拉普拉斯分布恰好滿足這樣的特征。因此可在平滑約束提供一個穩定解的基礎上，再加上微分拉普拉斯約束項，從而更好地刻畫地層邊界。因此需要定義兩個先驗權重系數λs、λb，平衡反演結果的噪聲、光滑度以及邊界刻畫能力。此時目標函數為

式中最后一項就是微分拉普拉斯約束一維表達式，其中D 是一階差分矩陣

可使用高斯—牛頓算法求解該目標函數，即

式中：H 為海森矩陣；k為迭代次數；梯度為

式中第一項表示損失函數對模型的導數，其中gk＝?G（mk）／?m，即對應正演問題的Fréchet導數矩陣；第二、第三項為平滑約束以及塊狀約束對模型的導數。如果問題的非線性程度不嚴重或者初始模型接近真實解時，海森矩陣H 的二階項可以忽略，則近似為

式中：E 是單位矩陣；Q 表示塊狀約束的對模型二階導數，即

在實際資料反演中，一般情況下通過井旁道反演測試，人為選取最佳權重系數λs、λb。

1.4 速度、密度參數提取

根據廣義聲波阻抗定義（式（13）），垂直入射和非垂直入射阻抗為

則有

式中vˉ 為低頻背景速度。

2 數值模擬測試

2.1 正演模擬分析

為突出解析解正演方法的優勢，首先使用厚層模型驗證傳播效應對一次反射波的影響。圖4a為對應的厚層阻抗模型，利用主頻為30Hz的Ricker子波合成零炮檢距地震記錄，對比常規阻抗正演方法與彈性聲波解析解方法結果（圖4b）。由于常規正演方法，只計算一次波，未考慮多次波以及透射損失，可以發現，透射損失的影響主要體現在振幅衰減上，隨著深度增加，透射損失將越來越嚴重，400ms處也存在少許多次波。在400ms處加入一組薄互層（圖4c）突出多次波，合成記錄如圖4d所示。可以看出，多次波（藍色方框所示）以及調諧效應都相當嚴重，400ms以下的一次反射均受到多次波干涉。由于薄層存在，透射效應加劇。當彈性參數差異較大時，傳播效應強烈影響一次反射波，深部反射受到強烈“污染”，并且當前地震處理技術無法完全消除這些影響。因此采用具有理論完備的正演方法，計算一次波的同時考慮層間多次波以及透射損失效應，有利于解決薄層問題。

圖4 常規阻抗正演方法（黑線）與聲波解析解方法（紅線）模擬的厚層和薄互層模型零炮檢距記錄對比

圖5 聲波介質反射系數隨入射角的變化曲線

計算不同速度變化率及密度變化率為40%的聲波界面反射系數隨入射角的變化曲線（圖5），可見：隨著入射角增大，速度變化導致的界面反射增強；速度差異越大，“AVO”效應越明顯；但40%的密度變化率并不會產生AVO 效應，即地震記錄對密度變化不敏感，這也是常常將介質假設為常密度的原因。圖6左為利用薄互層模型（圖4c）合成的聲波角度道集，其中子波主頻為30Hz，角度范圍為0～33°，角度間隔為3°。沿圖6左紅色虛線處提取的均方根振幅曲線如圖6 中所示，其變化特征與圖5一致。將不同角度的記錄與零入射角道相減，如圖6右所示，可見：當角度較小時，“AVO”效應可以忽略；當角度較大時則不能忽略。利用“AVO”效應可以提高速度信息的準確度，同時密度可根據阻抗和速度獲取，間接提高密度參數的精度。因此，通過近、遠部分疊加剖面振幅信息獲取較為準確的速度、密度信息是可行的。

圖6 薄互層模型聲波角度道集分析

2.2 模型數據反演測試

為驗證本文構建的反演算法的有效性，利用薄互層模型進行反演測試。將聲波解析解模擬結果作為輸入，假定其為尚未消除傳播效應影響的真實地震數據。分別使用常規和本文方法進行零入射角數據聲阻抗反演，結果如圖7a和圖7b所示。由于傳播效應以及調諧效應使得地震記錄振幅異常，常規聲阻抗反演結果誤差大，對分界面的刻畫不夠清晰。如圖7a黑色橢圓所示，由于層間反射無法消除，常規聲阻抗反演過程中，將其當作有效的一次反射，產生錯誤的反演結果。恰好相反，本文方法獲取了一個更為準確且邊界刻畫清晰的阻抗估計結果（圖7b）。與真實阻抗相比，二者殘差的歸一化均方根誤差為0.02。圖7c、圖7d為入射角為20°時常規和本文方法的廣義聲阻抗反演結果，與零角度阻抗反演結果相似。與之不同的是，阻抗在數值上變大。進一步提取的速度密度曲線如圖7e、圖7f所示，與真實速度、密度模型有較好的一致性。

圖7 模型數據反演結果

2.3 噪聲測試

利用阻抗井數據（圖9 黑線所示）合成地震記錄，并對地震數據加入噪聲測試算法的魯棒性。該模型數據在500ms以下阻抗參數的成層性較好，可用于驗證塊狀約束的有效性。圖8為加噪前、后地震記錄對比（噪聲與地震記錄最大絕對值之比，An／s＝36%）。圖9a、圖9b分別為相同條件下（相同初始模型、迭代次數等）高斯平滑約束以及塊狀約束的聲阻抗反演結果。明顯可見，平滑約束的反演結果分辨率低于塊狀約束的反演結果；塊狀約束對邊界刻畫更為清晰（箭頭所示）。為更好說明該問題，在不同噪聲水平（An／s）情況下，不加約束、平滑約束以及塊狀約束的情況下進行反演測試，計算反演結果與真實結果的相對誤差如圖10所示。當噪聲水平為零時，算法是穩定的，因此相對誤差幾乎為零；隨著噪聲水平的增加，相對誤差均會增大；不加約束時的相對誤差較大，高斯約束和塊狀約束都能獲取較為穩定的結果，但塊狀約束的相對誤差小于高斯約束。因此本文反演方法在高噪聲情況下依然穩定可靠。

圖8 加噪前（紅色）、后（黑色）地震記錄對比

圖9 不同約束條件的含噪數據反演結果

圖10 不同約束條件反演結果相對誤差隨信噪比變化曲線

3 實際數據測試

實際數據近炮檢距部分疊加剖面如圖11a所示，勘探目標是2300～2500ms附近的濁積扇砂體（紅色矩形框所示）。基于拾取的層位插值得到的初始阻抗模型，如圖11b所示。圖12為利用統計原理提取的子波。圖13為本文方法和常規方法反演的聲阻抗結果，其中圖上黑色曲線為根據測井數據計算的結果。圖14為實際地震記錄與聲波反演結果正演記錄的對比，二者之間殘差較小。圖15為測井數據（黑色，經Backus平均處理）與井旁地震數據反演結果（道號為100）對比，可以發現，本文方法反演結果（紅色）與測井數據匹配度更髙，與Backus平均后測井數據相關系數高達92.5%，而常規反演結果（藍色）與之相關系數僅為84.6%。與常規反演方法相比，本文方法提高了薄層的數值精度，且邊界刻畫更為清晰，分辨率更高（圖13、圖15紅色箭頭所示）。圖16a為遠炮檢距部分疊加剖面，由于復雜的地下傳播效應影響，遠炮檢距疊加剖面分辨率低于近炮檢距疊加剖面，導致對應的廣義聲阻抗反演結果分辨率（圖16b）也較低。應用Gardner公式作為低頻先驗信息，根據近、遠炮檢距反演結果提取穩定的速度和密度結果如圖17所示，與地震記錄有良好的一致性，但分辨率高于地震記錄，尤其是密度剖面。

圖11 近炮檢距部分疊加剖面（a）及初始阻抗模型（b）

圖12 基于統計理論估計的子波

圖13 本文方法（a）和常規方法（b）的阻抗反演結果

圖14 實際地震記錄（上）、反演結果的正演地震記錄（中）及其殘差（下）

圖15 兩種方法反演結果與測井曲線的對比

圖16 遠炮檢距部分疊加剖面（a）和本文方法的廣義阻抗反演結果（b）

圖17 提取的速度（a）和密度剖面（b）

4 結論與討論

本文基于一維聲波方程解析解，提出非線性聲阻抗反演方法。該方法充分考慮了透射損失、多次波的影響，是對常規聲阻抗反演方法的補充和發展。借鑒于圖形銳化技術，本文引入塊狀約束獲取高分辨率的穩定解。模型以及實際資料反演結果均表明，新方法可以提供更準確且分辨率更高的反演結果，且邊界刻畫清晰。利用廣義聲波阻抗反演結果，可進一步提取速度與密度信息，也可將其作為一種地震屬性用于流體識別。

本文反演方法不需要去除層間多次波，但需去除與表層相關的多次波，考慮橫向變化的球面擴散補償技術是必須的處理流程。

雖然本文推導的遞歸算法比遞歸矩陣正演算法快10倍左右，但執行單道反演仍然需數秒時間，可通過三個途徑提高計算效率：首先，通過井旁道反演選取最優反演迭代次數；其次由于常規正演方法是解析解的零階近似，可將常規反演結果作為初始模型減少迭代次數，加快收斂；最后，該方法為單道反演算法，易于并行計算，可通過Open MP 實現并行加速。

基于廣義聲波阻抗反演方法具有理論完備性，考慮了全波場信息，可將其推廣到彈性波阻抗反演和相應彈性參數的提取。