基于獨(dú)立分量分析的地震信號(hào)盲源分離及應(yīng)用

孟會(huì)杰，蘇 勤，曾華會(huì)，徐興榮，劉 桓，張小美

（中國石油勘探開發(fā)研究院西北分院，蘭州 730020）

0 引言

實(shí)際采集的地震資料受野外采集條件的影響，通常包含嚴(yán)重的噪聲干擾，資料信噪比較低，成像精度差。因此，地震資料的去噪優(yōu)化處理對(duì)油氣勘探至關(guān)重要［1］。提高地震資料信噪比是整個(gè)資料處理流程的關(guān)鍵，尤其是針對(duì)目前的復(fù)雜油氣藏勘探，必須保證一定的信噪比。然而，由于實(shí)際采集的地震資料包含的噪聲類型復(fù)雜多樣，因此，在地震資料處理過程中尋找有效的去噪方法至關(guān)重要。

許多專家學(xué)者根據(jù)噪聲和有效信號(hào)的特點(diǎn)，提出了不同的噪聲去除方法來對(duì)地震資料進(jìn)行處理。例如，基于域變換和濾波理論的噪聲壓制方法，包括t-x域、f-x濾波、f-k濾波、K-L變換等，這類方法主要是通過將待處理信號(hào)轉(zhuǎn)換到某種域空間來消除噪聲［2-5］。還有基于稀疏變換的去噪方法，如小波變換、Randon 變換、Curvelet 變換等［6-8］，這些方法主要是通過將地震數(shù)據(jù)轉(zhuǎn)換為稀疏域，然后對(duì)稀疏系數(shù)作閾值處理，再將處理后的稀疏系數(shù)轉(zhuǎn)換到時(shí)空域達(dá)到去噪的目的［8-13］。通常這種變換為了達(dá)到所需要的去噪效果，在采用閾值函數(shù)對(duì)稀疏系數(shù)進(jìn)行處理時(shí)，一般都是濾除小系數(shù)，保留少數(shù)較大的稀疏系數(shù)，即可反映原始數(shù)據(jù)的主要特征，這類方法也經(jīng)常被用于實(shí)際資料處理中。然而，利用此類算法進(jìn)行處理時(shí)，確定合適的閾值或者濾波器算子對(duì)去噪效果具有較大的影響，如果閾值較小或者濾波器算子較短，信噪分離程度低，處理后的資料中仍會(huì)存在較多噪聲；反之，如果完全將噪聲從資料中去除，難免也會(huì)傷及有效信號(hào)。

Herault 等［14］在一次計(jì)算機(jī)神經(jīng)網(wǎng)絡(luò)大會(huì)上首次提出盲源分離方法。獨(dú)立分量分析（ICA）算法是隨盲源信號(hào)分析理論發(fā)展而來的一種獨(dú)立源信號(hào)提取算法，該概念最早是由Comon［15］提出。ICA 算法的基本思想是假設(shè)源信號(hào)相互獨(dú)立，根據(jù)混合信號(hào)的高階統(tǒng)計(jì)量信息，可以提取混合信號(hào)的各獨(dú)立分量，并最終得到相互獨(dú)立的源信號(hào)的最佳近似值。后經(jīng)不斷改進(jìn)，該算法在各領(lǐng)域得到廣泛應(yīng)用。實(shí)際采集的野外地震資料，其有效信號(hào)和隨機(jī)噪聲是由不同信號(hào)源產(chǎn)生的，滿足相互獨(dú)立的特點(diǎn)。因此，可以利用ICA 算法對(duì)地震資料隨機(jī)噪聲進(jìn)行壓制。國內(nèi)許多專家學(xué)者對(duì)ICA 算法進(jìn)行了改進(jìn)，并應(yīng)用于地震資料噪聲壓制中，取得了不錯(cuò)的效果。劉喜武等［16-17］利用ICA 方法，在地震信息處理中進(jìn)行了初步的應(yīng)用研究，發(fā)現(xiàn)該方法在地震信號(hào)特征提取與信號(hào)分離方面具有應(yīng)用前景；呂文彪等［18］對(duì)ICA 方法進(jìn)行了改進(jìn)，利用兩部特征值分解法去除噪聲；李大衛(wèi)等［19］為提高分離信號(hào)的質(zhì)量，引入模擬退火算法，大幅度提高了地震資料的信噪比；張念等［20］對(duì)ICA 算法數(shù)學(xué)模型進(jìn)行求解分析，并利用FastICA 算法對(duì)仿真地震信號(hào)進(jìn)行了去噪測(cè)試；張銀雪等［21］提出一種步長(zhǎng)自適應(yīng)ICA 算法，實(shí)現(xiàn)了對(duì)疊前地震信號(hào)的去噪；王維強(qiáng)等［22］利用EMD 對(duì)信號(hào)自適應(yīng)分解的優(yōu)點(diǎn)，結(jié)合ICA 算法，實(shí)現(xiàn)了隨機(jī)噪聲和有效信號(hào)的分離；袁星虎等［23］將五階收斂速度的牛頓迭代式引進(jìn)FastICA 算法，實(shí)現(xiàn)了對(duì)混合地震信號(hào)的分離；逯宇佳等［24］將ICA 算法與動(dòng)態(tài)時(shí)間規(guī)整（DTW）算法相結(jié)合，實(shí)現(xiàn)了對(duì)非平緩地層中隨機(jī)噪聲的壓制。

標(biāo)準(zhǔn)的ICA 模型要求觀測(cè)信號(hào)數(shù)大于等于源信號(hào)數(shù)，且傳統(tǒng)的ICA 算法是基于相同道或者相似道數(shù)據(jù)的多次觀測(cè)，而實(shí)際地震資料通常是由成千上萬的道記錄組成，直接利用ICA 算法進(jìn)行去噪，效果不一定理想。實(shí)際地震資料道數(shù)較多，每一道地震數(shù)據(jù)可以看成是一個(gè)獨(dú)立的一維信號(hào)，如果針對(duì)每一道進(jìn)行單獨(dú)處理，則觀測(cè)信號(hào)數(shù)小于源信號(hào)數(shù)，標(biāo)準(zhǔn)ICA 模型轉(zhuǎn)化為一個(gè)欠定問題的求解，也就是從一個(gè)觀測(cè)信號(hào)中分離出多個(gè)源信號(hào)。相空間重構(gòu)是隨混沌理論發(fā)展而來一種信號(hào)重構(gòu)技術(shù)，是非線性動(dòng)力學(xué)去噪方法的理論基礎(chǔ)，被廣泛應(yīng)用于低維信號(hào)到多維信號(hào)的重構(gòu)。根據(jù)Takens 嵌入定理［25］，一維信號(hào)經(jīng)過重構(gòu)后可以轉(zhuǎn)換到等價(jià)的多維相空間中，并可反映原始信號(hào)的系統(tǒng)動(dòng)力學(xué)特征，其關(guān)鍵在于相空間重構(gòu)參數(shù)（嵌入維數(shù)和時(shí)間延遲）的選取，準(zhǔn)確地重構(gòu)參數(shù)可以保證重構(gòu)后的相空間能充分反映原始信號(hào)的動(dòng)力學(xué)特征。重構(gòu)參數(shù)選取的方法也比較多，在實(shí)際數(shù)據(jù)處理中，可以通過一個(gè)固定值來確定2 個(gè)參數(shù)，該固定值被稱為時(shí)間窗，可以表示為：tw=(m-1)τ，其中m為嵌入維數(shù)，τ為時(shí)間延遲，通過一個(gè)參數(shù)的變化來確定另一個(gè)參數(shù)，從而確定最佳的2 個(gè)參數(shù)組合。然而，這2 個(gè)參數(shù)的選取不應(yīng)該是相互獨(dú)立的，而應(yīng)該與時(shí)間窗有關(guān)，且隨著時(shí)間窗的變化而變化［26］。1999 年Kim 等［27］提出了一種C-C 算法，該算法利用和嵌入維數(shù)有關(guān)的關(guān)聯(lián)積分，并對(duì)其進(jìn)行求解，可以同時(shí)確定這2 個(gè)重構(gòu)參數(shù)，同時(shí)，該方法計(jì)算量相對(duì)較小，而且獲得的結(jié)果也更為準(zhǔn)確。

本文提出一種基于相空間重構(gòu)的地震信號(hào)單通道盲源分離算法，來對(duì)地震資料進(jìn)行處理。首先，利用C-C 算法確定信號(hào)的重構(gòu)參數(shù)；然后，利用相空間重構(gòu)技術(shù)將一維信號(hào)進(jìn)行重構(gòu)；最后，對(duì)重構(gòu)后的高維信號(hào)利用ICA 算法結(jié)合數(shù)據(jù)自身的高階統(tǒng)計(jì)特性，對(duì)有效信號(hào)和噪聲進(jìn)行分離，進(jìn)而提取有效信號(hào)記錄，從而達(dá)到保真去噪的目的，為后期精細(xì)解釋提供質(zhì)量更高的處理結(jié)果。。

1 方法原理

1.1 相空間重構(gòu)

野外采集的地震信號(hào)一般是由有效信號(hào)和噪聲2 個(gè)部分組成，因此，每一道地震信號(hào)也都可以看作是一個(gè)由噪聲和有效信號(hào)組成的一維時(shí)間序列。針對(duì)每一道地震信號(hào)，利用ICA 算法進(jìn)行處理是一個(gè)欠定問題（觀測(cè)信號(hào)數(shù)小于源信號(hào)數(shù)），所以要對(duì)一維信號(hào)進(jìn)行重構(gòu)。本文采用相空間重構(gòu)技術(shù)，將一維信號(hào)重構(gòu)到高維的相空間中，相空間重構(gòu)的過程可簡(jiǎn)單描述為：假設(shè)時(shí)間序列為x(ti),i=1,2,…,n,n為采樣點(diǎn)數(shù)；令時(shí)間延遲為τ，嵌入維數(shù)為m，根據(jù)Takens 嵌入定理，則可以重構(gòu)一個(gè)m維相空間：

式中：X(·)表示相空間中的相點(diǎn)；N是相空間中的相點(diǎn)數(shù)，N=n-(m-1)τ。

1.2 重構(gòu)參數(shù)選取

重構(gòu)參數(shù)對(duì)相空間重構(gòu)結(jié)果具有很大的影響，進(jìn)而直接影響到最終去噪的效果。為了確定準(zhǔn)確的重構(gòu)參數(shù)，本文采用C-C 算法［27］對(duì)重構(gòu)參數(shù)進(jìn)行計(jì)算，計(jì)算過程可描述為

式中：C為定義的關(guān)聯(lián)積分表達(dá)式；Xi為重構(gòu)的相空間中第i個(gè)相點(diǎn)。

定義如下函數(shù)：

對(duì)時(shí)間序列進(jìn)行不同劃分來確定合適的時(shí)間延遲，具體的劃分方法如下：實(shí)際觀測(cè)的某一時(shí)間序列x={x1,x2,…,xn}，當(dāng)時(shí)間序列被劃分為一個(gè)子序列，即t=1 時(shí)，劃分后的時(shí)間序列就是該時(shí)間序列本身。

當(dāng)t=2 時(shí)，原時(shí)間序列可被劃分為{x1,x3,…,xn-1}和{x2,x4,…,xn} 2 個(gè)子序列，此時(shí)有

當(dāng)時(shí)間序列被劃分為t個(gè)子序列時(shí)，則有

當(dāng)n→∞時(shí)，

定義r對(duì)應(yīng)最大值和最小值的差量：

σ為時(shí)間序列的標(biāo)準(zhǔn)差，根據(jù)式（10），可以同時(shí)確定時(shí)間序列的重構(gòu)參數(shù)，式（10-Ⅱ）的第一個(gè)極小值點(diǎn)的對(duì)應(yīng)值就是所要選取的時(shí)間延遲值，根據(jù)式（10-Ⅲ）的最小值點(diǎn)可以確定對(duì)應(yīng)的時(shí)間窗口，根據(jù)選取的時(shí)間延遲和確定的時(shí)間窗，利用公式tw=(m-1)τ，可以得到合適的嵌入維數(shù)值，進(jìn)而將低維時(shí)間序列重構(gòu)到高維相空間中。

1.3 獨(dú)立分量分析

獨(dú)立分量分析（ICA）算法是隨盲源分析理論發(fā)展而來的一種獨(dú)立源信號(hào)提取算法。該算法要求觀測(cè)信號(hào)數(shù)，也就是待處理的信號(hào)數(shù)要大于源信號(hào)的個(gè)數(shù)，且源信號(hào)需相互獨(dú)立，服從非高斯性分布。地震資料中的有效信號(hào)和隨機(jī)噪聲通常是相互獨(dú)立的，另外，地震子波通常也是非高斯分布的。根據(jù)測(cè)井?dāng)?shù)據(jù)的統(tǒng)計(jì)結(jié)果可以得出，有效波的反射系數(shù)是超高斯分布的［28］，地震信號(hào)是地層反射系數(shù)與地震子波褶積的結(jié)果。因此，地震信號(hào)恰好滿足這一特點(diǎn)。所以，可以利用該算法對(duì)地震資料進(jìn)行處理。

ICA 算法的去噪過程可以簡(jiǎn)單描述如下：

假設(shè)有n個(gè)獨(dú)立的源信號(hào)可以表示為：S=[s1,s2,…,sn]，通過一個(gè)混合矩陣A，A=[a1,a2,…,an]與S進(jìn)行混合，可以得到混合后的觀測(cè)信號(hào)X，其表示為

式（11）就是一個(gè)標(biāo)準(zhǔn)的ICA 模型，實(shí)際情況下，觀測(cè)信號(hào)是已知的，混合矩陣和源信號(hào)是未知的。假設(shè)混合矩陣是一個(gè)可逆矩陣，那么，ICA 算法的本質(zhì)就是求解一個(gè)分離矩陣W，也就是混合矩陣的逆矩陣，來實(shí)現(xiàn)對(duì)觀測(cè)信號(hào)X的分離，得到分離信號(hào)

因此ICA 算法的核心就是分離矩陣的求取，在源信號(hào)和混合矩陣未知的情況下，首先對(duì)觀測(cè)信號(hào)進(jìn)行預(yù)處理包括白化和中心化，中心化的目的是使得每一個(gè)觀測(cè)信號(hào)的均值為零，白化是指對(duì)經(jīng)過中心化處理后的多維信號(hào)，通過一個(gè)變換矩陣進(jìn)行線性變換處理，變?yōu)槊總€(gè)分量具有單位方差且互不相關(guān)的信號(hào)，對(duì)信號(hào)進(jìn)行白化處理的目的是使分離算法收斂速度更快并且更穩(wěn)定。

在求解分離矩陣時(shí)，選取負(fù)熵作為非高斯性的度量方式，可以表示為

式中：z是經(jīng)過白化和中心化預(yù)處理后的觀測(cè)信號(hào)；zgauss為與矩陣z具有相同協(xié)方差矩陣的高斯隨機(jī)向量；H（·）為隨機(jī)變量的微分熵；pz（·）為隨機(jī)變量的概率。

由于其不易求取，因此，式（13）和式（14）可分別變換為

v是零均值、單位方差的高斯變量，G（·）的取值類型比較多，式（16）僅是其中的2 種，對(duì)不同類型的信號(hào)取值不同，對(duì)式（15）進(jìn)行推導(dǎo)，可得到ICA算法的迭代式［20］為

式中：g和g′分別是G（·）的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)；w為分離矩陣W的行向量。

設(shè)定一個(gè)迭代終止條件和迭代次數(shù)，通過更新迭代求解得到所有的w，即得到最終要求的分離矩陣W。圖1 是ICA 算法的處理流程圖。

圖1 ICA 算法流程圖Fig.1 Workflow of ICA algorithm

1.4 基于相空間重構(gòu)的地震信號(hào)盲源分離算法應(yīng)用

利用本文方法對(duì)實(shí)際地震資料進(jìn)行處理，整個(gè)處理過程可簡(jiǎn)單描述如下：

（1）針對(duì)每一道地震信號(hào)，利用C-C 算法計(jì)算信號(hào)的重構(gòu)參數(shù)——嵌入維數(shù)和時(shí)間延遲；

（2）利用計(jì)算的重構(gòu)參數(shù)，根據(jù)Takens 嵌入定理，利用相空間重構(gòu)技術(shù)對(duì)地震信號(hào)進(jìn)行重構(gòu)，將一維信號(hào)重構(gòu)到高維相空間中，從而滿足ICA 算法的處理?xiàng)l件；

（3）對(duì)重構(gòu)后的數(shù)據(jù)分別進(jìn)行中心化和白化預(yù)處理；

（4）對(duì)預(yù)處理后的數(shù)據(jù)，利用ICA 算法進(jìn)行信噪分離處理，最終得到分離后的噪聲信號(hào)和有效信號(hào)，達(dá)到保真去噪的目的。

2 合成數(shù)據(jù)測(cè)試

為了驗(yàn)證本文方法的可行性，利用該方法分別對(duì)合成數(shù)據(jù)和實(shí)際數(shù)據(jù)進(jìn)行處理，合成數(shù)據(jù)分為2個(gè)，一個(gè)簡(jiǎn)單的一維信號(hào)和一個(gè)簡(jiǎn)單的二維線性模型，一維合成信號(hào)［圖2（a）］子波采用雷克子波，合成信號(hào)的采樣點(diǎn)數(shù)為2 000 個(gè)，采樣間隔為1 ms。對(duì)原始合成信號(hào)添加隨機(jī)噪聲［圖2（b）］，得到含噪的地震記錄［圖2（c）］；然后，利用本文方法對(duì)含噪信號(hào)進(jìn)行去噪處理，分別得到去噪后的有效信號(hào)［圖2（d）］和去除的噪聲［圖2（e）］。

圖2 單道合成數(shù)據(jù)測(cè)試Fig.2 Single channel synthetic data test

為得到較好的去噪效果，需要確定合適的重構(gòu)參數(shù)，利用C-C 算法對(duì)單道含噪信號(hào)計(jì)算重構(gòu)參數(shù)（圖3）。圖3（a）第一個(gè)極小值點(diǎn)對(duì)應(yīng)的τ值為計(jì)算的時(shí)間延遲τ=5；通過圖3（b）確定時(shí)間窗τw=15；利用時(shí)間延遲、嵌入維數(shù)以及時(shí)間窗三者之間的關(guān)系：tw=(m-1)τ，計(jì)算得到嵌入維數(shù)m=4。利用計(jì)算的重構(gòu)參數(shù)，根據(jù)Takes 嵌入定理，可以對(duì)原信號(hào)進(jìn)行相空間重構(gòu)（圖4）。

圖3 C-C 算法計(jì)算結(jié)果Fig.3 Calculation results of C-C algorithm

圖4 相空間重構(gòu)結(jié)果Fig.4 Result of phase space reconstruction

本文又利用二維線性模型進(jìn)行去噪測(cè)試，二維線性模型如圖5（a）所示，對(duì)該模型數(shù)據(jù)添加一定程度的噪聲［圖5（b）］，然后，利用本文方法對(duì)含噪數(shù)據(jù)進(jìn)行去噪處理，信噪分離結(jié)果分別如圖5（c）和圖5（d）所示。從分離的有效信號(hào)和噪聲結(jié)果上可以看到，隨機(jī)噪聲得到了有效去除，有效信號(hào)更容易識(shí)別，資料的信噪比以及分辨率也得到了大幅度的提高。

圖5 模擬數(shù)據(jù)去噪結(jié)果Fig.5 Denoising results of simulated data

3 實(shí)際資料處理

在實(shí)際資料處理中，隨機(jī)噪聲的去除方法也比較多。大多數(shù)去噪方法都是在資料具有一定的信噪比基礎(chǔ)上進(jìn)行隨機(jī)噪聲壓制，而且在處理過程中，由于容易傷及有效信號(hào)，所以，針對(duì)隨機(jī)噪聲的壓制一般都比較慎重，尤其是當(dāng)資料信噪比較低時(shí)，有效地去除隨機(jī)噪聲更為困難。

本文選取四川盆地某地區(qū)三維資料的部分疊加剖面，從剖面整體來看，淺層噪聲比較發(fā)育，部分噪聲信號(hào)能量強(qiáng)［圖6（a）］。利用f-xRNA 隨機(jī)噪聲壓制方法對(duì)該資料進(jìn)行處理后，疊加剖面上部分噪聲得到了有效壓制，信噪比較之前也得到了一定的改善［圖6（b）］，但是，還是可以明顯看到部分異常噪聲。使用本文方法進(jìn)行處理后的效果相比RNA要更好一些，噪聲去除相對(duì)更徹底，尤其是針對(duì)有效信號(hào)比較弱的淺層區(qū)域，經(jīng)過處理后有效信號(hào)更加突出，資料信噪比、分辨率更高［圖6（c）］。此外，從去除的噪聲剖面［圖6（d）］上來看，利用本文方法在消除噪聲的同時(shí)，對(duì)有效信號(hào)幾乎沒有損傷。

圖6 實(shí)際資料處理結(jié)果Fig.6 Processing results of actual data

4 結(jié)論

（1）單通道盲源信號(hào)分離是一種欠定問題求解，已知信息較少，需要估計(jì)的未知因素比較多，本文基于相空間重構(gòu)理論，可以將低維信號(hào)重構(gòu)到高維相空間，進(jìn)而滿足多通道盲源分離算法的要求，易于求解。

（2）重構(gòu)參數(shù)的選取對(duì)相空間重構(gòu)具有重要影響，進(jìn)而影響最終的信噪分離效果，本文采用C-C算法，根據(jù)信號(hào)自身的特征，可以準(zhǔn)確地計(jì)算時(shí)間延遲和嵌入維數(shù)的值，進(jìn)而為相空間重構(gòu)提供準(zhǔn)確合適的重構(gòu)參數(shù)。

（3）根據(jù)重構(gòu)相空間中噪聲和有效信號(hào)重構(gòu)軌道的幾何特征差異，利用ICA 算法，結(jié)合數(shù)據(jù)自身的高階統(tǒng)計(jì)特性，可以有效的實(shí)現(xiàn)信噪分離，合成數(shù)據(jù)和實(shí)際數(shù)據(jù)處理結(jié)果證明，該方法可以有效分離源信號(hào)，達(dá)到保真去噪的目的。