基于反射波動方程的地震反射數據波形成像

陳 生 昌(浙江大學地球科學學院，浙江杭州 310027)

1 概論

利用地震波可獲取地下地質體的信息。地質體是由塊狀體和層狀體組成，在地震勘探中這些塊狀體和層狀體產生散射數據和反射數據。利用地震數據對這些地質體的反演成像有三個不同層次： ①反演地質體的物性參數，即所謂的地震數據全波形反演[1，2]，如果僅考慮一次波數據且初始模型較合適，則該全波形反演就退化為波形線性反演[3]； ②反演地質體的邊界物性參數(如反射率或反射系數)； ③反演地質體的邊界位置信息，即所謂的地震數據偏移(構造)成像。從由①反演出的地質體物性參數可得到地質體的邊界物性參數和邊界位置信息； 由②反演出的地質體邊界物性參數可得到地質體的邊界位置信息，但不能得到地質體的物性參數分布； 由③反演出的地質體邊界位置信息不能得到地質體的物性參數和邊界物性參數。

地震數據偏移成像是油氣勘探開發中獲取地下地質體三維構造(邊界)圖像的主要方法技術[4-7]，為構造型油氣藏的勘探開發作出了巨大的貢獻。隨著油氣勘探開發目標的日趨復雜化和精細化，如致密砂巖油氣藏、頁巖油氣藏和火成巖油氣藏等的勘探開發，人們希望把地震數據的第三層次的構造成像提升到第一層次的全波形反演或第二層次的邊界物性參數反演，所得到的成像結果不僅能滿足復雜構造成像的要求，還要能反映地下巖性的變化，以滿足復雜油氣藏勘探開發中地震數據巖性處理解釋的需要。

現代地震數據偏移成像方法發源于上世紀70年代初期[8-11]，其目的是利用數學處理的方法把在地表觀測到的地震反射波場送回到產生它的位置上去，以得到反映地下反射體位置信息的偏移成像波場[12-17]。經過四十余年的發展，偏移成像方法已從疊后走向了疊前，從二維走向了三維，從時間偏移成像走向了深度偏移成像，從各向同性介質走向了各向異性介質，從標量聲波方程走向了矢量彈性波方程[18，19]，從不需迭代的常規偏移成像走向了真振幅偏移成像和基于迭代的最小二乘偏移成像[20-27]，并嘗試通過角度域的真振幅偏移或最小二乘偏移實現速度和波阻抗反演[28，29]。

地震反射波的傳播包含入射波傳播、反射、反射波傳播三個過程，在形式上可用Berkhout[30]提出的WRW概念模型表示。當前的偏移方法都是基于該反射波傳播概念模型和Claerbout“時間一致性”成像原理或“爆炸反射面”原理發展起來的。因此偏移成像方法的實現顯式或隱式地包含兩步： 一是基于波動方程的波場外推； 二是基于“時間一致性”成像原理或“爆炸反射面”原理的成像。所謂的顯式偏移成像方法是先進行波場外推然后再進行成像的偏移成像方法，如基于波動方程的單程波方程偏移成像方法[31-37]和雙程波方程逆時偏移成像方法[38-41]；而所謂的隱式偏移成像方法是把波場外推和成像合二為一的偏移成像方法，如基于射線追蹤的Kirchhoff積分偏移成像方法[42-44]和束偏移成像方法[45-47]。

在地震反射數據的逆時偏移中，假定偏移速度模型的運動學特征準確，基于波動方程的震源波場正向(順時)外推是對炮點激發的入射波場的模擬，可以得到旅行時和振幅都正確的入射波場，而基于波動方程的記錄波場逆時外推是對地震反射數據中反射傳播的近似逆(伴隨)傳播，這種伴隨傳播可消除反射傳播中的旅行時，但不能恢復反射傳播中的振幅變化。因此，由記錄波場逆時外推重建出的地下反射面上的反射波場相比于反射面上的真實反射波場，雖然有相同的旅行時，但振幅不同，而且分辨率低。Claerbout的“時間一致性”成像原理是對反射面處入射波場與反射波場間旅行時關系的定量描述，該成像原理在地震波傳播的旅行時方面是滿足波動方程的，但它不是由波動方程推導出來的。由基于“時間一致性”成像原理建立的成像公式，不論是反褶積形式的還是互相關形式的，所得到的偏移成像結果僅是指示地下反射面的存在位置，和反射面上的物性參數不存在數學物理意義上的聯系，因為成像公式所表示的入射波場與反射波場間的關系是不滿足波動方程的。基于Kirchhoff積分的疊前偏移成像方法(其對應的偏移成像計算式隱含了互相關成像公式)所得到的偏移成像結果也僅是指示反射面的存在位置和反射面上的物理參數也不存在數學物理意義上的聯系，因為Kirchhoff積分僅是表示波場的傳播，而沒有考慮波場的反射。當前的地震反射數據偏移成像方法雖然波場外推或旅行時計算和振幅變化的計算是基于波動方程的，但由于其成像公式僅是由地震波的旅行時關系建立的，而不是由反射波傳播所滿足的波動方程推導出來的，所以本文認為當前的偏移成像方法不是一種完全的波動方程方法。

地震數據觀測系統的不規則、不完整和反射面上覆介質物性變化會引起地震波的照明不均勻。逆時偏移中的地震波逆時外推既不能有效地消除地震波照明不均勻對反射波場重建的影響，也不能恢復地震波傳播中的幾何擴散效應，并由此造成偏移成像陰影和降低偏移成像分辨率。此外，互相關成像公式也會造成偏移成像陰影和偏移成像分辨率降低。為了消除逆時外推和互相關成像公式產生的偏移成像陰影和提高偏移成像結果的分辨率，人們提出了基于迭代求解方式的最小二乘偏移方法[48-52]。在當前的地震數據偏移中不區分散射和反射，當前的最小二乘偏移方法以一階Born近似的散射波傳播方程為基礎，得到的是速度相對擾動變化，而不是對反射面位置的成像。此外，散射波方程僅適合地下介質速度的局部小擾動體產生的散射波，而對于地下介質中大尺寸(相對于地震波主波長)的速度異常體(層)所產生的地震反射波用散射波方程進行描述是不合適的。因此，本文認為當前反射數據最小二乘偏移成像方法所基于的正演方程不合適，求得的偏移成像結果也不是對反射面位置的成像。從地震數據偏移和反演的觀點，當前的所謂最小二乘偏移僅是求取速度相對擾動的一種線性反散射方法，而不是一種以反射面位置為目標的偏移方法。本文認為基于當前的以獲取反射面位置信息為目標的偏移成像方法是不能構建所謂的最小二乘偏移成像方法的，因為“時間一致性”原理成像公式所得到的反射面位置是一個幾何變量，而不是一個物性參數變量，無法用于地震反射數據的反偏移，即目前還無法建立一個基于反射面位置信息的反射波傳播方程，適合反射數據偏移成像的反射波動方程是建立反射數據最小二乘偏移的基礎。

為了用波動理論描述地震數據偏移成像方法，Beylkin[53]和Bleistein[54]基于散射波動方程，利用地震反演中的線性反散射理論、擬微分算子理論和Fourier積分算子理論及其微局部分析，建立了一套嚴謹的完全基于散射波動方程的散射偏移成像方法，該方法對于局部弱散射體產生的散射地震數據具有很好的適應性，可實現對散射體位置的成像。對于反射體產生的地震反射數據的偏移成像，Bleistein等[15]仿照基于散射波動方程的散射數據偏移成像方法，利用Kirchhoff近似獲得基于反射系數的反射波場的Kirchhoff積分形式表達式(即Kirchhoff積分形式的反射波場正演計算表達式)，然后以此正演表達式為基礎，提出了改進的Kirchhoff積分偏移成像方法——真振幅Kirchhoff積分偏移成像方法，即對反射系數的真振幅Kirchhoff積分偏移成像方法。但是在Kirchhoff近似中對于反射面上的反射波場與入射波場間的關系表達仍然采用了類似Claerbout“時間一致性”成像公式中的表達，即ur=Rui(ur表示反射波場，R表示所定義界面的反射系數，ui表示入射波場)。為了獲得波動方程(單程波方程和雙程波方程)的真振幅偏移成像方法，一些研究者[55-57]根據上述的散射偏移成像方法和真振幅Kirchhoff積分偏移成像方法對Claerbout“時間一致性”成像公式進行修改，得到所謂的真振幅成像公式，并由此得到的真振幅偏移成像方法具有與真振幅Kirchhoff積分偏移成像方法相類似的真振幅特性。筆者認為僅對成像公式進行修改，對于建立波動方程真振幅偏移成像方法只是治標不治本，而應先建立與地震偏移成像對應的反射數據正演問題，再利用反演的方法構建出真正的反射數據波動方程真振幅偏移成像方法。

根據上述認識，筆者認為無論是當前以反射面位置為目標的反射數據偏移成像方法(包括基于微分方程求解的波動方程方法和基于射線追蹤的Kirchhoff積分方法和束偏移成像方法)還是真振幅偏移成像方法，以及以反射體的速度相對擾動為目標的地震反射數據線性反演方法(當前的最小二乘偏移成像方法)所面臨的最大問題是缺乏描述反射波傳播的反射波動方程，致使不能建立基于反射波動方程的反射數據反演成像公式。在地震數據的反演成像中，如果所建立的反演成像方法不是完全基于描述波場傳播的波動方程的，則該反演成像方法就不能完全利用地震數據的波形信息。由于成像公式的緣故，當前的地震數據偏移成像方法不是一種完全的波動方程方法，不能完全利用地震數據的波形信息。正演問題是進行反演的基礎，為了建立完全基于反射波動方程的反射數據成像方法，需推導控制反射波傳播的反射波動方程，然后在反演理論指導下建立可完全利用地震數據波形信息對反射波動方程中模型參數進行成像的波形成像方法。本文的波形成像包括波形線性反演和波形偏移。在有關波現象和地下介質物理參數空間變化的高頻近似條件下，利用地震波傳播的擾動理論研究描述地震反射波傳播的反射波動方程，推導出基于阻抗相對擾動的反射波動方程和基于反射率變化的反射波動方程。在僅考慮一次反射波的假定下，分別利用基于阻抗相對擾動的一次反射波動方程和基于反射率的一次反射波動方程推導反演阻抗相對擾動的波形線性反演方法和反演反射率的波形偏移方法，因為反射率不僅是一種物性參數，而且其分布還指示了反射面位置，所以把反演反射率的波形線性反演稱為波形偏移。在波形線性反演和波形偏移方法中都包含了利用反射波傳播算子的伴隨算子的近似方法和利用反射波傳播算子的最小二乘逆算子的最小二乘方法。最后把所提出的反射數據波形成像方法應用于合成數據，以驗證方法的正確性和有效性。

2 反射波動方程

為了構建基于反射波傳播方程的反射數據波形線性反演與波形偏移方法，必須構建適合反射數據波形線性反演與波形偏移的反射波動方程，即需要建立基于物性參數相對擾動的反射波動方程和基于反射面物性參數—反射率變化的反射波動方程。散射理論及其散射波動方程是研究地震波傳播的基本方法理論[58]，在有關波現象和地下介質物理性質空間變化的高頻近似條件下，產生散射波的散射體可近似為反射體，相應的散射波退化為反射波。地震勘探中的地震反射波在地下運動的物理過程在直觀上可表述為，地表震源激發的入射波場向下傳播，與反射體作用產生反射波場，反射波場向上傳播到地表被檢波器接收。因此地震反射波的地下運動物理過程由波的反射和傳播組成。

2.1 地震波的傳播與反射

在應力與應變線性近似假定下，震源激發的地震波在地下的運動可用下述非線性微分算子方程描述

L(m)u=f

(1)

式中：L為與地下模型參數m有關的二階偏微分算子，也稱為波動算子；u為表示地震波在地下運動狀態的狀態變量，也稱為地震波場；f為激發地震波的震源函數。式(1)也稱為非齊次波動方程，或輻射方程。如果式(1)的右端項等于零(無源項)，則得到齊次波動方程，即

L(m)u=0

(2)

式(1)、式(2)能描述震源和地下模型m變化產生的各種波，如直達波、散射波、反射波、透射波和轉換波(對于彈性波動方程)等，因此它們也稱為通用波動方程，是地震波場模擬、地震數據全波形反演和波場外推的基礎方程。在描述地震波場運動狀態的通用波動方程(式(1))中，地震反射波的傳播與反射是交織在一起的，這是因為式(1)中的模型參數m不僅控制了波場傳播，也控制了地震波的反射。為了清楚地描述地震反射波的運動狀態，根據對地震反射波傳播的直觀認識，希望有描述反射波場傳播的微分方程和描述反射波場反射的微分方程，也就是希望能在模型m中將傳播與反射分開，得到控制傳播的傳播模型mp和控制反射的反射模型mr，在形式上把m表示為m=mp+mr。

為了得到對傳播模型mp和反射模型mr的認識，將地震波傳播的散射理論應用于式(1)。即令m=mb+δm，其中mb表示背景模型，δm為模型擾動； 相應的波場有u=ub+δu，其中δu為波場擾動，ub表示背景模型下的波場。于是有

L(mb)ub=f

(3)

將上述模型和波場擾動表達代入式(1)，有

L(mb)δu=δL(mb,δm)u

(4)

式中：L(mb)為背景模型下的波動算子； 擾動波場δu也稱為散射波場； δL(mb,δm)為由背景模型和擾動模型綜合產生的擾動算子，也稱為散射算子，有δL(mb,δm)=L(mb+δm)-L(mb)；u為與模型m對應的全波波場。式(4)也稱為散射波動方程。在一階Born近似下，對于一次散射波，可寫為

L(mb)δu=δL(mb,δm)ub

(5)

在散射理論的Born近似中，對于δm不僅要求在數值上相對mb要小很多，即|δm|?mb，而且在空間上要求是局部分布。由于地震波場是一種具有一定頻帶范圍的波動數據，因此δm的空間局部分布要求與地震波場的波長(主頻波長)大小有關。通常認為，如果δm的空間變化特征尺度小于等于主頻波長，則δm可視為局部分布的散射體，產生散射波場。在高頻情況下，相對于短波長，如果δm的空間局部分布要求難以滿足，因而認為δm相對于短波長的地震波在空間上具有一定延續度。因此，在高頻近似條件下，δm可視為反射體，散射波場退化為反射波場，相應的散射波動方程(式(4)、式(5))就退化為反射波動方程。

在光學中，所謂的反射都是鏡面反射，是由鏡面兩側介質的速度變化產生的[59]。在地震勘探中，地下介質主要由塊狀體和層狀體組成。層狀介質產生地震反射波，塊狀體的空間尺寸與地震波的波長關系決定它是產生散射波的散射體還是產生反射波的反射體。本文所討論的地震反射是高頻近似條件下對地震散射的近似，即在地震波場的波長比較短的情況時，δm的空間局部分布對于地震波可視為反射體，δm的空間變化在地震波場的波長尺度下可視為緩慢變化甚至常數，它與入射波場的作用產生反射波場。

根據上述認識，在高頻近似條件下，假定背景模型mb是光滑變化，即式(3)中波場ub不包含反射波場，則地面觀測到的地震反射都是由δm產生。與mb對應的波動算子L(mb)控制和描述了地震波場的傳播，稱之為傳播算子； 與δm有關的擾動算子δL(mb,δm)控制和描述了地震波的反射，稱之為反射算子。因此光滑的背景模型mb就是所謂的地震波傳播模型，δm和mb共同組成了所謂的反射模型。

2.2 基于阻抗相對擾動的反射波動方程

為了體現簡潔和一般性，令波動方程式(1)中的模型參數m為速度v(x)和密度ρ(x)，則式(1)可表述如下

u(x,t;xs)=f(t)δ(x-xs)

(6)

(7)

(8)

=f(t)δ(x-xs)

(9)

如果不考慮密度變化，則聲波動方程式(6)就簡化為標量波動方程，即

(10)

(11)

(12)

(13)

2.3 基于反射率變化的反射波動方程

(14a)

根據反射率的定義，反射率r(x,σ)也有下述定義式，即

(14b)

上述定義的反射率是一個與反射半開角σ有關的變量，是介質聲波阻抗空間變化的反映。

把式(7)、式(8)變換到頻率域，有頻率域的基于聲波阻抗相對擾動的全反射波和一次反射波的聲波反射波動方程，即

(15)

(16)

把反射率定義式(14)代入方程式(15)、式(16)，可得到頻率域的基于反射率變化的全反射波和一次反射波的聲波反射波動方程，即

(17)

(18)

把式(17)、式(18)變換到時間域，有

(19)

(20)

式(19)、式(20)分別稱為時間域的基于反射率變化的全反射波和一次反射波的聲波反射波動方程。

對于標量波動方程式(10)，為了得到基于反射面上反射率變化的反射波動方程，定義地下反射率為速度相對擾動vr(x,σ)沿入射波傳播方向I的方向導數，即

(21)

利用定義式(21)，由基于速度相對擾動的標量反射波動方程式(11)、式(12)，通過與上述相類似的推導，可得到時間域的基于反射率變化的全反射波和一次反射波的標量反射波動方程，即

(22)

(23)

上述的反射波動方程都是在高頻近似條件下推導出來的，反射波動方程中的波場是高頻近似下的局部平面波場。

3 反射數據的波形成像

在高頻近似條件下得到的全反射波和一次反射波的反射波動(式(7)和式(8)、式(11)和式(12)、式(19)和式(20)、式(22)和式(23))是進行全反射波和一次反射波模擬和反演成像的基礎方程。對于一次反射波的模擬，給定震源函數、光滑的背景模型和阻抗(速度)相對擾動或反射面上的反射率，利用推導出的一次反射波動方程(式(8)、式(12)、式(20)和式(23))就可實現對一次反射波場的數值模擬。對于一次反射波的反演成像，給定地震一次反射數據、震源函數和光滑的背景模型，利用推導出的一次反射波動方程(式(8)、式(12)、式(20)和式(23))就可實現對阻抗(速度)相對擾動或反射面上反射率的反演成像。

基于反射波動方程構建的地震反射數據反演成像方法可有效利用反射波形信息，正如基于波動方程構建的地震數據全波形反演方法一般可有效地利用地震數據的波形信息。一次反射波動方程是關于阻抗(速度)相對擾動或反射率的線性方程。根據一次波反射波動方程，利用地震一次反射數據對阻抗(速度)相對擾動或反射率進行成像是線性反演問題求解。將這種完全根據反射波動方程可利用反射數據波形信息對阻抗(速度)相對擾動和反射率進行反演的方法分別稱為地震反射數據的阻抗(速度)相對擾動波形線性反演方法和反射率波形偏移方法。依據波形線性反演中具體求解方式的不同，對于波形線性反演又分為基于反射波傳播算子的伴隨算子的近似反演和基于反射波傳播算子的最小二乘逆算子的最小二乘反演；同樣，對于波形偏移也分為基于反射波傳播算子的伴隨算子的波形偏移和基于反射波傳播算子的最小二乘逆算子的最小二乘波形偏移。由于反射波傳播算子的時空變化特征，反射波傳播算子不是一個酉算子，其伴隨算子不等于其逆算子，因此基于伴隨算子的波形線性反演與波形偏移實質是對阻抗(速度)相對擾動和反射率的近似反演。不同于當前基于Claerbout成像原理的地震反射數據偏移方法，地震反射數據波形偏移是一種完全波動方程方法，它不僅可實現對反射面位置的成像，也可對反射面的物性參數反射率進行近似估計，最小二乘波形偏移是波形偏移的發展，是對反射面反射率的最小二乘反演。

3.1 阻抗相對擾動的近似反演與最小二乘反演

對于式(8)所表示的基于阻抗相對擾動的一次反射波動方程，借助背景模型下的Green函數可寫為下述積分形式

(24)

式中：“*”代表時間褶積；g(x,t;y)為背景模型下的Green函數，有

=δ(t)δ(x-y)

(25)

把式(24)變換到頻率域，有

(26)

令

(27)

其中

(28)

且F(ω)為震源函數的頻譜。

把式(27)代入式(26)，并根據陳生昌等[60]提出的時間二階積分波場的全波形反演方法，式(26)可寫為

(29)

對于線性積分方程(式(29))可寫為下述矩陣方程形式

u=Wm

(30)

m=W-gu=(W*W)-1W*u

(31)

式中：W-g代表W的最小二乘廣義逆矩陣；W*為W的伴隨矩陣。式(31)的運算在實質上是對波場數據u的最小二乘逆傳播。

(32)

式中Re表示取實部。

考慮到式(31)計算的穩定性和逆矩陣計算的復雜度，可用伴隨矩陣W*代替式(31)中的最小二乘廣義逆矩陣W-g，有

ma=W*u

(33)

式(33)的運算在實質上是對波場數據u的伴隨傳播(逆時傳播)。

再考慮式(32)中波場相除運算的不穩定性問題，并結合式(33)，把式(32)改寫為

(34)

式中上標“*”代表復共軛運算。用式(34)代替式(32)，雖然可以保證計算的穩定性和減小計算的復雜度，但會造成反演結果的保真性變差、分辨率降低。

式(28)、式(33)和式(34)所表示的運算構成了對阻抗相對變化Ir(x,σ)的近似反演。把對Ir(x,σ)的近似反演運算轉換到時間域，就有：

(1)地下入射波場的計算(對應式(28)所表示的運算)

=f(t)δ(x-xs)

(35)

(2)地下時間二次積分反射波場的伴隨(逆時)重建(對應式(33)所表示的運算)

(36)

(3)對阻抗相對變化Ir(x,σ)的近似反演(對應式(34)所表示的運算)

(37)

(38)

利用基于速度相對擾動的標量波反射波動方程(式(12))和相類似的公式推導，可得到對速度相對變化vr(x,σ)的近似反演，其相應的計算公式為如下。

(1)地下入射波場的計算

(39)

(2)地下時間二次積分反射波場的伴隨(逆時)重建

(40)

(3)對速度相對變化vr(x,σ)的近似反演

(41a)

或寫成

(41b)

式(41a)、式(41b)得到的速度相對變化近似反演結果可以用來估計速度擾動。

由于式(33)、式(34)中的近似，致使上述的波形偏移方法會受到地震數據觀測系統、地下介質變化、地震波幾何擴散以及近似反演公式等對成像結果的分辨率、信噪比和保真性產生的不利影響，所得到的結果對是對阻抗(速度)相對變化的近似估計。

顯然，式(28)、式(31)和式(32)雖然在理論上給出了利用最小二乘反演阻抗相對變化的計算公式，可以得到一個高分辨率、高信噪比和高保真性的阻抗(速度)相對變化估計，但存在一些數值計算方面的難題： 一是式(31)右端中的逆矩陣(W*W)-1的大規模計算問題和它的不穩定性問題； 二是式(32)右端中的波場相除運算也會出現計算的不穩定問題。

(42)

Δuik=ui-Wmk

(43)

Δmik=W*Δuik

(44)

(45)

利用迭代所得到的有關阻抗相對變化的反演，即稱之為基于波場最小二乘逆傳播的反演，也稱為阻抗相對擾動的最小二乘反演。

(46)

Δuik=ui-Wmk

(47)

Δmik=W*Δuik

(48)

(49)

利用迭代所得到的有關速度相對變化的反演，稱之為基于波場最小二乘逆傳播的反演，也稱為速度相對擾動的最小二乘反演。

通過迭代方式的最小二乘反演得到的相對波阻抗(速度)擾動相對于近似反演得到的相對波阻抗(速度)擾動具有保真性好、分辨率高，但相應的計算量會大量增加。

3.2 反射率的波形偏移與最小二乘波形偏移

根據基于反射率變化的一次波反射波動方程(式(20))，利用地震一次反射數據，不僅可對反射面位置進行成像，即實現對反射面位置的偏移成像，還能對反射面的物性參數反射率進行估計。對于式(20)所表示的一次反射波動方程借助Green函數可寫為積分形式，有

(50)

采用與地震反射數據的阻抗相對擾動近似反演相類似的公式推導，得到下列基于波場伴隨傳播的對反射率r(x,σ)進行近似估計的波形偏移時間域計算公式。

(1)地下入射波場的計算

=f(t)δ(x-xs)

(51)

(2)地下反射波場的伴隨(逆時)重建

(52)

(3)對反射率r(x,σ)的近似估計

(53)

式(52)中的dr(x,t;xs)為共炮道集地震反射數據；式(53)中的波場ur(x,t;xs)的時間一階導數對獲得高分辨率的r(x,σ)波形偏移結果有幫助。

根據標量波一次反射波動方程(式(23))，利用波場的伴隨傳播運算和相應的公式推導，得到下列對r(x,σ)進行近似估計的波形偏移時間域計算公式。

(1)地下入射波場的計算

(54)

(2)地下反射波場的伴隨(逆時)重建

(55)

(3)對反射率r(x,σ)的近似估計

(56)

上述導出的波形偏移成像方法不同于當前基于Claerbout成像原理的波動方程逆時偏移成像方法，它是一種基于反射波動方程的反演方法，不需要應用所謂的偏移成像原理及其成像公式，能夠在波動方程意義下真實地應用地震數據的波形信息，是一種真正的波動方程偏移成像方法。

(57)

Δuk=u-Wmk

(58)

Δmk=W*Δuk

(59)

rk+1(x,σ)=rk(x,σ)+αρb(x)vb(x)×

(60)

上述迭代所得到的有關反射率的波形偏移，可稱之為基于波場最小二乘逆傳播的波形偏移，也稱為反射率的最小二乘波形偏移。

(61)

Δuk=u-Wmk

(62)

Δmk=W*Δuk

(63)

rk+1(x,σ)=rk(x,σ)+αvb(x)×

(64)

迭代方式的最小二乘波形偏移得到的反射率相對于波形偏移得到的反射率具有保真性好、分辨率高的優點，但相應地計算量會大量增加。當前的最小二乘逆時偏移方法得到的是阻抗(速度)的相對擾動，而不是逆時偏移方法所得到的反射面位置，對于缺少低頻信息和高頻信息的帶限地震數據，最小二乘逆時偏移得到的是阻抗(速度)的相對擾動不能清晰地刻畫反射面的位置，即當前的最小二乘逆時偏移結果與逆時偏移結果是不同的。最小二乘逆時偏移不是對反射面位置的成像，因此嚴格地說最小二乘逆時偏移不是偏移，而本文提出的最小二乘波形偏移與波形偏移都是對反射面反射率的成像。

根據地下非均勻體的體積大小和波阻抗變化快慢與地震波長之間的相對關系，把非均勻體劃分為產生散射波的散射體和產生反射波的反射體，然后在高頻近似條件下推導出兩種不同參數的反射波方程，即基于波阻抗相對擾動的反射波方程和基于反射率的反射波方程，并把它們分別作為反射數據的正演方程。本文認為對于散射波和反射波應有不同的散射波方程和反射波方程，相應地對于散射數據和反射數據也應有不同的散射數據波形成像方法和反射數據波形成像方法[63]，但當前的偏移成像中不區分散射和反射。為了更精細地實現對地下目標的成像，趨向認為在偏移成像中對于散射波和反射波應有不同的方法。本文波形成像方法的目標是產生反射數據的反射體和反射面。如果波形成像的目標是產生散射數據的散射體，就需要應用基于散射波方程的散射數據波形成像方法。本文認為在地震數據的成像中，為了得到高分辨、高保真的成像結果，對于散射數據和反射數據應該區別對待，使用不同的方法。根據相關經驗，基于反射波方程的反射數據波形成像方法對于散射體產生的散射數據也能進行成像，只是成像效果不如基于散射波方程的散射數據波形成像方法的。同樣，基于散射波方程的散射數據波形成像方法對于反射體產生的反射數據也能進行成像，只是成像效果不如基于反射波方程的反射數據波形成像方法的。

4 數值試驗

為了驗證上文所提出的地震反射數據波形成像方法的正確性和有效性，采用不考慮密度變化的標量波合成地震反射數據進行了速度相對擾動的近似反演、反射率的波形偏移和最小二乘波形偏移試驗。

4.1 速度相對擾動的近似反演試驗

為了驗證所提出的速度相對擾動近似反演方法的正確性和有效性，本文構建了如圖1a所示的字母模型，該模型包含2層，上部低速層速度為2500m/s，下部高速度層速度為3000m/s，在上部的低速層中含有6個速度為2750m/s的高速字母體。利用標量波動方程(式(10))的有限差分方法進行合成地震數據模擬。在本試驗總共模擬了501炮，炮間距為40m，為中間放炮兩邊接收，每炮1001道，道間距為20m，時間采樣間隔是4ms，時間采樣點數為1250，震源為主頻20Hz的Ricker子波。這些字母體相對于主頻波長可視為反射體，因此合成得到的地震數據是反射波數據。圖1b為用于近似反演試驗的光滑速度模型，它是通過對圖1a的速度模型進行光滑處理得到的。圖1c為圖1a所示的真實速度模型與圖1b所示的光滑速度模型之間的速度擾動。圖1d為利用本文提出的速度相對擾動近似反演方法得到的速度相對擾動近似反演結果。對比圖1c、圖1d可看出，近似反演得到的速度擾動結果相對于理論速度擾動雖然存在保真性和分辨率方面的不足，但是反演的速度擾動與理論速度擾動有相同的變化形態。該試驗結果驗證了本文所提出的速度相對擾動近似反演方法的正確性和有效性。

圖1 速度相對擾動的近似反演效果(a)字母模型的速度模型； (b)光滑的背景速度模型； (c)圖a與圖b之間的速度擾動； (d)近似反演得到的速度擾動

4.2 反射率的波形偏移試驗

為了驗證所提出的反射率波形偏移方法的正確性和有效性，對兩套合成數據進行了波形偏移試驗，第1個是上述字母模型數據的試驗，第2個是Sigbee2A模型數據的試驗。

4.2.1 字母模型數據的試驗

利用上節得到的字母模型的反射數據和圖1b所示的光滑背景速度模型，分別應用常規逆時偏移成像方法和本文提出的波形偏移成像方法對字母模型進行偏移成像。圖2a為常規逆時偏移的結果，圖2b為波形偏移的結果。圖3為圖1c所示字母模型真實速度與背景速度之間速度擾動的垂向反射率。對比圖2a、圖2b和圖1a模型可看出，逆時偏移方法和本文的波形偏移方法雖然都能對字母體邊界和速度分界面進行很好的成像，但逆時偏移得到的結果有90°的相位差，即分界面處于波峰與波谷的過渡帶上，圖2a逆時偏移結果的分辨率也不如圖2b波形偏移結果的高。對比圖2b和圖3可看出，波形偏移方法可以很好地對字母體邊界和速度分界面上的反射率進行成像，波形偏移得到的反射率與圖3中的垂向理論反射率有相同的極性和相位。逆時偏移和波形偏移都在淺部反射面上方出現了低頻噪聲，試驗中雖然用Laplace濾波方法進行了壓制，但還有較明顯殘余。由于波形成像方法與逆時偏移相比多了對震源波場的時間一階導數，因此其低頻噪聲相對于逆時偏移更強，為了減弱低頻噪聲對波形成像結果多做一次Laplace濾波，所以造成成像能量的相對減弱，使直立90°邊界成像看起來不如逆時偏移的結果。深部出現的噪聲主要是用于逆時偏移和波形偏移的光滑速度模型與真實速度模型的偏差引起的。由于伴隨傳播和反射率近似估計公式(式(56))，致使反射率波形偏移結果的分辨率和保真性受到不利影響。

圖2 常規逆時偏移(a)和波形偏移(b)結果

圖3 模型速度擾動的理論垂向反射率

4.2.2 Sigbee2A模型數據的試驗

Sigbee2A模型含有高速鹽體、精細的反射面。圖4a是用于偏移成像試驗的背景速度模型。模型尺寸為3201×1201個網格點，水平和垂直方向的網格間距均為7.62m，該模擬數據的觀測系統為海上拖纜觀測系統。該地震數據含有500炮，炮間距為45.72m，拖纜長度為7932.42m，檢波點間距為22.86m，每一炮的道數不均，最少為57道，最多為348道，時間采樣間隔為8ms，時間采樣點數為1500。

分別應用常規逆時偏移成像方法和本文提出的波形偏移成像方法對Sigbee2A模型數據進行偏移成像。圖4b和圖4c分別為常規逆時偏移成像結果和波形偏移成像結果。對比圖4b和圖4c可見，雖然兩種偏移成像方法都能對Sigbee2A模型進行很好的成像，但本文提出的波形偏移的結果相對于逆時偏移的結果有更高的分辨率，因為波形偏移結果中反射界面的同相軸相對更細。本文的方法雖然是針對反射體(面)產生的反射數據進行成像，但對于鹽體邊界上的角點和模型中的點散射體產生的繞射波和散射波也能進行很好的成像，而且所得到的成像結果也好于常規的逆時偏移。為了獲得更精細的成像結果，建議若成像目標是產生反射數據的反射體(面)，就采用基于反射波動方程建立的反射數據波形成像方法； 若成像目標是產生散射數據的散射體，就采用基于反射波動方程建立的反射數據波形成像方法。

圖4 針對Sigbee2A背景速度模型(a)的常規逆時偏移(b)和波形偏移(c)結果

圖5 部分Sigsbee2a模型的波形偏移結果、最小二乘波形偏移結果和理論垂向反射率

(a)部分Sigsbee2a模型的速度分布； (b)用于偏移的光滑速度模型； (c)波形偏移結果； (d)35次迭代的最小二乘波形偏移結果； (e)模型速度擾動的理論垂向反射率

上述兩個試驗的結果驗證了本文所提出的波形偏移成像方法的正確性和有效性，以及該方法相對于常規逆時偏移成像方法的優勢。

4.3 反射率的最小二乘波形偏移試驗

為了驗證反射率最小二乘波形偏移方法的正確性和有效性，截取Sigsbee2a速度模型的一部分進行試驗，如圖5a所示。利用標量波動方程式(10)的有限差分方法進行合成地震數據模擬。在本試驗總共模擬了300炮，炮點范圍是0～6km，炮間距為20m，為中間放炮兩邊接收，每炮601道，道間距為10m，時間采樣間隔是0.5ms，時間采樣點數為3000，震源為主頻18Hz的Ricker子波。圖5b為用于最小二乘波形偏移試驗的光滑速度模型，它同樣是通過對圖5a的速度模型進行光滑處理得到的。

在試驗中，首先應用提出的反射波動方程波形偏移方法對合成的炮道集地震數據進行波形偏移成像，圖5c為利用波形偏移方法得到的偏移成像結果。由圖5c可看出，波形偏移方法對該模型可獲得很好的偏移成像結果。圖5d為應用本文提出的反射波動方程最小二乘波形偏移成像方法對部分Sigsbee2a模型的偏移成像結果，共進行了35次迭代。對比圖5c與圖5d可看出，最小二乘波形偏移結果的分辨率相對于波形偏移結果的有明顯地提高，而且層位面與斷層面的成像更清晰。與圖5e所示的部分Sigsbee2a模型的垂向理論反射率對比可看出圖5d的最小二乘波形偏移結果的振幅保真性相對圖5c的波形偏移結果也有明顯的提高。

上述模型數據的數值試驗結果驗證了本文所提出的最小二乘波形偏移成像方法的正確性和有效性。由于缺乏實際數據和研究條件的限制，本文所提出的方法沒有進行實際數據的測試。根據前文的計算公式， 認為本文所提出的方法具有很強的實用性，因為本文所提出的方法計算過程與當前的逆時偏移和最小二乘逆時偏移相同，所不同的只是成像公式。

5 結論

根據地震數據所滿足的波動方程進行反演成像的求解可充分利用地震數據的波形信息。基于阻抗相對擾動的反射波動方程和基于反射率的反射波動方程分別是阻抗相對擾動反演和反射率波形偏移的正演問題。地震一次反射數據的阻抗相對擾動反演和反射率的波形偏移都屬于波形線性反演。利用反射波傳播算子的伴隨算子所構建的近似反演和波形偏移方法雖然可對阻抗相對擾動和反射率進行成像，但保幅性、分辨率和信噪比存在不足；利用反射波傳播算子的最小二乘算子所構建的最小二乘反演和最小二乘波形偏移方法雖然計算量增大(至少增加一個量級的計算量)，但成像結果的保幅性、分辨率和信噪比會有明顯提高。本文所提出的波形偏移成像方法是一種真正的波動方程方法，不需要應用成像原理及其成像公式。本文的波形成像方法不僅可以實現地下介質參數構造成像，還能進行物性成像，是對當前地震數據偏移成像方法的升級換代。合成數據的數值試驗結果驗證了本文所提出的波形成像方法的正確性和有效性。本文所提出的波形偏移方法與當前的逆時偏移方法的計算過程一致，不會增加計算量，適用條件也一樣。本文從數學物理正反演的角度對地震反射數據的偏移成像進行了重新定義和推導，使得到的偏移成像結果具有明確的數學物理意義，可為后續偏移成像結果的進一步應用(如基于角度域成像道集的AVA分析和巖性反演)提供基礎。

參考文獻

[1] Lailly P．The seismic inverse problem as a sequence of before stack migrations．Conference Expanded Abstracts on Inverse Scattering，Theory and Application，Society for Industrial and Applied Mathematics，1983，206-220.

[2] Tarantola A．Inversion of seismic reflection data in the acoustic approximation．Geophysics，1984，49(8)：1259-1266.

[3] Tarantola A．Linearized inversion of seismic reflection data.Geophysical Prospecting,1984,32(4)：998-1015.

[4] Gray S，Etgen J，Dellinger J et al．Seismic migration problems and solutions．Geophysics，2001，66(5)：1622-1640.

[5] Yilmaz O．Seismic Data Analysis．SEG，2001.

[6] Etgen J，Gray S，Zhang Y．An overview of depth imaging in exploration geophysics．Geophysics，2009，74(6)：WCA5-WCA17.

[7] Leveille J，Jones I，Zhou Z et al．Subsalt imaging for exploration，production，and development：A review．Geophysics，2011，76(5)：WB3-WB20.

[8] Claerbout J．Toward a unified theory of reflector mapping．Geophysics，1971，36(3)：467-481.

[9] Claerbout J．Fundamentals of Geophysical Data Processing．McGraw-Hill Book Co Inc，1976.

[10] Claerbout J．Imaging of the Earth’s Interior．Blackwell Scientific Publication，1985.

[11] Loewenthal D，Lu L，Roberson R et al．The wave equation applied to migration．Geophysical Prospecting，1976，24(2)：380-399.

[12] 馬在田．地震成像技術：有限差分偏移．北京：石油工業出版社，1989.

[13] Robein E．Seismic imaging：A review of the techniques，their principles，merits and limitations．EAGE，2010.

[14] Stolt R，Benson A．Seismic Migration：Theory and Practice．Geophysical Press，London，1986.

[15] Bleistein N，Cohen J，Stockwell J．Mathematics of Multi-dimensional Seismic Imaging，Migration，and Inversion． Springer， New York，2001.

[16] 李振春等編著．地震疊前成像理論與方法．山東東營：中國石油大學出版社，2011.

[17] 李振春．地震偏移成像技術研究現狀與發展趨勢．石油地球物理勘探，2014，49(1)：1-21．

Li Zhenchun．Research status and development trends for seismic migration technology．OGP，2014，49(1)：1-21.

[18] Du Q， Zhu Y， Ba J．Polarity reversal correction for elastic reverse time migration．Geophysics，2012，77(2)：S31-S41.

[19] Chang W，McMechan G．A 3-D elastic prestack reverse-time depth migration．Geophysics，1994，59(4)：597-609.

[20] Zhang Y，Duan L and Xie Y．A stable and practical implementation of least-squares reverse time migration．Geophysics，2014，80(1)：V23-V31.

[21] Yao G，Jakubowicz H．Least-squares reverse time migration．SEG Technical Program Expanded Abstracts，2012，31：1-5.

[22] Yao G，Jakubowicz H．Non-linear least-squares reverse time migration．SEG Technical Program Expanded Abstracts，2012，31：1-5.

[23] Dong S，Cai J，Guo M et al．Least-squares reverse time migration：towards true amplitude imaging and impoving the resolution．SEG Technical Program Expanded Abstracts，2012，31：1-5.

[24] Gray S．True-amplitude seismic migration：A comparison of three approaches．Geophysics，1997，62(3)：926-936.

[25] 郭振波，李振春．最小平方逆時偏移真振幅成像．石油地球物理勘探，2014，49(1)：113-120．

Guo Zhenbo，Li Zhenchun．True amplitude imaging based on least-squares reverse time migration．OGP，2014，49(1)：113-120.

[26] 黃建平，孫隕松，李振春等．一種基于分頻編碼的最小二乘裂步偏移方法．石油地球物理勘探，2014，49(4)：702-707.

Huang Jianping，Sun Yunsong，Li Zhenchun et al．Least-squares split step migration based on frequency division encoding．OGP，2014，49(4)：702-709.

[27] 周華敏，陳生昌，任浩然等．基于照明補償的單程波最小二乘偏移．地球物理學報，2014，57(8)：2644-2655．

Zhou Huamin，Chen Shengchang，Ren Haoran et al．One way wave equation least squares migration based on illumination compensation．Chinese Journal of Geophysics，2014，57(8)：2644-2655.

[28] Zhang Y，Ratcliffe A，Duan L．Velocity and impedance inversion based on true amplitude reverse time migration．SEG Technical Program Expanded Abstracts， 2013，32： 949-953.

[29] Duan L，Zhang Y and Peng L．Band-limited impedance perturbation inversion using cross-correlative least-squares RTM．SEG Technical Program Expanded Abstracts，2014，33：3720-3725.

[30] Berkhout A．Seismic Migration：Imaging of Acoustic Energy by Wavefield Extrapolation，Part A：Theoretical Aspects(2nd Ed).Elsevier，1982.

[31] Gazdag J．Wave equation migration with the phase-shift method．Geophysics，1978， 43(7)：1342-1351.

[32] Berkhout A．Steep-dip finite-difference migration.Geophysics，1979，44(1)：196-213.

[33] Huang L，Wu R．Prestack depth migration with acoustic screen propagators．SEG Technical Program Expanded Abstracts，1996，15：415-418.

[34] Ristow D，Ruhl T．Fourier finite-difference migration．Geophysics，1994，59(12)：1882-1893.

[35] Chen J，Liu H．Two kinds of separable approxima-tions for the one-way wave operator．Geophysics，2006，71(1)：T1-T5.

[36] 陳生昌，馬在田，Wu Rushan．波動方程偏移成像陰影的照明補償．地球物理學報，2007，50(3)：844-850．

Chen Shengchang，Ma Zaitian，Wu Rushan．Illumination compensation for wave equation migration shadow．Chinese Journal of Geophysics，2007，50(3)：844-850.

[37] 劉少勇，韓冰凱，顧漢明等．帶限局部平面波傳播算子與偏移方法．石油地球物理勘探，2017，52(5)：948-955．

Liu Shaoyong，Han Bingkai，Gu Hanming et al．Band-limited local plane wave propagator and migration．OGP，2017，52(5)：948-955.

[38] Hemon C．Wave equations and models．Geophysical Prospecting，1978，26(4)：790-821.

[39] McMechan G．Migration by extrapolation of time-dependent boundary values．Geophysical Prospecting，1983，31(3)：413-420.

[40] Baysal E，Kosloff D and Sherwood J．Reverse time migration．Geophysics，1983，48(11)：1514-1524.

[41] Whitmore D．Iterative depth migration by backward time propagation．SEG Technical Program Expanded Abstracts，1983，2：382-385.

[42] French W．Two-dimensional and three-dimensional migration of model-experiment reflection．Geophy-sics，1974，39(1)：265-277.

[43] Schneider W．Integral formulation for migration in two and three dimensions．Geophysics，1978，43(1)：49-76.

[44] Sun J．Limited aperture migration．Geophysics，2000，65(2)：584-595.

[45] Cerveny V，Popov M and Psencilk I．Computation of wave fields in inhomogeneous media Gaussian beam approach．Geophysical Journal of the Royal Astronomical Society，1982，70(1)：109-128.

[46] Costa C，Raz S and Kosloff D．Gaussian beam migration．SEG Technical Program Expanded Abstracts，1989，8：1169-1171.

[47] Hill N．Prestack Gaussian-beam depth migration．Geo-physics，2001，66(4)：1240-1250.

[48] Tang Y，Biondi B．Least-squares migration/inversion of blended data．SEG Technical Program Expanded Abstracts，2009，28：2859-2863.

[49] Ren H，Wu R and Wang H．Wave equation least square imaging using the local angular Hessian for amplitude correction．Geophysical Prospecting，2011，59(4)：651-661.

[50] 黃建平，李振春，孔雪等．碳酸鹽巖裂縫型儲層最小二乘偏移成像方法研究．地球物理學報，2013，56(5)：1716-1725．

Huang Jianping，Li Zhenchun，Kong Xue et al．A study of least squares migration method for fractured type carbonate reservoir．Chinese Journal of Geophy-sics，2013，56(5)：1716-1725.

[51] 李慶洋，黃建平，李振春等．優化的多震源最小二乘逆時偏移．石油地球物理勘探，2016，51(2)：334-341．

Li Qingyang，Huang Jianping，Li Zhenchun et al．Optimized multi-source least-squares reverse time migration．OGP，2016，51(2)：334-341.

[52] 李娜．基于Huber范數的多震源最小二乘逆時偏移．石油地球物理勘探，2017，52(5)：941-947．

Li Na．Multi-source least-squares reverse time migration based on Huber norm．OGP，2017，52(5)：941-947.

[53] Beylkin G．Imaging of discontinuities in the inverse scattering problem by the inversion of a causal Radon transform．Journal Mathematics Physics，1985，26(1)：99-108.

[54] Bleistein N．On the imaging of reflectors in the earth．Geophysics，1987，52(7)：931-942.

[55] Zhang Y，Zhang G，Bleistein N．Theory of true amplitude one-way wave equation and true amplitude common-shot migration．SEG Technical Program Expanded Abstracts，2003，22：925-928.

[56] Xu S，Zhang Y，Tang B．3D angle gathers from reverse time migration．Geophysics，2011，76(2)：S77-S92.

[57] Duprat V，Baina R．Prestack true amplitude imaging condition．SEG Technical Program Expanded Abstracts，2015，34：4075-4079.

[58] Devaney J．Mathematical Foundations of Imaging，Tomography and Wavefield Inversion．Cambridge University Press，2012.

[59] Goodman J．Introduction to Fourier Optics(3rd Ed).Roberts and Company Publishers Inc，2005.

[60] 陳生昌，陳國新．時間二階積分波場的全波形反演．地球物理學報，2016，59(10)：3765-3776．

Chen Shengchang，Chen Guoxin．Full waveform inversion of the second-order time integral wavefield．Chinese Journal of Geophysics，2016，59(10)：3765-3776.

[61] Sava P， Fomel S．Angle-domain common-image ga-thers by wavefield continuation methods．Geophy-sics，2003，68(3)：1065-1074.

[62] Kirsch A．An Introduction to the Mathematical Theory of Inverse Problems(2nd Ed).Springer，2011.

[63] 陳生昌，周華敏．再論地震數據偏移成像．地球物理學報，2016，59(2)：643-65．

Chen Shengchang，Zhou Huamin．Re-exploration to migration of seismic data．Chinese Journal of Geophysics，2016，59(2)：643-654．