張文征 唐 杰* 劉英昌 孟 濤 陳學(xué)國

（①中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院，山東青島266850；②中國石化勝利油田勘探開發(fā)研究院，山東東營257015）

0 引言

實際地震數(shù)據(jù)由于采集儀器或者環(huán)境等因素的影響，往往存在非平穩(wěn)隨機噪聲［1］，在地震數(shù)據(jù)中以脈沖形式出現(xiàn)，具有較強的能量，主要包括猝發(fā)脈沖、異常擾動、聲波干擾等［2］，很難用人工方法完全剔除野值和串狀干擾［3］。如果中、深層地震資料中存在強能量干擾，則疊前時間偏移剖面會出現(xiàn)不同程度的“畫弧”現(xiàn)象［4］?，F(xiàn)代地震資料處理的目標(biāo)是既能衰減噪聲、提高地震資料的信噪比，又能保留和增強有效不連續(xù)地震反射信息及儲層流體信息［5］。常規(guī)濾波方法常常放大了噪聲的影響，同時噪聲的存在也限制了分辨率的提升［6］，并“平滑”了地震數(shù)據(jù)中的不連續(xù)信息。通過自動識別地震記錄中的強能量干擾，確定噪聲出現(xiàn)的空間位置，根據(jù)定義的閾值和衰減系數(shù)用適當(dāng)?shù)姆绞綁褐品瞧椒€(wěn)隨機噪聲［7］。在動校正后或同相軸相對平直的地震剖面上，中值濾波器具有壓制強能量干擾的能力［8-9］，但很難選取合適的閾值算子兼顧保護有效信號和壓制噪聲［10］。由于異常噪聲既不連續(xù)也不相關(guān)，采用局部中值簡單地代替噪聲位置處的數(shù)據(jù)，容易壓制有效信號，尤其在同相軸彎曲和間斷處［11］。為此，Boashash等［12］提出了一種瞬時頻率估計算法——時頻峰值濾波算法，基于時頻分析理論壓制非平穩(wěn)隨機噪聲，將含噪信號調(diào)制為解析信號的瞬時頻率信號，然后取解析信號時頻分布的峰值估計有效信號，經(jīng)時頻峰值濾波后信號得到增強，噪聲被壓制。林紅波等［13］提出了基于絕對級差統(tǒng)計量（ROAD）的徑向時頻峰值濾波非平穩(wěn)隨機噪聲壓制方法，結(jié)合局部時頻峰值濾波和徑向時頻峰值濾波壓制地震非平穩(wěn)隨機噪聲。

近年來機器學(xué)習(xí)算法發(fā)展迅速，深度學(xué)習(xí)算法在數(shù)據(jù)去噪領(lǐng)域的應(yīng)用也越來越廣泛。Burger等［14］推出了具有三層網(wǎng)絡(luò)的多層感知機模型（MLP），雖然網(wǎng)絡(luò)結(jié)構(gòu)不深，但是處理效果較好。Schmidt等［15］、Chen等［16］分別提出聯(lián)合收縮方法（CSF）、反應(yīng)擴散模型（TNRD），這兩種方法的信噪比提高和邊緣保護效果超過了塊匹配方法（BM3D）。Mao等［17］利用機器學(xué)習(xí)算法實現(xiàn)了地震數(shù)據(jù)處理、屬性提取、斷層預(yù)測等方面的云地震數(shù)據(jù)分析平臺。Zhang 等［18］受迭代收縮閾值算法（ISTA）的啟發(fā)，提出了一種新的結(jié)構(gòu)化深度網(wǎng)絡(luò)（Iterative Shrinkage-Thresholding Algorithm Network，ISTA-net），能快速、準(zhǔn)確地實現(xiàn)數(shù)據(jù)的去噪和壓縮感知重建。

本文利用深度學(xué)習(xí)算法中的監(jiān)督數(shù)據(jù)驅(qū)動算法——迭代 啟 發(fā) 網(wǎng) 絡(luò)（Iterative scheme inspired network，IIN）壓制非平穩(wěn)隨機噪聲［19］，該網(wǎng)絡(luò)結(jié)構(gòu)簡單、緊湊。IIN 由交替方向乘子算法（Alternating Direction of Method of Multipliers，ADMM）的迭代過程推導(dǎo)而來，利用L1范數(shù)優(yōu)化變分模型。在訓(xùn)練階段，通過增加一個新的輔助變量，將目標(biāo)函數(shù)的極值轉(zhuǎn)化為增廣拉格朗日格式，使用L-BFGS算法判別、訓(xùn)練所有網(wǎng)絡(luò)參數(shù)，最終得到最優(yōu)去噪模型參數(shù)。在實驗階段，首先利用模型數(shù)據(jù)和實際數(shù)據(jù)建立訓(xùn)練數(shù)據(jù)庫，使用訓(xùn)練階段得到的去噪模型去噪。實驗結(jié)果表明，與現(xiàn)有的壓制非平穩(wěn)隨機噪聲方法相比，文中方法具有明顯的優(yōu)勢。

1 基本理論與技術(shù)流程

1.1 迭代收縮—閾值啟發(fā)網(wǎng)絡(luò)算法

噪聲消除的目的是從采集到的含噪數(shù)據(jù)y中恢復(fù)無噪數(shù)據(jù)。傳統(tǒng)壓縮感知（Compressive Sensing，CS）方法通過求解以下優(yōu)化問題恢復(fù)無噪數(shù)據(jù)，其目標(biāo)函數(shù)為［20］

式中：x為無噪數(shù)據(jù)的估計值；Φ 為變化矩陣，將x變換到不同域；Ψ 為非線性變換函數(shù)；λ 為正則化參數(shù)。

ISTA-net網(wǎng)絡(luò)算法結(jié)合CS算法，網(wǎng)絡(luò)由固定的層數(shù)組成，每層中一次迭代采用傳統(tǒng)ISTA 算法，該網(wǎng)絡(luò)結(jié)合了ISTA 和CS的優(yōu)勢［21］。該網(wǎng)絡(luò)通過迭代兩個更新步驟

估計去噪 數(shù) 據(jù)［22］。式 中：k 為 迭 代 次 數(shù)；s（k）為 步長，根據(jù)之前的步長參數(shù)由學(xué)習(xí)確定；T 為軟閾值算子的閾值；F（x）為參數(shù)可學(xué)習(xí)的非線性變換函數(shù)。

1.2 IIN 算法

非平穩(wěn)隨機噪聲主要表現(xiàn)為脈沖噪聲，以不規(guī)則脈沖或噪聲“尖峰”的形式隨機出現(xiàn)在地震數(shù)據(jù)中，持續(xù)時間短，具有譜密度較寬和振幅相對較高的特點［23］。由于L2范數(shù)對局部異常值的檢測不穩(wěn)定，在壓制非平穩(wěn)隨機噪聲時會限制算法的處理效果。因此本文采用的IIN 算法在迭代收縮—閾值啟發(fā)網(wǎng)絡(luò)算法的基礎(chǔ)上，采用L1范數(shù)優(yōu)化噪聲壓制方程，增強對局部強能量噪聲的敏感度［24］，結(jié)合ADMM 網(wǎng)絡(luò)（ADMM-net）［25］，通過優(yōu)化問題

得到去噪數(shù)據(jù)。式中l(wèi)∈［1，2，…，L］為濾波器序號，L 為單層濾波器的總數(shù)。濾波器D 的參數(shù)由網(wǎng)絡(luò)進行訓(xùn)練。

引入變量t＝x－y，z＝Dx，則式（3）的增廣拉格朗日函數(shù)為［26］

式中：ε、ρ為權(quán)重系數(shù)；α、ω 為拉格朗日乘數(shù)矩陣。根據(jù)ADMM 算法，引入變量p＝ω／ε，βl＝αl／ρl，則式（4）分解為6個問題

式中：n為層數(shù)；S（·）為非線性收縮函數(shù)；I為各元素為1的矩陣；Hnl為第n 層中第l 個濾波器矩陣；τ、η為可學(xué)習(xí)的權(quán)重參數(shù)，控制學(xué)習(xí)速率，在網(wǎng)絡(luò)訓(xùn)練中逐層修正。

1.3 技術(shù)流程

在IIN 網(wǎng)絡(luò)算法中，每層有6個計算節(jié)點（對應(yīng)式（5）），具體步驟如下。

（1）利用式（5a）重構(gòu)輸入的地震數(shù)據(jù)。若該層是第一層（n＝1）時，式（5a）變?yōu)?/p>

（2）對重構(gòu)數(shù)據(jù)x（n）進行濾波處理。

（3）根據(jù)式（5c）對濾波數(shù)據(jù)進行非線性變換，其中S［·］為分段函數(shù)，其形式為

由于分段線性函數(shù)擬合性程度高，可以根據(jù)數(shù)據(jù)的特點，學(xué)習(xí)更加合適的非線性變換函數(shù)。其中

（5）對重構(gòu)數(shù)據(jù)按步驟（3）進行非線性變換，如果n＝1，即為第1層時，令p（n－1）＝0。

（6）根據(jù)式（5f）更新乘數(shù)p（n）。

圖1為第k層正向傳播流程。

在網(wǎng)絡(luò)訓(xùn)練過程中，通過計算估計數(shù)據(jù)與輸入數(shù)據(jù)之間的誤差作為損失函數(shù)

式中Θ 表示上述流程中各個參數(shù)的集合。在反向傳播過程中計算損失函數(shù)E（Θ）對參數(shù)集Θ 的梯度，利用L-BFGS算法，通過最小化式（8）的取值，優(yōu)化去噪數(shù)據(jù)（y，Θ）的參數(shù)集Θ。通過不斷循環(huán)反向傳播過程，使E（Θ）的值不再下降、保持穩(wěn)定為止，以保證在去噪的同時不壓制有效信號。

圖2為第k層反向傳播流程。通過反向傳播計算得到去噪網(wǎng)絡(luò)參數(shù)模型，利用該模型對輸入地震數(shù)據(jù)進行正向傳播便可得到去噪后的地震數(shù)據(jù)。

圖1 第k層正向傳播流程

圖2 第k層反向傳播流程

2 理論模型試算

為了驗證上述去噪算法的可行性和有效性，分別對模型數(shù)據(jù)和實際地震數(shù)據(jù)進行去噪實驗。通過比較不同層數(shù)的網(wǎng)絡(luò)處理效果設(shè)置網(wǎng)絡(luò)總層數(shù)，過多的層容易導(dǎo)致過擬合，過少的層數(shù)不能有效去除噪聲。綜合去噪效果和計算時間設(shè)置單層濾波器總數(shù)目。根據(jù)經(jīng)驗設(shè)置NC值，關(guān)系到分段線性函數(shù)的擬合程度。合成地震數(shù)據(jù)添加了非平穩(wěn)隨機噪聲，生成含噪數(shù)據(jù)作為輸入數(shù)據(jù)，無噪數(shù)據(jù)作為輸出數(shù)據(jù)；單炮疊前地震數(shù)據(jù)作為實際數(shù)據(jù)，存在大量的非平穩(wěn)隨機噪聲，采用空變中值濾波和時變中值濾波去噪得到去噪后數(shù)據(jù)，再選擇去噪效果理想的數(shù)據(jù)作為數(shù)據(jù)庫數(shù)據(jù)。通過水平翻轉(zhuǎn)和90°旋轉(zhuǎn)將訓(xùn)練數(shù)據(jù)集擴充，以增加數(shù)據(jù)的多樣性。然后對上述數(shù)據(jù)訓(xùn)練得到去噪模型。圖3 為訓(xùn)練庫部分?jǐn)?shù)據(jù)。

圖4為損失函數(shù)曲線與峰值信噪比（PSNR）曲線。由圖可見：①在第6次迭代之前，損失函數(shù)值呈明顯下降趨勢，當(dāng)?shù)螖?shù)大于6后，損失函數(shù)值收斂于小值，說明去噪誤差達(dá)到最小值，并且其變化趨于平穩(wěn)，證明訓(xùn)練模型去噪結(jié)果穩(wěn)定（圖4a）。②前幾次迭代生成的去噪模型的去噪能力不強，經(jīng)過幾次迭代之后，去噪能力趨于穩(wěn)定，尤其在迭代次數(shù)大于7 之后，去噪數(shù)據(jù)的PSNR 較高且穩(wěn)定（圖4b）。由于采用簡單、緊湊的網(wǎng)絡(luò)結(jié)構(gòu)，因此曲線收斂速度較快，表明文中方法用相對較小的數(shù)據(jù)集和較少的訓(xùn)練時間得到了穩(wěn)定的去噪模型。

根據(jù)模型中有效信號的能量在合成數(shù)據(jù)（圖5a）加入非平穩(wěn)隨機噪聲得到含噪數(shù)據(jù)（圖5b），再利用訓(xùn)練好的網(wǎng)絡(luò)模型去噪得到去噪結(jié)果（圖5c）?？梢?，經(jīng)過去噪后噪聲被壓制，去除的噪聲（圖5d）中基本不存在有效信號，提高了信噪比。為了對比對疊后模型的去噪效果，分別利用本文方法與傳統(tǒng)中值濾波方法對含噪模型（圖6a）去噪，結(jié)果表明：本文方法的去噪結(jié)果（圖6b）沒有殘留噪聲，去除的噪聲（圖6c）中不存在有效信號；中值濾波方法的閾值設(shè)置影響了去噪結(jié)果（圖6d），部分有效信號被壓制，如果閾值設(shè)置太大，則會殘留噪聲，很難選取合適的閾值有效地壓制噪聲。

圖3 訓(xùn)練庫部分?jǐn)?shù)據(jù)

圖4 損失函數(shù)曲線（a）與峰值信噪比（PSNR）曲線（b）

圖5 數(shù)據(jù)模型的去噪效果

圖6 對疊后模型的去噪效果

圖7為圖6第10道數(shù)據(jù)去噪前、后對比。由圖可見：由于本文方法噪聲壓制方程的保真項采用L1范數(shù)，可更敏感地檢測地震數(shù)據(jù)中的非平穩(wěn)隨機噪聲并進行壓制（圖7a）；中值濾波去噪結(jié)果中有效信號也被壓制，去噪效果不理想（圖7b）。

圖7 圖6第10道數(shù)據(jù)去噪前、后對比

3 實際地震數(shù)據(jù)的應(yīng)用

分別采用時頻峰值濾波方法和本文方法對實際地震數(shù)據(jù)去噪并分析去噪結(jié)果。圖8為F 區(qū)實際地震數(shù)據(jù)去噪結(jié)果。由圖可見：86道實際地震記錄中有幾道含有較強的非平穩(wěn)隨機噪聲，干擾有效信號并影響弱有效信號的拾?。▓D8a）；時頻峰值濾波對圖8a的去噪結(jié)果中殘留噪聲（圖8b）；經(jīng)本文方法對圖8a去噪后有效壓制了噪聲，沒有破壞有效信號且很好地保護了弱有效信號（圖8d）。

單點高密度地震資料信號頻帶寬、波場信息豐富，反映了地下真實的反射信息，但同時也存在大量的噪聲，因此出現(xiàn)信噪比較低、有效信號與噪聲混雜的問題，不易發(fā)揮高密度資料的潛在優(yōu)勢。為了檢驗本文方法對單點高密度地震數(shù)據(jù)的去噪效果，在訓(xùn)練庫中增加相應(yīng)的單點高密度數(shù)據(jù)，并選取G 區(qū)的單點高密度單炮道集進行去噪實驗。圖9 為G區(qū)單點高密度地震數(shù)據(jù)去噪結(jié)果。由圖可見：經(jīng)本文方法對單點高密度地震數(shù)據(jù)（圖9a）去噪后，在壓制非平穩(wěn)隨機噪聲的同時，保護了有效信號（圖9d箭頭處）；時頻峰值濾波方法對圖9a的去噪結(jié)果中殘留噪聲（圖9b箭頭處）。圖10 為圖9 的局部放大。由圖可見：時頻峰值濾波處理結(jié)果中殘留噪聲（圖10b）；本文方法去噪結(jié)果的同相軸清晰、連續(xù)，有效壓制了噪聲（圖10c）。

圖8 F區(qū)實際地震數(shù)據(jù)去噪結(jié)果

圖9 G 區(qū)單點高密度地震數(shù)據(jù)去噪結(jié)果

理論模型試算和實際地震數(shù)據(jù)去噪結(jié)果表明，本文方法有效減小了信號畸變，較好地壓制了各種非平穩(wěn)隨機干擾，且對實際數(shù)據(jù)具有較強適應(yīng)性。

圖10 圖9的局部放大

4 結(jié)束語

針對地震數(shù)據(jù)中的非平穩(wěn)隨機噪聲，本文提出了采用基于迭代啟發(fā)網(wǎng)絡(luò)算法的去噪方法，采用深度學(xué)習(xí)算法，利用L1范數(shù)對目標(biāo)函數(shù)進行優(yōu)化，可更敏感地檢測并壓制地震數(shù)據(jù)中的非平穩(wěn)隨機噪聲，同時較好地保護有效信號。理論模型及實際資料的去噪結(jié)果表明：

（1）利用迭代啟發(fā)網(wǎng)絡(luò)（IIN）算法，通過構(gòu)建包含不同形態(tài)有效信號同相軸的數(shù)據(jù)，并結(jié)合實際地震數(shù)據(jù)建立訓(xùn)練數(shù)據(jù)庫，由訓(xùn)練得到的去噪模型根據(jù)有效信號的特征，在去噪的同時可保留同相軸的形狀特征；采用的迭代網(wǎng)絡(luò)簡單、緊湊，加快了網(wǎng)絡(luò)的收斂速度，用相對較小的數(shù)據(jù)集和較少的訓(xùn)練時間快速訓(xùn)練去噪模型，去噪效果較好。

（2）本文方法具有較強的適應(yīng)性，有效壓制了常規(guī)地震數(shù)據(jù)中的非平穩(wěn)隨機噪聲。對于單點高密度地震數(shù)據(jù)，較傳統(tǒng)的時頻峰值濾波算法，可有效減小信號畸變，較好地壓制各種非平穩(wěn)隨機干擾。