龍生海，鄒時(shí)林，王勝平

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

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

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

1 理論方法

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

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

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

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

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

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

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

2 處理過程

2.1 質(zhì)量控制

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

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

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

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

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

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

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

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

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

3 試驗(yàn)分析

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

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

4 結(jié) 論

本文通過實(shí)驗(yàn)分析得出如下結(jié)論：

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

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

3）地磁日變數(shù)據(jù)質(zhì)量控制是地磁日變數(shù)據(jù)處理的一個(gè)重要環(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.