基于數(shù)字信號處理技術(shù)的地磁日變數(shù)據(jù)處理方法研究

龍生海，鄒時林，王勝平

（東華理工大學(xué) 測繪工程學(xué)院，江西 南昌 330013）

地磁是地球的一個重要的地球物理特征。地磁場不是恒定不變的，而是隨時間緩慢變化的。研究這種變化的時空分布規(guī)律對于了解地球內(nèi)部物質(zhì)的性質(zhì)和運動具有重要意義。近年來地磁在水下無源導(dǎo)航以及地震預(yù)報方面的研究也取得了突破性的進(jìn)展［1－3］，已經(jīng)成為地震監(jiān)測的重要手段。地磁臺站是最基本的地磁監(jiān)測手段。我國的地磁臺網(wǎng)發(fā)展迅速，1Hz采樣率的數(shù)據(jù)已經(jīng)成為現(xiàn)實［4－5］。如何對如此高密度的地磁時間序列數(shù)據(jù)進(jìn)行快速準(zhǔn)確的分析計算成為研究的重點，包含準(zhǔn)確變化特征的地磁日變數(shù)據(jù)是地震科學(xué)研究的根本保障。

正是基于此背景，本文提出基于數(shù)字信號處理技術(shù)的地磁日變處理技術(shù)，并對整個數(shù)據(jù)處理過程進(jìn)行深入分析。

1 理論方法

傅里葉變化是數(shù)字信號處理技術(shù)中最經(jīng)典的算法，對于離散的地磁觀測序列，離散傅里葉變換實現(xiàn)了頻域離散化，在數(shù)字信號處理技術(shù)中起著及其重要的作用，它可以直接用來分析信號的頻譜、計算濾波器頻率響應(yīng)，以及實現(xiàn)信號通過線性系統(tǒng)的卷積運算等。該方法既可以大大減少運算時間，在選取合適的截止頻率的情況下又可以實現(xiàn)有用信息的最大保留。

傅里葉變化可以看做是傅里葉級數(shù)的連續(xù)形式。傅里葉級數(shù)把定義在［－π，π］上的信號分解為頻率為整數(shù)倍關(guān)系的諧波分量組合，相應(yīng)的傅里葉變換將一個無限時寬的信號分解為頻率為λ的一系列頻率分量，其中λ可以是任意實數(shù)（甚至是復(fù)數(shù)）。如果x（t）是定義在整個實軸上的實值或復(fù)值函數(shù)，則其傅里葉變換可由式（1）給出。

若對任意參數(shù)f，上述積分都存在，則式（1）確定了一個函數(shù)X（f），稱為X（f）的傅里葉變換，如果已知X（f），則利用如下的傅里葉逆變換，還可以復(fù)原x（t）。

若x（t）和X（f）同時滿足式（1）和式（2），則稱它們是一個傅里葉變換對，記為x（t）?X（f），通常X（f）是一個復(fù)函數(shù)，因此可以寫成以下2部分：

式中：R（f）和I（f）分別是X（f）的實部和虛部。將式（3）表示成指數(shù)形式常常是方便的。

快速傅氏變換（FFT），是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實等特性，對離散傅里葉變換的算法進(jìn)行改進(jìn)獲得的。它對傅氏變換的理論并沒有新的發(fā)現(xiàn)，但是對于在計算機(jī)系統(tǒng)或者說數(shù)字系統(tǒng)中應(yīng)用離散傅里葉變換是很大進(jìn)步。

該算法可以實現(xiàn)離散信號從時間域到頻率域的轉(zhuǎn)換，在得到信號的頻率之后，可以根據(jù)實際需要刪除或保留自己需要的頻率段，將所保留的數(shù)據(jù)再進(jìn)行傅里葉逆變換。最終實現(xiàn)對離散地磁數(shù)據(jù)的平滑處理。

2 處理過程

2.1 質(zhì)量控制

各類地磁變化具有不同的相態(tài)和時空分布特征，并且常常疊加在一起，起伏變化有時比較平緩，有時比較劇烈。有時就形成了一種復(fù)雜的干擾。同時地磁數(shù)據(jù)對外界的干擾十分的敏感。如果測站周圍有任何磁性體都會對原始地磁數(shù)據(jù)序列產(chǎn)生很大的干擾，從而引起地磁數(shù)據(jù)的異常突變。因此在數(shù)據(jù)處理之前必須首先對原始觀測數(shù)據(jù)進(jìn)行相應(yīng)的質(zhì)量控制［6－8］。

本文所述方法采用抗差卡爾曼濾波與人工交互式方法相結(jié)合的辦法實現(xiàn)了對地磁日變數(shù)據(jù)的質(zhì)量控制。

2.2 基于FFT的地磁數(shù)據(jù)平滑方法研究

目前地磁數(shù)據(jù)的采樣率已經(jīng)可以到達(dá)1Hz的頻率。如此高采樣率的地磁數(shù)據(jù)雖然可以更好反應(yīng)地磁變化，但是由于各種原因引起的小的噪聲，使數(shù)據(jù)呈現(xiàn)一種鋸齒狀的變化特征。具有不同時變特征的數(shù)據(jù)有不同的用處，因此研究可以方便為不同用戶快速提供不同特征數(shù)據(jù)的方法顯得至關(guān)重要。本文基于第一部分提出的FFT數(shù)字信號處理技術(shù)，展開對地磁日變數(shù)據(jù)的平滑處理。

地磁日變數(shù)據(jù)的平滑處理一直是地磁數(shù)據(jù)處理中一個重要的環(huán)節(jié)。利用FFT算法可以將地磁日變觀測數(shù)據(jù)從時間域轉(zhuǎn)換到頻率域，再根據(jù)相應(yīng)的截止提取地磁日變序列中的長周期信號，最后將提取出的信號通過FFT逆變換轉(zhuǎn)換到時間域。該方法也可以實現(xiàn)地磁數(shù)據(jù)的平滑。但是其中截止周期的確定是至關(guān)重要的一個環(huán)節(jié)。

本文通過研究表明，選擇不同的截止周期可以為用戶提供不同特征的數(shù)據(jù)，可以滿足不同的需求。圖1～4為采樣間隔為1s的原始數(shù)據(jù)分別在周期為100s、300s以及500s所對應(yīng)的結(jié)果圖，圖中細(xì)線代表原始數(shù)據(jù)，粗線代表計算結(jié)果。

由圖1～4可以看出，采用FFT算法，選取不同的周期提取數(shù)據(jù)的特征各不相同。圖1選用周期為100s，提取數(shù)據(jù)很好地反映了數(shù)據(jù)細(xì)微的變化，同時也將可能是噪聲的誤信號帶入其中。該數(shù)據(jù)可以用于物探、地震、磁暴預(yù)測等對數(shù)據(jù)細(xì)節(jié)要求很高的鄰域，并不是數(shù)據(jù)平滑的最佳截止周期。圖2選用周期為300s，提取數(shù)據(jù)很好地實現(xiàn)了數(shù)據(jù)的平滑，與圖3相比在細(xì)節(jié)部分對數(shù)據(jù)的反映更加真實。圖4選用周期為1 000s，該數(shù)據(jù)反映了地磁變化趨勢，通過選用合適的濾波門限可以用于地磁日變數(shù)據(jù)的消噪，但是其也丟失了過多的地磁數(shù)據(jù)時變特征，不能很好反映整個地磁序列的細(xì)部特征。因此，建議確定300s為地磁日變數(shù)據(jù)平滑的截止周期。同時通過上述研究，分別得出了幾種數(shù)據(jù)截止周期（見表1）。

表1 1Hz采樣率原始數(shù)據(jù)不同的截止周期選擇 s

圖1 周期為100s時的計算結(jié)果

3 試驗分析

通常地磁數(shù)據(jù)的平滑多選用多項式擬合逼近的辦法來實現(xiàn)。如果曲線的階數(shù)選擇的過小，擬合效果不好；如果曲線的階數(shù)過高，雖然數(shù)據(jù)點上看到效果好，數(shù)據(jù)點之間會出現(xiàn)有數(shù)據(jù)振蕩的問題且計算速度很慢效率較低，同時可調(diào)節(jié)性過小。為驗證本文方法的準(zhǔn)確性與可靠性以及與傳統(tǒng)方法相比較具有的優(yōu)缺點。本文采用中國西部某地磁日變站2008年某一天的地磁日變原始觀測數(shù)據(jù)作為試驗數(shù)據(jù)。基于文中所述FFT方法及截止周期的確定，選取300s的截止周期，將數(shù)據(jù)進(jìn)行平滑，并與5階多項式模型進(jìn)行比較（見圖5、圖6）。

通過圖5、圖6可以看出，二者在圖像上都實現(xiàn)了與原始數(shù)據(jù)的逼近，均可以對原始地磁序列的平滑處理，但是FFT算法可以更真實地反映數(shù)據(jù)的變化趨勢，與原始數(shù)據(jù)的趨勢吻合的更好。多項式方法的擬合結(jié)果雖然也實現(xiàn)了趨勢的反映，但是其對圖像的一些細(xì)部特征未能很好地表達(dá)，由圖5還可以發(fā)現(xiàn)傳統(tǒng)算法有一定的時延效應(yīng)。平滑曲線與原始趨勢相比有一定的滯后。

4 結(jié) 論

本文通過實驗分析得出如下結(jié)論：

1）基于FFT數(shù)字信號處理技術(shù)的地磁日變處理方法是正確和可靠的。用于1Hz采樣率的地磁日變數(shù)據(jù)平滑的截止周期建議選擇300s。同時利用該法還可以根據(jù)自己的需要對數(shù)據(jù)進(jìn)行細(xì)部特征提取及消噪。

2）傳統(tǒng)的多項式擬合逼近的辦法可以實現(xiàn)地磁日變數(shù)據(jù)的平滑，但是其反映數(shù)據(jù)細(xì)部特征的能力有限，同時還存在一定的滯后性，性能略低于FFT算法。

3）地磁日變數(shù)據(jù)質(zhì)量控制是地磁日變數(shù)據(jù)處理的一個重要環(huán)節(jié)，在地磁日變數(shù)據(jù)處理之前，必須進(jìn)行原始數(shù)據(jù)的質(zhì)量控制。

［1］Wang Tanwen.Research on Geomagnetic Field Model.Recent Developments In World Seismology［J］.Vol.4.

［2］AN Zhen－chang，XU Yuan－fang，WANG Yue－h(huán)ua.Derivation and analysis of the main geomagnetic field models in China for 1950－1980.Chinese Journal of Geophysics［J］.1991，34（5）.

［3］AN Zhen－chang.Review of Geomagnetic Surveys，Geomagnetic Charts and Geomagnetic Field Models in China.Chinese Journal of Geophysics［J］.45：189－196.

［4］Gao Jintian，An Zhenchang，Gu Zuowen，Han Wei，Zhan Zhijia，Yao Tongqi.Selections of the Geomagnetic Normal Field and Calculations of the Geomagnetic Anomalous Field.Chinese Journal of Geophysics.48（1）：51－62.

［5］Gu Zuowen，An Zhenchang，Gao Jintian，Zhan Zhijia，Yao Tongqi，Han Wei，Chen Bin.Computation and Analysis of the Geomagnetic Field Model In China and ITS Adjacent Area for 2003.Acta Seismologica Sinica［J］.28（2）：141－149.

［6］Behrooz Kamgar－Parsi，et al.Toward an auto－mated system for a correctly registered bathym etric ch art.IEEE Journal of Oceanic Engineering［J］.1989，14（4）：314－325.

［7］Zhang Z.Iterative point matching for registration of freeform curves and surfaces.International Journal of Computer Vision［J］.1994，13（2）：116－152.

［8］John J Leonard，Andrew Bennett.Autonimous Underwater Vehicle Navigation.MIT Marine Robotics Labotics Laboratory Technical Memorandum［J］.1998，（1）：1－11.