本發(fā)明涉及一種地震數(shù)據(jù)處理去噪的方法，尤其是工頻干擾嚴(yán)重的地震記錄去噪方法，利用地震記錄各道的平均能量值識(shí)別干擾道以及利用主成分分析技術(shù)壓制工頻干擾的方法。

背景技術(shù)：
：

地震勘探資料采集過(guò)程中往往受到工頻干擾的影響，采集到的地震記錄中?；煊泄ゎl噪聲，當(dāng)工頻干擾很強(qiáng)，且貫穿于一道地震數(shù)據(jù)的始終時(shí)，該地震道常表現(xiàn)為壞道，引起地震資料品質(zhì)下降。通過(guò)對(duì)大量壞道數(shù)據(jù)的分析，很多地震道并不是源于地面附近的強(qiáng)干擾源，往往來(lái)自地震檢波器本身的問(wèn)題，與地震記錄的時(shí)間深度無(wú)關(guān)，當(dāng)前對(duì)壞道的處理是剔除壞道，這種處理降低了地震數(shù)據(jù)質(zhì)量，而本專利提供的技術(shù)，能夠有效識(shí)別由于工頻干擾產(chǎn)生的壞道，這為進(jìn)一步壓制工頻干擾提供了可能。針對(duì)工頻干擾當(dāng)前還沒(méi)有非常有效的解決方法。目前主要的壓制方法有陷波法、工頻回歸相減抑制法，自適應(yīng)濾波法三類。文章“地震資料中工頻干擾的自動(dòng)識(shí)別與壓制”提出陷波法濾除工業(yè)干擾，但該方法同時(shí)也會(huì)濾除50Hz頻率的地震信號(hào)；文章“強(qiáng)工頻干擾波的提取與消除方法”提出時(shí)間域內(nèi)運(yùn)用共軛梯度算法提取單頻干擾，但該方忽略了實(shí)際地震資料中工頻頻率的可變性，研究表明在不同區(qū)域，不同時(shí)間，工頻干擾頻率是變化的；文章“壓制地震資料中工業(yè)電干擾的余弦逼近法的改進(jìn)及應(yīng)用”提出基于余弦逼近法濾除工業(yè)干擾，但該方法無(wú)法實(shí)現(xiàn)自動(dòng)識(shí)別干擾道；文章“在檢波點(diǎn)域分離50Hz工業(yè)干擾”提出炮集轉(zhuǎn)換濾除工業(yè)干擾，但該方法效率低并無(wú)法自動(dòng)識(shí)別干擾道；文章“基于維納濾波的50Hz工業(yè)干擾去噪方法及應(yīng)用”提出維納濾波器來(lái)調(diào)整參考道的振幅和相位的匹配來(lái)預(yù)測(cè)單頻噪聲，但該方法不適應(yīng)于平穩(wěn)性差的工頻干擾；文章“基于獨(dú)立分量分析的工頻干擾消除技術(shù)”提出從混合信號(hào)中分離出相互獨(dú)立的各個(gè)信號(hào)分量達(dá)到消除工業(yè)噪聲的目的，但該方法依靠經(jīng)驗(yàn)性，缺乏嚴(yán)格的理論依據(jù)；

CN104570118A公開的“一種基于雙因素的自動(dòng)識(shí)別與去除工業(yè)干擾的方法”提出正余弦加權(quán)逼近法處理工業(yè)干擾，但該方法忽略了工頻的頻率，相位和幅值的不固定性；

CN103630935A公開的“去除地震數(shù)據(jù)中交流電干擾信號(hào)的方法”提出頻率域提取工頻頻率的方法，但該方法沒(méi)有考慮工頻諧波的存在；

CN101907726A公開的“一種自動(dòng)識(shí)別和消除地震勘探工業(yè)電干擾的方法”提出基于自相關(guān)理論確定工頻干擾，但該方法自動(dòng)識(shí)別工頻噪聲時(shí)依賴于初至?xí)r間之前地震資料質(zhì)量。

技術(shù)實(shí)現(xiàn)要素：
：

本發(fā)明的目的就在于針對(duì)上述現(xiàn)有技術(shù)的不足，提供一種基于主成分分析的地震資料工頻干擾自動(dòng)識(shí)別與壓制方法。

本發(fā)明的主要思想是：地震采集過(guò)程中地震資料往往會(huì)受到工頻干擾的影響，從而造成地震記錄中的大量壞道或強(qiáng)干擾道，嚴(yán)重影響地震勘探數(shù)據(jù)的質(zhì)量，本發(fā)明是通過(guò)自動(dòng)識(shí)別受工頻干擾的地震道，然后通過(guò)主成分分析技術(shù)分離出干擾，再對(duì)重構(gòu)后的數(shù)據(jù)進(jìn)行保幅，實(shí)現(xiàn)了地震資料工頻干擾的自動(dòng)識(shí)別與壓制。

本發(fā)明是通過(guò)以下技術(shù)方案實(shí)現(xiàn)的：

基于主成分分析的地震資料工頻干擾自動(dòng)識(shí)別與壓制方法，包括以下步驟：

a、針對(duì)待處理的原始地震記錄，選取記錄時(shí)間位于中下段的時(shí)窗T_W，計(jì)算該地震記錄在T_W內(nèi)的各道平均能量值，如公式

其中E_i為時(shí)窗內(nèi)地震記錄第i道平均能量值，x_i(t_j)為時(shí)窗內(nèi)地震記錄第i道的第j個(gè)采樣點(diǎn)值，i＝1,2…n，n為地震記錄總道數(shù)，j＝1,2…k，k為時(shí)窗內(nèi)地震記錄的采樣點(diǎn)數(shù)，時(shí)窗T_W的長(zhǎng)度建議為總記錄長(zhǎng)度的到時(shí)間內(nèi)；

b、定義第l道特征因子根據(jù)λ_l的取值識(shí)別干擾道，E_l為地震記錄第l道在T_W內(nèi)的平均能量值，q為T_W內(nèi)具有最小平均能量值的地震道號(hào)，其中l(wèi)＝1,2,…n，識(shí)別干擾道具體步驟如下：

b1、若對(duì)所有道均存在λ_l＜5，表明該地震記錄無(wú)顯著工頻干擾，無(wú)需后續(xù)去噪過(guò)程，否則執(zhí)行步驟b2；

b2、若λ_l≥5，表明地震記錄有顯著工頻干擾，記第l道為干擾道，記為U_l；

c、提取與U_l相鄰且包括U_l的連續(xù)5道記錄組成子地震干擾道集，記為X，X中每一道為一列，求出X每列的均值，并用X中的每個(gè)元素減去該列的均值，得到X處理后的子地震道集，記為X^*；

d、計(jì)算X^*的協(xié)方差矩陣Γ，如公式

其中b為X^*的列數(shù)，X^*T為X^*的轉(zhuǎn)置矩陣，“·”表示矩陣乘法；

e、利用奇異值分解法計(jì)算出協(xié)方差矩陣Γ的特征值矩陣Λ和特征向量矩陣R，則存在公式

Γ＝R·Λ·R^T (3)

其中Λ為特征值從大到小排列的對(duì)角矩陣，R每一列為對(duì)應(yīng)特征值的特征向量，R^T為轉(zhuǎn)置矩陣，滿足R^T·R＝R·R^T＝E,E為單位矩陣；

f、X經(jīng)R^T線性映射，得到主成分矩陣，如公式

Φ＝R^T·X (4)

其中Φ為主成分矩陣；

g、將Φ第一行置0，得到Φ'，令

X₁^*＝R·Φ' (5)

其中X₁^*為重構(gòu)矩陣；

h、提取X₁^*中對(duì)應(yīng)干擾道的數(shù)據(jù)，記為U_l'，對(duì)U_l'進(jìn)行保幅處理，如公式

其中U_l^*為保幅處理后干擾道信號(hào)，A₁為去噪后干擾道信號(hào)振幅絕對(duì)值的平均值，A₂為與U_l相鄰兩道數(shù)據(jù)的振幅絕對(duì)值平均值之和的平均。

有益效果：經(jīng)試驗(yàn)，本發(fā)明公開的一種基于主成分分析的地震資料工頻干擾的自動(dòng)識(shí)別與壓制方法能夠?qū)崿F(xiàn)在地震勘探中壓制工頻干擾源帶來(lái)的噪聲,并且能夠自動(dòng)識(shí)別干擾道，該算法能夠修復(fù)工頻干擾引起的壞道，有效提高了地震資料質(zhì)量，降低了采集成本。

附圖說(shuō)明：

圖1實(shí)際地震記錄

圖2壓制工頻干擾后的地震記錄

具體實(shí)施方式：

下面結(jié)合附圖和實(shí)施例對(duì)本發(fā)明做進(jìn)一步的詳細(xì)說(shuō)明：

基于主成分分析的地震資料工頻干擾自動(dòng)識(shí)別與壓制方法，包括以下步驟：

a、采用實(shí)際單炮記錄，本例T_W＝700ms～800ms，計(jì)算T_W內(nèi)的E_i，如公式

其中E_i為地震記錄第i道平均能量值，x_i(t_j)為地震記錄第i道的第j個(gè)采樣點(diǎn)元素，本例中i＝1,2…n，n＝48，j＝1,2…k，k＝100；

b、定義第l道特征因子根據(jù)λ_l的取值識(shí)別干擾道，E_l為地震記錄第l道在T_W內(nèi)的平均能量值，本例中q為T_W內(nèi)的第一道，E₁＝0.0042，其中l(wèi)＝1,2,…n，識(shí)別干擾道具體步驟如下：

b1、若對(duì)所有道均存在λ_l＜5，表明該地震記錄無(wú)顯著工頻干擾，無(wú)需后續(xù)去噪過(guò)程，本例中滿足λ_l＜5的道為第1～7道，第9～42道，第44～48道，否則執(zhí)行步驟b2；

b2、若λ_l≥5，表明地震記錄有顯著工頻干擾，記第l道為干擾道，記為U_l，本例中滿足λ_l≥5的道為第8和第43道，記為U₈和U₄₃，E₈＝0.078，E₄₃＝0.23，λ₈＝18.57，λ₄₃＝54.76；

c、以U₈為例，提取第6，7，8，9，10道組成子地震干擾道集，記為X，X中每一道為一列，求出X每列的均值，并用X中的每個(gè)元素減去該列的均值，得到X處理后的子地震道集，記為X^*；

d、計(jì)算X^*的協(xié)方差矩陣Γ，如公式

其中b為X^*的列數(shù)，X^*T為X^*的轉(zhuǎn)置矩陣，“·”表示矩陣乘法，式中b＝5；

e、利用奇異值分解法計(jì)算出協(xié)方差矩陣Γ的特征值矩陣Λ和特征向量矩陣R，則存在公式

Γ＝R·Λ·R^T (3)

其中Λ為特征值從大到小排列的對(duì)角矩陣，R每一列為對(duì)應(yīng)特征值的特征向量，R^T為轉(zhuǎn)置矩陣，滿足R^T·R＝R·R^T＝E,E為單位矩陣；

f、X經(jīng)R^T線性映射，得到主成分矩陣，如公式

Φ＝R^T·X (4)

其中Φ為主成分矩陣；

g、將Φ第一行置0，得到Φ'，令

X₁^*＝R·Φ' (5)

其中X₁^*為重構(gòu)矩陣；

h、提取X₁^*中對(duì)應(yīng)干擾道的數(shù)據(jù)，記為U_l'，對(duì)U_l'進(jìn)行保幅處理，如公式

其中U_l^*為保幅處理后干擾道信號(hào)，A₁為去噪后干擾道信號(hào)振幅絕對(duì)值的平均值，A₂為與U_l相鄰兩道數(shù)據(jù)的振幅絕對(duì)值平均值之和的平均,本例中A₂＝0.0689,A₁＝0.0057。