抑制地下流體觀測數(shù)據(jù)干擾的期望最大化算法?

周 洋 王 俊

(1.中國地震局地震研究所(地震大地測量重點(diǎn)實(shí)驗(yàn)室)，湖北 武漢 430071；2.湖北省地震局，湖北 武漢 430071；3.中南民族大學(xué)計(jì)算機(jī)科學(xué)學(xué)院，湖北 武漢 430071)

湖北區(qū)域地球物理臺網(wǎng)目前共計(jì)23 個(gè)臺站83套儀器，其中地下流體儀器41 套，SWY－II 型數(shù)字式水位儀9 套、SZW－1A 及SZW－II 數(shù)字式溫度計(jì)15套、RTP－II 型氣溫氣壓雨量綜合觀測儀3 套、WYY－1型氣象三要素觀測儀14 套[1]。地下流體儀器種類繁多、在網(wǎng)運(yùn)行數(shù)量也是在各觀測手段(重力、形變、電磁、流體)中居首，同時(shí)受干擾因素眾多且難以區(qū)分。故尋求一種應(yīng)用前景更為廣闊的數(shù)據(jù)干擾抑制技術(shù)顯得尤為重要。

近年來興起的小波分析方法以其良好的時(shí)頻分析能力迅速成為非平穩(wěn)信號處理的有力工具，基于小波分析的去噪方法更是大量涌現(xiàn)，且被證明其具有傳統(tǒng)傅里葉變換去噪方法所不能比擬的優(yōu)越性[2－7]。去除地震信號中隨機(jī)噪聲技術(shù)雖然已經(jīng)有了很大發(fā)展，但仍然不能滿足人們對地震數(shù)據(jù)去噪效果的要求，因此改進(jìn)現(xiàn)有去噪技術(shù)，尋求更加精確、有效、快速、便捷的去噪算法，是一項(xiàng)十分有意義、尚待研究的課題[8－12]。目前國內(nèi)在數(shù)字信號處理如降噪方面的技術(shù)已日趨成熟，如文獻(xiàn)[13－16]介紹了基于時(shí)空域數(shù)字信號處理的地震資料降噪方法，而文獻(xiàn)[17－21]則研究了基于深度學(xué)習(xí)的自適應(yīng)學(xué)習(xí)去噪方法。

本文提出一種基于神經(jīng)網(wǎng)絡(luò)的期望最大化算法，用于地震地下流體觀測數(shù)據(jù)噪聲抑制。該算法較其他小波算法及神經(jīng)網(wǎng)絡(luò)算法去噪在保留原始數(shù)據(jù)信號能量前提下，更能快速、準(zhǔn)確定位并抑制噪聲干擾。與常規(guī)幅度閾值相比，所得到的概率閾值更適合于數(shù)據(jù)。它通過概率閾值為用戶提供了量化是否將大幅度異常視為噪聲的置信度的可能性。該方法顯示出比優(yōu)化的常規(guī)方法稍好的性能，且參數(shù)測試和變化少得多。

1 期望最大化算法

1.1 簡介

地震數(shù)據(jù)總是會(huì)因各種噪聲而損壞，從而降低了數(shù)據(jù)質(zhì)量。因此，噪聲衰減是有助于解釋地震數(shù)據(jù)處理的重要步驟。可以通過基于模型的信號處理來分離信號和噪聲。

該方法假定信號和噪聲具有可以由特定數(shù)學(xué)模型捕獲的不同特性。該模型可以基于諸如波動(dòng)理論之類的物理學(xué)定律，該定律被用于許多分析方法中。或者，模型可以利用確定性特征，例如不同的跡線偏移或統(tǒng)計(jì)屬性(例如獨(dú)立性)以區(qū)分信號和噪聲。

我們關(guān)注高振幅噪聲衰減的問題。考慮一個(gè)簡單的噪聲模型，其中所有大于特定閾值的數(shù)據(jù)樣本都可能是噪聲。閾值可以以多種形式指定，但通常用閾值因子乘以某些數(shù)據(jù)統(tǒng)計(jì)量(例如，分析窗口內(nèi)計(jì)算的平均值，中位數(shù)或均方根值)來表示。閾值因子和數(shù)據(jù)統(tǒng)計(jì)信息是用戶定義的參數(shù)。通常由于其簡單性而調(diào)用此模型。但是，其主要缺點(diǎn)是閾值的頻繁變化，導(dǎo)致整個(gè)數(shù)據(jù)噪聲功率的變化。

我們提出了一種自動(dòng)閾值確定技術(shù)，用于檢測選擇域中的大振幅噪聲。以頻率偏移域中的衰減噪聲衰減為例，但結(jié)果可以推廣到其他類型的噪聲，例如衍射的多個(gè)噪聲和尖峰，其中在某個(gè)時(shí)間或時(shí)間偏移中會(huì)出現(xiàn)大幅度采樣窗口被檢測到，然后通過插值刪除。我們的技術(shù)與常規(guī)方法相比，它具有更好的數(shù)據(jù)適應(yīng)性。

首先，我們簡要回顧一下涌浪噪聲的問題。然后，我們概述了自動(dòng)閾值確定的新方法。最后，我們結(jié)合湖北地震臺地下流體數(shù)據(jù)給出一些真實(shí)的數(shù)據(jù)示例。

1.2 降噪衰減

涌浪噪聲是由惡劣的天氣條件引起的，并且是獲取地震數(shù)據(jù)時(shí)經(jīng)常遇到的問題。它對地震數(shù)據(jù)質(zhì)量有不利影響，甚至可能導(dǎo)致暫停采集。它的特點(diǎn)是振幅大，頻率低。

討論了產(chǎn)生涌浪噪聲的幾種可能機(jī)制，得出結(jié)論，對于現(xiàn)代的泡沫填充拖纜，最可能的原因是強(qiáng)烈的海浪或拖纜表面上的動(dòng)態(tài)壓力變化引起海洋垂直運(yùn)動(dòng)，從而引起靜水壓力波動(dòng)，這是由于拖纜周圍存在湍流層而導(dǎo)致的。

抑制涌浪噪聲的常規(guī)技術(shù)首先計(jì)算滑動(dòng)窗口內(nèi)數(shù)據(jù)的幅度或功率譜。窗口內(nèi)給定頻率上超過某個(gè)閾值的所有頻譜值均被視為噪聲。然后對這些噪聲樣本進(jìn)行衰減或內(nèi)插。但是，在此過程中必須定義一個(gè)閾值。

對于給定的頻率f，我們用rk＝｜Dk(f)｜2定義Sn＝{r1，r2，…，rn}，其中Dk(f)是數(shù)據(jù)窗口中第k條跡線的傅里葉變換。因此，Sn代表給定頻率下的一組跡線的功率譜值。閾值通常計(jì)算為

1.3 離群值檢測及閾值計(jì)算

本節(jié)介紹了一種通過檢測異常值來確定適當(dāng)閾值的自動(dòng)方法。

考慮一個(gè)數(shù)據(jù)集Sn＝{r1，r2，…，rn}，其中樣本rk被假定為獨(dú)立且均勻分布的，由概率密度函數(shù)(PDF)g(r)生成。在應(yīng)用統(tǒng)計(jì)中，在集合Sn中查找異常樣本稱為異常值檢測。目的是找到那些與其余數(shù)據(jù)顯示不同統(tǒng)計(jì)特性的樣本。在我們的問題中，離群值(噪聲＋信號)的幅度大于常規(guī)數(shù)據(jù)(僅信號)的幅度。因此，異常值與常規(guī)數(shù)據(jù)的總體不同之處在于，可以通過使用PDF 建模來捕獲一些明顯的統(tǒng)計(jì)量。

便于標(biāo)記，離群點(diǎn)的PDF 用p(r｜θ1)表示，常規(guī)數(shù)據(jù)的PDF 用p(r｜θ0)表示。為簡單起見，假設(shè)兩個(gè)分布都屬于同一系列的參數(shù)化PDF，但具有不同的參數(shù)值，即θ1≠θ0。讓標(biāo)量表示集合Sn中離群值的分?jǐn)?shù)。該參數(shù)的統(tǒng)計(jì)含義是:從Sn隨機(jī)抽取的樣本ri是離群值的先驗(yàn)概率(即前數(shù)據(jù)建模)。現(xiàn)在，數(shù)據(jù)PDF 包含兩種模型的混合結(jié)構(gòu):

可以從Sn估算方程(2)中定義的模型的所有參數(shù)。估算完成后，可以使用貝葉斯規(guī)則計(jì)算給定離群值樣本r＝ri的后驗(yàn)概率(即后數(shù)據(jù)建模)。在貝葉斯規(guī)則中，假設(shè)事件B已發(fā)生，則A和B是兩個(gè)事件，而P(A｜B)是事件A發(fā)生的概率。它表示兩個(gè)事件的條件似然，例如P(A｜B)P(B)＝P(B｜A)P(A)。這里，r是離群值。

如果離群值ri的后驗(yàn)概率大于概率閾值β，即

例如，對于β＝0.5，式(4)表示僅當(dāng)模型確定至少有50%是異常值時(shí)才選擇異常值。式(4)中概率閾值β的選擇比式(1)中任意閾值因子α客觀得多，并且反映了我們可能需要將任何數(shù)據(jù)分類為異常值的統(tǒng)計(jì)置信度。

計(jì)算閾值時(shí)，涉及三個(gè)步驟:選擇一個(gè)PDF，識別算法并解釋結(jié)果。

為了使用式(3)給出的檢測標(biāo)準(zhǔn)，需要選擇式(2)中p(r｜θ)的參數(shù)形式。當(dāng)r增加時(shí)，p(r｜θ)的形式應(yīng)提供不同的衰減率。需要這樣做，以確保與常規(guī)數(shù)據(jù)相比，離群值人群中大振幅值的可能性更大。使用傅里葉變換的實(shí)部和虛部是具有等方差的零均值高斯隨機(jī)變量的近似值，可以證明功率譜樣本的PDF 具有指數(shù)形式。

指數(shù)分布具有一個(gè)參數(shù)，θ＝λ，它是分布的平均值:

然后，式(2)中的模型變?yōu)?/p>

用λ1>λ0表示一個(gè)事實(shí)，即異常值的分布更可能產(chǎn)生較大的振幅。

在使用方程式進(jìn)行推論之前，人們可能會(huì)質(zhì)疑模型在式(6)中的擬合優(yōu)度。這個(gè)問題的答案很重要，因?yàn)樗枰獙υ僭O(shè)進(jìn)行統(tǒng)計(jì)檢驗(yàn):數(shù)據(jù)遵循式(6)中的模型。數(shù)學(xué)模型有助于我們理解和描述現(xiàn)實(shí)世界，并且公認(rèn)的是，真正的數(shù)據(jù)生成模型比任何可能的假設(shè)模型都要復(fù)雜得多。因此，模型的正確性問題尚未正式解決，但本文稍后將討論模型敏感性測試。

式(6)中的模型具有三個(gè)未知參數(shù)(λ0，λ1，ε)，將從Sn進(jìn)行估計(jì)。我們建議在以下優(yōu)化問題中使用最大似然估計(jì)器(MLE)，以獲得其所需的統(tǒng)計(jì)屬性，例如一致性和效率:

此問題沒有封閉形式的解決方案，因此參數(shù)使用迭代程序估算值。此過程屬于一類稱為期望最大化(EM)算法的技術(shù)，用于在數(shù)據(jù)丟失或不完整時(shí)查找模型參數(shù)的MLE。迭代EM 算法由以下步驟給出:

(2)如下更新參數(shù)

在該特定應(yīng)用中，EM 算法的收斂速度很快，并且在很大程度上與初始條件無關(guān)。設(shè)置初始條件，使得等于[ε(0)n]個(gè)最大振幅的平均值等于其余數(shù)據(jù)的平均值。

在將未知參數(shù)λ0，λ1和ε替換為其估計(jì)值的同時(shí)，結(jié)合式(3)、式(4)和式(6)會(huì)得出基于幅度的閾值檢測標(biāo)準(zhǔn)。幅度r是概率為β(如果r>rexp)的離群值，其中exp 表示表達(dá)式中的指數(shù)。

式(12)中給出的閾值考慮了數(shù)據(jù)的統(tǒng)計(jì)信息(通過式(6)中模型參數(shù)的估計(jì))和β值的置信度要求。當(dāng)log(β/(1－β))＝0，β＝0.5，用戶對閾值的影響最小。這與以下事實(shí)相吻合:對于要做出的任何二元決策，每個(gè)概率為0.5 的結(jié)果表示信息量不足的情況下，該概率的先驗(yàn)知識不會(huì)影響決策。因此，對于建議的概率閾值，β＝0.5 將是更好的默認(rèn)值。

對式(12)中定義的新閾值標(biāo)準(zhǔn)的分析揭示了幾點(diǎn)。首先，當(dāng)用戶增加值β時(shí)，接受異常值所需的統(tǒng)計(jì)置信度就會(huì)增加，從而會(huì)增加閾值水平。因此，β具有與常規(guī)方法的閾值因子α類似的作用。其次，在統(tǒng)計(jì)模型下，較小的表示存在少量異常值，并且增加閾值水平以反映該異常值。另一方面，當(dāng)時(shí)，模型會(huì)推斷出大多數(shù)數(shù)據(jù)都是異常值，因此閾值水平降低了。最后，當(dāng)時(shí)，方程(6)中兩個(gè)混合模型的假設(shè)沒有數(shù)據(jù)支撐。因此，的值沒有統(tǒng)計(jì)意義。在這種情況下，我們假定數(shù)據(jù)僅由常規(guī)樣本組成，因?yàn)榕懦藬?shù)據(jù)僅為異常值樣本的情況。閾值會(huì)自動(dòng)設(shè)置為較大的值，因此不會(huì)選擇異常值。

為了研究本技術(shù)對所選統(tǒng)計(jì)模型的敏感性，我們考慮在數(shù)據(jù)建模中使用瑞利分布而不是指數(shù)分布。

一旦通過式(12)中的閾值標(biāo)準(zhǔn)識別出噪聲幅度，就可以將它們從Sn中或從相鄰樣本中刪除(例如，通過設(shè)置ri＝0)。但是，我們通常用恒定因子對噪聲樣本進(jìn)行重新縮放以選擇更保守的選項(xiàng)，使得噪聲樣本的新均值等于指數(shù)模型的常規(guī)數(shù)據(jù)的均值，以及噪聲樣本的新均方根值等于瑞利模型的常規(guī)數(shù)據(jù)的均方根值。

2 實(shí)驗(yàn)及結(jié)果分析

本文在上述模型算法的基礎(chǔ)上，結(jié)合湖北省地震地下流體數(shù)據(jù)實(shí)例進(jìn)行仿真實(shí)驗(yàn)，并與傳統(tǒng)神經(jīng)網(wǎng)絡(luò)及小波去噪算法比較。

2.1 信號去噪

選取湖北省房縣三海村臺、鐘祥馬嶺臺、秭歸臺水溫儀觀測數(shù)據(jù)做仿真實(shí)驗(yàn)，時(shí)間跨度200 d。如圖1～圖3。

圖1 房縣三海村臺水溫?cái)?shù)據(jù)

圖2 鐘祥馬嶺臺水溫?cái)?shù)據(jù)

圖3 秭歸臺水溫?cái)?shù)據(jù)

從圖中得到的結(jié)果來看，期望最大化算法對噪聲有很好的抑制作用，保留信號能量的同時(shí)，去噪效果明顯，為實(shí)際工程應(yīng)用及地震數(shù)據(jù)分析預(yù)報(bào)提供了一個(gè)新的研究思路。實(shí)現(xiàn)仿真時(shí)間段可根據(jù)實(shí)際需要實(shí)時(shí)調(diào)整，系統(tǒng)可移植性、靈活性強(qiáng)。

2.2 與傳統(tǒng)算法的比較

同傳統(tǒng)的數(shù)據(jù)處理方法相比，小波濾波方法具有獨(dú)特的優(yōu)勢，即能夠在去除噪聲的同時(shí)，很好保留信號的突變部分或圖像的邊緣。而本文提出的一種新的期望最大化算法則比傳統(tǒng)的小波濾波算法更具優(yōu)越性。

以2020 年5 月荊門臺1 測點(diǎn)水溫?cái)?shù)據(jù)為例，給出數(shù)據(jù)信號在4 種不同的去噪算法處理下的結(jié)果，如圖4。

由圖4 可知本文提出的算法較傳統(tǒng)神經(jīng)網(wǎng)絡(luò)及小波算法(如模極大值法、空域相關(guān)法、小波閾值法)去噪相比，效果更理想。表1 為上述幾種算法的檢測性能比較，表2 給出幾種不同算法的定性比較。

圖4 幾種去噪方法比較

表1 幾種算法檢測的結(jié)果比較

表2 幾種濾波方法的定性比較

從SNR 和RMSE 的定義可知，SNR 值越大，檢測效果越好；RMSE 反之。能量成分比越小，檢測效果越好，標(biāo)準(zhǔn)殘差反之。而局部線性模型法很好地保留了信號發(fā)展初期的高頻特性，最大限度地反映了原信號本身的性質(zhì)，且性能參數(shù)優(yōu)于其他幾種算法。

3 結(jié)論

本文提出一種基于神經(jīng)網(wǎng)絡(luò)的期望最大化算法，用于地震地下流體觀測數(shù)據(jù)干擾抑制的新方法。結(jié)合湖北地震臺站地下流體數(shù)據(jù)實(shí)例做仿真實(shí)驗(yàn)，并與其他神經(jīng)網(wǎng)絡(luò)及小波分析方法做出比較。得出以下結(jié)論:

(1)使用統(tǒng)計(jì)建模可以自動(dòng)檢測大幅度的噪聲。最終的統(tǒng)計(jì)檢測方法等效于優(yōu)化的常規(guī)方法，但參數(shù)測試和變化少得多。噪聲衰減測試顯示，在整個(gè)數(shù)據(jù)范圍內(nèi)各種變化的噪聲水平下，結(jié)果都是一致的。該方法是通用的，可以推廣到其他類型的異常高振幅噪聲抑制，例如尖峰和衍射倍頻噪聲抑制。

(2)該算法對地下流體觀測數(shù)據(jù)信號消噪效果好，能有效抑制各種因素干擾，為臺站數(shù)據(jù)分析預(yù)報(bào)提供一種新的研究方向和方法，為儀器維護(hù)人員提供一種新的處理手段和途徑。

(3)該算法不依賴于具體的對象的數(shù)學(xué)模型，可移植性強(qiáng)。該算法可推廣應(yīng)用至其他觀測手段的數(shù)據(jù)干擾抑制，如形變、重力、電磁等。