溫志平， 方江雄， 劉 軍， 鄭成龍

(東華理工大學(xué) 地球物理與測控技術(shù)學(xué)院，南昌 330013)

0 引言

地震記錄初至波(P波)走時(shí)拾取,是地震層析正演計(jì)算和靜校正技術(shù)中一項(xiàng)基礎(chǔ)而又重要的工作，對(duì)震源定位、油藏描述、探明地質(zhì)缺陷分布等具有重要意義。初至波走時(shí)的手工拾取已不符合現(xiàn)今大規(guī)模地震勘探數(shù)據(jù)的需求，一般的初至拾取方法能滿足于地震記錄背景較平靜且初至波形簡單的地震資料，但對(duì)于具有較強(qiáng)背景噪聲干擾，初至波形變化較大，各種波相互干擾的地震記錄，往往難以獲得準(zhǔn)確的初至?xí)r間[1]。因此，在運(yùn)用自動(dòng)拾取算法前對(duì)地震隨機(jī)噪聲進(jìn)行壓制處理，是提高地震初至波走時(shí)拾取精度的有效手段。

初至拾取方法主要有三類:①基于地震記錄瞬時(shí)特征的方法[2-4](如極值法、能量比值法和AIC法等)，當(dāng)?shù)卣鹩涗浶旁氡容^低時(shí)，單一方法的使用導(dǎo)致拾取誤差較大；②基于地震記錄整體特征的方法[5-8](如相關(guān)法、線性最小平方預(yù)測法等)，這類方法對(duì)地震隨機(jī)噪聲一定的抑制作用，但受地震道集相關(guān)性等因素影響，在處理地表情況復(fù)雜的地震資料時(shí)效果不理想；③綜合地震多維信息的人工智能方法[9-12](如神經(jīng)網(wǎng)絡(luò)法和分形維數(shù)法等)，該類方法參數(shù)設(shè)置復(fù)雜，拾取效率較低。

地震隨機(jī)噪聲壓制一直是國內(nèi)、外眾多學(xué)者研究的對(duì)象，針對(duì)不同類型的噪聲干擾，提出了不同的噪聲壓制方法。當(dāng)前，針對(duì)地震隨機(jī)噪聲的壓制方法主要有：中值濾波[13]、f-x域預(yù)測濾波[14]、獨(dú)立成分分析[15]、小波變換[16]、經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition, EMD)方法[17-19]等。

完備集合經(jīng)驗(yàn)?zāi)B(tài)分解(Complete Ensemble Empirical Mode Decomposition，CEEMD)，避免了EMD方法中模態(tài)混疊的問題，有效改善集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition，EEMD)方法中加入隨機(jī)白噪聲對(duì)原始信號(hào)造成污染的缺點(diǎn)。地震信號(hào)經(jīng)CEEMD分解得到各級(jí)IMF分量，采用有效指標(biāo)函數(shù)(Effective Index Function，EIF)確定有效信號(hào)和噪聲的能量分界點(diǎn)，對(duì)部分高頻IMF分量再次進(jìn)行小波閾值去噪，并保持低頻IMF分量不變；最后將去噪后的高頻IMF分量和低頻IMF分量以及殘余量進(jìn)行重構(gòu)，即得到去噪后的地震信號(hào)。采用改進(jìn)滑動(dòng)時(shí)窗能量比聯(lián)合AIC準(zhǔn)則法，對(duì)去噪后的地震初至波走時(shí)進(jìn)行拾取。該方法首先根據(jù)改進(jìn)滑動(dòng)時(shí)窗能量比法大致確定初至波時(shí)刻，然后在該時(shí)刻前后取一時(shí)窗，使用AIC方法計(jì)算AR模型在該時(shí)窗范圍內(nèi)的局部極小值，即得到精確的初至波走時(shí)。

1 CEEMD聯(lián)合小波閾值地震噪聲壓制

1.1 CEEMD完備集合經(jīng)驗(yàn)?zāi)B(tài)分解

CEEMD基于EMD和EEMD改進(jìn)而來，其具體步驟如下：

1)向原信號(hào)中添加正負(fù)白噪聲對(duì)，得到兩組模態(tài)集合M1、M2：

(1)

2)EMD分解第i次添加白噪聲后的信號(hào)，得到對(duì)應(yīng)的模態(tài)分量IMFij。

3)最終分解結(jié)果為各個(gè)IMF分量整體取平均：

(2)

CEEMD通過在原地震信號(hào)中自適應(yīng)加入正負(fù)成對(duì)白噪聲的方法，得到成對(duì)的i組j個(gè)分量，再進(jìn)行加總平均將高頻隨機(jī)噪聲或間斷脈沖信號(hào)提出來，以達(dá)到消除模態(tài)混疊的目的，幾乎不存在重構(gòu)誤差，在地震隨機(jī)壓制方面具有較好的應(yīng)用效果。

筆者將CEEMD與小波閾值法相結(jié)合，將CEEMD分解出的部分高階IMF分量再次進(jìn)行小波閾值去噪，然后將去噪后的部分高頻IMF分量與低頻IMF分量以及殘余分量進(jìn)行重構(gòu)。該方法可有效規(guī)避CEEMD方法中直接舍棄高頻IMF分量導(dǎo)致其中有效信息損失的弊端。

1.2 小波閾值去噪

小波閾值去噪首先將原始地震信號(hào)轉(zhuǎn)換至小波域，再做閾值處理，壓制含隨機(jī)噪聲的小波系數(shù)。最后重構(gòu)小波系數(shù)，即得到去噪后的地震信號(hào)。閾值處理有硬閾值和軟閾值法兩種，筆者采用后者，軟閾值函數(shù)的表達(dá)式為：

(3)

1.3 CEEMD聯(lián)合小波閾值地震噪聲壓制

地震信號(hào)經(jīng)CEEMD分解得到各級(jí)IMF分量，采用EIF有效指標(biāo)函數(shù)確定有效信號(hào)和噪聲的能量分界點(diǎn)，對(duì)部分高頻IMF分量再次進(jìn)行小波閾值去噪，并保持低頻IMF分量不變；最后將去噪后的高頻IMF分量和低頻IMF分量以及殘余量進(jìn)行重構(gòu)，即得到去噪后的地震信號(hào)。EIF有效指標(biāo)函數(shù)如式(4)所示。

(4)

式中：N代表信號(hào)長度；IMF(t)為原信號(hào)分解出的各級(jí)模態(tài)分量；EIF越小就表示IMF(t)分量越近似于原始信號(hào)，EIF的較大值被認(rèn)為是噪聲。

CEEMD聯(lián)合小波閾值去噪重構(gòu)即得到去噪后的信號(hào)為式(5)。

rn(t)

(5)

式中：n為EMD分解得到的IMF分量級(jí)數(shù)；m為需要進(jìn)行小波閾值去噪的高頻IMF分量數(shù)目；IMFi(t)為需要進(jìn)行小波閾值去噪的高頻IMF分量；IMFi(t)為不需要處理的低頻IMF分量；rn(t)為EMD分解后的殘余分量。

2 改進(jìn)方法

2.1 改進(jìn)滑動(dòng)時(shí)窗能量比法

滑動(dòng)時(shí)窗能量比法是常用的初至拾取方法，其原理是根據(jù)地震記錄初至波到達(dá)時(shí)，該時(shí)刻的前后時(shí)窗能量比值R(t)達(dá)到最大值，拾取R(t)最大值對(duì)應(yīng)的時(shí)間即為初至?xí)r間。在前、后時(shí)窗能量中加入穩(wěn)定因子，增強(qiáng)其穩(wěn)定性，改進(jìn)公式可表示為式(6)。

(6)

2.2 AIC信息準(zhǔn)則

持續(xù)時(shí)間序列可以截?cái)酁槎鄠€(gè)局部平穩(wěn)部分，其各部分均可表示為自回歸(Auto Recursive，AR)過程。當(dāng)模型階數(shù)固定時(shí)，AIC函數(shù)的最小值位置，可充當(dāng)兩段時(shí)間序列的分界點(diǎn)。對(duì)于一個(gè)給定長度為N的時(shí)間序列[xn]，用自回歸模型可將其表示為式(7)。

(7)

定義第k點(diǎn)為有效地震信號(hào)和隨機(jī)噪聲分量的最佳分界點(diǎn)，原始地震信號(hào)在第k點(diǎn)被截?cái)酁閮啥?，其過程采用式(8)計(jì)算AIC值。

(8)

式中：σ1和σ2則分別表示隨機(jī)噪聲和有效地震信號(hào)分量的方差；C為常數(shù)。

在地震初至波來臨時(shí)波形發(fā)生較大變化，隨機(jī)噪聲和地震信號(hào)具有較大的數(shù)學(xué)統(tǒng)計(jì)差異，兩者擬合度最差，由式(8)求得AIC值最小。因此，計(jì)算地震波的AIC極小值可用于地震初至?xí)r刻的確定。

2.3 改進(jìn)滑動(dòng)時(shí)窗能量比聯(lián)合AIC準(zhǔn)則拾取初至波走時(shí)

單一的初至拾取方法往往無法完全適應(yīng)各種類型的實(shí)際地震信號(hào)，導(dǎo)致某些初至波走時(shí)的拾取錯(cuò)誤。因此，采用改進(jìn)滑動(dòng)時(shí)窗能量比聯(lián)合AIC信息準(zhǔn)則方法，對(duì)去噪后的地震波進(jìn)行初至走時(shí)拾取。該方法根據(jù)改進(jìn)滑動(dòng)時(shí)窗能量比法獲取大致的初至?xí)r刻，在該時(shí)刻前后取一時(shí)窗，使用AIC方法計(jì)算AR模型在該時(shí)窗范圍內(nèi)的局部極小值，即得到精確的初至波走時(shí)。該聯(lián)合方法可有效緩解能量比方法中時(shí)窗大小的設(shè)定對(duì)拾取精度的影響。只在局部范圍內(nèi)計(jì)算AIC值，可有效避免在真實(shí)初至?xí)r刻來臨之前，出現(xiàn)多個(gè)局部AIC極小值，從而導(dǎo)致初至走時(shí)拾取出錯(cuò)的問題。

3 應(yīng)用實(shí)例

3.1 CEEMD聯(lián)合小波閾值地震噪聲壓制應(yīng)用效果

基于Ricker子波人工合成一組單同相軸地震數(shù)據(jù)，采樣點(diǎn)數(shù)為1 000，道數(shù)為40，波速為3 000 m/s，道間距為50 m，采樣頻率為1 000 Hz。

從圖1(c)中可以發(fā)現(xiàn)，EMD去噪后的記錄可以有效地將淹沒在噪聲中的同相軸恢復(fù)出來，具有一定去噪能力，但一些較小的突變?nèi)源嬖?，模態(tài)混疊導(dǎo)致EMD去噪過程不徹底；從圖1(d)可以看出，基于CEEMD小波閾值算法的去噪效果要優(yōu)于EMD，有效地平滑了EMD方法中的高頻毛刺，使得去噪后的結(jié)果與純凈信號(hào)更為接近。圖2為合成單同相軸地震信號(hào)采用不同去噪方法處理后各道SNR分布對(duì)比。表1為各方法去噪后信噪比SNR及均方根誤差RMES的數(shù)值對(duì)比，由圖2、表1可知,本文聯(lián)合去噪方法可有效去除地震信號(hào)中的隨機(jī)噪聲干擾，提升地震信號(hào)信噪比，降低去噪信號(hào)的重構(gòu)誤差。

圖3為不同信噪比的三分量實(shí)際地震資料，采用本文聯(lián)合噪聲壓制方法和小波閾值法去噪后的效果對(duì)比。本文聯(lián)合方法中CEEMD采用的參數(shù)為ceemd(s,0.2,100,10)，小波閾值法中選取sym8小波基，采用軟閾值法。圖3各小圖中：第一行為原始地震信號(hào)，第二行為CEEMD聯(lián)合小波閾值法的去噪效果，第三行為小波閾值法的去噪效果。由圖3實(shí)際去噪效果可知本文CEEMD聯(lián)合小波閾值法對(duì)不同信噪比地震數(shù)據(jù)均具有較好的去噪效果，小波閾值法對(duì)高信噪比地震數(shù)據(jù)去噪效果較好，在低信噪比情況下仍有較多毛刺及干擾噪聲未完全去除，去噪效果不佳。

表2為小波閾值法、EMD法、EEMD法和CEEMD聯(lián)合小波閾值法四種方法去噪后SNR及RMES變化對(duì)比。由表2可知，采用本文CEEMD聯(lián)合小波閾值法去噪后的地震信號(hào)信噪比提升最為明顯，且均方根誤差最小，驗(yàn)證了圖3中去噪效果。

3.2 改進(jìn)滑動(dòng)時(shí)窗能量比聯(lián)合AIC準(zhǔn)則拾取初至波走時(shí)應(yīng)用效果

選取1 000組不同信噪比的合成數(shù)據(jù)，采用不同方法進(jìn)行初至拾取對(duì)比，圖4為各信噪比數(shù)據(jù)對(duì)應(yīng)不同拾取方法拾取誤差分布圖。由圖4可知，當(dāng)信噪比大于22 dB時(shí)，本文改進(jìn)滑動(dòng)時(shí)窗能量比聯(lián)合AIC信息準(zhǔn)則法、STA/LTA、AR-AIC三種方法的初至拾取誤差均較??；當(dāng)信噪比在10 dB～22 dB時(shí)本文聯(lián)合初至拾取方法依然具有較高的拾取精度，STA/LTA和AR-AIC的拾取誤差呈增大趨勢；當(dāng)信噪比低于10 dB時(shí)，本文初至拾取方法的拾取誤差仍然保持在可接受的范圍內(nèi)，而其余兩種拾取方法的拾取誤差已出現(xiàn)較大的拾取誤差，將對(duì)后續(xù)層析計(jì)算造成較大的影響。因此，本文聯(lián)合初至拾取方法對(duì)不同信噪比的地震數(shù)據(jù)均具有較高的拾取精度。

圖1 合成地震信號(hào)去噪效果對(duì)比Fig.1 Comparison of the effect of synthetic seismic waves denoising

圖2 合成地震信號(hào)各道去噪后SNR分布對(duì)比Fig.2 Comparison of SNR distribution of synthetic seismic waves denoising

表1 合成地震數(shù)據(jù)各方法去噪效果對(duì)比Tab. 1 Comparison of denoising effects of synthetic seismic data

圖3 不同信噪比實(shí)際地震數(shù)據(jù)各方法去噪效果對(duì)比Fig.3 Comparison of different SNR of real seismic data denoising(a)SNR=15 dB；(b)SNR=6 dB；(c)SNR=2 dB

圖4 各信噪比數(shù)據(jù)對(duì)應(yīng)不同拾取方法拾取誤差分布圖Fig.4 The different SNR data corresponding to different first-break picking methods of pick up error profile

表2 不同信噪比實(shí)際地震數(shù)據(jù)各方法去噪后SNR及RMES變化對(duì)比

表3 對(duì)應(yīng)拾取誤差參數(shù)對(duì)比Tab.3 Comparison of first-break picking error parameter

圖5 各去噪方法與本文聯(lián)合初至拾取方法組合后拾取效果對(duì)比Fig.5 Comparison of the picking effect of each denoising method combined with the first-break picking method

將40道合成單同相軸地震數(shù)據(jù)經(jīng)不同方法去噪后的地震道集，采用改進(jìn)滑動(dòng)時(shí)窗能量比聯(lián)合AIC準(zhǔn)則法拾取初至波走時(shí)，驗(yàn)證本文拾取方法的有效性。圖5為各去噪方法與本文聯(lián)合初至拾取方法組合后對(duì)應(yīng)40道拾取效果對(duì)比。表3為對(duì)應(yīng)拾取均方根誤差RMES及其在相應(yīng)絕對(duì)拾取誤差范圍內(nèi)的占比情況。由圖5和表3數(shù)據(jù)可知，采用本文CEEME聯(lián)合小波閾值法去除地震隨機(jī)噪聲后，再采用本文改進(jìn)滑動(dòng)時(shí)窗能量比聯(lián)合AIC信息準(zhǔn)則對(duì)去噪后地震數(shù)據(jù)進(jìn)行初至拾取，其拾取誤差較小，各道初至?xí)r間分布與純凈信號(hào)經(jīng)人工拾取相當(dāng)，拾取誤差明顯低于未去噪直接采用本文聯(lián)合初至拾取方法和EMD去噪后采用本文聯(lián)合初至拾取方法的組合模式。本文聯(lián)合去噪方法與本文聯(lián)合初至拾取方法的組合，在50 ms內(nèi)的拾取誤差范圍內(nèi)，全部準(zhǔn)確拾??；在10 ms拾取誤差范圍內(nèi)，達(dá)95%準(zhǔn)確率；在5 ms拾取誤差范圍內(nèi)，達(dá)87%準(zhǔn)確率。明顯優(yōu)于其他組合模式。驗(yàn)證了本文噪聲壓制方法和本文初至拾取方法的有效性。

圖6 四組實(shí)際地震單炮記錄Fig.6 Four groups of actual seismic records

圖6為四組實(shí)際地震單炮道集記錄。圖7為對(duì)應(yīng)四組單炮記錄采用不同初至拾取方法的拾取結(jié)果分布圖。由圖7可知本文聯(lián)合去噪方法與本文聯(lián)合初至拾取方法的組合拾取結(jié)果與人工拾取相當(dāng)，相比其他去噪和拾取方法組合具有較低的拾取誤差；在圖7中，由于壞道或干擾道的存在使得其余初至拾取組合產(chǎn)生明顯的拾取錯(cuò)誤，本文聯(lián)合去噪方法與本文聯(lián)合初至拾取方法組合有效避免了上述錯(cuò)誤的發(fā)生，拾取精度最高。圖8為對(duì)應(yīng)拾取精度ROC曲線圖，其中ROC曲線右下方的面積被定義為AUC值，其值越接近“1”，表示該算法準(zhǔn)確性越高。在本文的初至拾取任務(wù)，在給定拾取誤差ξ的情況下可轉(zhuǎn)化為二分類問題(即正確拾取與錯(cuò)誤拾取)。因此拾取結(jié)果可以劃分成兩類：正類(Positive)和負(fù)類(Negative)，則會(huì)出現(xiàn)四種情況，即真正類(Tree Positive，TP)、假正類(False Positive，F(xiàn)P)、假負(fù)類(False Negative，F(xiàn)N)及真負(fù)類(True Negative，TN)。ROC曲線的橫坐標(biāo)負(fù)正類率(FPR)及縱坐標(biāo)真正類率(TPR)的計(jì)算公式如下：

(9)

(10)

拾取誤差ξ的大小將對(duì)ROC曲線的形狀有影響，但對(duì)AUC值的大小次序沒有影響。因此，ROC曲線可消除人為界定參數(shù)影響，對(duì)模型算法本身效果做出評(píng)價(jià)。

四組實(shí)際地震記錄中采用本文聯(lián)合去噪方法與本文聯(lián)合初至拾取方法組合的ROC曲線AUC平均值達(dá)0.96，明顯高于其余組合。綜合以上分析可知，本文聯(lián)合去噪方法與本文聯(lián)合初至拾取方法的組合模式，可取得最好的初至拾取效果。

圖7 四組實(shí)際地震記錄初至拾取效果對(duì)比Fig.7 Comparison of the first-break picking effect

4 結(jié)論與認(rèn)識(shí)

圖8 四組實(shí)際地震記錄初至拾取ROC曲線圖Fig.8 Comparison of the first-break picking ROC curve(a)ROC曲線；(b)ROC曲線；(c)ROC曲線；(d)ROC曲線

作者采用CEEMD聯(lián)合小波閾值法有效壓制地震隨機(jī)噪聲，CEEMD方法避免了EMD方法中模態(tài)混疊的問題，有效改善EEMD方法中加入隨機(jī)白噪聲對(duì)原始信號(hào)造成污染的缺點(diǎn)。地震信號(hào)經(jīng)CEEMD分解得到各級(jí)IMF分量，采用EIF有效指標(biāo)函數(shù)確定有效信號(hào)和噪聲的能量分界點(diǎn)，對(duì)部分高頻IMF分量再次進(jìn)行小波閾值去噪，并保持低頻IMF分量不變；最后將去噪后的高頻IMF分量和低頻IMF分量以及殘余量進(jìn)行重構(gòu)，即得到去噪后的地震信號(hào)。該方法有效改善傳統(tǒng)方法中直接舍棄部分IMF高頻分量，造成其中有效信息損失，進(jìn)而導(dǎo)致重構(gòu)后信號(hào)失真的問題。

本文采用改進(jìn)滑動(dòng)時(shí)窗能量比聯(lián)合AIC準(zhǔn)則拾取去噪后地震信號(hào)初至波走時(shí)，該方法首先根據(jù)改進(jìn)滑動(dòng)時(shí)窗能量比法獲取大致的初至?xí)r刻，然后在該時(shí)刻前后取一時(shí)窗，使用AIC方法計(jì)算AR模型在該時(shí)窗范圍內(nèi)的局部極小值，即得到精確的初至波走時(shí)。該聯(lián)合方法可有效緩解能量比方法中時(shí)窗大小的設(shè)定對(duì)拾取精度的影響；只在局部范圍內(nèi)計(jì)算AIC值，可有效避免在真實(shí)初至?xí)r刻來臨之前，出現(xiàn)多個(gè)AIC局部極小值，從而導(dǎo)致初至走時(shí)拾取出錯(cuò)的問題。

地震信號(hào)去噪及其初至走時(shí)拾取方法中，多方法的巧妙融合可有效避免單一方法的局限性，充分發(fā)揮兩種方法的優(yōu)點(diǎn)，有效改善地震噪聲壓制的效果和提高初至波走時(shí)的拾取精度，為后續(xù)地震層析的準(zhǔn)確性創(chuàng)造條件。