二階同步提取變換的沉積旋回界面定位與追蹤識(shí)別

李雪英, 王 鑫

(1.東北石油大學(xué) 地球科學(xué)學(xué)院, 黑龍江 大慶 163318; 2.黑龍江省油氣藏形成機(jī)理與資源評(píng)價(jià)重點(diǎn)實(shí)驗(yàn)室, 黑龍江 大慶 163318)

0 引 言

沉積旋回研究是層序地層學(xué)中的一項(xiàng)重要工作,諸多學(xué)者已經(jīng)對(duì)沉積旋回的規(guī)律做了大量研究[1-4]。地球物理觀測(cè)數(shù)據(jù)反映了地層的巖性和物性,包含了大量與沉積旋回有關(guān)的信息[5]。通常情況下,由現(xiàn)場(chǎng)測(cè)井?dāng)?shù)據(jù)進(jìn)行沉積旋回的劃分,顯然,這種方法對(duì)于井點(diǎn)附近的劃分結(jié)果準(zhǔn)確度高,但是對(duì)于無井點(diǎn)地區(qū)的分析,僅能通過已有地質(zhì)規(guī)律推導(dǎo),這樣導(dǎo)致了該地段沉積旋回劃分結(jié)果不盡人意。地震資料對(duì)于地下信息收集比較全面,也沒有井點(diǎn)的制約,這種特性有利于更準(zhǔn)確地研究地層變化。在地下各類地層介質(zhì)中,泥巖與砂巖的占比比較大,也比較常見。而且通過砂巖泥巖重疊形式和厚度變化進(jìn)行劃分,可以作為重要指標(biāo)指導(dǎo)判別沉積體系以及海侵-海退過程,也是使用時(shí)頻分析方法判定沉積旋回模式的重要條件。時(shí)頻分析方法判定沉積旋回模式的依據(jù)為:薄層時(shí)頻響應(yīng)機(jī)理中的升頻降幅作用[6];旋回性薄互層中的小層厚度遞增或遞減變化,會(huì)相應(yīng)地引起瞬時(shí)頻譜峰值頻率逐漸減小或增大。

通過時(shí)頻分析方法識(shí)別沉積旋回的最早研究可以追溯到Бродов等提出的СФИ技術(shù)[7]。此后,加窗傅立葉變換[8]、小波變換[9-10]、廣義S變換[11]和希爾伯特-黃變換[12]等已被開發(fā)用于確定沉積旋回。利用時(shí)頻分析方法確定沉積旋回模式的過程中,越高的時(shí)頻分辨率越有利于旋回模式的判別。為了實(shí)現(xiàn)理想的時(shí)頻分析目標(biāo),一些先進(jìn)的時(shí)頻后處理方法成為近年來研究的重點(diǎn)。它主要包括RS[13]、SST[14]以及SET[15]等。經(jīng)過推導(dǎo)與實(shí)驗(yàn)證明,二階同步提取變換方法可以產(chǎn)生比RS、SST和SET方法能量聚焦更好的時(shí)頻分析結(jié)果[16-17],這一特點(diǎn)非常符合沉積旋回判別的需求。

綜上所述,筆者建立多套沉積旋回模型,以模擬地下各種情況下的組合形式,運(yùn)用二階同步提取變換對(duì)該模型的地震道數(shù)據(jù)進(jìn)行時(shí)頻分析處理,探究沉積旋回模型中,各單元頂界面附近的時(shí)頻特征,給出沉積旋回頂界面的識(shí)別與追蹤方法。

1 二階同步提取變換原理

二階同步提取變換(SET_2)是基于短時(shí)Fourier變換的時(shí)頻分析后處理方法,其是對(duì)同步提取變換進(jìn)行改良的算法;其主要思想是:首先對(duì)信號(hào)進(jìn)行短時(shí)傅里葉變換分析,得到信號(hào)的時(shí)頻譜,之后利用求取的二階瞬時(shí)頻率和delta函數(shù),求得二階同步提取算子(SEO_2),最后將算子與短時(shí)傅里葉變換結(jié)果相乘,求得二階同步提取變換時(shí)頻譜[20]。該思路的主要公式如下。

短時(shí)傅里葉變換公式為

式中:t——時(shí)間;

x——積分的自變量;

g(x-t)——高斯窗函數(shù);

y(x)——輸入信號(hào);

Ge(t,ω)——短時(shí)傅里葉變換結(jié)果。

二階估算頻率求解公式[18]為

將delta函數(shù)與二階估算頻率結(jié)合,得到二階同步提取算子,則二階同步提取變換的表達(dá)式為

Te(t,ω)=Ge(t,ω)×δ(ω-ω1(t,ω)),

式中,δ——二階同步提取算子。

二階同步提取算子表達(dá)式為

2 時(shí)頻特征與旋回界面定位方法

設(shè)計(jì)多套不同的沉積旋回單元模型,在泥巖背景下進(jìn)行波動(dòng)方程正演模擬,各模型的旋回單元厚度范圍hr、厚度遞變梯度hg、分割旋回單元的隔層厚度hc等參數(shù),如表1所示。地震子波采用零相位的雷克子波,子波主頻為39 Hz;通過正演模擬獲取單炮紀(jì)錄,并在單炮記錄中提取零偏移距的地震道數(shù)據(jù)進(jìn)行之后的研究。

表1 旋回單元組合模型參數(shù)

經(jīng)過正演模擬之后各模型的時(shí)域波形及短時(shí)傅里葉變換(STFT)、二階同步提取變換(SET2)結(jié)果,如圖1～3所示。當(dāng)模型中各旋回單元的厚度大,而且旋回單元中各小層厚度變化劇烈時(shí),時(shí)域波形振幅在較厚小層處比較大,且波形較為稀疏,通過這個(gè)特征可以劃分旋回模式及其厚度。當(dāng)模型中各旋回單元厚度較小時(shí)(由模型c所示,其小層厚度均為入射波的1/8波長以下,均小于地震分辨率[19]),通過各模型的時(shí)域波形無法識(shí)別旋回特征。當(dāng)模型中個(gè)旋回單元厚度較大,但是各小層厚度變化不大時(shí)(模型a),各時(shí)域波形在厚度較大的小層處振幅大,但是從整體來看,波形稀疏程度趨于一致,僅在旋回模型的頂界面附近出現(xiàn)強(qiáng)反射尖峰。

圖1 模型a、b時(shí)頻特征Fig.1 Time-frequency characterization of models a and b

圖2 模型c、d時(shí)頻特征Fig.2 Time-frequency characterization of models c and d

圖3 模型e、f時(shí)頻特征Fig.3 Time-frequency characterization of models e and f

綜合分析短時(shí)傅里葉變換結(jié)果和二階同步提取變換結(jié)果可知,二階同步提取變換得到的時(shí)頻譜能量更加聚焦,具有較高的時(shí)頻分辨率。從圖1～3可以看出,二階同步提取變換頻譜圖中的每個(gè)頻譜變化單元都對(duì)應(yīng)著一個(gè)沉積旋回單元,而且對(duì)于不同旋回單元厚度和不同旋回單元厚度變化梯度,二階同步提取變換都可以清晰反映出時(shí)頻譜中高頻到低頻的變化特征,從這個(gè)方面可以對(duì)旋回單元的模式進(jìn)行判別。短時(shí)傅里葉變換頻譜圖中在旋回單元較厚時(shí)頻譜變化能夠?qū)?yīng)沉積旋回單元,然而當(dāng)旋回單元較薄時(shí),無法準(zhǔn)確判斷旋回單元。兩種方法得到的時(shí)頻譜都表現(xiàn)出了向模型中小層層厚增大的方向時(shí)移的特點(diǎn)。

通過時(shí)頻分析結(jié)果分析總結(jié)沉積旋回單元頂界面的特征,當(dāng)模型中沉積旋回單元的總厚度較大、各小層的厚度變化比較明顯時(shí),二階同步提取變換的高頻極值對(duì)應(yīng)著旋回單元的頂界面(模型d、e、f中的SET),當(dāng)厚度變化梯度比較小時(shí),則需要結(jié)合二階同步提取算子進(jìn)行輔助判別,可以看到在頂界面附近,二階同步提取算子呈現(xiàn)出高頻突變現(xiàn)象,其拱形形態(tài)的尖峰對(duì)應(yīng)著沉積旋回單元的頂界面。

薄層具有升頻降幅的作用[6],在圖1表示的沉積旋回模型中,由于薄層的升頻作用,可以看出頻率最高的位置總是指向薄層。薄層降幅作用導(dǎo)致薄層產(chǎn)生的高頻信號(hào)能量弱,二階同步提取變換對(duì)于弱信號(hào)具有較強(qiáng)的提取能力,可以很好地提取薄層產(chǎn)生的高頻成分。由模型CSEO2可以看出,在旋回單元厚度階躍的位置,厚層產(chǎn)生的強(qiáng)能量的低頻信息及薄層產(chǎn)生的弱能量的高頻信息均得到了清晰的刻畫。在識(shí)別這類旋回單元時(shí),需要結(jié)合二階同步于是取變換的時(shí)頻譜和二階同步提取算子的結(jié)果進(jìn)行判斷,最終得出頂界面的位置。

綜上分析可知,在二階同步提取變換得到的時(shí)頻譜中,沉積旋回單元的頂界面位置對(duì)應(yīng)著時(shí)頻譜的高頻極值部分,可以根據(jù)這一性質(zhì)對(duì)沉積旋回單元進(jìn)行判別與劃分,進(jìn)而預(yù)測(cè)沉積旋回的厚度。

3 沉積旋回的噪聲魯棒性的判別

在實(shí)際的地震信號(hào)中,必定會(huì)有噪聲的干擾,為了使文中提出的方法具有更廣泛的適用性,探究二階同步提取變換方法判斷沉積旋回單元的抗噪聲能力,選取一種前文建立的地質(zhì)模型,加入一定能量的高斯白噪聲使得該模型能夠更好的模擬一般情況下底層情況,使得經(jīng)驗(yàn)更具有一般性。

由前文可知,時(shí)域波形隨著噪聲能量的增大失去了本身的特征。對(duì)于二階同步提取算子計(jì)算結(jié)果來說,該結(jié)果受噪聲影響較大,但是通過綜合比對(duì)二階同步提取變換時(shí)頻譜可知,算子計(jì)算結(jié)果仍然保留了頂界面處的高頻突變現(xiàn)象,呈現(xiàn)為規(guī)則的拱形形態(tài),并且該形態(tài)與噪聲的頻率分量有明顯的差別,如圖4所示,拱形形態(tài)在70 Hz左右。而從二階同步提取變換的時(shí)頻譜來看,加入各強(qiáng)度的噪聲之后,噪聲能量在時(shí)頻譜中有一定的響應(yīng),但對(duì)于反映沉積旋回模型中各單元小層厚度變化的頻譜能量并沒有很大的干擾,如圖5所示。

圖4 二階同步提取算子抗噪能力試驗(yàn)Fig.4 Noise immunity test of second-order synchronous extraction operator

圖5 二階同步提取變換抗噪能力試驗(yàn)Fig.5 Second-order simultaneous extraction transform noise immunity test

而實(shí)際地震探測(cè)中,噪聲能量遠(yuǎn)不及實(shí)驗(yàn)中所加最大能量。綜上所述,通過綜合二階同步提取變換以及二階同步提取算子的計(jì)算結(jié)果,能夠在噪聲的干擾下較為準(zhǔn)確的識(shí)別沉積旋回單元頂界面位置。

4 沉積旋回界面追蹤方法

綜合上述分析,沉積旋回單元頂界面在二階同步提取變換時(shí)頻譜中,必然表現(xiàn)出高頻極值的特點(diǎn),則可以通過這個(gè)高頻極值來判別沉積旋回單元的頂界面位置。為了驗(yàn)證這個(gè)結(jié)論的合理性以及探究頂界面追蹤方法,設(shè)計(jì)一個(gè)由三個(gè)正旋回單元組合的地質(zhì)模型(圖6),地質(zhì)模型部分參數(shù),如表2所示。

圖6 三個(gè)正旋回單元組合厚度模型Fig.6 Combined thickness model of three positive cycle elements

表2 旋回單元組合地質(zhì)模型部分參數(shù)

利用軟件制作該模型的合成地震記錄剖面(圖7),然后使用二階同步提取變換對(duì)正演結(jié)果的各地震道數(shù)據(jù)進(jìn)行處理。

圖7 由圖6模型正演的地震剖面局部放大Fig.7 Local enlarged view of seismic profile modeled by model in figure 6

4.1 局限性分析

由圖7可知,時(shí)域波形可以反映出沉積旋回單元的大致位置,當(dāng)夾層厚度較大時(shí),沉積旋回單元能夠顯示明顯的分離界面,從時(shí)域波形的特征來看,能夠明顯看出沉積旋回單元的位置,但是單元界面與波形的峰谷對(duì)應(yīng)程度不高;當(dāng)夾隔層厚度較小時(shí),沉積旋回單元界面的時(shí)域波形的特征不明顯,無法準(zhǔn)確判斷各沉積旋回單元位置;另外時(shí)域波形易受噪聲等因素的影響,很難指定精準(zhǔn)的判別標(biāo)準(zhǔn)。由此,應(yīng)通過時(shí)頻分析結(jié)果結(jié)合時(shí)域波形綜合討論。

由圖7地震道1～130可見,沉積旋回單元厚度不小于波長的一倍時(shí):從時(shí)域波形來看,小層厚度較大的位置波形比較稀疏,反之小層厚度較小的位置波形比較密集。如圖8中地震道波形98～131所示。時(shí)頻譜中變化趨勢(shì)清晰反映了沉積旋回單元的厚度變化趨勢(shì),并且各單元的頂界面均對(duì)應(yīng)頻率的極大值處,從時(shí)頻譜中可準(zhǔn)確判斷沉積旋回單元頂界面位置及沉積旋回模式。

圖8 圖7剖面不同地震道的二階同步提取變換結(jié)果及自動(dòng)追蹤位置與實(shí)際位置對(duì)比Fig.8 Second order synchronous extraction and transformation results of different seismic tracks in profile in fig.7 and comparison between automatic tracking position and actual position

如圖7中地震道130～170可見,沉積旋回單元厚度在一半到一倍波長之間時(shí):從時(shí)域波形來看,沉積旋回單元的波形特征仍比較清晰,但在圖8中對(duì)應(yīng)的時(shí)頻譜中,頂界面的對(duì)應(yīng)效果不佳,結(jié)合圖6中的二階同步提取算子計(jì)算結(jié)果則能夠很好的矯正出頂界面的位置。如圖7中地震道170～200所示,當(dāng)沉積旋回單元厚度不大于一半波長時(shí):時(shí)域波形趨于對(duì)稱,頂?shù)捉缑娓浇嬖诿黠@的強(qiáng)能量振幅;時(shí)頻譜無法反應(yīng)旋回單元頻率變化趨勢(shì),但薄層位置依然對(duì)應(yīng)高頻極值,此時(shí),應(yīng)先根據(jù)時(shí)域波形限定旋回單元頂界面識(shí)別范圍,進(jìn)而追蹤時(shí)頻譜頻率極大值位置定位旋回單元頂界面。

對(duì)于夾層來說,其在模型中起到了分割各沉積旋回單元的作用,由于夾層的存在且夾層厚度對(duì)于界面附近的頻譜有一定的影響,需要注意夾層對(duì)于追蹤結(jié)果的干擾。當(dāng)沉積旋回單元的厚度大于波長的一半時(shí),夾層對(duì)于判別旋回模式及識(shí)別頂界面沒有影響。當(dāng)沉積旋回單元厚度小于波長的一半并且夾層厚度大于波長的一半時(shí),沉積旋回單元可以通過時(shí)域波形分辨;當(dāng)夾層厚度小于波長的一半時(shí),則需要直接從時(shí)頻譜中進(jìn)行判別。

4.2 沉積旋回單元頂界面追蹤實(shí)驗(yàn)

根據(jù)上文總結(jié)的結(jié)論,制定出沉積旋回單元頂界面的追蹤方法,如圖9所示。分析時(shí)域波形的頻率范圍,給定閾值進(jìn)行掃描;根據(jù)掃描范圍內(nèi)的時(shí)頻譜頻率極大值的位置,對(duì)沉積旋回頂界面進(jìn)行初次識(shí)別追蹤;將追蹤結(jié)果結(jié)合地震剖面信息進(jìn)行分析,對(duì)給定誤差閾值的追蹤結(jié)果結(jié)合二階同步提取算計(jì)計(jì)算結(jié)果進(jìn)行二次矯正,最終得出較為準(zhǔn)確的沉積旋回單元頂界面位置。

圖9 沉積旋回頂界面識(shí)別方法流程Fig.9 Sedimentation spinning back top interface identification method process

文中使用前文所建立模型進(jìn)行實(shí)驗(yàn)驗(yàn)證。首先根據(jù)時(shí)域波形指定掃面范圍為時(shí)頻譜中大于最大振幅30%的值。首次識(shí)別結(jié)果在圖8中以虛線劃出,在圖10b的地震剖面上以實(shí)現(xiàn)劃出,可以看出識(shí)別結(jié)果能夠較為準(zhǔn)確的反映出沉積旋回單元的頂界面。但也可以在地震剖面圖上看到由大幅度突變點(diǎn),結(jié)合二階同步提取算子計(jì)算結(jié)果(圖10)對(duì)這些突變點(diǎn)附近的地震道頂界面進(jìn)行二次矯正(圖10中虛線所示),第二次矯正結(jié)果在圖11b中以點(diǎn)劃線表示,可以看出經(jīng)過二次矯正之后,突變點(diǎn)附近的識(shí)別結(jié)果更貼近于真實(shí)情況。

圖10 圖7剖面二階同步提取算子計(jì)算結(jié)果及位置對(duì)比Fig.10 Calculation results of second order synchronous extraction operators in profile in fig.7 and comparison of tracking positions

受噪聲干擾下的合成地震剖面情況,如圖11a所示。實(shí)驗(yàn)中加入的噪聲強(qiáng)度為10 dB。從圖中的識(shí)別追蹤結(jié)果來看,受噪聲干擾,追蹤結(jié)果圍繞真實(shí)界面有一定的波動(dòng)情況,但仍然能較清晰的看出界面位置,經(jīng)過二次矯正之后,識(shí)別結(jié)果更加貼近真實(shí)界面位置,相對(duì)誤差基本能夠保持在30%以內(nèi),追蹤結(jié)果與實(shí)際位置基本吻合,如圖11a中點(diǎn)劃線所示。

綜上,文中提出的方法具有較好的應(yīng)用價(jià)值。

5 結(jié) 論

(1)當(dāng)沉積旋回單元的厚度大于波長的一半時(shí),夾層對(duì)于判別旋回模式及識(shí)別頂界面沒有影響。當(dāng)沉積旋回單元厚度小于波長的一半并且夾層厚度大于波長的一半時(shí),沉積旋回單元可以通過時(shí)域波形分辨;當(dāng)夾層厚度小于波長的一半時(shí),則需要直接從時(shí)頻譜中進(jìn)行判別。

(2)二階同步提取變換經(jīng)由同步提取變換改進(jìn)而得,該方法具有更高的時(shí)頻分辨率,具有良好的抗干擾性能,魯棒性更好。能夠更好的提取弱能量信號(hào),對(duì)于沉積旋回單元頂界面的薄層信號(hào)具有比一般時(shí)頻分析方法更好的識(shí)別能力,提升了頂界面的定位精度。

(3)二階同步提取變換計(jì)算的時(shí)頻譜在沉積旋回單元頂界面位置出現(xiàn)穩(wěn)定的頻率極值,或出現(xiàn)一個(gè)具有穩(wěn)定窄帶拱形特征的頻率極值;當(dāng)二階同步提取變換計(jì)算的時(shí)頻譜在旋回單元頂界面頻率極值能量弱無法識(shí)別時(shí),結(jié)合二階同步提取算子計(jì)算結(jié)果進(jìn)行定位調(diào)整將提升定位精度。

(4)根據(jù)沉積旋回單元頂界面與同步提取變換結(jié)果頻率極值之間的穩(wěn)定關(guān)系,建立旋回單元頂界面追蹤方法;通過模型結(jié)果分析,討論方法的局限性,結(jié)合二階同步提取算子計(jì)算結(jié)果建立人機(jī)交互的反饋調(diào)節(jié)機(jī)制,通過沉積旋回模型頂界面追蹤實(shí)驗(yàn)證明了方法的可行性與實(shí)用性。