基于數據自適應加權的疊前深度偏移成像方法

吳成梁,王華忠,胡江濤,馬建波

(1.波現象與智能反演成像研究組(WPI),同濟大學海洋與地球科學學院,上海200092;2.油氣藏與開發工程國家重點實驗室,成都理工大學,四川成都610059;3.中國石油化工股份有限公司中原油田分公司物探研究院,河南濮陽457000)

隨著“兩寬一高”地震數據采集技術的發展和進步,油氣地震勘探逐步進入了高分辨率、高保真的反演成像階段,同時對偏移成像技術要求也越來越高。Bayes估計理論提供了地震波反演成像的基本思想。無論是去噪、偏移成像,還是層析建模等方法技術,都可在(也應該在)Bayes反演框架下討論,因為地震數據總受隨機噪聲的干擾,且實測數據總是不完全。

根據Bayes估計理論,假設給定的正演算子能夠預測觀測數據,在預測噪聲和反演解具有先驗統計分布特征(一般為高斯分布或廣義高斯分布)基礎上,定義一定范數意義下的預測誤差,尋找預測誤差最小的反演解。根據統計分布特征,一般數據誤差采用L2范數,模型約束采用L1或Cauchy范數等。若預測噪聲為高斯分布,此時的Bayes估計可以在最小二乘意義下實現;進一步地,根據正算子的線性與否,選擇用線性最小二乘方法或非線性最小二乘方法求解反演問題。

反問題的本質是Bayes框架下的后驗概率密度最大化的決策或估計問題。更具實際計算可行性的使代價函數最小的優化問題的解是反問題的本質。再進一步可落實到算子投影上,TARANTOLA[1]定義了參數空間及其對偶空間、數據空間及其對偶空間,原空間與對偶空間由協方差算子聯系在一起,由伴隨算子實現數據空間向模型空間的映射,得到模型參數的估計結果。從原空間到對偶空間的投影過程為數據和模型的特征表達過程。模型協方差算子包含地下介質的結構特征信息,體現了模型不同點之間的聯系。數據協方差算子包含了觀測數據中相干信號的統計特征,高維數據的投影過程等價于數據的特征提取過程,通過對數據進行預處理,可以提取數據中的相關特征。然而在實際應用中很難得到關于數據和模型的協方差(逆)算子,也沒有具體的解析公式,如何合理利用數據和模型的協方差(逆)算子成為了研究的重點。

在Bayes反演框架下,TARANTOLA[1]全面深入地分析了地震波全波形反演(FWI)的理論問題。理論上,FWI能夠獲得地下全波數帶參數場分布。然而經典的FWI是強非線性反演問題,實際地震數據的地震波FWI受多種因素(比如:缺失低頻長偏移距數據、子波未知、正演算子無法模擬所有波現象、數據噪聲類型未知等)的制約,因而會出現收斂到局部極值或者不收斂的情況。解決此類問題的基本辦法是將這個強非線性問題轉化為一系列線性問題或一個較凸的問題求解[2]。最小二乘逆時偏移(LS_RTM)[3-5]是FWI中典型的線性化子問題,是在背景速度固定的情況下,根據模擬數據與觀測數據的波形差異估計地下介質的高波數擾動。從反演的角度理解,常規的疊前深度偏移成像(PSDM)可以認為是最小二乘偏移成像的第1次迭代結果。在FWI和LS_RTM中,學者們對模型參數的正則化[6-7]進行了深入探討,提出了很多有效的方法,促進了成像質量的提高。但是,他們對數據正則化的關注不充分,如何將數據質量的信息融合在成像過程中以進一步提高成像質量,是一個值得重視的問題。

本文首先從概率論的角度討論了Bayes框架下的地震波反演成像問題,說明了觀測噪聲為高斯白噪的情況下,Bayes估計可以在最小二乘意義下實現;分析了吉洪諾夫(Tikhonov)正則化與Bayes估計的等價性;重點分析了連接原空間與對偶空間的協方差算子的作用,證明了數據特征表達和模型特征表達的必要性;從加權最小二乘誤差泛函角度分析了數據協方差(逆)算子的作用,初步探討了數據協方差(逆)算子在偏移成像中的應用;提出了采用傾角掃描和動態時間規整算法來估計數據加權系數;最后,理論和實際數據的數值實驗結果表明數據協方差(逆)算子能夠有效提高成像質量。

1 Bayes框架下的地震波反演成像

地震波反演成像建立在Bayes框架下估計理論基礎之上。整個地震數據處理與反演成像過程(去噪、偏移成像、速度分析、層析建模等)都可認為在此框架下實現[8-9]。

實際地震數據是一個隨機過程,實測到的地震數據則認為是該隨機過程的一次具體實現。因此,可以將觀測的地震數據視為一組隨機信號,其具有均值、方差等統計特性,滿足一定的概率分布。同樣地,反演的地下模型參數也被認為是隨機的。假設有模型空間M和數據空間D,在其定義的聯合空間M×D上,有同時反映觀測數據d以及模型參數m信息的聯合先驗概率密度ρ(d,m),同樣有反映聯系觀測數據與模型參數物理規律的概率密度Θ(d,m)。利用Bayes公式,可得出聯合空間M×D上的后驗概率密度σ(d,m):

(1)

其中,k為歸一化常數,μ(d,m)為聯合空間中關于觀測數據與模型參數的均勻概率密度函數。均勻概率密度函數與坐標系統有關,表示特定數據和模型參數化下(選定的坐標系統)概率空間中單位區域內的概率密度函數,代表了對模型空間與數據空間的參數化認識。

進一步地,關于模型空間的后驗概率密度σM(m)以及關于數據空間的后驗概率密度σD(d)可由邊緣概率密度公式計算得到:

假設數據空間與模型空間為兩個獨立空間,聯合空間中的均勻概率密度μ(d,m)為:

(4)

其中,μD(d)和μM(m)分別表示關于觀測數據和模型參數的均勻概率密度函數。(4)式表示模型空間中的模型參數化與數據空間中的數據參數化無關。

假設觀測數據與模型參數之間的先驗信息無關,公式(1)中的聯合先驗概率密度ρ(d,m)變為:

(5)

其中,ρD(d)和ρM(m)分別表示關于觀測數據和模型參數的先驗概率密度函數。

地球物理反演(包括其它基于數據的反演)的本質都是從數據中學習到有用的信息并據此得到客觀認識。實測數據和要估計參數之間存在一定的關系是反演的基礎,若二者沒有關系或關系極弱,則不可能從實測數據中得到要估計的參數信息。據此,聯合概率密度Θ(d,m)可用關于模型參數的均勻概率密度函數及其對應的條件概率密度函數來表示:

(6)

其中,θ(d|m)為在給定的模型參數m的條件下能夠預測觀測數據d的條件概率密度函數,表示要反演的參數與實測數據之間存在的關系。均勻概率密度函數μM(m)的引入表示模型參數的選擇與地震數據的理論預測無關。

因此,關于模型空間的后驗概率密度σM(m)變為:

(7)

(8)

此時的估計結果是用后驗概率密度函數得到的加權平均結果。只有當后驗概率密度函數是單峰的,譬如高斯或廣義高斯,此期望值才有意義。

(9)

方差越大,參數估計的精度在概率意義上越低。

顯然,從概率論的角度看,利用公式(8)進行參數估計是合理的。但是目前缺乏有效的方法計算高維模型參數的聯合后驗概率密度函數。(8)式定義的積分也需要巨大的計算量。為此,將求后驗概率密度函數的期望值轉為求取使得后驗概率密度最大化的模型參數估計結果作為反演解:

(10)

由于公式(7)中的后驗概率密度函數較為復雜,無法直接求解該式對應的最優化問題。需要進一步引入先驗信息。

首先假設模型空間M和數據空間D均為線性空間,而且模型參數和觀測數據在M和D中的采樣沒有先驗認識,則均勻概率密度函數為:

(11)

假設聯系觀測數據與模型參數之間的物理規律不確定性可以忽略,則:

(12)

假設關于模型參數m的先驗信息是具有高斯分布的概率密度函數,則:

(13)

假設數據空間中數據d與實測數據dobs之差滿足高斯分布:

(14)

融合上述假設后,模型空間的后驗概率密度σM(m)變為:

其中,S(m)為誤差泛函。由單調性可知,后驗概率密度函數最大化等價于誤差泛函最小化:

(16)

進一步地,誤差泛函有以下形式:

(17)

基于Bayes估計理論可以從概率統計的角度理解反演成像,對反演結果進行精確評價。數據協方差(逆)陣和模型協方差(逆)陣攜帶了數據和模型參數的不確定性信息,據此可以更方便地進行正則化處理。另外,只有當觀測數據與模型參數之間存在弱非線性關系,后驗概率密度函數為高斯或廣義高斯分布時,后驗概率密度最大化以及據此發展出的代價函數最小化的方法才有可能用于實際計算。線性的或非線性的最小二乘方法是最基本的算法。

2 數據協方差與模型協方差

從以上分析可以看出,(15)式和(17)式中的數據協方差CD與模型協方差CM在模型參數的優化求解中具有特別重要的意義,關于數據和模型參數的先驗信息很多都蘊含在這兩個協方差算子中。所謂關于數據的正則化和關于模型參數的正則化處理主要在于如何利用好這兩個協方差矩陣。

首先,從算子投影的角度,引入數據和模型參數的原空間以及它們各自的對偶空間。假設M和D為線性空間,M*和D*分別為其對應的對偶空間,m∈M,d∈D,m*∈M*,d*∈D*。令L為從模型空間映射到數據空間的線性算子,即:

(18)

LT定義為從D*映射到M*的線性算子,存在如下的對偶內積關系:

(19)

假設在線性空間Z中存在一個加權算子W,該算子是一個線性、對稱、正定的映射算子,它將線性空間Z中的元素z映射到其對偶空間Z*中。線性空間Z中任意兩個向量u,z的標量積定義為:

(20)

由加權算子W建立了原空間與對偶空間的映射,其中,W算子的逆算子C=W-1稱為協方差算子。令z與z*分別為原空間Z與對偶空間Z*的兩個元素,則有:

(21)

因此,可以說協方差矩陣將一個線性空間投影到其對偶空間中。而轉置算子定義了從數據對偶空間向模型參數對偶空間的映射關系。標量積的定義引入了協方差算子,建立起原空間的元素到對偶空間元素的映射關系,反之亦然。

從映射的角度看線性反演,核心是建立起原數據空間向原模型參數空間的映射關系。這需要引入伴隨算子的概念。L的伴隨算子L*定義為從原數據空間D映射到原模型參數空間M的線性算子,且對所有的元素m∈M和d∈D滿足:

(22)

經過標量積和轉置算子的運算可得到:

(23)

因此伴隨算子為:

(24)

最終從原數據空間D映射到原模型參數空間M的映射過程(圖1)為:

(25)

據此可以看到,數據和模型的協方差(逆)算子連接起了原空間與對偶空間的橋梁。

接著從線性反演的角度來進一步理解協方差算子。假設聯系數據和模型之間的正問題是線性、確定性的,此時正問題(12)式變成公式(18),則對應的Bayes估計后驗概率密度為:

圖1 原空間與對偶空間關系示意圖解

(26)

誤差泛函為:

(27)

對誤差泛函求導:

(28)

對應的法方程為:

(29)

法方程的解為:

(30)

已知協方差陣是對稱正定陣,可以對其進行正交分解。數據協方差陣對數據向量進行變換可以認為是對數據進行了特征表達。模型參數對應的協方差陣具有同樣的作用。因此,可以認為對偶空間是特征數據和特征模型參數各自形成的線性空間。特征數據和特征模型參數之間的關系更緊密、或者更線性、或更弱非線性。轉置算子建立起它們之間一對一的投影關系。這對提高反演結果的穩定性非常有益。這是反演問題中正則化處理最本質的內涵。

分別對模型和數據的協方差陣做正交分解:

(31)

其中,U和V分別為左奇異向量和右奇異向量,Σ為奇異值組成的對角矩陣。

整理(31)式,得到特征表達后的解為:

(32)

分析(32)式,得到以下認識:

1)Vd算子是將數據向量d投影到正交坐標系中進行特征提取,相當于對數據進行預處理,提取數據中的相干信號;在高維數據中,可認為是對數據進行方向譜分解,等價于提取局部平面波過程。

2)VL算子是對傳播算子進行相應處理,VL對應局部平面波傳播算子。(VL)T是對預條件處理后的局部平面波殘差(或數據)進行反投影的算子。

3)Umprior算子是將模型參數向量mprior投影到對應的正交坐標系中,等價于在模型參數場中提取線性結構信息。

對數據實施預條件處理,提出特征波[2]的概念,同時配合相應的模型參數特征表達,是反演成像中更為根本的正則化思想。據此可以將一個強非線性問題提成凸問題,或將一個很病態問題提成臨近良態問題,降低反問題求解的難度,提高反演成像的精度。

3 數據協方差在反演(偏移)成像中的應用

在地震波反演成像中,人們對模型參數的正則化進行了深入的研究,提出了很多有效的方法,但是對數據的正則化考慮較少。根據前面討論的基本概念,我們引入數據正則化方法,提出了自適應數據變化的加權疊前深度偏移成像方法。

3.1 加權最小二乘偏移

忽略模型的正則化項,關于數據匹配項的最小二乘誤差泛函可寫成:

(33)

其中,e表示觀測數據與模擬數據的誤差(殘差)向量。在Bayes反演框架下,誤差向量(寫成分量形式為:e=(e1,e2,…,eN)為隨機變量,該隨機變量的協方差矩陣可表示為:

(34)

當誤差變量ei(i=1,2,…,N)各分量之間高斯統計無關時,Cee變成Σ,更理想地,變成σ2I。(33)式考慮的就是σ2=1時高斯白噪聲情況下的反演成像。根據線性高斯情況下Bayes估計理論假設,對于一般的正問題預測及實際噪聲情況,可認為Cee滿足經典的正態分布。當反演解趨于(等于)真解時,數據殘差滿足高斯白噪分布。此時,誤差向量e中沒有任何有效的信號成分。若正問題不能模擬所有信號成分,或噪聲不符合高斯白噪,導致Cee偏離理想的期望情況時,應引入協方差逆算子懲罰掉數據中的相關信號成分以滿足上述誤差為高斯白噪的假設。因此需要在地震波成像中將(33)式修改為加權最小二乘定義的誤差泛函。

加權最小二乘誤差泛函可表示為:

(35)

其中,Wd為數據加權函數,表示對數據殘差中的不同成分懲罰程度不同,體現了數據殘差中不同成分的可靠性。依據上述討論,最理想的加權函數應取為誤差向量e的協方差矩陣的逆矩陣:

(36)

此時,加權最小二乘偏移成像對應的法方程的解為:

(37)

對應的迭代公式可寫為:

(38)

其中,k為迭代次數,μk為迭代步長。

然而在實際應用中,很難直接通過數據計算出完整的協方差陣及其逆矩陣。數據協方差(逆)陣是與觀測系統有關的矩陣,一般需要根據其物理含義賦予不同的假設形式。

3.2 Beam偏移

(39)

若假設加權系數為對角陣:

(40)

此時數據殘差中各個分量之間不相關且權重不同,權系數依據數據中的噪聲水平或數據的可靠性選擇合適的權重。對于噪聲強的數據給予小的權重,可靠性強的數據,給予大的權重。

鑒于勘探地震中殘差向量存在空間相關性(這是地震波場的空間連續性決定的),可以用殘差向量相鄰道的數據相關性來構造加權函數。直接在時間空間域考察(殘差)數據的空間相關性比較困難,一般地,需要在變換域按照特征數據來確定合適的權系數。

根據地震波在地下介質中的傳播特征,可以選擇用波束合成(Beam Forming)或線性Radon變換的方式構建權函數。該思想利用局部平面波波前到達不同檢波器的時差,對觀測數據進行時移疊加,使局部平面波相干疊加而其它信號相互抵消,從而去除地震數據中的某些噪聲。這就是勘探地震中Beam偏移的思想,可看作將數據向量d投影到正交坐標系中進行特征提取的過程。比如Gaussian Beam偏移[10-12]是在中心射線的傍軸近似展開,采用局部高斯函數計算方式改善焦散、陰影區和臨界區振幅計算不準確問題。SUN等[13]采用在局部地震道內進行束疊加后再偏移的方法,提高了成像結果的信噪比和計算效率。

實測地震記錄中若存在波動方程無法描述的波現象,將會導致預測誤差不符合高斯分布假設。一種有效的做法是將正演算子所能夠描述的波現象從觀測記錄中剝離出來,然后用分離后的波現象進行參數反演[9]。

3.3 自適應數據加權偏移

對比公式(37)與公式(39)可以看出,數據協方差矩陣主要作用在常規偏移成像結果和Hessian矩陣中。LS_RTM是在Bayes反演框架下一種有效的估計模型參數高波數擾動的保真成像技術。理想的LS_RTM需要估計精確的Hessian逆算子并作用到常規的成像結果上。然而由于Hessian矩陣過于龐大,無法直接估計精確的Hessian逆算子。即使可以采用各種迭代法近似求取Hessian逆算子,仍受限于觀測數據不完備、正算子不準確以及地震子波未知等多種因素影響,需要付出巨大的計算代價才能收斂,甚至出現不收斂現象。因此本文忽略Hessian矩陣,僅在常規的疊前深度偏移成像中考慮應用數據協方差。此時加權最小二乘法方程的解退化為:

(41)

疊前深度偏移成像中誤差向量e是觀測的地表數據減去背景波場(直達波、透射波等)后的一次波(及散射波)。其中包含相干的波場成分,不是隨機噪聲,更不是白噪聲。此時引入加權函數提高成像質量的物理解釋是:We將數據進行了特征表達,各特征成分之間不相關,高斯假設下相互獨立,滿足了線性最小二乘反演對殘差向量的假設,可以得到無偏和方差最小的反演解,這是高質量反演成像的目標。

考慮到數據協方差矩陣逆算子的作用,直接用Cee求逆形成加權系數矩陣不可行。如何合理選擇加權算子W是重要的研究內容。利用數據中蘊含的空間相關信息近似地確定加權函數W是基于數據自適應加權的PSDM方法要處理的核心問題。最基本的思想是測量數據中蘊含的相關信息,據此構造加權算子W。

基于上述理論分析,假設加權系數矩陣為對角矩陣,加權系數可由觀測波場中的相干成分近似計算。采用傾角掃描和動態時間規整算法確定加權系數,其基本流程為:

1) 在局部的共炮(CS)/共中心點(CMP)道集中提取局部波場;

2) 在選擇的局部波場中,采用傾角掃描的方式獲取不同道之間的線性相位擾動;

3) 采用動態時間規整算法校正非線性相位擾動;

4) 在拉平后的道集上利用相似系數譜方法估計信息中的相干成分,從而估計出加權系數;

5) 重復步驟1)到步驟4),直到估計出完整的加權系數;

6) 采用公式(41)進行偏移成像。

關于具體的傾角掃描和動態時間規整算法可參考文獻[14]。

加權算子的施加方式有:①直接改造疊前數據;②作用于Kirchhoff積分公式中的加權系數上;③改造波場外推成像方法中的相關成像條件。本文僅實現了第一種方法,后兩種方法在今后的研究中逐步展開。

3.4 數值實驗

3.4.1 模擬數據

采用圖2所示單層理論模型模擬數據對本文方法進行測試。圖2a為采用有限差分法模擬的單炮記錄,共601道,道間隔10m,時間采樣點數為5000,時間采樣率為0.5ms;圖2b為加入隨機噪聲后的單炮記錄,其中信噪比為2;圖2c為采用本文方法計算的加權系數,由圖2c可以看到,有效信號(相干信號)存在的地方,加權系數接近于1,而在非相干信號區域,加權系數基本為0;圖2d為采用數據加權處理后得到的單炮記錄,由圖2d可以看出,噪聲得到了有效消除。圖3a和圖3b分別給出了數據加權處理前、后的含噪單炮記錄局部放大結果。由圖3可以看出,原始單炮記錄中包含能量較強的隨機干擾噪聲,經過數據加權處理后,隨機干擾噪聲得到了消除,較為完整地保留了有效信號。

圖2 單層理論模型模擬數據測試結果a 原始單炮記錄; b 含噪單炮記錄; c 采用本文方法計算的加權系數; d 數據加權處理后的單炮記錄

3.4.2 實際地震數據

采用某工區二維炮集數據測試本文方法的有效性。圖4a為原始單炮記錄,可以看到,原始單炮記錄中含有較多的相干噪聲,深層同相軸被噪聲掩蓋,無法識別;采用本文方法處理后的單炮記錄如圖4b所示,可以看出,噪聲得到有效去除,無論是淺層還是深層,同相軸連續性更好。對應的偏移成像結果如圖5所示。從圖5可以清晰地看到,采用本文方法處理后的偏移剖面成像質量得到明顯提升,同相軸更加清晰連續,成像剖面上噪聲明顯減弱。

采用某工區三維CMP道集數據測試本文方法的效果。圖6a為原始的CMP道集;圖6b為數據加權處理后的CMP道集。由圖6可以看出,淺層能量強的同相軸表現出較強的相似性,在處理后的道集中得到增強,而在沒有相干信號的地方,原始道集中含有較多的噪聲,處理后的道集中噪聲得到有效去除。圖7為某條測線采用本文方法處理前、后的偏移成像結果。對比圖7a和圖7b發現,采用本文方法處理后,成像質量得到明顯改善。

圖4 采用本文方法處理前(a)、后(b)的單炮記錄

圖5 采用本文方法處理前(a)、后(b)的偏移成像剖面

圖6 數據加權處理前(a)、后(b)的CMP道集

圖7 某測線采用本文方法處理前(a)、后(b)的偏移成像結果

4 結論與討論

“兩寬一高”地震勘探的目的是進行更精細的油藏描述。“兩寬一高”地震數據采集提供了更完善的地震數據,但是地震波成像技術的發展相對滯后,不能充分挖掘出所采集的巨量數據中蘊含的與儲層有關的信息。地震波成像應該逐漸地奠定在Bayes估計理論基礎上,從信息綜合的角度提高地震波成像的質量。FWI和LS_RTM充分關注了模型參數的正則化,促進了成像質量的提高。但是,對數據正則化的關注有所欠缺,將數據中蘊含的先驗信息融合在成像過程中進一步提高成像質量是值得重視的。

為此,本文首先在Bayes框架下討論了地震波反演成像問題。在地震波正問題(波場和數據間的關系)是弱線性的假設下,反演成像的核心問題是基于后驗概率密度最大化原則的參數估計問題。若觀測噪聲(嚴格地講是預測誤差)為高斯白噪,模型參數的概率分布為高斯分布情況下,最小二乘估計可以在最小二乘意義下實現。但是,實際數據并不滿足上述假設,即預測誤差滿足有色高斯分布、甚至不滿足高斯分布。此時,引入數據協方差(逆)矩陣是提高反演成像質量所必須的。數據和模型的協方差(逆)矩陣連接了各自的原空間與對偶空間,通過協方差矩陣進行數據和模型參數的正交投影,相當于分別對數據和模型參數進行了特征表達。由此可以在數據預條件的基礎上提取特征波場,并對模型參數進行特征表達,建立起更為緊密的特征波場與特征參數之間的聯系,構建更凸的反問題。這是進行數據自適應加權疊前深度偏移的數學物理理論基礎。

本文從最小二乘誤差泛函角度分析了數據協方差(逆)算子的作用,指出在地震波反演成像中加權最小二乘引入的必要性。在假設數據協方差矩陣為對角陣的基礎上,針對疊前深度偏移成像,提出了一種自適應數據加權的數據預條件處理方法。該方法通過傾角掃描和動態時間規整方法自動計算數據加權系數,將其融合在偏移成像過程中。理論和實際資料測試結果驗證了該方法的有效性。

“兩寬一高”地震勘探是一項綜合性技術,采集技術還有巨大的發揮空間,反演成像技術的差距更大。基于Bayes理論,充分地融合關于數據的先驗信息和關于模型參數的先驗信息是提高成像質量的必由之路。地震波反演成像中對數據不確定引起的成像結果的不確定性研究還處于初始階段,本文的研究工作僅僅是較為初步的探索。

致謝:感謝中石油勘探開發研究院及西北分院、中海油研究院和湛江分公司、中國石油化工股份有限公司石油物探技術研究院和勝利油田分公司對波現象與智能反演成像研究組(WPI)研究工作的資助與支持。