劉財(cái)， 王博,2， 劉洋*

1 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 130026 2 石家莊經(jīng)濟(jì)學(xué)院， 石家莊 050031

?

基于Duffing振子混沌系統(tǒng)的地震速度分析方法

劉財(cái)1， 王博1,2， 劉洋1*

1 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 130026 2 石家莊經(jīng)濟(jì)學(xué)院， 石家莊 050031

強(qiáng)隨機(jī)噪聲干擾是導(dǎo)致地震勘探資料低信噪比的主要原因，如何在強(qiáng)隨機(jī)噪聲干擾下獲取有效的信息是值得關(guān)注的問(wèn)題.Duffing振子混沌系統(tǒng)是一個(gè)非線性的動(dòng)力學(xué)系統(tǒng)，其對(duì)強(qiáng)隨機(jī)噪聲具有免疫能力，而對(duì)特定的周期性信號(hào)具有敏感性.本文提出一種基于Duffing振子混沌系統(tǒng)的速度分析方法.對(duì)CMP道集按照時(shí)距曲線關(guān)系進(jìn)行移動(dòng)窗口截取，將所截取的信號(hào)構(gòu)建為待測(cè)信號(hào)加入Duffing振子混沌系統(tǒng)，通過(guò)相圖網(wǎng)格分割方法(GPM)判斷系統(tǒng)狀態(tài)的改變，從而在強(qiáng)隨機(jī)噪聲背景下獲得高分辨率的速度譜.理論模型和實(shí)際資料的處理結(jié)果表明，與傳統(tǒng)的水平疊加速度分析方法相比，本方法能夠在強(qiáng)隨機(jī)噪聲背景下獲得更準(zhǔn)確的速度分析結(jié)果.

混沌系統(tǒng)； Duffing振子； 強(qiáng)隨機(jī)噪聲； 涌浪噪聲； 速度分析

1 引言

地震勘探資料中的隨機(jī)噪聲是影響地震勘探資料信噪比的主要原因之一.隨機(jī)噪聲的產(chǎn)生通常與接收、激發(fā)以及儀器本身等因素有關(guān).近年來(lái)，國(guó)內(nèi)外學(xué)者就隨機(jī)噪聲的壓制以及有效地震信號(hào)的提取與重建進(jìn)行了大量的研究并取得了一定的進(jìn)展(劉洋等，2009，2011；林紅波等，2011；白蘭淑等，2014).在海上拖纜地震勘探中，涌浪噪聲是一種常見的強(qiáng)隨機(jī)噪聲干擾.涌浪噪聲具有低頻、強(qiáng)能量的特點(diǎn)，在原始地震記錄上會(huì)形成強(qiáng)烈的低頻噪聲背景(郭建卿等，2007).國(guó)內(nèi)外一些研究人員對(duì)涌浪噪聲的產(chǎn)生機(jī)制與壓制方法進(jìn)行了相關(guān)研究.通常認(rèn)為檢波器掛上異物或電纜平衡不好而造成的干擾波是涌浪噪聲的主要產(chǎn)生原因.Parrish(2005)研究發(fā)現(xiàn)拖纜弦波能夠?yàn)橛坷嗽肼暤奶攸c(diǎn)提供合理的解釋.Elboth等(2009)從流體動(dòng)力學(xué)角度分析了涌浪噪聲的產(chǎn)生原因，認(rèn)為與天氣條件相關(guān)的水流擾動(dòng)與拖纜之間的相互作用是產(chǎn)生涌浪噪聲的原因，并提出了相應(yīng)的壓制方法.宋家文等(2001)采用地震道隨機(jī)重排聯(lián)合EM(Expectation Maximization，最大期望)算法對(duì)高振幅涌浪噪聲進(jìn)行衰減壓制.徐善輝(2012)使用HHT(Hilbert-Huang Transform，希爾伯特黃變換)技術(shù)進(jìn)行時(shí)頻域?yàn)V波及低切方法壓制涌浪噪聲.在強(qiáng)隨機(jī)噪聲的干擾下，由于地震同相軸淹沒(méi)在噪聲中，傳統(tǒng)的水平疊加速度分析方法會(huì)產(chǎn)生不準(zhǔn)確的速度拾取結(jié)果.近年來(lái)，有研究者提出了一些新的速度分析方法.Fomel(2009)提出了基于AB相似性的速度分析方法，與傳統(tǒng)速度分析方法相比，該方法受振幅變化影響更小.Blias(2009)針對(duì)上覆地層速度橫向變化的情況，提出一系列技術(shù)方法以提高疊加速度的求取精度.Luo和Hale(2012)采用加權(quán)相似的方法，提高了基于相似性方法速度分析的分辨率.而針對(duì)強(qiáng)隨機(jī)噪聲干擾的地震速度分析方法研究具有重要意義，如何獲取強(qiáng)隨機(jī)噪聲干擾下的速度信息是一個(gè)值得關(guān)注的問(wèn)題.

在信號(hào)檢測(cè)領(lǐng)域中，國(guó)內(nèi)外很多相關(guān)學(xué)者對(duì)含有噪聲情況下的有效信號(hào)提取與檢測(cè)進(jìn)行了研究.曹開田和楊震(2010)提出了基于隨機(jī)矩陣?yán)碚?Tulino and Verdú，2004)的頻譜感知方法，該方法對(duì)噪聲的不確定性具有較強(qiáng)的魯棒性.而在基于非線性方法的弱信號(hào)檢測(cè)方法中，對(duì)混沌系統(tǒng)的研究受到了廣泛的關(guān)注.Brix和Pipenberg(1992)提出將混沌理論應(yīng)用于弱信號(hào)檢測(cè)，通過(guò)實(shí)驗(yàn)從隨機(jī)高斯噪聲背景下檢測(cè)出了有效的信號(hào).Haykin和Li (1995)利用人工神經(jīng)網(wǎng)絡(luò)的方法實(shí)現(xiàn)了對(duì)混沌噪聲下的目標(biāo)弱信號(hào)檢測(cè).王冠宇等(1997)利用Duffing振子對(duì)強(qiáng)噪聲背景下的弱信號(hào)檢測(cè)進(jìn)行了相關(guān)研究.Wang等(1999)利用振子混沌系統(tǒng)實(shí)現(xiàn)了在白噪聲背景下低信噪比的正弦信號(hào)的檢測(cè).聶春燕和石要武(2001)提出了將互相關(guān)與Duffing振子相結(jié)合進(jìn)行信號(hào)檢測(cè)的方法.李月等(2005)針對(duì)地震勘探資料湮沒(méi)在隨機(jī)噪聲中的微弱同相軸問(wèn)題，提出了基于混沌理論的混沌振子檢測(cè)算法.李月等(2006)構(gòu)建了由兩個(gè)Duffing方程耦合確定的混沌系統(tǒng)，并通過(guò)仿真實(shí)驗(yàn)，表明此類混沌系統(tǒng)具有更好的抗噪聲能力和更穩(wěn)定的周期相態(tài).Li等(2009)將Duffing振子混沌系統(tǒng)用于地震弱信號(hào)檢測(cè)，在共炮點(diǎn)道集上實(shí)現(xiàn)了強(qiáng)噪聲背景下的地震同相軸識(shí)別.

本文提出了一種基于Duffing振子混沌系統(tǒng)的地震速度分析方法.以不同速度對(duì)CMP道集數(shù)據(jù)進(jìn)行移動(dòng)窗口截取，并構(gòu)建成待測(cè)信號(hào)，將構(gòu)建的待測(cè)信號(hào)加入到臨界狀態(tài)的Duffing振子混沌系統(tǒng)中，利用系統(tǒng)狀態(tài)的改變來(lái)判斷待測(cè)信號(hào)是否具有特定的周期性，從而獲取準(zhǔn)確的動(dòng)校正速度.同時(shí)提出一種基于混沌系統(tǒng)相圖的網(wǎng)格分割法(GPM)對(duì)混沌系統(tǒng)的狀態(tài)進(jìn)行判斷，定量地反映混沌系統(tǒng)的狀態(tài)，為自動(dòng)拾取地震速度信息提供有效的方法.理論模型與實(shí)際數(shù)據(jù)驗(yàn)證的結(jié)果表明，與傳統(tǒng)的水平疊加速度分析方法相比，本文所提出的方法能夠在強(qiáng)隨機(jī)噪聲背景下獲得更準(zhǔn)確的速度分析結(jié)果.

2 理論基礎(chǔ)

2.1 Duffing方程回顧

Duffing方程是由Duffing(1918)提出的一個(gè)用以描述軟彈簧弱阻尼振動(dòng)的非線性方程.經(jīng)典Holmes型Duffing方程的形式是

(1)

式中，γcos(t)為混沌系統(tǒng)的周期性策動(dòng)力(參考信號(hào))，k為阻尼系數(shù)，一般令k=0.5，ax+bx3為非線性恢復(fù)力.取a=-1，b=1時(shí)，公式(1)的狀態(tài)方程形式為

(2)

用4階龍格庫(kù)塔(Runge-Kutta)法對(duì)此狀態(tài)方程進(jìn)行求解，龍格庫(kù)塔(Runge-Kutta)法的步長(zhǎng)取信號(hào)的采樣步長(zhǎng)，求得每一時(shí)間采樣點(diǎn)處的x(n+1)和y(n+1)，輸出x(n+1)和y(n+1)得到系統(tǒng)的相圖.首先，固定k=0.5，研究γ對(duì)系統(tǒng)狀態(tài)的影響.當(dāng)γ從0逐漸增大時(shí)，系統(tǒng)狀態(tài)隨γ的變化呈現(xiàn)出規(guī)律性的變化：經(jīng)歷平衡點(diǎn)、同宿軌道、分叉狀態(tài)、混沌狀態(tài)和大尺度周期狀態(tài).圖1為Duffing系統(tǒng)在混沌狀態(tài)時(shí)和大尺度周期狀態(tài)時(shí)的相圖.從圖1中可以看出，系統(tǒng)從混沌狀態(tài)到大尺度周期狀態(tài)，相圖產(chǎn)生了非常明顯的變化.

2.2 相圖網(wǎng)格分割法(GPM)

目前，對(duì)于Duffing振子混沌系統(tǒng)狀態(tài)的判斷依據(jù)主要有以下幾種方法：對(duì)相圖的直接觀察、Lyapunov指數(shù)(Wolf et al.，1985)、Melnikov法(Wiggins，1988)、Poincare截面(Dubois et al.，1982)等.其中，對(duì)相圖的直接觀察最為簡(jiǎn)單直接，但是效率比較低，無(wú)法使用到實(shí)際的地震速度分析中，且難以實(shí)現(xiàn)自動(dòng)化.而諸如Lyapunov指數(shù)等方法又存在計(jì)算效率的問(wèn)題，過(guò)高的計(jì)算成本使其不適合應(yīng)用在地震速度分析上.本文提出了對(duì)系統(tǒng)相圖進(jìn)行網(wǎng)格分割的判斷方法，將二維的相圖轉(zhuǎn)化成為一維的參數(shù)分析，實(shí)現(xiàn)對(duì)系統(tǒng)狀態(tài)的定量判斷.在提高計(jì)算效率的同時(shí)，為速度的掃描與自動(dòng)拾取提供了條件.

圖1 Duffing系統(tǒng)在不同狀態(tài)下的相圖(a)γ=0.824，混沌狀態(tài)；(b)γ=0.828，大尺度周期狀態(tài).Fig.1 The phase plane diagram of different system states(a) γ=0.824, chaotic state; (b) γ=0.828, large-scale periodic state.

相圖網(wǎng)格分割方法的原理如圖2所示，將Duffing系統(tǒng)的相圖按一定尺寸的網(wǎng)格(如正方形網(wǎng)格)進(jìn)行平均分割，例如每個(gè)網(wǎng)格的邊長(zhǎng)為gx.對(duì)每一網(wǎng)格進(jìn)行判斷賦值，若網(wǎng)格中有相軌跡經(jīng)過(guò)，對(duì)其賦權(quán)值為1；而對(duì)沒(méi)有相軌跡經(jīng)過(guò)的網(wǎng)格，則賦權(quán)值為0.對(duì)賦值后的全部網(wǎng)格權(quán)值進(jìn)行求和統(tǒng)計(jì)，得出系統(tǒng)狀態(tài)判斷參數(shù)p.從圖2可以看出，在網(wǎng)格大小gx固定時(shí)，大尺度周期狀態(tài)的系統(tǒng)狀態(tài)判斷參數(shù)p要明顯小于混沌狀態(tài)時(shí)的系統(tǒng)狀態(tài)判斷參數(shù).

利用Duffing振子混沌系統(tǒng)檢測(cè)微弱特征信號(hào)，主要根據(jù)混沌系統(tǒng)對(duì)特定小信號(hào)的敏感性以及對(duì)噪聲“免疫”的特點(diǎn)進(jìn)行檢測(cè).將Duffing方程中的參數(shù)γ設(shè)置在臨界值γc附近，使系統(tǒng)處于變化的邊緣.根據(jù)混沌學(xué)理論，這時(shí)噪聲對(duì)系統(tǒng)的影響很小，不會(huì)引起系統(tǒng)相態(tài)的變化，而微弱特征信號(hào)對(duì)系統(tǒng)狀態(tài)的改變起著決定性作用.當(dāng)待測(cè)信號(hào)加入系統(tǒng)后，觀察系統(tǒng)是否從混沌狀態(tài)變化到穩(wěn)定的大尺度周期狀態(tài)，以此判斷待測(cè)信號(hào)中是否含有特定的周期信號(hào).因此，需要預(yù)先確定系統(tǒng)臨界狀態(tài)的臨界值γc的大小.通過(guò)取網(wǎng)格分割法的網(wǎng)格大小gx，使混沌狀態(tài)時(shí)的系統(tǒng)判斷參數(shù)pc與大尺度周期狀態(tài)時(shí)的系統(tǒng)判斷參數(shù)pp差別最大.經(jīng)過(guò)試驗(yàn)對(duì)比，選擇gx=0.02作為正方形網(wǎng)格參數(shù).當(dāng)系統(tǒng)無(wú)外加信號(hào)輸入時(shí)，γ從0.750逐漸增加至0.900的過(guò)程中，隨著系統(tǒng)由混沌狀態(tài)變?yōu)榇蟪叨戎芷跔顟B(tài)，系統(tǒng)判斷參數(shù)p有出現(xiàn)階躍躍變，如圖3a所示.進(jìn)一步通過(guò)對(duì)系統(tǒng)狀態(tài)判斷參數(shù)(圖3a)求導(dǎo)，得到系統(tǒng)的臨界狀態(tài)閾值γc為0.826(圖3b).

2.3 用于地震速度分析的Duffing振子混沌系統(tǒng)

根據(jù)公式(1)，標(biāo)準(zhǔn)的Duffing方程形式是

(3)

為了實(shí)現(xiàn)對(duì)不同頻率信號(hào)進(jìn)行檢測(cè)，對(duì)公式(3)進(jìn)行坐標(biāo)變換.令t=ωτ，則x(t)=x(ωτ)=xτ(τ)，可得

(4)

(5)

代入公式(3)得：

(6)

寫成狀態(tài)方程形式為：

(7)

與公式(3)進(jìn)行比較后發(fā)現(xiàn)，x和y都變成了原來(lái)的ω倍，系統(tǒng)以原來(lái)ω倍的速率進(jìn)行運(yùn)動(dòng)，但系統(tǒng)的動(dòng)力學(xué)性質(zhì)并未發(fā)生變化.這樣，通過(guò)選取不同的ω

圖2 相圖網(wǎng)格分割法原理示意圖(a)混沌狀態(tài)的相圖；(b)大尺度周期狀態(tài)的相圖.Fig.2 The illustration of grid partition method(a) The phase plane diagram of chaotic state; (b) The phase plane diagram of large-scale periodic state.

值可以實(shí)現(xiàn)對(duì)不同頻率信號(hào)的檢測(cè).為了表述方便，令τ=t，并消去角標(biāo)，將外加信號(hào)加入系統(tǒng)后，公式(6)變?yōu)?/p>

=γcos(ωt)+ξR(t)，

(8)

其中，R(t)為系統(tǒng)外加信號(hào)，這里R(t)=s(t)+n(t)，s(t)為待測(cè)信號(hào)，n(t)為噪聲，ξ為外置信號(hào)的振幅可調(diào)參數(shù).由于待檢測(cè)信號(hào)可以看作是系統(tǒng)內(nèi)置周期策動(dòng)力的補(bǔ)充，因此，這種“補(bǔ)充”是否適當(dāng)，還應(yīng)當(dāng)考慮周期性策動(dòng)力與待檢測(cè)信號(hào)之間的相位關(guān)系.若待測(cè)信號(hào)為R(t)=cos(100t+0.5π)，令ξ=0.01，將公式(8)對(duì)應(yīng)的Duffing振子混沌系統(tǒng)調(diào)至臨界混沌狀態(tài).此時(shí)，系統(tǒng)的周期性策動(dòng)力與待測(cè)信號(hào)的相位差為0.5π.將待測(cè)信號(hào)加入Duffing振子混沌系統(tǒng)中，系統(tǒng)仍處于混沌狀態(tài)，系統(tǒng)狀態(tài)沒(méi)有改變，如圖4所示.

為解決該問(wèn)題，加入相位參數(shù)φ，Duffing振子混沌系統(tǒng)變?yōu)椋?/p>

=γcos(ωt+φ)+ξR(t).

(9)

令φ=0.5π，此時(shí)，系統(tǒng)的周期性策動(dòng)力與待測(cè)信號(hào)間的相位差為0.從時(shí)域圖上看，在消除了相位差之后，待測(cè)信號(hào)對(duì)處于閾值狀態(tài)的周期性策動(dòng)力形成了很好的“補(bǔ)充”，使系統(tǒng)進(jìn)入大尺度周期狀態(tài)，如圖5所示.因此，加入相位參數(shù)后，能夠提高Duffing振子混沌系統(tǒng)的檢測(cè)精度.

2.4 基于Duffing振子混沌系統(tǒng)的速度分析方法

首先，建立一個(gè)理論無(wú)噪聲共中心點(diǎn)(CMP)合成地震記錄，如圖6a所示.時(shí)間采樣4 ms，記錄時(shí)間4 s，80個(gè)記錄道，最小炮檢距為零，道間距50 m.

圖3 混沌系統(tǒng)臨界值曲線(a)狀態(tài)判斷參數(shù)p隨γ變化的關(guān)系曲線；(b)狀態(tài)判斷參數(shù)導(dǎo)數(shù)曲線(γc=0.826).Fig.3 Determine the critical value of chaotic system(a) The judgment parameter p varies as the γ increases; (b) The derivative of judgment parameter (γc=0.826).

圖4 系統(tǒng)周期性策動(dòng)力與待測(cè)信號(hào)存在相位差時(shí)混沌系統(tǒng)狀態(tài)圖(a)時(shí)域圖(實(shí)線為系統(tǒng)周期性策動(dòng)力，虛線為外加待測(cè)信號(hào))；(b)系統(tǒng)相圖.Fig.4 The phase plane diagram of chaotic system with the phase difference between the periodic force and input signal(a) Time-domain plot (the solid line is periodic force of system; the dash line is input signal); (b) the phase plane diagram.

圖5 消除系統(tǒng)周期性策動(dòng)力與待測(cè)信號(hào)間的相位差后，系統(tǒng)進(jìn)入大尺度周期狀態(tài)(a)時(shí)域圖(實(shí)線為系統(tǒng)周期性策動(dòng)力，虛線為外加待測(cè)信號(hào))；(b)系統(tǒng)相圖.Fig.5 The phase plane diagram of chaotic system without the phase difference between the periodic force and input signal(a) Time-domain plot (the solid line is periodic force of system; the dash line is input signal); (b) the phase plane diagram.

圖6 理論共中心點(diǎn)合成地震記錄(a)無(wú)噪聲合成地震記錄；(b)加入強(qiáng)隨機(jī)噪聲后的合成地震記錄.Fig.6 Synthetic CMP record(a) Noise-free synthetic CMP record; (b) Synthetic CMP record with strong random noise.

其中，在t0=1.2 s和2.6 s處，兩條同相軸所對(duì)應(yīng)的均方根速度Vr ms分別為1800 m·s-1和2500 m·s-1.合成地震記錄采用的子波是雷克子波，主頻fr=25 Hz.在如圖6a的人工合成共中心點(diǎn)(CMP)地震記錄中，加入高斯白噪聲，信噪比約為-16.14 dB(圖6b).從圖6b中可以看出，原本在t0=1.2 s和t0=2.6 s處的兩條地震同相軸已經(jīng)完全湮沒(méi)在噪聲之中，肉眼很難分辨.

對(duì)共中心點(diǎn)地震記錄進(jìn)行移動(dòng)窗口截取.對(duì)于每一個(gè)自激自收時(shí)間t0，用不同速度在各道上截取窗口大小為ws的信號(hào).在同一掃描速度下，將每道所截取到的信號(hào)首尾相接組成待檢測(cè)信號(hào).截取窗口中心點(diǎn)滿足雙曲線方程

(10)

其中，t0取值范圍為0～4000ms，x為炮檢距，V為掃描速度.選擇截?cái)啻翱诖笮s時(shí)，應(yīng)當(dāng)保證在窗口范圍內(nèi)能完整將子波的波形包含其中.同時(shí)，截取窗口大小ws與系統(tǒng)周期性策動(dòng)角頻率ω之間的關(guān)系應(yīng)滿足：

(11)

式中，n為正整數(shù)，即截?cái)啻翱诖笮s應(yīng)為系統(tǒng)周期性策動(dòng)力周期的整數(shù)倍.對(duì)模型所選用的主頻為25 Hz的雷克子波，本文選用的截?cái)啻翱诖笮s為100 ms.

按公式(10)的關(guān)系在CMP道集上進(jìn)行雙曲型移動(dòng)窗口截取.在t0=1.2 s(或2.6 s)處，若分別按1200 m·s-1(或2600 m·s-1)的速度正確進(jìn)行移動(dòng)窗口截取，所截取到的待測(cè)信號(hào)如圖7a所示，此時(shí)所截取到的待測(cè)信號(hào)為持續(xù)時(shí)間為8 s的雷克子波序列.而在t0=1.2 s(或2.6 s)處，當(dāng)掃描截取速度不正確時(shí)，所截取的待測(cè)信號(hào)則無(wú)法包含完整的80個(gè)雷克子波(圖7b).在正確的截取速度下，對(duì)如圖6b所示的混合強(qiáng)隨機(jī)噪聲后的CMP道集進(jìn)行移動(dòng)窗口截取，從時(shí)域圖上看，待檢測(cè)的雷克子波序列已經(jīng)完全湮沒(méi)在噪聲中(圖7c).此時(shí)，疊加噪聲后的信號(hào)信噪比為

(12)

其中，PS為信號(hào)的平均功率，PN為噪聲的平均功率.

在t0=1.2 s處，將掃描速度正確時(shí)所截取的待測(cè)信號(hào)(圖7a)作為待測(cè)信號(hào)R(t)，加入公式(9)所描述的Duffing振子混沌系統(tǒng)中.為使檢測(cè)更加準(zhǔn)確，調(diào)整γ稍小于γc使系統(tǒng)處于臨界混沌狀態(tài)邊緣，即令γ=0.824，周期策動(dòng)力角頻率ω=20×2π，周期策動(dòng)力相位φ=π/2，外加信號(hào)振幅可調(diào)參數(shù)ξ=0.02.加入待測(cè)信號(hào)R(t)后，原本處于臨界混沌狀態(tài)的系統(tǒng)變?yōu)榇蟪叨戎芷跔顟B(tài)(圖8a).掃描速度不正確時(shí)所截取的待測(cè)信號(hào)(圖7b)加入到上述臨界狀態(tài)的Duffing振子混沌系統(tǒng)中，由于待測(cè)信號(hào)不具有周期性，系統(tǒng)狀態(tài)依然處于混沌狀態(tài)(圖8b).將混合了噪聲后的信號(hào)加入臨界Duffing振子混沌系統(tǒng)中進(jìn)行檢測(cè).從相圖(圖8c)上看，與無(wú)噪聲時(shí)的情況(圖8a)相比，雖然噪聲的加入使相圖的輪廓變的有些“粗糙”，但系統(tǒng)的狀態(tài)依然是大尺度周期狀態(tài)，噪聲并未對(duì)系統(tǒng)狀態(tài)的改變產(chǎn)生影響.

圖7 檢測(cè)信號(hào)對(duì)比圖(a)掃描速度正確時(shí)的待測(cè)信號(hào)；(b)掃描速度不正確時(shí)的待測(cè)信號(hào)；(c)掃描速度正確時(shí)的含噪聲信號(hào).Fig.7 The comparison diagram of reconstructed signal(a) The reconstructed signal intercepted by correct velocity; (b) The reconstructed signal intercepted by incorrect velocity; (c) The reconstructed signal with noise intercepted by correct velocity.

圖8 加入待測(cè)信號(hào)后的系統(tǒng)相圖(a)圖7a對(duì)應(yīng)的待測(cè)信號(hào)；(b)圖7b對(duì)應(yīng)的信號(hào)；(c)圖7c對(duì)應(yīng)的待測(cè)信號(hào).Fig.8 The phase plane diagrams when added the different signals(a) The signal corresponding to Fig.7a; (b) The signal corresponding to Fig.7b; (c) The signal corresponding to Fig.7c.

3 理論模型分析

對(duì)疊加強(qiáng)隨機(jī)噪聲的地震記錄(圖6b)分別進(jìn)行水平疊加速度分析和基于Duffing振子混沌系統(tǒng)的速度分析.兩種分析方法的速度掃描步長(zhǎng)均選取為25 m·s-1，速度掃描范圍為1000～4000 m·s-1.水平疊加速度分析結(jié)果如圖9a所示.Duffing振子混沌系統(tǒng)的參數(shù)在每一個(gè)t0時(shí)間均為γ=0.824，周期策動(dòng)力角頻率ω=20×2π，周期策動(dòng)力相位φ=π/2，外加信號(hào)振幅可調(diào)參數(shù)ξ=0.02.對(duì)CMP地震記錄進(jìn)行窗口截?cái)鄴呙?，在各自激自收時(shí)間t0下，將每一截?cái)鄴呙杷俣人厝〉男盘?hào)作為待測(cè)信號(hào)R(t)加入系統(tǒng)中，采用相圖網(wǎng)格分割方法對(duì)系統(tǒng)狀態(tài)進(jìn)行判斷，網(wǎng)格大小gx=0.2，將系統(tǒng)判斷參數(shù)p輸出，得到速度譜(圖9b).從水平疊加速度分析的結(jié)果(圖9a)可以看到，由于強(qiáng)隨機(jī)噪聲的存在，使得對(duì)應(yīng)的速度譜幾乎沒(méi)有可分辨的能量團(tuán)存在，因此無(wú)法識(shí)別出同相軸的自激自收時(shí)間t0和相應(yīng)的RMS速度.由于噪聲的存在，在速度譜的低速區(qū)(1000～1500 m·s-1)產(chǎn)生強(qiáng)烈干擾，會(huì)導(dǎo)致自動(dòng)速度拾取錯(cuò)誤.圖9b為基于Duffing振子混沌系統(tǒng)的速度分析方法所獲得的速度譜.在t0=1.2 s，V=1800 m·s-1和t0=2.6 s，V=2500 m·s-1處，由于對(duì)應(yīng)合成地震記錄(圖6a)中的兩條同相軸，此時(shí)加入Duffing振子混沌系統(tǒng)的待測(cè)信號(hào)為混合了隨機(jī)噪聲的周期雷克子波序列，系統(tǒng)由混沌狀態(tài)轉(zhuǎn)變?yōu)榇蟪叨戎芷跔顟B(tài).因此系統(tǒng)狀態(tài)判斷參數(shù)p在這兩處為局部極小值，分別為162和145.與水平疊加速度分析相比，基于Duffing振子混沌系統(tǒng)的速度分析方法能夠獲得高分辨率的速度譜，且表現(xiàn)出很強(qiáng)的抗噪能力.

圖9 速度分析方法對(duì)比(a)水平疊加速度分析；(b)基于Duffing振子混沌系統(tǒng)的速度分析方法.Fig.9 Comparison on two velocity analysis methods(a) Standard velocity analysis spectrum; (b) Duffing chaotic system velocity analysis spectrum.

圖10 實(shí)際數(shù)據(jù)測(cè)試(a)某海上實(shí)際CMP數(shù)據(jù)體；(b)基于Duffing振子混沌系統(tǒng)速度分析方法的自動(dòng)速度拾取結(jié)果.Fig.10 Field data test(a) A marine dataset with strong swell noise; (b) Automatic picked velocities using Duffing chaotic system.

圖11 實(shí)際數(shù)據(jù)速度分析對(duì)比圖(a)傳統(tǒng)水平疊加速度分析；(b)基于Duffing振子混沌系統(tǒng)速度分析方法.Fig.11 Comparison on the velocity analysis results from field data(a) Standard velocity analysis result; (b) Velocity analysis result using Duffing chaotic system.

圖12 水平疊加結(jié)果對(duì)比圖(a)傳統(tǒng)水平疊加速度分析；(b)基于Duffing振子混沌系統(tǒng)速度分析方法.Fig.12 Comparison on stacking results(a) Standard velocity analysis; (b) Velocity analysis using Duffing chaotic system.

4 實(shí)際數(shù)據(jù)測(cè)試

在實(shí)際資料處理中，選取海上某地區(qū)CMP數(shù)據(jù)體，數(shù)據(jù)中包含較強(qiáng)的涌浪噪聲，如圖10a所示.該記錄共200個(gè)CMP道集.圖中三個(gè)截面分別為：左上截面為1.6 s位置的時(shí)間切片，左下截面為炮檢距為0.8 km時(shí)的共炮檢距剖面，右下為1.2933 km處的共中心點(diǎn)道集.涌浪噪聲導(dǎo)致地震資料的信噪比很低.分別采用傳統(tǒng)水平疊加方法和基于Duffing振子混沌系統(tǒng)的方法進(jìn)行速度分析，圖10b為采用Duffing振子混沌系統(tǒng)所獲得的RMS速度拾取結(jié)果.這里，Duffing振子混沌系統(tǒng)方法的參數(shù)為：截?cái)啻翱诖笮s取100 ms，γ=0.823，周期策動(dòng)力角頻率ω=20×2π，周期策動(dòng)力相位φ=0，外加信號(hào)振幅可調(diào)參數(shù)ξ=0.02，相圖網(wǎng)格分割方法(GPM)的網(wǎng)格大小gx=0.2.為進(jìn)行比較，選取1.2933 km處的CMP道集，采用兩種方法分別進(jìn)行速度分析，傳統(tǒng)水平疊加速度分析所獲得的速度譜如圖11a所示.由于涌浪噪聲的干擾，傳統(tǒng)速度分析方法所獲得的速度譜質(zhì)量較差，導(dǎo)致自動(dòng)速度拾取結(jié)果不正確.圖11b為采用基于Duffing振子混沌系統(tǒng)進(jìn)行速度分析所獲得的速度譜，與傳統(tǒng)方法相比，基于Duffing振子混沌系統(tǒng)的速度分析方法受涌浪噪聲的影響更小，自動(dòng)速度拾取結(jié)果更加準(zhǔn)確.需要說(shuō)明的是，由于在實(shí)際數(shù)據(jù)中地震子波隨偏移距的變化，由截?cái)啻翱趻呙杷鶚?gòu)成的待測(cè)信號(hào)的周期性將受到諸如遠(yuǎn)偏移距子波能量衰減、AVO效應(yīng)等因素的影響.待測(cè)信號(hào)的非周期性增加會(huì)造成基于Duffing振子混沌系統(tǒng)速度分析方法的分辨率降低.對(duì)采用傳統(tǒng)水平疊加方法和基于Duffing振子混沌系統(tǒng)的方法所獲得的速度結(jié)果進(jìn)行水平疊加，結(jié)果如圖12所示.其中，圖12b為采用基于Duffing振子混沌系統(tǒng)速度分析所得到的疊加剖面，與水平疊加速度分析結(jié)果所得到的疊加剖面(圖12a)相比，本文提出的方法能夠獲取更高信噪比的疊加結(jié)果.

5 結(jié)論

本文提出一種基于Duffing振子混沌系統(tǒng)速度分析方法.由于Duffing振子混沌系統(tǒng)對(duì)于強(qiáng)隨機(jī)噪聲具有“免疫”特性，使得該方法能夠在強(qiáng)噪聲背景下取得比常規(guī)速度分析方法更為準(zhǔn)確的速度分析結(jié)果.主要內(nèi)容包括：對(duì)CMP道集進(jìn)行移動(dòng)窗口截取，將截?cái)嗨玫牡卣鹱硬?gòu)建成Duffing振子混沌系統(tǒng)所能夠檢測(cè)的周期性信號(hào)，然后對(duì)加入信號(hào)后的系統(tǒng)狀態(tài)進(jìn)行判斷.本文提出相圖網(wǎng)格分割法對(duì)混沌系統(tǒng)的狀態(tài)進(jìn)行判斷，能夠獲得高效、穩(wěn)定且定量的系統(tǒng)狀態(tài)判斷結(jié)果，進(jìn)而可將系統(tǒng)狀態(tài)判斷結(jié)果用于自動(dòng)獲取速度信息.通過(guò)對(duì)理論模型和實(shí)際地震記錄進(jìn)行處理，驗(yàn)證了該速度分析方法能夠在強(qiáng)隨機(jī)噪聲背景下取得比常規(guī)速度分析方法更準(zhǔn)確的速度信息，進(jìn)而得到高信噪比的水平疊加結(jié)果.

Bai L S, Liu Y K, Lu H Y, et al. 2014. Curvelet-domain joint iterative seismic data reconstruction based on compressed sensing.ChineseJ.Geophys. (in Chinese), 57(9): 2937-2945, doi: 10.6038/cjg20140919. Birx D L, Pipenberg S J. 1992. Chaotic oscillators and complex mapping feed forward networks (CMFFNS) for signal detection in noisy environments. ∥ International Joint Conference on Neural Networks, 1992. IJCNN. Baltimore, MD: IEEE, 2: 881-888. Blias E. 2009. Stacking velocities in the presence of overburden velocity anomalies.GeophysicalProspecting, 57(3): 323-341.

Cao K T, Yang Z. 2010. DET cooperative spectrum sensing algorithm based on random matrix theory.JournalofElectronics&InformationTechnology(in Chinese), 32(1): 129-134.

Dubois M, Bergé P, Croquette V. 1982. Study of non steady convective regimes using Poincaré sections.JournaldePhysiqueLettres, 43(9): 295-298. Duffing G. 1918. Erzwungene Schwingungen bei ver?nderlicher Eigenfrequenz und ihre technische Bedeutung. Braunschweig, F. Vieweg & Sohn.

Elboth T, Reif B A, Andreassen ?. 2009. Flow and swell noise in marine seismic data.Geophysics, 74(2): Q17-Q25.

Fomel S. 2009. Velocity analysis using AB semblance.GeophysicalProspecting, 57(3): 311-321.

Guo J Q, Zhang B Q, Zhao Y S, et al. 2007. Case analysis of offshore seismic data processing in Mexico Bay.OilGeophysicalProspecting(in Chinese), 42(Suppl.): 132-135.Haykin S, Li X B. 1995. Detection of signals in chaos.ProceedingsoftheIEEE, 83(1): 95-122.

Li Y, Yang B J, Zhao X P, et al. 2005. An algorithm of chaotic vibrator to detect weak events in seismic prospecting records.ChineseJ.Geophys. (in Chinese), 48(6): 1428-1433, doi: 10.3321/j.issn:0001-5733.2005.06.028.

Li Y, Lu P, Yang B J, et al. 2006. Applying a special kind of two coupled Duffing oscillator system to detect periodic signals under the background of strong colored noise.ActaPhysicaSinica(in Chinese), 55(4): 1672-1677.

Li Y, Yang B J, Badal J, et al. 2009. Chaotic system detection of weak seismic signals.GeophysicalJournalInternational, 178(3): 1493-1522.

Lin H B, Li Y, Xu X C. 2011. Segmenting time-frequency peak filtering method to attenuation of seismic random noise.ChineseJ.Geophys. (in Chinese), 54(5): 1358-1366, doi: 10.3969/j.issn.0001-5733.2011.05.025.

Liu Y, Fomel S, Liu C, et al. 2009. High-order seislet transform and its application of random noise attenuation.ChineseJ.Geophys. (in Chinese), 52(8): 2142-2151, doi: 10.3969/j.issn.0001-5733.2009.08.024.

Liu Y, Wang D, Liu C, et al. 2011. Weighted median filter based on local correlation and its application to poststack random noise attenuation.ChineseJ.Geophys. (in Chinese), 54(2): 358-367, doi: 10.3969/j.issn.0001-5733.2011.02.012.

Luo S, Hale D. 2012. Velocity analysis using weighted semblance.Geophysics, 77(2): U15-U22.

Nie C Y, Shi Y W. 2001. The research of weak signal detection based on cross-correlation and chaos theory.ChineseJournalofScientificInstrument(in Chinese), 22(1): 32-35.

Parrish J F. 2005. Streamer string waves and swell noise. ∥ 2005 SEG Annual Meeting. Society of Exploration Geophysicists. 72-75.

Song J W, Chen X H, Liu G C. 2011. High-amplitude swell noise attenuation using seismic trace recording and EM algorithm (in Chinese). ∥ 27th Annual General Assembly of the Chinese Geophysical Society. 667.

Tulino A M, Verdú S. 2004. Random Matrix Theory and Wireless Communications. Hanover, USA: Now Publisher Inc.

Wang G Y, Chen D J, Lin J Y, et al. 1999. The application of chaotic oscillators to weak signal detection.IEEETransactionsonIndustrialElectronics, 46(2): 440-444. Wang G Y, Tao G L, Chen X, et al. 1997. Application of chaotic oscillators in signal detection in the strong noise background.ChineseJournalofScientificInstrument(in Chinese), 18(2): 209-212.

Wiggins S. 1988. Global Bifurcations and Chaos. Berlin: Springer. Wolf A, Swift J B, Swinney H L, et al. 1985. Determining Lyapunov exponents from a time series.PhysicaD:NonlinearPhenomena, 16(3): 285-317. Xu S H. 2012. The study of seismic noise attenuation and thin-layer prediction technique based on Hilbert-Huang transform [Ph. D. thesis] (in Chinese). Changchun: Jilin University.

附中文參考文獻(xiàn)

白蘭淑, 劉伊克, 盧回憶等. 2014. 基于壓縮感知的Curvelet域聯(lián)合迭代地震數(shù)據(jù)重建. 地球物理學(xué)報(bào), 57(9): 2937-2945, doi: 10.6038/cjg20140919.

曹開田, 楊震. 2010. 基于隨機(jī)矩陣?yán)碚摰腄ET合作頻譜感知算法. 電子與信息學(xué)報(bào), 32(1): 129-134.

郭建卿, 張保慶, 趙貽水等. 2007. 墨西哥灣海洋地震數(shù)據(jù)處理案例分析. 石油地球物理勘探, 42(增刊): 132-135.

李月, 楊寶俊, 趙雪平等. 2005. 檢測(cè)地震勘探微弱同相軸的混沌振子算法. 地球物理學(xué)報(bào), 48(6): 1428-1433, doi: 10.3321/j.issn:0001-5733.2005.06.028.

李月, 路鵬, 楊寶俊等. 2006. 用一類特定的雙耦合Duffing振子系統(tǒng)檢測(cè)強(qiáng)色噪聲背景中的周期信號(hào). 物理學(xué)報(bào), 55(4): 1672-1677.

林紅波, 李月, 徐學(xué)純. 2011. 壓制地震勘探隨機(jī)噪聲的分段時(shí)頻峰值濾波方法. 地球物理學(xué)報(bào), 54(5): 1358-1366, doi: 10.3969/j.issn.0001-5733.2011.05.025.

劉洋, Fomel S, 劉財(cái)?shù)? 2009. 高階seislet變換及其在隨機(jī)噪聲消除中的應(yīng)用. 地球物理學(xué)報(bào), 52(8): 2142-2151, doi: 10.3969/j.issn.0001-5733.2009.08.024.

劉洋, 王典, 劉財(cái)?shù)? 2011. 局部相關(guān)加權(quán)中值濾波技術(shù)及其在疊后隨機(jī)噪聲衰減中的應(yīng)用. 地球物理學(xué)報(bào), 54(2): 358-367, doi: 10.3969/j.issn.0001-5733.2011.02.012.

聶春燕, 石要武. 2001. 基于互相關(guān)檢測(cè)和混沌理論的弱信號(hào)檢測(cè)方法研究. 儀器儀表學(xué)報(bào), 22(1): 32-35.

宋家文, 陳小宏, 劉國(guó)昌. 2001. 地震道隨機(jī)重排聯(lián)合EM算法衰減高振幅涌浪噪聲. ∥中國(guó)地球物理學(xué)會(huì)第二十七屆年會(huì)論文集. 667.

王冠宇, 陶國(guó)良, 陳行等. 1997. 混沌振子在強(qiáng)噪聲背景信號(hào)檢測(cè)中的應(yīng)用. 儀器儀表學(xué)報(bào), 18(2): 209-212.

徐善輝. 2012. 基于 Hilbert-Huang變換的地震噪聲衰減與薄層預(yù)測(cè)技術(shù)研究[博士論文]. 長(zhǎng)春: 吉林大學(xué).

(本文編輯 何燕)

Seismic velocity analysis based on the Duffing oscillator chaotic system

LIU Cai1, WANG Bo1,2, LIU Yang1*

1CollegeofGeo-explorationScienceandTechnology，JilinUniversity，Changchun130026，China2ShijiazhuangUniversityofEconomics，Shijiazhuang050031，China

The strong random noise is the main reason for low signal-to-noise ratio of seismic data. It is a much concerned issue how to extract useful information from data under a strong random noise background. The Duffing oscillator chaotic system is a nonlinear dynamic one that is immune to noise, but very sensitive to particular periodic signals. This paper presents a seismic velocity analysis method based on the Duffing oscillator chaotic system.In a chaotic system, one needs to make the system reach a critical state before it is used to analyze signal. We propose a new method called the grid partition method (GPM) to estimate the chaotic system state. The principle of GPM is dividing the phase plane diagram into many small square grids. If the phase trajectories pass any small grid, this grid will be assigned to 1; otherwise, the value of the small grid is 0. Summing all the values of small grids will obtain the grid judgment parameter, which provides a stable criterion for seismic velocity analysis. A moving window intercept method following the time-distance curve is used to reconstruct the signal from the CMP gathers. Then we add the reconstructed signal to the Duffing oscillator chaotic system and use the grid partition method to judge the state of system. Therefore, the high resolution velocity spectrum can be obtained by the proposed method from the data with a strong random noise background.We apply the proposed method to the synthetic model which is a CMP gather contains two events. Adding strong random noise to the synthetic model and two events are buried in noise, which cannot be identified. The signal-to-noise ratio (SNR) is -16.14 dB. For comparison, we use standard velocity scan and chaotic system with GPM to calculate velocity spectra of noisy data, respectively. Velocity analysis with stacking criterion fails in providing correct velocity trends. The results of the chaotic system with the GPM velocity analysis method shows higher resolution and random noise has less influence on the low velocity zone, which may cause the incorrect velocity picking. We use a marine field dataset with swell noise to further evaluate the proposed method. For standard velocity scan, the swell noise causes low quality of velocity panels, which leads to inaccurate velocity picking. The Duffing chaotic system is less affected by the strong swell noise and shows more accurate velocity picking results. The velocity analysis method using the Duffing chaotic system produces a stacking result with high SNR.Because the Duffing chaotic system is less sensitive to strong noise, so suitable for detecting weak seismic events. A new approach called grid partition method (GPM) can judge state of the Duffing oscillator system, which provides a stable criterion to detect accurate RMS velocities. To reconstruct the periodic seismic signal, we use an intercept window to select seismic wavelets from the CMP gathers. Where after, the Duffing chaotic system provides a better way for acquiring velocity information in the case strong random noise exists. Results of applying the proposed method to a synthetic example and field data show that this new method, compared with traditional velocity analysis, can obtain more accurate velocities in a strong random noise environment.

Chaotic system; Duffing oscillator; Strong random noise; Swell noise; Velocity analysis

10.6038/cjg20150620.

國(guó)家自然科學(xué)基金重點(diǎn)基金(41430322)，國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)課題(2013CB429805)，國(guó)家自然科學(xué)基金項(xiàng)目(41274119，41174080，41004041)，以及國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)重大項(xiàng)目(2012AA09A20103)資助.

劉財(cái)，男，教授，博士生導(dǎo)師.主要從事地震波場(chǎng)正反演理論、綜合地球物理等研究.E-mail:liucai@jlu.edu.cn

*通訊作者 劉洋，男，教授，博士生導(dǎo)師.主要從事地震數(shù)據(jù)處理、海洋電磁數(shù)據(jù)處理和地質(zhì)-地球物理綜合研究等工作. E-mail:yangliu1979@jlu.edu.cn

10.6038/cjg20150620

P631

2014-09-01，2015-01-18收修定稿

劉財(cái)，王博，劉洋.2015.基于Duffing振子混沌系統(tǒng)的地震速度分析方法.地球物理學(xué)報(bào)，58(6):2057-2068,

Liu C, Wang B, Liu Y. 2015. Seismic velocity analysis based on the Duffing oscillator chaotic system.ChineseJ.Geophys. (in Chinese),58(6):2057-2068,doi:10.6038/cjg20150620.