帶有橫向約束的全局優(yōu)化波阻抗反演方法及應(yīng)用

朱劍兵，高照奇，田亞軍，梁興城

(1.中國石油化工股份有限公司 勝利油田分公司物探研究院，山東 東營 257022；2.西安交通大學(xué) 信息與通信工程學(xué)院，陜西 西安 710049)

0 引言

地震反演在油氣藏勘探領(lǐng)域中具有重要作用，該方法能夠基于地震以及其他地球物理數(shù)據(jù)估計(jì)地下介質(zhì)參數(shù)[1-3]。地震反演往往是通過最小化描述觀測地震數(shù)據(jù)與合成地震數(shù)據(jù)誤差的目標(biāo)函數(shù)來實(shí)現(xiàn)的。地震波阻抗反演是一類地震反演問題，其特點(diǎn)是：①模型的非唯一性；②模型和數(shù)據(jù)的非線性關(guān)系，目標(biāo)函數(shù)存在多個(gè)局部極小值。

地震波阻抗反演問題的常用求解方法是局部優(yōu)化算法。這類方法從某一預(yù)設(shè)的初始模型出發(fā)，利用目標(biāo)函數(shù)的梯度信息更新模型參數(shù)，如共軛梯度法[4]，最速下降法[5]等。局部優(yōu)化方法的不足是優(yōu)化過程中可能陷入目標(biāo)函數(shù)的局部極小值。與局部優(yōu)化算法不同，全局優(yōu)化算法不依賴于目標(biāo)函數(shù)的梯度信息，具有跳出目標(biāo)函數(shù)局部極值的能力，因此此類方法越來越多地被應(yīng)用于求解非線性地震反演問題。在地球物理學(xué)中，應(yīng)用最廣泛的全局優(yōu)化類算法有模擬退火算法及其變體[6-7]，粒子群優(yōu)化算法[8-9]，遺傳算法[10-12]，差分進(jìn)化算法(differential evolution,DE)[13-15]等。

近年來，DE已經(jīng)成功地應(yīng)用于求解地震反問題[16]。但是，維數(shù)災(zāi)難問題阻礙了DE在高維模型空間反問題中的應(yīng)用。為了解決這個(gè)問題，Potter等[17]提出了協(xié)同進(jìn)化的策略將高維問題分解為多個(gè)低維問題進(jìn)行求解。隨后，Wang等[18]提出了協(xié)同變異差分進(jìn)化算法(cooperative coevolutionary DE,CCDE)，其將一個(gè)高維模型分解為多個(gè)低維子成分，并用“局部適應(yīng)度函數(shù)”來評(píng)價(jià)每個(gè)子成分的質(zhì)量，用來指導(dǎo)后續(xù)變異方向。盡管CCDE在求解高維反問題上比傳統(tǒng)DE具有明顯的優(yōu)勢，但其對(duì)搜索空間的范圍敏感。Gao等[19]提出了結(jié)合DE/rand/1變異策略和CCDE變異策略的多組變異差分進(jìn)化算法(multi-mutation DE,MMDE)。該方法對(duì)搜索空間的范圍不敏感，具有更快的收斂速度。

盡管上述全局優(yōu)化方法在地震波阻抗反演問題上表現(xiàn)出了優(yōu)越的性能，但這類方法均采用了逐道反演策略，未考慮模型的空間相關(guān)性，導(dǎo)致反演結(jié)果的橫向連續(xù)性差?？紤]到地下介質(zhì)參數(shù)是連續(xù)變化的，兩條相鄰地震道之間存在一定的聯(lián)系。因此，在反演過程中，應(yīng)考慮當(dāng)前地震道和相鄰地震道之間的相似性。為了解決這一問題，研究者們提出了一些多道反演策略來保證地下結(jié)構(gòu)的橫向連續(xù)性。Gholami[20]提出了帶有TV正則化的多道地震波阻抗反演方法，Zhang等[21]提出了一種基于各向異性全變差(ATV)正則化的多道反演方法。盡管這些假設(shè)先驗(yàn)信息的正則化方法可以提高反演模型的橫向連續(xù)性，但當(dāng)先驗(yàn)信息與實(shí)際情況不符時(shí)，這種方法可能會(huì)失效。最近，Cheng等[22]提出了一種受兩種范數(shù)約束的多道反演方法。此外，還針對(duì)橫向連續(xù)性問題引入了卡爾曼濾波器。然而，上述解決橫向連續(xù)性問題的方法與全局優(yōu)化算法結(jié)合存在難度。

本文提出了一種帶有橫向約束的多組變異差分進(jìn)化地震波阻抗反演方法。通過將旁道最優(yōu)解信息融入到當(dāng)前地震道的搜索空間中，有效地約束了當(dāng)前道波阻抗反演的搜索空間范圍，實(shí)現(xiàn)了逐道反演搜索空間的精確初始化。通過將本文所提出方法與傳統(tǒng)多組變異差分進(jìn)化算法進(jìn)行對(duì)比，實(shí)驗(yàn)結(jié)果表明，本文所提出方法不僅收斂速度更快，而且反演得到的波阻抗剖面橫向連續(xù)性更優(yōu)。此外，本文所提出方法被應(yīng)用于勝利油田某區(qū)塊波阻抗反演中，反演結(jié)果和測井資料具有很好的一致性，有效地刻畫了儲(chǔ)層砂巖厚度。

1 帶有橫向約束的多組變異差分進(jìn)化算法

1.1 多組變異差分進(jìn)化算法

多組變異差分進(jìn)化算法(MMDE)是一種全局優(yōu)化算法，其算法流程如圖 1所示。MMDE主要包含初始化、變異、交叉和選擇4個(gè)部分，下面將詳細(xì)介紹這幾部分內(nèi)容。

圖1 MMDE算法流程

1.1.1 初始化

假設(shè)MMDE擬求取一個(gè)目標(biāo)函數(shù)f(x)的全局極小值，其中x是一個(gè)D維的參數(shù)向量x={x1,x2,…,xD}。MMDE首先隨機(jī)初始化一個(gè)種群，種群中有N個(gè)個(gè)體(候選解)，每一個(gè)初始種群中的個(gè)體可以表示為：

在搜索空間內(nèi)均勻隨機(jī)初始化得到的。假定搜索空間的下限為xmin={x1,min,x2,min,…,xD,min}，上限為xmax={x1,max,x2,max,…,xD,max}，則初始群體中第i個(gè)個(gè)體的第j維可以表示為：

(2)

其中，randi,j(0,1)是均勻分布于0和1之間的隨機(jī)數(shù)。

1.1.2 變異

在群體進(jìn)化中，變異操作用于提高候選解的質(zhì)量。在MMDE中，種群中的每個(gè)個(gè)體都通過結(jié)合DE/rand/1和CCDE[18]的多組變異策略進(jìn)行更新。具體過程如下：

,

(3)

,

(5)

1.1.3 交叉

交叉操作用于實(shí)現(xiàn)種群中個(gè)體的基因交換從而提高種群的多樣性。具體來說，交叉操作通過將當(dāng)前候選解與其對(duì)應(yīng)的變異向量混合，以構(gòu)成一個(gè)試探向量。對(duì)于當(dāng)前種群中的第i個(gè)候選解，試探向量的第j個(gè)分量可以表示為：

(6)

1.1.4 選擇

選擇操作用于確定候選解還是其對(duì)應(yīng)的試探向量將被保留進(jìn)入下一次迭代，這一操作是通過比較兩者的目標(biāo)函數(shù)值而實(shí)現(xiàn)的：

(7)

經(jīng)過多次迭代，最終群體中目標(biāo)函數(shù)值最小的個(gè)體即為優(yōu)化問題的解。

1.2 帶有橫向約束的模型搜索空間設(shè)定方法

如圖 2a所示，傳統(tǒng)MMDE求解地震波阻抗反演問題時(shí)，每一道阻抗反演的搜索空間是通過在該道低頻阻抗的基礎(chǔ)上加減一個(gè)常數(shù)而設(shè)定的，可用如下公式表達(dá)：

,

(8)

考慮到地下介質(zhì)具有橫向連續(xù)性，相鄰兩道波阻抗之間應(yīng)具有很強(qiáng)的空間相關(guān)性。為充分利用這一先驗(yàn)信息，我們提出一種新的帶有橫向約束的模型搜索空間設(shè)定方法。如圖 2b所示，這一新方法利用旁道波阻抗反演結(jié)果約束本道波阻抗反演的搜索空間范圍，具體可表示如下：

a—傳統(tǒng)MMDE進(jìn)行地震波阻抗反演時(shí)采用的模型搜索空間； b—本文所提出的帶有橫向約束的模型搜索空間

(9)

1.3 基于帶有橫向約束的全局優(yōu)化地震波阻抗反演

地震波阻抗反演問題中，地震道d和波阻抗模型z可通過正演算子G(·)聯(lián)系起來，如下式所示：

d=G(z)

,

(10)

具體來說，地震道d可以表示為地震子波與反射系數(shù)序列的褶積：

d=w?r

,

(11)

其中符號(hào)“?”表示褶積運(yùn)算。反射系數(shù)序列與波阻抗的關(guān)系可以表示為：

(12)

其中，ri、zi、ρi和vi分別表示反射系數(shù)、阻抗、密度和縱波速度的第i個(gè)分量。

基于上述正演模型，地震波阻抗反演可描述為一個(gè)D維優(yōu)化問題，該問題的目標(biāo)函數(shù)可表示為：

,

(13)

其中，dobs和G(z)分別表示觀測地震數(shù)據(jù)與合成地震數(shù)據(jù)。本文采用上述所提出的帶有橫向約束的多組變異差分進(jìn)化算法優(yōu)化式(13)所示的目標(biāo)函數(shù)。

地震波阻抗反演中低頻信息的引入對(duì)得到絕對(duì)波阻抗參數(shù)至關(guān)重要。在本文所提出方法中，每一個(gè)地震道的低頻阻抗信息是在設(shè)定算法搜索空間范圍時(shí)引入的。如式(9)所示，每一道的搜索空間范圍除與旁道最優(yōu)解相關(guān)外，還和本道的低頻阻抗相關(guān)。在搜索空間范圍內(nèi)，算法通過迭代最小化式(13)所示的目標(biāo)函數(shù)逐漸恢復(fù)波阻抗的中高頻信息。

2 基于Marmousi II模型的算例分析

本節(jié)將所提出的理論方法應(yīng)用于基于Marmousi II模型的實(shí)驗(yàn)算例并與傳統(tǒng)方法進(jìn)行對(duì)比。實(shí)驗(yàn)中使用的阻抗模型如圖3a所示，該模型是通過對(duì)原始深度域Marmousi II模型進(jìn)行深時(shí)轉(zhuǎn)換及重采樣得到的。該模型橫向共有1 134道，每道長度為200 ms，采樣間隔為2 ms。通過對(duì)真實(shí)阻抗模型每一道使用巴特沃茲低通濾波器提取其低頻分量，構(gòu)建了如圖3b所示的低頻阻抗模型作為反演的初始模型。實(shí)驗(yàn)中，我們采用主頻為40 Hz的雷克子波生成地震數(shù)據(jù)，反演時(shí)假設(shè)子波已知。實(shí)驗(yàn)中分別采用所提出新方法與傳統(tǒng)多組變異差分進(jìn)化算法反演該模型的波阻抗，為保證對(duì)比的公平性，實(shí)驗(yàn)中兩種對(duì)比方法除模型搜索空間初始化方式不同外，其它反演參數(shù)保持一致。實(shí)驗(yàn)結(jié)果如圖 4、圖5和表 1所示。

a—真實(shí)阻抗模型；b—反演使用的初始阻抗模型

圖4 對(duì)比所提出新方法與傳統(tǒng)多組變異差分進(jìn)化算法的收斂曲線

a—所提出新方法反演的阻抗模型；b—多組變異差分進(jìn)化算法反演的阻抗模型；c、d—分別是a和b在CDP100～300處的阻抗模型放大對(duì)比

表1 定量化對(duì)比不同反演方法得到阻抗模型的質(zhì)量

圖4所示為兩種反演方法的收斂曲線對(duì)比。在全部1134道數(shù)據(jù)中，我們選擇第100道數(shù)據(jù)的反演過程進(jìn)行對(duì)比。如圖 4中藍(lán)色曲線所示，傳統(tǒng)多組變異差分進(jìn)化算法在經(jīng)過100次迭代后目標(biāo)函數(shù)值從7.3936下降到0.4822，這說明多組變異差分進(jìn)化算法可以在這一反演問題中收斂。與多組變異差分進(jìn)化算法相比，本文所提出的理論方法的收斂速度更快，該方法經(jīng)過100次迭代后目標(biāo)函數(shù)值從6.6910下降到5.7949 e-5。此外，該方法僅需要12次迭代目標(biāo)函數(shù)值就可以下降到0.4177，這一數(shù)值比多組變異差分進(jìn)化算法100次迭代后的目標(biāo)函數(shù)值還低?；谏鲜鰧?shí)驗(yàn)結(jié)果，可以看出所提出新方法具有顯著優(yōu)于傳統(tǒng)多組變異差分進(jìn)化算法的收斂速度，可在相同迭代次數(shù)的前提下得到更優(yōu)的反演結(jié)果。

圖5所示為兩種反演方法構(gòu)建的波阻抗模型。圖5a所示為所提出新方法反演得到的波阻抗模型，該模型相對(duì)于真實(shí)模型的信噪比(signal-to-noise ratio，SNR)為36.0658 dB、結(jié)構(gòu)相似性指數(shù)為0.9332。這些數(shù)據(jù)說明了該反演阻抗模型在阻抗結(jié)構(gòu)以及數(shù)值上均和真實(shí)模型有很好的一致性。此外，圖5c所示為所提出方法反演的波阻抗模型在CDP100～300處的放大圖。從圖中可見，反演阻抗模型具有良好的橫向連續(xù)性。與所提出理論方法不同，圖5b所示的傳統(tǒng)多組變異差分進(jìn)化算法反演的波阻抗模型精度不足。該模型相對(duì)于真實(shí)模型的信噪比SNR為33.2153 dB、結(jié)構(gòu)相似性指數(shù)為0.8682，均顯著低于所提出新方法的反演結(jié)果。此外，如圖5d所示的反演結(jié)果表明傳統(tǒng)多組變異差分進(jìn)化算法的反演結(jié)果橫向連續(xù)性差。具體來說，該反演阻抗模型充滿噪聲，真實(shí)模型中橫向連續(xù)的地質(zhì)結(jié)構(gòu)在該模型中呈不連續(xù)狀態(tài)。

綜上所述，基于Marmousi II模型的實(shí)驗(yàn)算例驗(yàn)證了所提出理論方法的有效性，并驗(yàn)證了其相比于傳統(tǒng)方法在效率和反演精度上的優(yōu)勢。

3 實(shí)際資料應(yīng)用

本節(jié)將上述理論方法應(yīng)用于勝利油田某工區(qū)實(shí)際三維地震資料。實(shí)驗(yàn)中采用整個(gè)三維工區(qū)的部分?jǐn)?shù)據(jù)，這部分?jǐn)?shù)據(jù)的范圍是Inline 7080～7600，Crossline 3180～3450，時(shí)間方向1400～2450 ms。該數(shù)據(jù)的時(shí)間采樣間隔為2 ms。目的層位于三角洲前緣，發(fā)育大量濁積巖砂體，具有良好的油氣顯示。圖6a 所示為目標(biāo)區(qū)連井地震剖面。目的層的主力儲(chǔ)層砂體發(fā)育在T1、T2和T3層位之間，為典型的薄互層結(jié)構(gòu)，砂層厚度薄單層在2～10 m，組合厚度4～20 m之間。此外，目的層砂體橫向相變快，僅憑地震反射數(shù)據(jù)難以落實(shí)砂體展布邊界。從地震波形上看(圖6a)，由于薄互層干涉效應(yīng)導(dǎo)致內(nèi)部同相軸能量弱，橫向連續(xù)性差，難以追蹤。

采用本文提出的阻抗反演方法對(duì)該工區(qū)三維地震數(shù)據(jù)進(jìn)行處理。這里，反演使用的子波是通過統(tǒng)計(jì)的方法從地震數(shù)據(jù)中提取。圖6b為采用本文所提出方法計(jì)算的波阻抗連井剖面。根據(jù)測井解釋成果，我們知道目標(biāo)層段砂巖表現(xiàn)為阻抗高值，泥巖表現(xiàn)為阻抗低值。井點(diǎn)位置黑色曲線展示的是利用測井?dāng)?shù)據(jù)計(jì)算的波阻抗曲線。我們也將測井巖性解釋結(jié)果覆蓋到井點(diǎn)處。從巖性解釋結(jié)果中不難看出，T1～T2層段發(fā)育厚度非常薄的砂體，在地震剖面上表現(xiàn)為不連續(xù)的同相軸(圖6a)，難以分辨薄砂體的組合關(guān)系與橫向分布范圍。利用本文所提出方法計(jì)算的波阻抗結(jié)果雖然無法刻畫T1～T2層段之間5 m以下的薄層砂體(圖6b中紅色箭頭指示)，但對(duì)于厚度較大的砂體則能清晰地展示砂體橫向展布范圍，并且與井?dāng)?shù)據(jù)計(jì)算的波阻抗有良好的匹配關(guān)系(圖6b中藍(lán)色箭頭指示)。T2～T3層段內(nèi)部發(fā)育組合厚度較厚的砂體，采用本文方法計(jì)算的波阻抗能夠很好地表征其內(nèi)部的砂體展布規(guī)律，且在井點(diǎn)處能很好地匹配測井計(jì)算結(jié)果。圖7為本文方法計(jì)算結(jié)果的三維展示。為了更好的對(duì)比效果，我們沿圖6b所示的T2_1層位做沿層切片，得到圖8a所示的沿層阻抗。進(jìn)一步與同層位井插值得到的砂地比(圖8b)進(jìn)行對(duì)比，不難看出高阻抗區(qū)域?qū)?yīng)高砂地比，低阻抗區(qū)域?qū)?yīng)低砂地比。這進(jìn)一步驗(yàn)證了本文方法的有效性。

a—地震剖面；b—波阻抗剖面；井點(diǎn)位置處黑色曲線為利用測井?dāng)?shù)據(jù)計(jì)算的波阻抗

a—?jiǎng)倮吞锬彻^(qū)三維疊后地震資料；b—所提出理論方法反演的波阻抗數(shù)據(jù)體

a—反演波阻抗模型沿T2_1層位的切片；b—井插值砂地比剖面

4 結(jié)論及討論

1) 帶有橫向約束的模型搜索空間設(shè)定方法在本文的全局優(yōu)化波阻抗反演中起到了保證反演結(jié)果橫向連續(xù)性和提升反演效率的效果。

2) 帶有橫向約束的全局優(yōu)化波阻抗反演方法進(jìn)行了實(shí)際應(yīng)用，反演的波阻抗能夠很好地表征儲(chǔ)層內(nèi)部的砂體展布規(guī)律，且在井點(diǎn)處與測井結(jié)果匹配度高。驗(yàn)證了本研究的實(shí)用性。

3) 進(jìn)一步提升全局優(yōu)化波阻抗反演方法的效率，是下一步研究的方向。