基于單程波算子的全波場最小二乘偏移方法

葉月明，鐘世超，范國章，任浩然，王兆旗

（1. 中國石油杭州地質研究院，浙江杭州 310023；2. 北京理工大學長三角研究院，浙江嘉興 314019；3. 浙江大學先進技術研究院， 浙江杭州 310027）

0 引言

在海洋勘探中，由于海水面和海底兩個強波阻抗界面的存在， 海洋地震數據普遍發育能量強的多次反射波。在常規的地震資料處理中， 多次波在地震剖面上表現為構造假象，且與真實構造難以有效區分，這影響了解釋人員對地下真實情況的判斷，尤其是微幅度構造下的儲層更易受多次波假象干擾。所以長期以來多次波一直被視為噪聲，需盡可能地從地震數據中壓制，以免在地震資料解釋中造成誤解［1-6］。事實上，多次波與一次波一樣也是反射信號，它攜帶了豐富的地層信息。在逆掩構造和高速鹽丘等復雜地質背景下，一次波經常由于照明不足出現成像陰影區，而多次波卻能通過多次反射“照射”到一次波難以觸及的地方。此外，多次波較小的入射角使其具有更高的垂向分辨率［7-8］。一次波雖然是地震數據的主體，但多次波所攜帶的信息對于地震成像的幫助也是不容忽視的。完整的地震成像應該是全波場（一次波、表面多次波和層間多次波）的貢獻。為獲得全面準確的地下信息，全波場成像是地震資料處理技術發展的必然趨勢。

野外采集的陸地資料相對復雜，一般包含直達波、面波、折射波、繞射波、一次波和多次波等，所以陸地資料的全波場成像非常困難。海洋地震資料不受近地表的影響，具有較高的信噪比，發育的波場主要是一次波和多次波，因此就全波場成像研究而言，海洋地震資料具有優越的數據基礎。對于海洋深水地震資料，根據實現方法的不同，利用多次波成像方法主要有四類。

第一類是高階多次波轉化為一次波后的偏移成像。該類方法的主要思想是將高階多次波降階為一次反射波，再利用現有的一次波成像技術。Wapenaar等［9］基于地震波干涉理論，將共炮點道集內的各道互相關生成虛擬炮記錄，在相同位置疊加后將多次反射波轉變為一次反射波，再按照一次波偏移算法成像［10］。利用聚焦變換策略［11］，將原始地震記錄的反射波能量聚焦到零時間點，可實現高階多次波向低階多次波的轉變［12］。單國健［13］通過互相關技術將多次波轉化為準一次波，再以一次波地表接收點為準震源對準一次波進行偏移成像。

第二類方法是全波場的直接偏移成像。這類技術與面炮偏移相似［14］，將多個檢波器接收到的多次反射波認為是多個震源同時激發的一次反射響應。常規的單程波偏移和逆時偏移均可修改為多次波直接偏移成像。Berkhout 等［15］闡述了以含多次波數據為震源、表面多次波為記錄的多次波成像理論，并指出成像會受不同程度干涉噪聲影響［16］。Muijs等［17］以分離出的上行波場為震源、下行波為記錄，利用單程波延拓和反褶積成像條件獲得準確成像的同時壓制了部分干擾噪聲。Lu 等［18］在三維SEAM 模型上驗證了多次波對鹽下照明的優勢。葉月明等［19］將基于單程波算子的多次波成像技術應用到海洋深水實際資料處理中，照明均衡度與垂向分辨率均有顯著改善。隨著計算機技術的進步，全波方程逆時偏移技術得到了快速發展。Liu 等［20］修改了逆時偏移的上、下行波場，提出了多次波逆時偏移（RTM）技術，類似的策略也成功應用到海底地震儀（OBS）數據［21］。以子波和地震記錄為震源、含多次波數據為記錄可實現一次波與多次波的聯合成像，Ye 等［22］分別對一次波和多次波成像，再進行匹配疊加以避開子波估算。正向延拓一次波或低階多次波與反向延拓高階多次波產生的相干噪聲是多次波偏移的固有難題，Wang等［23］通過曲率差異在局部角度域成像道集上衰減干涉噪聲。基于波場分解和立體成像的方法也可壓制相干噪聲［24-26］。分離混疊在一起的多階多次波，形成多個單階多次波后再分別成像也可避免假象的產生，利用類似的方法可對層間多次波成像［27］。以表面多次波為上行波場，按照一次波的偏移方法可預測出表面多次波假象，再通過匹配相減在成像域壓制震源與多次波產生的強能量地震假象［28］。Liu等［29］延伸了單階多次波成像方法，提出了基于相位編碼的多階多次波RTM 技術，在提高計算效率的同時顯著壓制了噪聲。

第三類方法是基于反演理論的全波場成像方法。基于全波場的最小二乘偏移（LSM）技術有三種實現策略：①基于旅行時計算的射線追蹤方法，其關鍵是多次波地下傳播旅行時的求取。Brown 等［30］提出了在CMP 域的全波場最小二乘Kirchhoff 偏移，并利用懲罰系數壓制干涉噪聲。②修改上、下行波場的LSM，這種方法與第二類“全波場的直接偏移成像”的區別在于通過反演壓制干涉噪聲。Zhang等［31］以預測出的表面多次波為上行波場、包含多次波的記錄為下行波場實現了多次波的LSM；Wong 等［32］將該方法擴展到海底電纜（OBS）數據的成像，增強了地下照明度，并證明不同階次表面多次波單獨成像同樣可以壓制干涉噪聲。以子波震源和預測出的表面多次波為下行波場、以表面多次波為上行波場能夠實現一次波與表面多次波的LSM，但若不能準確估計子波震源，同樣會產生較強干擾噪聲［33］。Lu 等［34］利用單程波算子進行一次波和多次波LSM，應用構造導向正則化條件保障反演的穩定性，提升了墨西哥灣寬方位海洋實際資料的成像質量，在增強照明的同時還減弱了采集腳印影響。③基于單程波算子的全波場成像技術，能夠利用包含表面多次波和層間多次波的全波場偏移成像，利用反射系數和地震速度正演，通過最小化模擬數據和偏移數據間的殘差來更新成像結果［35］。Verschuur［36］對全波場成像方法的機遇和挑戰做了系統的闡述。Davydenko 等［37］提出了基于RTM 框架的全波場成像方法，利用二次源項定義了波場關系，通過反射率對多次散射建模，克服了單程波偏移算子對成像角度的限制。

第四類方法是Marchenko 成像。Marchenko 方法可用于多次波預測和多次波成像。基于逆散射理論，Broggini等［38］提出了地震超越干涉法（BSI），把波場聚焦到一維介質的內部，形成Marchenko 方法的雛形。Wapenaar 等［39］結合BSI 格林函數重構以及基于互相關函數和多維反褶積的成像，將這種成像方法命名為 Marchenko 成像，并給出了完整的理論推導。Slob 等［40］基于 Marchenko 方程，利用層間多次波對地震反射界面成像。Broggini 等［41］通過數據驅動聚焦完成了格林函數重構，并通過多維反褶積研究了層間多次波的數據成像。Singh 等［42］拓展了Marchenko 方程，同時計算了一次波、自由表面多次波和層間多次波的格林函數，實現了基于Marchenko 方程的全波場成像。Singh 等［43］系統研究了自由表面多次波的Marchenko 成像，并指出該方法能夠面向局部目標成像。靳中原等［44］提出了基于低頻信息補償的數據驅動Marchenko 成像，并討論了含有自由表面多次波的地震數據在Marchenko 成像中的應用方法。

本文研究的基于單程波傳播算子的全波場最小二乘偏移方法隸屬上文所述第三類，在DELPHI工業組織的全波場偏移（FWM）技術基礎上開展。首先，利用地層上、下界面反射系數和背景速度，考慮地層的透射和反射效應，并基于惠更斯二次源的思想，推導了基于單程波算子閉循環延拓的全波場正演模擬算子；其次，在二范數意義下求解全波場最小二乘偏移的誤差泛函和梯度項表達式，構建基于反演框架下的全波場LSM；最后，針對透鏡體與加入水層的Marmousi 模型開展偏移測試分析，驗證了本文方法的有效性。多次反射波信息的利用顯著提高了成像品質，與此同時，多輪次反演迭代也壓制了復雜波場產生的相干假象。該方法的研究拓寬了地震成像的手段，尤其在多次波發育的海域地震資料處理中具有實用和推廣價值。

1 基于單程波算子的全波場正演模擬

地球物理方法中，正問題是反問題的研究基礎，單程波是對全波方程近似而分裂出的一種簡化形式。它相對于射線追蹤技術能夠適應速度場更強的橫向變化，相對于全波算子又具有更高的計算效率，并且具有靈活控制波場上、下行傳播方向的特點，本節將闡述如何基于單程波算子進行全波場的正演模擬。

1.1 常規單程波傳播算子

在一次波偏移方法中，多次反射波被認為是噪聲。單程波偏移方法將一次波場激發與接收的過程分解為兩個單程波的傳播過程，即入射波場的下行傳播和反射波場的上行傳播。基于Berkhout 等［15］的“WRW”模型表述方法，從地表起始深度z0傳播至地下深度zm的下行波場在頻率域可以表示為

式中：m=1，2，…，M，其中M表示最深層序號；zM指最大深度；W+(zm，z0)是從z0傳播至zm的下行波傳播算子矩陣；是震源波場，其中下標j代表炮的序號。

從zm傳播至z0的上行波場在頻率域表示為

式中：W-(zm，zn)是從zn傳播至zm的上行傳播算子矩陣；R∪(zn，zn)是zn上界面的反射系數矩陣或算子。式（1）和式（2）分別表示了一次波在地下的激發與接收過程。圖1a是下行波的傳播示意圖，下行波場只有來自上層入射波場。圖1b所示的是上行波場，主要是來自層界面下所有反射波場上行傳播后的疊加（圖1b中實線所示）。按照逐層遞推的方式，傳播矩陣算子W?展開后表示為

圖1 下行波場（a）、上行波場（b）一次波傳播示意圖

圖2 上（a）、下（b）行波場界面處的連續性分析

1.2 考慮透射效應的單程波傳播算子

基于單程波算子的傳統一次波偏移方法存在一個明顯缺陷，即其反射界面處的波場是非連續的。反射界面zn之上的總波場包含了入射波場和入射波場在界面上發生反射后的波場而界面之下的波場只包含了該界面的入射波場在層界面之下的波場包含了上行入射波場和入射波場在界面下反射后的波場而界面上的總波場只有入射波場單程波偏移中的波場延拓特點不符合波場連續性，在波場傳播過程中并未考慮波場的透射效應，因此，雖然基于單程波算子的偏移方法能適應速度場的強烈橫向變化，但在振幅保真上是欠缺的。通過引入波場的透射算子能夠修正透射效應的影響，地層上、下界面透射算子可以表示為

式中：I是單位矩陣；T-、T+是上、下行波透射算子。在界面兩側橫波速度差異較小的條件下，滿足δT+(zn，zn)=R∪(zn，zn)，δT-(zn，zn)=R∩(zn，zn)，將透射算子加權到延拓的波場中，相當于考慮了波場在界面處的連續性，也保證了波場的能量守恒（入射波場的能量等于透射波場和反射波場的能量和），考慮透射效應后的上、下行單程波延拓可以表示為

比較式（6）與式（3）可見，修正后的傳播算子包含波場經過界面的透射效應保證了波場在界面處的連續性。

1.3 考慮透射效應和多次反射的單程波傳播算子

全波場傳播算子不僅考慮了波場傳播過程中的透射，還需要考慮多次反射效應（包括表面相關多次波和層間多次波）。借鑒表面多次波預測方法中的反饋模型思路，通過單程波場在z0與zM間的多輪延拓，能夠模擬出多次反射波場，每一輪的波場延拓能夠正演模擬出一個階次的多次反射波，在z0處注入下行波場能夠模擬出表面相關多次波，在地下每個位置處的波場可以認為是惠更斯二次震源，分別沿上行和下行傳播，從而模擬出層間多次反射波。因此，考慮多次反射波的單程波延拓算子的波場和上、下行波場［34］分別表示為

式中：下行波延拓算子（第二式）與一次波偏移中的下行波延拓算子（式（5）第二式）相比，多了一項基于反射系數R?的反射波場，該項就是產生表面多次波或層間多次波的二次場源。上行波延拓算子（式（7）第一式）與一次波偏移中的上行波延拓算子（式（5）第一式）相比，多了一項zM處的上行波場。將式（7）中的透射傳播算子展開后，下行波場和上行波場可以分別表示為

式中：z-和z+分別代表第m層的上界面和下界面；代表了地下網格點位置處雙向二次源，分別表示為

式（9）如果忽略在zn處的波場轉換，那么有它包含了反射項和透射項，對應關系如圖3所示的地下網格點波場分布。其中，處的波場包含了下行波處的反射波場和上行波從透射至的透射波場處波場包含了上行波在處的反射波場和下行波從透射至的透射波場數值相等、方向相反，保證了波場在zn層網格點的連續性。式（8）就是考慮多次反射波和透射效應的單程波傳播算子，而常規基于一次反射波的單程波偏移方法中，假設了下界面反射系數R∩為零（沒有多次波）和不存在透射算子δT±（忽略了波場在界面的傳播影響）。

圖3 地下網格點波場散射關系

1.4 基于反饋模型逐層遞推的全波場正演模擬

單程波傳播算子的波場延拓方式是以一定的深度步長逐層延拓，如圖4所示。對于離開zm-1處的下行波場包含了處的入射波場和散射波場在下行延拓算子的作用下延拓至處，得到處下行入射波場

圖4 單程波場延拓示意圖

圖5 基于單程波算子的全波場模擬流程

2 全波場最小二乘偏移成像

LSM 是利用模擬數據（或反偏移記錄）與觀測數據之間的匹配程度來解決常規成像結果中振幅失真及采樣不規則等問題。基于1.4節介紹的全波場正演模擬技術，全波場LSM 在迭代偏移過程中考慮了表面和層間等多次反射波對成像的貢獻，可以通過如下的最優化反演問題來描述［32］

式中：ΔP是每一炮全波場正演模擬Pmod與實際記錄數據Pobs間的殘差；f(R)是正則化函數，用以壓制假象，提高成像分辨率。本研究使用的是基于L1范數的柯西約束正則化條件，通過在目標函數（成本函數）中增加L1范數約束，系數化反演成像，便于特征值的提取。該項可以表示為

式中：σ是控制稀疏性的標量因子，通常是反射系數的5%；zn和l分表代表深度方向和橫向方向；ε是控制正則化項影響的權系數項。利用閉循環的方式，通過梯度下降法計算目標函數關于反射系數的梯度項

式中：ΔP+(zm)、ΔP-(zm)分別表示上、下行波場逆向延拓的殘差；上標H 表述矩陣的轉置。所有向下反射波場的殘差逆向向上延拓，可以表示為

圖6 展示了全波場偏移成像中反射系數估計過程，正演模擬波場與逆向延拓的殘差波場相關后形成反射系數的梯度項。上層界面反射系數貢獻主要來自一次反射波和下行多次散射波場，而下界面的反射系數貢獻主要來自層間多次反射波。在構造成像的情況下，選取對角線元素上的梯度并沿頻率求和，也就是傳統的互相關成像條件

圖6 全波場偏移中的反射系數估計

反射系數更新為

式中α和β均是使目標函數最小化的梯度項步長。當橫波變化速度較小時，上、下界面反射系數的關系為R∪=-R∩，因此，只需要求解上界面反射系數即可。

相對于傳統一次波成像方法，全波場最小二乘的計算更加耗時，每增加一個階次的閉循環波場延拓，相當于增加一倍的偏移時間，對于包含N階次多次波的全波場成像，偏移計算是單程波偏移時間的N倍。

3 數值測試

本節通過含透鏡體的層狀模型和含水層的Marmousi模型來驗證基于單程波算子的FWM 成像方法的有效性。透鏡體模型利用全波場正演模擬方法可以實現不同階次的多次波的模擬及分離；對于含水層的Marmousi模型，使用的是空間2 階和時間4階的有限差分算子模擬全波場數據。通過兩個模型數值算例說明，全波場最小二乘偏移方法（FW-LSM）可以有效利用多次波信息進行成像，提高復雜構造地區的成像品質。

3.1 透鏡體模型測試

高速透鏡體與下伏地層間容易產生層間多次反射波場，適用于本文研究的全波場LSM 成像方法測試。速度模型如圖7a 所示，淺部平層是海底，海水速度1500 m/s，水底下的兩套平層間夾有速度為3200 m/s 的高速透鏡體，透鏡體下伏三套平層和一套背斜地層，高速巖丘構造和目標層間會產生較強的層間多次反射波場。圖7b 所示的是密度場，圖7c 是用于計算的反射率模型。該模型橫向4000 m，縱向深1200 m，正演模擬地震子波為主頻20 Hz 的Ricker子波，采樣間隔為4 ms，記錄時間為2 s。圖8a是正演模擬得到的一次反射波，①～⑦箭頭所指為第1 層～第7 層的一次反射。圖8b 和圖8c 分別是含一階表面多次波和含一階及二階表面多次波的正演模擬炮集，圖8d 和圖8e 分別是分離后的一、二階表面相關多次反射波。當含有較強的波阻抗界面時，多次波信號對一次反射波地震信號造成了強烈的干擾。圖9a 是常規一次波LSM，雖然所有阻抗界面都可以成像，但由于多次反射波無法正確歸位導致干涉噪聲影響成像質量，而且透鏡體下的成像照明不足，多次波假象影響嚴重（圖9a中的箭頭所示）。圖9b 是利用本文方法迭代12 次得到的全波場LSM 結果，可以看出， FWLSM 成像通過多次迭代過程能夠抑制干涉假象，與此同時，多次波信息的利用也提高了成像品質，鹽丘下目的層成像照明增強。除此之外，由于充分利用了層間多次波的能量，地震成像的垂向分辨率也有了一定的提高。

圖7 本文使用的速度（a）、密度（b）及反射率（c）模型

圖9 LSM（a）與FW-LSM（b）透鏡體模型偏移效果對比

3.2 含水層的Marmousi 模型測試

透鏡體模型驗證了本文方法對簡單模型的有效性，下面對含水層Marmousi模型測試FW-LSM 技術對復雜模型的適用性。圖10a 為正演模擬速度場，淺層為速度是1500 m/s的水層，模型的寬度和深度分別為4100 m 和1100 m，該模型的空間被離散為821×221 的網格，CDP 間距為5 m，縱向深度步長為5 m。圖10b為用于進行FW-LSM 的背景平滑速度模型，用于計算反射率的密度值被視為均值，故密度參數不對地震波場的傳播產生影響。正演模擬數據震源子波為主頻20 Hz 的Ricker 子波，采樣間隔為4 ms，501個采樣點，模擬信號記錄長度為2 s。炮點從最左邊開始布置，炮間距為50 m，總共84 炮。檢波點均勻布設，間隔為10 m，共411道。圖10c為正演模擬得到的單炮記錄。圖11a為基于真實速度結構得到的反射率模型，利用FW-LSM 成像第一次迭代偏移結果如圖11b所示，可見成像分辨率較低，強阻抗差異界面間的弱反射層不能很好地成像。圖11c 為第10 次迭代后的FW-LSM，成像效果顯著改善，與同樣10 次迭代的一次波LSM（圖11d）相比，整體上成像分辨率有了提升，箭頭所指處成像品質有顯著改善，圓形框內斷層的成像也更加清晰。由此可見，全波場信息的引入在提高整體成像分辨率的同時，也能增加弱反射地層的照明。雖然多次波對信噪比有所影響，但隨著迭代次數的提高，這種影響會逐步減弱。

圖10 Marmousi 模型與正演記錄

圖11 Marmousi 模型偏移

4 結論

反射地震學利用地層反射波信息研究地下結構特征，傳統地震資料處理方法受限于波場一次反射的假設條件，只能利用一次反射波信息成像，豐富的地下多次反射波信息通常被認為是噪聲而進行了壓制。為了能夠合理利用地下全波反射波場信息，本文提出了一種基于單程波傳播算子的全波場成像方法，優化了傳統單程波傳播算子，引入上、下層反射系數以及透射算子影響補償，實現了全波場正演模擬和最小二乘偏移，充分利用了一次波和多次波信息對地下構造成像，在增強照明度的同時提高了垂向分辨率。

透鏡體模型和Marmousi模型測試結果均驗證了本文方法的有效性，該技術也為多次波發育資料（海洋地震資料）的處理提供了新的思路，具有實用性和潛在的應用價值。

當然，受限于資料的品質，對于強面波發育的陸地資料和極低信噪比資料，全波場成像技術現階段并不適用。面向深海的地震資料具有強多次波能量和高信噪比較的特點，是全波場最小二乘偏移技術最具潛力的應用領域。