地磁秒數據補齊策略探討

栗 寧 李從芬 賀北芳

1 湖北省地震局，武漢市洪山側路48號，4300712 中國地震局地震研究所，武漢市洪山側路40號，4300713 十堰市人民防空辦公室，湖北省十堰市江漢南路76號，442000

地磁數據在觀測中通常會受儀器故障、標定、雷擊等因素的影響，導致記錄缺數，如何盡可能不失真地重構補齊這些數據，并保證數據的連續性和可用性是一個亟需解決的問題[1]。目前國內對地磁秒數據補齊方法的研究較少[1-2]，相關的補齊算法有均值替代法、最近鄰域替代法、多重填補法及回歸替代法[3]等。這些算法針對個別缺數情況有較好的優勢，但無法應用于缺數在分鐘以上的地磁秒數據中。

針對這一問題，本文通過研究湖北省所有地磁臺站的秒數據文件發現，相鄰臺站數據的相關性高，且數據呈周期性，可依據相鄰臺站數據之間的相關性提取共同特性，再跟據數據自身個性補齊細節，從而盡可能不失真地補齊數據。

1 補數原理及實現

為了不失真地補齊數據，需要做到以下幾點：1)選取與所缺數據相關性高的數據，并確定移植的抽形圖；2)對缺失的數據依據抽形圖進行補數；3)對補好的抽形圖進行細節修補；4)對修補完成的數據進行檢測。

1.1 數據選取及傅里葉分解合成

對于一組包含缺失值的地磁數據，可在數據庫中找到一個與其最相似的對象。不同問題導致的缺數可選用不同標準進行相似性判定，最常見的是使用相關系數矩陣法，以選取與缺失數據關聯性最強的數據作為抽形圖。

抽形即擬合，選取何種擬合算法將直接決定修補數據的精度。本文選取傅里葉擬合法：首先將相似度高的地磁數據與所缺數據進行傅里葉分解，對分解后的正余弦波參數作對比分析，然后將調整后的波形合成便可得到所缺地磁數據的抽形圖。

在實際工作中遇到的大多數非正弦周期函數都滿足Dirichlet條件[4]，其函數f(t)可表示為：

(1)

式中，a0、ak、bk為傅里葉系數。由于地磁數據基本為秒采樣序列，設為S(t)：

m≤[n/2]

(2)

式(2)可寫成矩陣的形式AX=S，其中,

X=[a0,a1,b1,a2,b2,…,am,bm]T

(3)

S=[S(t1),S(t2),…,S(tn)]T

(4)

A=

(5)

利用最小二乘法求得X=(ATA)-1ATS，即可求出傅里葉級數擬合中各項的傅里葉系數，根據2組數據之間波形的對比便可得到缺數數據的擬合圖形。

1.2 基于ARMA預測模型的細節補齊算法

地磁數據是按時間排列的序列，可作為時間序列進行分析，建立自回歸滑動平均模型(ARMA模型)。ARMA模型能依據觀測值進行預測，因此將完好的數據作為歷史觀測值，缺失數據作為預測值，利用ARMA模型進行預測，完好數據又可作為檢驗值對預測值進行檢驗。

時間序列zt的ARMA模型具有以下形式：

(6)

式中，φi、θi為模型參數，p、q為模型的階次。

利用ARMA模型進行預測的前提是輸入的時間序列是平穩的，若為非平穩序列則需要通過差分方法轉化為平穩序列后再進行ARMA模型擬合[5]。ARMA模型建模的預測流程見圖1。

圖1 ARMA模型建模流程Fig.1 Path of ARMA model

2 實驗結果

選取2021-03-04恩施臺FGM01地磁H分量數據。該測點當天受人為進洞線路改造影響，世界時09:23～10:01原始數據丟失，00:00～10:33數據突跳，噪聲增加并產生臺階。為避免數據突跳對后期預測產生干擾，利用格拉布期準則對異常值進行剔除，并對實驗數據作歸一化處理，以消除量綱。

在數據庫中尋找與恩施臺FGM01地磁H分量預處理數據相似度高的數據。本文使用MATLAB中corrcoef函數篩選相關系數大于98%的數據，其中恩施臺GM4地磁的相關系數為0.983 9，應城臺GM4地磁、FGM01地磁及M15地磁的相關系數分別為0.993 6、0.994 3及0.994 2。因此，本文選擇2021-03-04應城臺FGM01地磁作為抽形圖框架，對其作歸一化處理后再進行傅里葉分解。為分辨擬合曲線，將同時段擬合曲線減去0.2 nT。從圖2可以看出，擬合曲線有截尾現象，但并不影響其整體形態。如果是尾部缺數，則可選取后1 d的部分數據作為驗證數據，并不影響補數精度。

圖2 應城臺FGM01地磁原始曲線、傅里葉擬合曲線及分解波形Fig.2 Original curve, Fourier fitting curve and decomposition waveform of FGM01 in Yingcheng station

將應城臺FGM01地磁抽形圖移植到恩施臺FGM01地磁(圖3)上發現，二者并不重合。取恩施臺FGM01地磁后半段完整數據與同時段應城臺FGM01地磁數據進行傅里葉分解，兩者具有高度相似性，可認為部分與整體呈對應比例關系，即可得出缺失數據的波形參數。

圖3 應城臺FGM01地磁擬合曲線與恩施臺FGM01地磁缺數曲線對比Fig.3 Comparison of Yingcheng FGM01 fitting curve and Enshi FGM01 missing curve

本文取波形個數為300個，對應波形參數如表1所示，得到恩施臺整體波形參數，合成后便可得到修改后的恩施臺FGM01地磁抽形圖(圖4)。

表1 波形參數調整

圖4 調整后應城臺FGM01地磁擬合曲線與恩施臺FGM01地磁缺數曲線對比Fig.4 Comparison of Yingcheng FGM01 fitting curveand Enshi FGM01 missing curve after adjusting

首先對結果進行平穩性與隨機性檢驗。使用Eviews軟件對數據進行單位根檢驗，結果顯示，P值顯著小于0.05，可以確認該數據為平穩數據。再利用Eviews軟件作自相關和偏相關分析，結果如圖5所示，可以看出，該數據自相關系數拖尾，而偏相關系數3階截尾。借助Eviews軟件得到3階自回歸過程，對應的模型表達式[6]為：

圖5 原始數據相關性及ARMA參數估計Fig.5 Correlation diagram of original data andestimation of ARMA parameters

Δyt=0.565 307Δyt-1+0.172 354Δyt-2+

0.170 496Δyt-3+vt

(7)

為了對比分析，將缺數數據減小0.05 nT，其補齊效果如圖6所示。

圖6 ARMA模型修補效果Fig.6 ARMA model supplement diagram

將補齊后的差值與調整后的擬合結果疊加，便可得到完整曲線，如圖7所示，這樣補齊的數據既能保證地磁數據的曲線形態，也能保留其自身特征，可認為是符合精度要求的。從圖7還可看出，修補數據的干擾較完整數據大，這是因為該時段本就受人為干擾，修補的數據是符合實際結果的。為方便對比，將修補的數據減小0.2 nT。

圖7 修補數據曲線與缺數曲線對比Fig.7 Comparison chart of complemented curve and missing curve

從圖8可以看出，修補的數據既保證了原有數據的完整性，又保留了地磁數據的特性，驗證了本文策略的可靠性。

圖8 修補后恩施臺FGM01地磁數據與其他臺站秒數據對比Fig.8 Comparison of Enshi FGM01 data and other second data of other stations after supplement

為進一步驗證本文策略的適用性與可靠性，選取2021-06-15恩施臺FGM01地磁Z分量數據作為實驗對象，將中間1 h數據人為去掉，其原始數據(減小2 nT)與修補數據(減小5 nT)的對比及殘差如圖9所示。

圖9 原始數據與修補數據曲線對比及對應殘差Fig.9 Curve comparison and corresponding residuals between original data and patched data

由圖9可以看出，修補數據與原始數據在形態與細節上幾乎一致，與真實值之差基本在0.2 nT范圍內，屬于背景噪聲，可認為是準確的。但隨著時間的推移，修補數據與真實值之間的差異越來越大，說明該修補模型僅適用于對短期數據的修補，不適合大面積補數。

3 結 語

本文基于素描的思路，提出一種新的地磁秒數據修補策略，以相似度高的地磁數據作輪廓，結合ARMA模型進行細節填充，從而完成對地磁秒數據的修補。具體步驟如下：1)對地磁秒數據進行預處理，依據格拉布期準則剔除突跳點并進行歸一化處理，以消除量綱；2)從數據庫中尋找相似度高的地磁秒數據作為輪廓依托，抽取擬合曲線；3)對擬合曲線分解出的波形進行參數調整，利用調整后的曲線擬合缺省數據；4)將調整后的擬合曲線與缺數數據的差值平穩化，便于ARMA模型進行預測；5)用ARMA模型補齊缺失差值，然后疊加調整后的擬合數據，便可得到修補后的數據。

由于ARMA模型不具備長期預測精度，若地磁數據缺失1 d以上，預測得到的值將不準確，因此本文主要針對1 d內缺數作研究。實驗結果顯示，本文策略可獲得較滿意的修補效果，為地磁秒數據補齊策略提供了新的思路。