基于廣義S 變換時頻濾波的MT數據去噪

蔡劍華，肖永良

(1.湖南文理學院,洞庭湖生態經濟區建設與發展協同創新中心，湖南常德 415000；2.湖南財政經濟學院信息技術與管理學院，湖南長沙 410205 )

0 引言

在石油天然氣勘探中，大地電磁(MT)勘探已成為首選的綜合物探技術，特別是在地震勘探實施困難區，MT勘探在確定油氣儲層埋深、油氣勘探的有利遠景區和提供地球物理背景資料等方面比較有優勢，被廣泛應用于油氣田普查勘探和工程物理勘測(王家映，1997；Szarka，1998；何展翔，2000；賈進斗和孔繁恕，2001；宋小超等，2016；張文浩等，2019)。但隨著經濟的發展和人類活動的日益頻繁，以天然電磁場為源的MT探測方法，其采集的數據越發受到各類噪聲的干擾(Szarka，1998；Jing，et al.，2012；尚世貴，2014；蔡劍華和熊銳，2016；陶龍等，2019；張文浩等，2019)。尤其在礦集區，鉆采機械等大功率設備負荷的改變或開閉，會產生大幅度的噪聲，嚴重污染MT數據(嚴家斌等，2008；湯井田等，2012;王丹丹等，2020；陳炳錦等，2021；王基元等，2021)。大地電磁工作者一直以來就在尋求各種濾波技術來抑制此類噪聲，其中小波變換的方法受到大家的重視，基于小波變換的硬、軟閾值去噪方法被廣泛用于MT信號的處理(Trad and Travassos，1996；Li and Oh，2002；尚帥等,2015；李肅義等，2016)，但小波變換也存在一些缺陷，如小波基函數的確定、分解層數的選取都需要人為干預，且選擇不同，得到的處理結果也相差很大(楊學亭等，2014)。廣義S變換是近期發展起來的一種新的時頻分析技術，具有很高的時頻分析精度，并可無損重構(Stockwell et al.，1996)，在語音信號處理、機械故障診斷和地震勘探等領域應用廣泛(Pinnegar and Mansinha，2003；Simon et al.，2007；王小杰等，2012；王長江等，2013)。廣義S變換也引起了大地電磁工作者的關注，如景建恩提出了基于S變換時頻譜的大地電磁測深數據 Robust處理方法(景建恩等，2012)；陳海燕采用時頻比值法研究了廣義 S 變換在 MT 數據處理中的應用(陳海燕等，2012)。本文將廣義S變換時頻濾波方法應用于MT數據處理，研究受充放電干擾MT數據在S域的時頻分布特征，提出基于廣義S變換的時頻閾值濾波方法，在時頻域實現對MT數據干擾的去噪，以減少人為參與，并基于仿真和實測數據測試結果，討論其應用效果。

1 基于S變換時頻濾波的基本原理

1.1 廣義S變換及逆變換

隨機信號x(t)的S變換可定義為

(1)

(2)

GSTx(τ，f)也被稱為廣義S變換后信號x(t)的時頻分布(Simon et al.，2007)。

工程應用中，直接利用傅里葉反變換來實現廣義S變換的逆變換(Stockwell et al.，1996；Simon et al.，2007)。實現過程中先對GSTx(τ，f)的時間變量進行積分，得到傅里葉頻譜Y(f):

(3)

然后再進行傅立葉反變換，實現原始信號x(t)的高效重構:

(4)

1.2 時頻域閾值濾波

廣義S變換的時頻展開能力和可重構特性(Stockwell et al.，1996;王小杰等，2012；王長江等，2013)，為在時頻域對MT信號進行去噪提供了新途徑。首先采用廣義S變換對MT信號進行時頻展開，得到信號的時頻分布；進一步根據噪聲在時-頻面上的分布特征，確定時頻濾波因子對其進行閾值濾波，最后利用時頻反變換方法重構分離出消噪后的MT數據，流程如圖1所示。

圖1 基于廣義S變換的MT數據去噪流程圖

離散域時頻濾波的表達式為

y(n)=IGST[GST[x(n)]H(n,f)]

(5)

式中,GST和IGST分別為S變換和反變換，x(n)為去噪前的MT信號，y(n)為去噪后的MT信號，H(n,f)為時頻濾波算子，定義為(Donoho，1995；Duan and Xiong，2007；王長江等，2013；蔡劍華和肖曉，2015)。

(6)

2 仿真數據分析

用仿真信號對方法的有效性進行驗證。圖2a所示為蒙特卡洛方法產生的仿真MT信號。仿真數據的統計參數為表1中“原始信號”欄所示，其中能量為各頻率成分能量的總和，方差的定義如下(Donoho,1995)：

圖2 仿真信號的去噪

(7)

式中x(n)為信號序列，M為信號的均值，N表示數據的個數。鑒于充放電信號幅值與時間呈指數關系，模擬的“充放電干擾”是用指數函數產生的數據。圖2b為要加入的充放電干擾，在0.07s～0.31s和0.67s～0.91s兩個時間段模擬了兩段充放電干擾，圖2c為加噪后的數據。顯然，由于噪聲的存在，原始仿真信號幾乎被埋沒，完全看不到以前的隨機形態特征。加噪后數據的統計參數為表1中 “加噪后”欄所示，干擾使信號總能量和方差增加了近100倍，使得信號變得極不平穩。

表1 仿真信號去噪前后信號的表征參數

對仿真信號進行廣義S變換，得到去噪前、后的S域(圖3)。圖3a中，原始信號能量比較平穩，均勻地分布在各時頻點。加入噪聲后，圖3b所示時頻譜中信號的能量被提升，尤其低頻部的能量增加到原來的幾十倍，這與類充放電干擾頻率和頻譜范圍較小，主要影響低頻段數據的實際相吻合。采用本文提出的基于廣義S變換的時頻濾波方法，在時頻域內根據公式(5)和公式(6)對圖3b所示的時頻譜進行閾值去噪處理，得到去噪后的S域時頻譜如圖3c所示。可見，去噪后信號能量在整個時頻段內分布比較均勻。根據公式(3)和公式(4)，對去噪后的時頻譜實施反變換，重構去噪后的信號。去噪后的仿真信號如圖2d所示，統計參數為表1中“去噪后”欄所示。可見，信號的參數基本恢復到了加噪前的狀態,去噪后的數據與原數據高度逼近，噪聲得到了很好的抑制。

圖3 仿真信號去噪前后的S域時頻譜

3 實測數據應用

3.1 信號去噪

類充放電干擾是電磁測深作業中經常遇到的噪聲，特別是在作業區接近礦區或就在礦區內，此類電磁干擾成為研究區重要的干擾源，經常同時出現在電場和磁場信號中，往往相關性較好，振幅也很強，對低頻響應參數影響非常大。實測數據是來自某礦區測點的EH-4數據，MT數據的4個分量如圖4左圖所示，數據受噪聲的干擾很大。為清晰給出信號的時域特征，圖4右圖僅取了信號的前512個數據點，可見4個分量都受到了高能量的充放電類噪聲的干擾，在時域，MT信號的細節很難被觀察到。

圖4 實測大地電磁信號(左)和僅取前512點的數據(右)

現以磁場信號Hy分量為例，闡述干擾壓制的過程。先對Hy分量進行廣義S變換，得到S域的時頻譜圖，如圖5a所示。由于受到噪聲的干擾，低頻部分的能量被顯著抬升，使得在同一坐標尺下，其他頻帶的信號能量顯得很小。采用本文提出的S域時頻濾波方法對時頻譜進行閾值濾波處理，處理后的時頻譜如圖5b所示。可見，在圖5b中，低頻部分的能量得到了有效抑制，整個時頻譜內信號的能量分布勻稱。進一步計算Hy分量去噪前后的功率譜，如圖5c和圖5d所示。顯然，功率譜反映的情況與時頻譜表征的信號能量分布特征是一致的：去噪前，信號的能量較強，且在低頻段表現的尤為突出；去噪后，信號的總體能量減小，有效抑制了低頻段信號的能量。

圖5 Hy分量去噪前后的S域時頻譜和功率譜

去噪后實測大地電磁信號的時域波形如圖6所示。同時，表2給出了4個分量去噪前后的統計參數。對比去噪前后的信號波形和統計參數可知：數據均值明顯減小，方差也減小了2個數量級(如去噪前Hy分量的均值M=-236.65 nT，VAR=2.3×106；消噪后均值M=0.9199 nT，VAR=8.4×104)，信號變得平穩，充放電類噪聲得到有效抑制，大地電磁信號的形態特征得到顯現。

表2 MT信號去噪前后信號的表征參數

圖6 去噪后的實測大地電磁信號

3.2 響應參數的估計

分別用去噪前后的信號估算MT響應參數曲線ρxy，ρyx和Фxy，Фyx，并將處理的結果與小波軟閾值方法的處理結果進行對比(小波方法中的基函數采用“db5”小波，分解為6層)。圖7對比了去噪前后估算的MT響應參數曲線。進一步引入方差，用曲線平滑系數(Smoothing coefficient，SC)來定量評估去噪效果，平滑系數定義如下(Donoho,1995)：

(9)

式中x(i)為計算得到的各點響應參數，m(i)為小波方法和本文方法在i點計算所得響應參數的均值。曲線平滑系數越小表明曲線越平滑。表3給出了兩種方法計算的曲線ρxy的VAR和SC。由圖7和表3數據的對比可知，去噪前后計算的響應參數曲線整體趨勢是一致的，但去噪后MT資料的質量得到了明顯改善；相比小波軟閾值方法，ρxy電阻率曲線的方差VAR由6.365×104減小到5.761×104，平滑系數SC由33.21減小到26.97，去噪后響應參數曲線更為平穩，MT資料質量有更明顯的改善。可見本文提出的基于廣義S 變換時頻濾波的MT數據去噪方法去噪效果更好。

表3 兩種方法去噪后ρxy曲線參數對比

圖7 去噪前后大地電磁響應參數對比

4 結論

將廣義S變換應用于大地電磁信號中的充放電干擾抑制是有效的，通過廣義S變換，得到信號的時頻譜分布，在時頻域進行閾值去噪后，再通過其逆變換，重構去噪后的MT數據。在S域時頻譜內，抑制噪聲能量的同時，可有效保留有用信號的信息，去噪后大地電磁信號的特性得到了顯現，計算的響應曲線變得平穩，為后續工作提供了更為準確的資料。鑒于一條剖面有眾多測點，數據是海量的，進一步可選取一條剖面的數據進行本文方法的去噪處理，并對比視電阻率和相位擬斷面圖特征和反演結果。

[附中文參考文獻]

蔡劍華,肖曉.2015.基于小波自適應閾值去噪的MT信號處理方法[J].地球物理學進展,30(6):2433-2439.

蔡劍華,熊銳.2016.基于頻率切片小波變換的時頻分析與MT信號去噪[J].石油物探,55(6):904-912.

陳炳錦,鄭崔勇,袁波,余紅波,王偉.2021.磁法三維建模在陜西省龍王溝磁鐵礦勘探中的應用[J]，地質與勘探，57(1)：146-155.

陳海燕，景建恩，魏文博.2012.廣義S變換時頻域濾波在MT數據處理中的應用[J].現代地質，26(6)：1211-1216.

段錦升,熊曉燕.2007.小波變換在弱信號檢測中的應用研究[J].機械工程與自動化，(5)：96-97.

何展翔.2000.非地震勘探技術的進步與發展趨勢[J].石油地球物理勘探，35(3)：354-360.

賈進斗，孔繁恕.2001.不同油氣勘探階段的非地震勘探技術[J].石油地球物理勘探，36(4)：444-450，479.

景建恩,魏文博,陳海燕.2012.基于廣義S變換的大地電磁測深數據處理[J].地球物理學報,55(12):4015-4022.

李肅義,蔣善慶,王躍洋.2016.海洋可控源電磁數據中海水擾動噪聲的小波校正方法研究[J].石油物探,55(5):657-663.

尚帥,韓立國,胡瑋.2015.壓縮小波變換地震譜分解方法應用研究[J].石油物探,54(1):51-55.

尚世貴.2014.小波變換在磁異常分離中的應用―以安徽廬樅礦集區小包莊鐵礦床為例[J].物探與化探，38(6):1119-1123.

宋小超,李代榮,劉濤,段凱波,譚云基.2016.基于小波分析的高精度磁測在內蒙某多金屬礦普查中的應用[J].地球物理學進展，31(6):2649-2656.

湯井田,李晉,肖曉.2012.數學形態濾波與大地電磁噪聲壓制[J].地球物理學報,55(5):1784-1793.

陶龍,張莎莎,蘭學毅,湯正江,安明,嚴加永,王云云,郭冬,葉林,洪大軍.2019.1∶5萬重磁勘查在安徽宣城覆蓋區地質找礦中的應用探索[J].中國地質，46(4)：894-905.

王丹丹,趙松,張文浩,劉衛彬,陳中,劉偉偉.2020.松遼盆地外圍通化地區高精度重磁電特征及其構造格架[J].中國地質，47(4)：1056-1068.

王家映.1997.我國大地電磁測深研究新進展[J].地球物理學報，40 (增刊):206-216.

王長江,楊培杰,羅紅梅.2013.基于廣義S變換的時變分頻技術[J].石油物探,52(5)：489-494.

王小杰,印興耀,吳國忱.2012.基于S變換的吸收衰減技術在含氣儲層預測中的應用研究[J].石油物探,51(1)：37-42.

王基元，楊春海，符俊，張衛文，謝俊，李加斌.2021.滇西金廠河金銅多金屬礦集區成礦模式與綜合找礦模型[J].地質與勘探，57(2)：254-268.

嚴家斌，劉貴忠，柳建新.2008.小波變換在天然電磁場信號時間序列處理中的應用[J].地質與勘探，44(3)：75-78.

楊學亭,劉財,劉洋.2014.基于連續小波變換的時頻域地震波能量衰減補償[J].石油物探，53(5)：523-529

張文浩,王丹丹,李世臻,周新桂,張交東,劉衛彬,周錫明,王佩業.2019.重磁電勘探在黑龍江三江盆地石炭系-二疊系油氣地質調查中的應用[J].中國地質，46(1)：191-202.

張翔,王輝,汪遠志,劉建軍,陳江源.2016.高精度航磁資料在冶山鐵礦找礦潛力擴大中的應用[J].地質與勘探，52(6)：1138-1146.