基于非平穩濾波算子的成像域最小二乘偏移

郭 頌,王華忠,胡江濤

(1.波現象與智能反演成像研究組(WPI),同濟大學海洋與地球科學學院,上海200092;2.油氣藏地質及開發工程國家重點實驗室,成都理工大學,四川成都610059)

油氣勘探的最終目標是儲層描述。我國油氣地震勘探的對象正在由大尺度的構造油氣藏向小尺度的(薄互層)巖性油氣藏轉變,因此需要采用寬帶波阻抗成像結果進行巖性油氣藏儲層的精確描述[1]。基于散射和逆散射的全波形反演(Full Waveform Inversion,FWI)成像技術試圖利用地表觀測的疊前數據估計地下寬波數帶的參數場,基于此反演結果,直接進行儲層描述。然而,實際情況下,地震波場與地下介質之間的關系十分復雜,具有很強的非線性性,僅僅利用地表觀測的有限帶寬、有限孔徑和無假頻數據,FWI的實用化并不成功。當前,它僅僅能在低頻和長偏移距地表觀測數據下,得到較好的背景速度的估計[2]。另一方面,勘探地震中地下介質以層狀為主,可以將寬帶阻抗分為背景阻抗和反射系數兩部分分別進行線性反演。隨著油氣地震勘探技術逐漸向單點高密度、寬帶、寬方位地震勘探方向轉變,“兩寬一高”的地震數據為精確的背景(偏移)速度和寬帶反射系數估計提供了數據基礎。結合來自測井數據的密度建模,采用合理的深度域反演方法,將背景阻抗和寬帶反射系數融合成寬帶波阻抗成像結果,是進行精確油氣藏描述的重要技術方向[1]。

假設已經得到準確的背景速度,本文聚焦于采用最小二乘疊前深度偏移(LS-PSDM)估計地下保真、高分辨率的反射系數,為寬帶的波阻抗成像奠定基礎。LS-PSDM有數據域和成像域兩種求解手段,其本質目標都是估計Hessian矩陣的逆,作用在常規的成像剖面上,從而得到真振幅、高分辨率的反演成像結果。原則上,數據域與成像域LS-PSDM等價,但是兩者的具體實現過程不同,也會導致不同的求解代價及反演結果。

數據域求解LS-PSDM基于反偏移數據與觀測數據的匹配,利用數據殘差計算梯度進行迭代求解。每輪迭代都需要進行反偏移、偏移以及計算步長,至少是兩次常規偏移的計算量。通常需要多達10次的迭代才可以得到滿意的結果,這導致數據域求解LS-PSDM的計算量十分巨大。另外,由于實際數據中含有噪聲、背景速度不準、子波未知等因素會導致迭代收斂緩慢甚至不收斂,無法得到滿意的結果。實際情況下,數據域LS-PSDM由于其高昂的計算成本以及較慢的收斂速率,通常不能很好地應用于實際生產中。

相對于數據域的LS-PSDM,成像域LS-PSDM可以視為圖像去模糊問題,直接對成像剖面進行處理,避免了數據域迭代中的反偏移、偏移計算,極大地減少了計算量。成像域LS-PSDM的核心是顯式計算Hessian矩陣(逆)。由于全Hessian矩陣十分龐大,無法直接計算、存儲,所以需要引入近似計算策略使得成像域LS-PSDM得以應用。

本文先對數據域與成像域的LS-PSDM方法原理進行對比,再介紹一種基于非平穩濾波算子構造近似的Hessian矩陣逆的方法[3]。采用該方法得到的Hessian矩陣逆的近似可以直接作用在常規成像結果上,得到成像域LS-PSDM的反演結果,也可以作為預條件算子,加速數據域LS-PSDM的迭代收斂。數值試驗結果證明,利用求解非平穩濾波算子近似構造的Hessian矩陣的逆實現成像域的LS-PSDM,可以有效地均衡常規成像結果的振幅,并且在一定程度上提高了成像分辨率。

1 數據域與成像域的LS-PSDM

假設波傳播算子f能夠很好地模擬觀測數據,并且數據中的噪聲滿足高斯分布,基于Bayes反演理論的全波形反演(FWI)通過求解如下目標函數估計地下模型參數:

(1)

其中,dobs代表地表觀測數據,m為地下模型參數。實際情況下,波場與地下介質之間的關系十分復雜,具有很強的非線性性,因此FWI是一個強非線性問題,理論上需要采用非線性尋優算法(例如蒙特卡洛方法)來求解。但是,對于龐大的疊前地震數據,由于計算成本等因素的限制,當前并沒有可行的非線性尋優算法求解FWI問題。因此,基于不同的初始模型,FWI通常被近似為一系列線性問題來求解。LS-PSDM為其中一個線性化子問題,在已知地下模型參數光滑背景m0的基礎上,估計反射系數。一般地,LS-PSDM利用線性化的正算子(反偏移算子)L建立反射系數Δm與一階反射數據dpri間的線性關系:

(2)

此時,可以反演線性目標函數(公式(3)),估計地下反射系數:

(3)

對于線性目標函數(公式(3)),數學上采用梯度導引的方法求解。利用正算子L預測當前模型下的數據dcal,構造線性正算子的伴隨算子LT,將觀測數據與預測數據的殘差反投影到模型空間,便可以得到線性目標函數(公式(3))的梯度:

(4)

其中,dcal和LT需要滿足如下關系:

(5)

其中,“〈〉”代表內積。梯度是求解線性問題(公式(3))的核心,不同的迭代算法如最速下降法、共軛梯度法等都是基于梯度構造不同的迭代策略來求解線性問題的方法。然而本質上,梯度是由正算子L決定的,只要給定具體的正算子,梯度便可以用公式(4)計算得到。這說明,對于線性化的地球物理反演問題,線性正算子L的選擇決定了線性反演的成敗,如果選定的線性正算子具有較好的表達數據的能力,那么線性反演便具有較好的收斂性。反之,線性反演會收斂緩慢甚至不收斂,這對于求解諸如LS-PSDM的大規模地球物理線性反問題是不可接受的。因此,求解LS-PSDM問題的第一步,也是最重要的一步,就是構造一個能夠較好地表達一次反射波數據的線性正算子。

在確定了線性正算子的基礎上,LS-PSDM剩下的問題便是數值求解線性目標函數(公式(3))。一般有數據域和成像域兩種求解手段,其目標都是估計線性算子L的Hessian矩陣的逆,并將其作用于常規的成像剖面上,從而得到真振幅、高分辨率的反演成像結果。

數據域的LS-PSDM基于梯度公式(4),迭代地求解線性目標函數(公式(3))。當目標函數收斂到給定的閾值范圍內(或者達到給定的迭代次數)時,認為得到了數據域LS-PSDM的解。以最速下降法為例,對于每一步迭代,利用線性正算子計算當前模型下的反偏移數據(第一次迭代該步驟省略),將預測數據與觀測的一次反射波數據的殘差通過正算子的伴隨算子反投影(偏移),得到本次迭代的梯度。由此可以看出,數據域LS-PSDM的一次迭代,至少需要一次全空間的線性正演(用于計算數據殘差),一次全空間的偏移(用于計算梯度),再加上估計步長的計算量。一般地,數據域LS-PSDM至少需要10次迭代才能收斂到令人滿意的水平[4],而對于實際數據而言,由于線性正算子不能很好地模擬實際波場,導致數據域的迭代收斂通常很緩慢,甚至不收斂[5]。

成像域的LS-PSDM從另一個途徑求解目標函數(公式(3))。對于線性目標函數,當其梯度(公式(4))等于0時,可得法方程:

(6)

令H=LTL,代表線性算子L的Hessian矩陣,m1=LTdpri,代表常規的偏移成像結果。成像域LS-PSDM試圖直接求解法方程(6),估計Hessian矩陣的逆并將其作用到常規成像結果上,從而得到LS-PSDM的解:

(7)

相較于數據域的LS-PSDM需要在全空間進行線性正演、偏移,成像域LS-PSDM操作靈活,可以僅對感興趣的目標區域反演,進而節約計算成本[6]。但是,成像域LS-PSDM需要顯式構造Hessian矩陣,并且對其求逆,而Hessian矩陣是一個極大規模矩陣,其規模大小為模型大小的平方,以目前的計算能力,無法對Hessian矩陣直接求解及存儲,求逆更是十分困難。所以,尋求Hessian矩陣逆的快速近似構造方法具有非常重要的實際意義。為了克服這一問題,使成像域LS-PSDM實用化,我們提出了一種基于非平穩濾波算子近似Hessian矩陣逆的方法,大幅降低了計算量,并且可以得到穩定的成像域LS-PSDM近似解,成像質量相較常規成像結果有了較大幅度的提升。

2 基于非平穩濾波算子的成像域LS-PSDM

給定觀測的一階反射數據dpri,可得常規偏移成像結果m1:

(8)

利用線性正算子L,基于常規偏移成像結果m1進行反偏移,獲得反偏移數據d1=Lm1,利用偏移算子LT對反偏移數據d1再做一次偏移,獲得成像結果m2:

(9)

此時,m1與m2有如下關系:

(10)

由(10)式可知,m2為m1經過Hessian矩陣模糊作用后的成像結果,常規成像結果m1為地下真實反射系數經過Hessian矩陣模糊作用后的結果。由于m1與m2已知,可以構造非平穩匹配濾波算子F,使得:

(11)

其中,F≈(LTL)-1,可認為是Hessian算子逆的低秩近似。

為了求取匹配濾波算子F,可構造如下最小二乘問題:

(12)

對于非平穩濾波算子,原則上對于模型空間每一個點i,存在一組濾波系數。Fj,i代表對應模型第i個點的第j個濾波系數。由于非平穩濾波算子F系數的數目大于模型空間的樣點數,所以(12)式為一個欠定問題,需要額外的先驗信息(正則化)對濾波器系數進行約束。由于濾波算子F為Hessian算子逆的近似,應該包含模型的結構信息,同時,模型不同樣點處的濾波系數應該光滑變化,因此,本文采用具有保邊界能力的總變差(total variation,TV)正則化對非平穩濾波算子F的求取進行約束,以便得到穩定的濾波算子求解結果。由于F≈(LTL)-1,將求得的濾波算子F作用到常規成像結果m1上,便可得到成像域LS-PSDM的近似結果:

(13)

綜上所述,基于非平穩濾波算子的成像域LS-PSDM方法可總結如下:

1) 給定觀測數據dpri,進行一次常規偏移,得到成像結果m1;

2) 對成像結果m1進行反偏移,對反偏移數據再進行一次偏移,得到成像結果m2;

3) 利用m1與m2估計非平穩濾波算子F;

4) 將估計得到的濾波算子F作用到成像結果m1上,得到近似的成像域LS-PSDM反演結果。

上述方法僅需要進行兩次偏移,一次反偏移,以及求取非平穩濾波算子的計算量。相對于迭代求解數據域LS-PSDM,計算量大大降低,并且避免了數據域迭代收斂緩慢甚至不收斂的問題。對于不同的傳播、偏移算子,m1與m2已經包含了所有波傳播的物理規律、算子特性,所以求取非平穩濾波算子的方式一樣。因此,該方法具備較好的通用性,適用于各類復雜的波傳播、偏移算子(如彈性波算子、各向異性算子等等)。

3 TV正則化約束下非平穩濾波算子的構造

對于匹配濾波問題,有如公式(11)所示的線性關系。由于濾波算子F為Hessian矩陣逆的近似,應該包含模型的結構信息,并且模型相鄰樣點處的濾波算子系數應該光滑變化,因此利用TV正則化約束模型不同點對應濾波算子系數的變化。構造如下有約束的濾波算子估計問題:

(14)

(15)式為無約束凸優化問題,可采用Split-Bregman算法[8],將基于L1范數的總變差約束項與基于L2范數的數據逼近項解耦。求解(15)式的步驟如下。

首先進行變量替換:

(17)

(17)式顯然與(15)式等價。利用Bregman迭代算法求解(17)式,有:

現在問題集中在求解(18)式,可分解為如下兩步:

(21)

(22)

其中,(21)式是一個二次凸優化問題,可以利用梯度導引法(共軛梯度法)求解。利用廣義的shrinkage公式,(22)式有顯式的求解算法:

其中,

(25)

綜上所述,TV正則化約束下,非平穩匹配濾波器F的計算方法為:

Initialize:F0=0,bi,j=di,j=0

fori=1toN(求解無約束問題公式(15))

end

end

采用以上求解算法,可以穩健、高效地計算出非平穩濾波器F,即Hessian算子逆的近似。

4 數值試驗

圖1為局部Sigsbee模型。根據該模型背景(光滑)成分和觀測數據,首先進行常規偏移,得到偏移成像結果m1,如圖2所示。基于該成像結果進行反偏移,對反偏移數據再進行一次偏移成像,得到第2個成像結果m2,如圖3所示。

由圖2和圖3可以看出,兩個成像結果有較大差別,m2在深層及邊界處照明不足的區域幅值更弱,分辨率更低。根據前文介紹的方法利用m1和m2求解非平穩濾波算子F,可視為Hessian矩陣逆的低秩近似。對于模型每一點,選取的濾波器相當于一個11×11的二維褶積算子,計算得到的非平穩濾波器零延遲系數如圖4所示。從圖4可以看出,計算得到的非平穩濾波算子零延遲系數對模型深層以及兩邊照明較弱的區域有較大的補償,而且攜帶了模型的構造信息,這也正體現了Hessian矩陣逆對角線的作用。將求得的非平穩濾波器F作用到常規成像結果m1上,可以得到成像域LS-PSDM反演結果,如圖5所示。對比圖5與圖2可以看出,成像域LS-PSDM反演結果比常規偏移成像結果振幅更均衡,弱照明以及深層成像區域的振幅得到了充分的提高,使得成像結果更保真,分辨率也得到了提高。

圖1 局部Sigsbee模型

圖2 局部Sigsbee模型數據常規偏移成像結果

圖3 局部Sigsbee模型數據常規偏移后再進行反偏移,對反偏移數據再進行偏移后的成像結果

圖4 局部Sigsbee模型非平穩濾波算子零延遲系數

對地表500m,3750m以及7000m處成像域LS-PSDM反演結果、常規偏移成像結果和真實反射系數進行抽道(圖6到圖8)。由圖6到圖8可以明顯地看出,與常規偏移成像結果相比,成像域LS-PS-DM反演可以得到更保真的反演成像結果,分辨率也有所提高,尤其是在弱照明和深層區域。

圖5 局部Sigsbee模型成像域LS-PSDM反演結果

圖6 地表500m處成像域LS-PSDM反演結果、常規偏移成像結果和真實反射系數的抽道結果

圖7 地表3750m處成像域LS-PSDM反演結果、常規偏移成像結果和真實反射系數的抽道結果

圖8 地表7000m處成像域LS-PSDM反演結果、常規偏移成像結果和真實反射系數的抽道結果

采用如圖9所示的Marmousi速度模型進一步說明本文方法對成像結果的改善效果。Marmousi速度模型的常規偏移結果如圖10所示。對常規偏移結果再進行反偏移、偏移后的成像結果如圖11所示。利用m1,m2構造非平穩濾波器,對于模型每一點,選取一個11×11的二維褶積算子,計算得到的非平穩濾波算子零延遲系數如圖12所示。將求得的非平穩濾波算子F作用到常規偏移成像結果m1上,得到成像域LS-PSDM反演結果,如圖13所示。對比圖13與圖10可以看出,通過計算近似的Hessian算子逆,作用到常規偏移成像結果上,得到的成像域LS-PSDM結果比常規偏移成像結果振幅更均衡,弱照明區域的振幅得到了充分的提高,成像結果更保真。

圖9 Marmousi速度模型

圖10 Marmousi速度模型的常規偏移成像結果

圖11 Marmousi速度模型常規偏移后再進行反偏移,對反偏移數據再進行偏移后的成像結果

圖12 Marmousi速度模型非平穩濾波算子零延遲系數

圖13 Marmousi速度模型的成像域LS-PSDM反演結果

5 討論

油氣地震勘探的目標是識別與描述含油氣儲層。寬帶波阻抗成像更有利于巖性油氣藏的(定量)精細描述?！皟蓪捯桓摺钡卣饠祿杉呀洺蔀榈卣饠祿杉暮诵姆椒夹g,為寬帶波阻抗的估計提供了數據基礎。盡管理論上FWI具有直接采用疊前地震數據估計寬帶波阻抗的潛力,但是由于地下波場與介質之間的強非線性關系,對于實際問題,很難采用求解非線性的FWI問題得到理想的寬帶阻抗成像結果。由于地震勘探中地下介質以層狀介質為主,可以將寬帶阻抗分為背景阻抗和反射系數兩部分分別進行線性反演,進而將二者融合成寬帶波阻抗成像結果,這是進行精確油氣藏描述的重要技術方向。假設通過層析反演已經得到較為準確的背景速度,利用測井等信息進行背景密度建模,可以得到背景阻抗。高保真、高分辨率的反射系數由最小二乘偏移估計得到。利用合理的深度域反演方法,可以將背景阻抗以及高保真、高分辨率的反射系數的波數成分融合,并轉換為寬帶(絕對)波阻抗[9]。

本文討論了高保真、高分辨率反射系數估計方法。相對于只能正確地定位反射系數出現空間位置的常規偏移成像,最小二乘偏移試圖將Hessian算子的逆作用在常規偏移成像結果上,具有補償由于觀測系統不完備以及上覆介質復雜造成的地下照明不均勻,校正成像振幅,提高成像結果縱、橫向分辨率,消除成像假象等能力。因此,最小二乘偏移成像將是今后主流的成像技術。

但是,實際應用中,常規數據域迭代的最小二乘偏移成像計算量非常大,且迭代收斂速度緩慢,甚至不收斂。究其原因,可歸結為:給定的背景速度不準,地震子波未知,線性化的正問題不能很好地預測數據等。鑒于上述原因,LS-PSDM是一個病態性極強的反問題,正問題的不確定性以及解的非唯一性使得LS-PSDM目標函數的凸性很差,迭代求解收斂緩慢甚至不收斂,無法獲得理想的成像結果。為了克服數據域迭代求解LS-PSDM計算量過大以及不易收斂的問題,使得LS-PSDM更實用化,基于圖像去模糊(反褶積)的成像域LS-PSDM試圖從另一個途徑實現LS-PSDM的功能。成像域LS-PSDM的關鍵是估計線性算子的Hessian矩陣或者Hessian矩陣的逆,進而對成像結果進行去模糊或者直接將Hessian矩陣的逆作用在成像剖面上,提高成像質量。對于線性成像問題,數據域迭代求解LS-PSDM與成像域的圖像去模糊原理上等價,但相對于數據域迭代的LS-PSDM,成像域LS-PSDM操作更靈活,可以僅對感興趣的目標區域進行反演成像,而且一旦估計出Hessian算子的逆,便不需要進行迭代,也就回避了數據域迭代收斂緩慢和不收斂的問題,節約大量的計算時間,更實用。但是,Hessian矩陣是一個超大規模矩陣,對于實際三維工區,以目前計算機能力,直接對其計算和存儲基本不可能。因此,研究成像域LS-PSDM Hessian矩陣(逆)的近似計算具有很現實的意義。

6 結論

本文研究了成像域LS-PSDM Hessian矩陣逆的近似估計方法。利用常規成像剖面以及反偏移、偏移剖面構造非平穩濾波算子,以近似Hessian矩陣的逆。由于非平穩濾波算子的求解具有較強的欠定性,考慮到濾波器系數應具備成像結果的構造特征,本文在求解濾波算子過程中引入具有保留圖像邊界能力的TV正則化約束,以得到穩健、合理的濾波器求解結果。求得的非平穩濾波算子作為Hessian矩陣逆的近似直接作用到常規成像剖面上,即可得到成像域LS-PSDM反演結果。該方法極大地降低了LS-PSDM的計算量,并且計算穩定,適用于各類偏移算子?；诰植縎igsbee模型以及Marmousi模型的數值試驗結果表明,求得的非平穩濾波算子能夠較好地實現Hessian矩陣逆的功能,得到的成像域LS-PSDM結果相對于常規偏移成像結果振幅更均衡,弱照明區域的振幅得到了充分的提高,使得成像結果更保真,并且分辨率也有了一定提升。另外,非平穩濾波算子也可以作為預條件算子,應用于數據域LS-PSDM的迭代求解過程中,加速迭代收斂,以獲得更高分辨率、振幅更加保真的反射系數估計結果。求解LS-PSDM得到的高保真、高分辨率反射系數,為寬帶波阻抗成像提供了重要的數據基礎。

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