紀(jì)永禎,張渝悅,朱立華,林正良

(中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

地震反射系數(shù)反演在推斷地下構(gòu)造和展現(xiàn)地震剖面中更多細(xì)節(jié)等方面扮演著重要角色,也得到了廣泛研究[1-4]。傳統(tǒng)的提高分辨率手段,如反褶積和譜白化方法,只能通過增強或恢復(fù)地震主頻帶范圍內(nèi)的信息來獲得帶限的反射系數(shù),因而并不能從根本上拓寬地震數(shù)據(jù)的頻帶寬度以獲得脈沖反射系數(shù),難以準(zhǔn)確描述地層邊界和一些特殊地質(zhì)體,如薄儲層、透鏡體、單砂體以及地層巖性尖滅位置等[5-6]?；诜囱莸姆瓷湎禂?shù)估計方法則可以從帶限地震數(shù)據(jù)中重構(gòu)地震頻帶以外的反射系數(shù)的頻率成分,獲得寬頻帶的脈沖反射系數(shù)來提高地震數(shù)據(jù)分辨率和地震數(shù)據(jù)解釋的精度[7]。

隨機稀疏反射系數(shù)反演是一種反射系數(shù)反演方法,其思路是在地層反射系數(shù)稀疏的假設(shè)下,將反射系數(shù)的估計過程看作是兩個部分的綜合體:第一部分是一個非線性問題,即非零反射系數(shù)位置的確定問題;第二部分是一個線性問題,即確定位置后的非零反射系數(shù)的幅值求取問題[8]。這種反射系數(shù)反演思路可以被看成是最嚴(yán)格意義上的稀疏反射系數(shù)反演[9],并且,采用該思路的反演方法具有可以使用不同的初始種子進行多次運算以減少反演的不確定性的優(yōu)點[8]。獲得的脈沖反射系數(shù)可以有效拓展地震數(shù)據(jù)的頻譜,對噪聲和子波的敏感性較小,能夠較好地保護弱反射信號,結(jié)合寬頻帶子波等手段可獲得高分辨率地震數(shù)據(jù)[8-10]。但是,此類方法均在單道模型下實施,且存在以下兩個問題。①反演的預(yù)設(shè)參數(shù),包括正則化因子和非零反射系數(shù)個數(shù),通常根據(jù)經(jīng)驗人工選定,不合理的參數(shù)選擇通常會引發(fā)過擬合問題[11]。VELIS[9]詳細(xì)分析了非零反射系數(shù)個數(shù)的選取方法和非零反射系數(shù)個數(shù)與正則化因子的相互關(guān)系,發(fā)現(xiàn)非零反射系數(shù)的個數(shù)設(shè)定在很大程度上會影響反演的最終結(jié)果,不僅會影響非零反射系數(shù)的位置,也會影響非零反射系數(shù)的振幅估計,建議謹(jǐn)慎選取非零反射系數(shù)的個數(shù)預(yù)設(shè)值,并且在反演前進行大量實驗。WANG等[10]在頻率域引入有限更新率理論,改善了選取預(yù)設(shè)非零反射系數(shù)的穩(wěn)定性,但未進行統(tǒng)計性分析,證據(jù)略顯不足。SEN等[12]將預(yù)設(shè)非零反射系數(shù)個數(shù)當(dāng)作一個未知數(shù),在估計反射系數(shù)的同時估計非零反射系數(shù)個數(shù),但這種策略使得反演的未知數(shù)維度增加,成本增大。②在反演過程中,算法的收斂取決于地震數(shù)據(jù)和反演反射系數(shù)的正演模擬數(shù)據(jù)之間的匹配程度,具有較大幅值的反射系數(shù)會對反演的趨勢產(chǎn)生巨大影響。因此,振幅較大的異常值對此類反演的結(jié)果危害很大[11,13],而基于單道模型的隨機稀疏反射系數(shù)反演方法對異常值影響問題未給出合適的解決辦法。

為了解決上述兩個問題,假設(shè)反射系數(shù)稀疏和橫向連續(xù),引入相鄰道數(shù)據(jù)作為約束條件,提出了多道隨機稀疏反射系數(shù)反演模型,進而構(gòu)建了一種多道隨機稀疏反射系數(shù)反演方法。該方法通過對目標(biāo)地震道數(shù)據(jù)的非零反射系數(shù)的位置和傾角信息進行同時非線性搜索,在一定程度上解決了前文分析的第一個問題,有效提升了反演預(yù)設(shè)參數(shù)的選擇穩(wěn)定性。與此同時,由于采用了多道模型,可以有效緩解地震數(shù)據(jù)中存在的少數(shù)異常值對反演的影響,解決了第二個問題。作為單道隨機稀疏反射系數(shù)反演方法的改進,該方法保留了進行多次反演實現(xiàn)的能力和對估計出的反射系數(shù)的不確定性進行分析的能力。

1 方法原理

疊后地震數(shù)據(jù)可以描述為反射系數(shù)與地震子波的褶積[14]。假設(shè)子波在某一時間和深度范圍內(nèi)是穩(wěn)態(tài)的,那么褶積模型可以寫成如下形式:

s(l)(t)=w(t)*r(l)(t)+n(l)(t)

l=1,2,…,Lt=1,2,…,N

(1)

(2)

將公式(2)代入公式(1)中,可以獲得:

s(l)(t)=w(t)*r(l)(t)+n(l)(t)

(3)

(4)

(5)

其中,floor(x)代表向下取整函數(shù),即取不超過實數(shù)x的最大整數(shù)。公式(5)的作用是利用目標(biāo)地震道數(shù)據(jù)中非零反射系數(shù)的位置計算出相鄰地震道數(shù)據(jù)中非零反射系數(shù)的位置。公式(4)用矩陣形式表示為:

(6)

通過最小化地震數(shù)據(jù)與正演計算獲得數(shù)據(jù)的殘差進行最優(yōu)化估計,構(gòu)建反演目標(biāo)函數(shù):

(7)

Δl=…,-1,0,1

求解過程中通過非線性搜索非零反射系數(shù)的位置和斜率,并利用最小二乘方法獲得反射系數(shù)的估計,公式(7)對應(yīng)的是常規(guī)最小二乘形式:

(8)

由(8)式可以獲得:

(9)

由于地震數(shù)據(jù)受噪聲影響,為使反演結(jié)果更加穩(wěn)定,反演時需要增加正則化約束。本文采用了二次型正則化約束[8],則公式(9)改寫為:

(10)

與公式(9)對應(yīng)的、加入了正則化約束后的最終目標(biāo)函數(shù)可以寫為:

(11)

1) 選定預(yù)設(shè)非零反射系數(shù)的數(shù)值、同時反演道數(shù)、最大迭代次數(shù)和終止迭代閾值。

2) 初始化時間延遲和斜率,計算當(dāng)前時間延遲和斜率下的相鄰地震道數(shù)據(jù)的時間延遲。

3) 在模擬退火算法的每一次迭代中進行以下計算:①更新目標(biāo)地震道數(shù)據(jù)的時間延遲和斜率,計算更新的時間延遲和斜率下的相鄰地震道數(shù)據(jù)的時間延遲;②利用公式(10)更新A和振幅;③利用公式(11)計算Ja;④通過模擬退火算法中的接受準(zhǔn)則,決定是否接受更新后的時間延遲和斜率;⑤利用公式(7)計算J,如果J大于終止迭代閾值或迭代次數(shù)還未達(dá)到最大迭代次數(shù),回到步驟①。

4) 在步驟3)收斂后,利用時間延遲和振幅,根據(jù)公式(1)重構(gòu)反射系數(shù)序列。

2 模擬數(shù)據(jù)測試

困擾單道隨機稀疏反演方法的問題之一是在反演前需要預(yù)設(shè)非零反射系數(shù)的個數(shù),而預(yù)設(shè)的個數(shù)只能通過人為估計和預(yù)測來確定,并且需要進行反復(fù)試驗才可獲得該個數(shù)的大概估計。而這個參數(shù)的選取會直接影響反演結(jié)果,并且會影響正則化因子的作用效果。一個相對大的預(yù)設(shè)參數(shù)會降低反演的分辨率,并出現(xiàn)高頻的毛刺狀的偽反射系數(shù)。本文提出的多道隨機稀疏反射系數(shù)反演可以有效緩解這個問題。為了對多次獨立反演結(jié)果進行統(tǒng)計分析,排除其它干擾,本文設(shè)計了一個簡單的反射系數(shù)模型,該模型包含10道地震數(shù)據(jù),每一道中含有3個非零反射系數(shù),如圖1a所示。將反射系數(shù)與采樣間隔為1ms、主頻為45Hz的雷克子波褶積,并加入隨機噪聲,生成了信噪比為5的模擬數(shù)據(jù),如圖1b所示。隨機設(shè)定模擬退火算法中的初始種子值,進行了50次獨立的反演試驗,并統(tǒng)計分析每一次的反演結(jié)果中的非零反射系數(shù)總數(shù)。很顯然,在理想情況下,每一次的反演結(jié)果中非零反射系數(shù)的理論值應(yīng)該是3×10=30個。然而,由于地震數(shù)據(jù)中存在隨機噪聲干擾,且預(yù)設(shè)的非零反射系數(shù)個數(shù)大于真實值,使得反演結(jié)果中的非零反射系數(shù)個數(shù)必定大于理論值。圖2為50次獨立反演獲得的非零反射系數(shù)個數(shù)的統(tǒng)計結(jié)果。在不同的預(yù)設(shè)K條件下,對比單道隨機稀疏反射系數(shù)反演和多道隨機稀疏反射系數(shù)反演的結(jié)果發(fā)現(xiàn),多道隨機稀疏反射系數(shù)反演估計的反射系數(shù)個數(shù)更加接近理論值,證明多道隨機稀疏反射系數(shù)反演的壓制噪聲和抑制高頻的毛刺狀反射系數(shù)的能力比單道隨機稀疏反射系數(shù)反演的更強。

在圖1所示的模擬含噪數(shù)據(jù)中加入異常值干擾,得到另一種模擬數(shù)據(jù)(圖3a)。常規(guī)單道隨機稀疏反演方法力求用最少的反射系數(shù)個數(shù)達(dá)到最好的數(shù)據(jù)匹配效果,因此振幅較大的異常值往往會被視為有效值而保留下來,對反演結(jié)果造成負(fù)面影響,如圖3c所示。而圖3b中的多道隨機稀疏反演獲得的結(jié)果沒有受奇異值影響,較好地還原了模型反射系數(shù)特征,展示了本文方法可以降低異常值干擾的特點。圖4給出了一個稍復(fù)雜并含有界面尖滅的模擬數(shù)據(jù)試驗結(jié)果,選取的同時反演道個數(shù)為3,即Δl∈[-1,0,1],預(yù)設(shè)非零反射系數(shù)個數(shù)參數(shù)K=12,正則化參數(shù)λ=0.1。對比單道隨機稀疏反射系數(shù)反演和多道隨機稀疏反射系數(shù)反演結(jié)果中紅框區(qū)域可以看出,多道隨機稀疏反射系數(shù)反演結(jié)果保留了識別較薄界面的能力,且對噪聲的抵抗能力更強,反演結(jié)果更加穩(wěn)定。

圖1 反射系數(shù)模型(a)及含噪地震數(shù)據(jù)(b)

圖2 50次獨立反演獲得的非零反射系數(shù)個數(shù)的統(tǒng)計結(jié)果

圖3 含噪聲和異常值的地震數(shù)據(jù)和反演結(jié)果a 含噪聲和異常值的地震數(shù)據(jù); b 多道隨機稀疏反射系數(shù)反演結(jié)果; c 單道隨機稀疏反射系數(shù)反演結(jié)果

圖4 模擬數(shù)據(jù)反演結(jié)果a 反射系數(shù)模型; b 含噪聲地震數(shù)據(jù); c 多道隨機稀疏反射系數(shù)反演結(jié)果; d 單道隨機稀疏反射系數(shù)反演結(jié)果

3 實際數(shù)據(jù)應(yīng)用

實際數(shù)據(jù)來自川東北某工區(qū),該區(qū)地下構(gòu)造較為復(fù)雜,地層橫向變化快,地層構(gòu)造和起伏的描述對后續(xù)解釋工作十分重要。反演所用子波來自井震標(biāo)定結(jié)果。選取同時反演道個數(shù)為3,即Δl∈[-1,0,1],預(yù)設(shè)非零反射系數(shù)個數(shù)參數(shù)K=20,正則化參數(shù)λ=0.25。地震數(shù)據(jù)如圖5a所示。圖5b和圖5c分別為多道隨機稀疏反射系數(shù)反演和單道隨機稀疏反射系數(shù)反演得到的高分辨率地震剖面,獲取過程是利用寬頻帶子波與反演反射系數(shù)結(jié)果進行褶積。對比圖5b和圖5c可以看出,多道隨機稀疏反射系數(shù)反演結(jié)果中斷層和構(gòu)造的變化展現(xiàn)得更加清晰,振幅較強同相軸間的信息得以展示,且剖面背景干凈,同相軸疊置關(guān)系清晰,證實了多道隨機稀疏反射系數(shù)反演在體現(xiàn)界面細(xì)節(jié)、斷層和反演結(jié)果穩(wěn)定程度上均具有一定的優(yōu)勢。實際數(shù)據(jù)的應(yīng)用體現(xiàn)了本文方法在實際數(shù)據(jù)中的應(yīng)用潛力和優(yōu)勢及其對改善地震解釋精度的意義。

圖5 實際地震數(shù)據(jù)及不同方法的反演結(jié)果a 實際地震數(shù)據(jù); b 采用多道隨機稀疏反射系數(shù)反演得到的高分辨率地震剖面; c 采用單道隨機稀疏反射系數(shù)反演得到的高分辨率地震剖面

4 結(jié)論

本文提出了一種多道隨機稀疏反射系數(shù)反演方法,通過同時非線性搜索反射系數(shù)位置和斜率,將多道反射系數(shù)褶積模型和單道隨機反射系數(shù)反演方法相結(jié)合,該方法可以用于地震數(shù)據(jù)的高分辨率處理和反演。多道隨機稀疏反射系數(shù)反演方法保留了單道隨機稀疏反演方法識別薄層的能力和進行不確定性分析的優(yōu)勢,并提高了選取預(yù)設(shè)參數(shù)的穩(wěn)定性和抵抗奇異值噪聲的能力。模擬數(shù)據(jù)測試結(jié)果證明了本文方法的有效性和對單道隨機稀疏反射系數(shù)反演的改進,其中的統(tǒng)計性試驗證明了方法對參數(shù)選擇穩(wěn)定性的提高。實際地震數(shù)據(jù)的應(yīng)用揭示了本文方法提高地震資料分辨率和解釋效果的潛力。本文方法仍存在進一步研究的空間,如何加入測井資料約束是下一步研究的方向。