地磁秒數(shù)據(jù)補(bǔ)齊策略探討

栗 寧 李從芬 賀北芳

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

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

針對這一問題，本文通過研究湖北省所有地磁臺站的秒數(shù)據(jù)文件發(fā)現(xiàn)，相鄰臺站數(shù)據(jù)的相關(guān)性高，且數(shù)據(jù)呈周期性，可依據(jù)相鄰臺站數(shù)據(jù)之間的相關(guān)性提取共同特性，再跟據(jù)數(shù)據(jù)自身個(gè)性補(bǔ)齊細(xì)節(jié)，從而盡可能不失真地補(bǔ)齊數(shù)據(jù)。

1 補(bǔ)數(shù)原理及實(shí)現(xiàn)

為了不失真地補(bǔ)齊數(shù)據(jù)，需要做到以下幾點(diǎn)：1)選取與所缺數(shù)據(jù)相關(guān)性高的數(shù)據(jù)，并確定移植的抽形圖；2)對缺失的數(shù)據(jù)依據(jù)抽形圖進(jìn)行補(bǔ)數(shù)；3)對補(bǔ)好的抽形圖進(jìn)行細(xì)節(jié)修補(bǔ)；4)對修補(bǔ)完成的數(shù)據(jù)進(jìn)行檢測。

1.1 數(shù)據(jù)選取及傅里葉分解合成

對于一組包含缺失值的地磁數(shù)據(jù)，可在數(shù)據(jù)庫中找到一個(gè)與其最相似的對象。不同問題導(dǎo)致的缺數(shù)可選用不同標(biāo)準(zhǔn)進(jìn)行相似性判定，最常見的是使用相關(guān)系數(shù)矩陣法，以選取與缺失數(shù)據(jù)關(guān)聯(lián)性最強(qiáng)的數(shù)據(jù)作為抽形圖。

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

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

(1)

式中，a0、ak、bk為傅里葉系數(shù)。由于地磁數(shù)據(jù)基本為秒采樣序列，設(shè)為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，即可求出傅里葉級數(shù)擬合中各項(xiàng)的傅里葉系數(shù)，根據(jù)2組數(shù)據(jù)之間波形的對比便可得到缺數(shù)數(shù)據(jù)的擬合圖形。

1.2 基于ARMA預(yù)測模型的細(xì)節(jié)補(bǔ)齊算法

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

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

(6)

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

利用ARMA模型進(jìn)行預(yù)測的前提是輸入的時(shí)間序列是平穩(wěn)的，若為非平穩(wěn)序列則需要通過差分方法轉(zhuǎn)化為平穩(wěn)序列后再進(jìn)行ARMA模型擬合[5]。ARMA模型建模的預(yù)測流程見圖1。

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

2 實(shí)驗(yàn)結(jié)果

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

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

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

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

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

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

表1 波形參數(shù)調(diào)整

圖4 調(diào)整后應(yīng)城臺FGM01地磁擬合曲線與恩施臺FGM01地磁缺數(shù)曲線對比Fig.4 Comparison of Yingcheng FGM01 fitting curveand Enshi FGM01 missing curve after adjusting

首先對結(jié)果進(jìn)行平穩(wěn)性與隨機(jī)性檢驗(yàn)。使用Eviews軟件對數(shù)據(jù)進(jìn)行單位根檢驗(yàn)，結(jié)果顯示，P值顯著小于0.05，可以確認(rèn)該數(shù)據(jù)為平穩(wěn)數(shù)據(jù)。再利用Eviews軟件作自相關(guān)和偏相關(guān)分析，結(jié)果如圖5所示，可以看出，該數(shù)據(jù)自相關(guān)系數(shù)拖尾，而偏相關(guān)系數(shù)3階截尾。借助Eviews軟件得到3階自回歸過程，對應(yīng)的模型表達(dá)式[6]為：

圖5 原始數(shù)據(jù)相關(guān)性及ARMA參數(shù)估計(jì)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)

為了對比分析，將缺數(shù)數(shù)據(jù)減小0.05 nT，其補(bǔ)齊效果如圖6所示。

圖6 ARMA模型修補(bǔ)效果Fig.6 ARMA model supplement diagram

將補(bǔ)齊后的差值與調(diào)整后的擬合結(jié)果疊加，便可得到完整曲線，如圖7所示，這樣補(bǔ)齊的數(shù)據(jù)既能保證地磁數(shù)據(jù)的曲線形態(tài)，也能保留其自身特征，可認(rèn)為是符合精度要求的。從圖7還可看出，修補(bǔ)數(shù)據(jù)的干擾較完整數(shù)據(jù)大，這是因?yàn)樵摃r(shí)段本就受人為干擾，修補(bǔ)的數(shù)據(jù)是符合實(shí)際結(jié)果的。為方便對比，將修補(bǔ)的數(shù)據(jù)減小0.2 nT。

圖7 修補(bǔ)數(shù)據(jù)曲線與缺數(shù)曲線對比Fig.7 Comparison chart of complemented curve and missing curve

從圖8可以看出，修補(bǔ)的數(shù)據(jù)既保證了原有數(shù)據(jù)的完整性，又保留了地磁數(shù)據(jù)的特性，驗(yàn)證了本文策略的可靠性。

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

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

圖9 原始數(shù)據(jù)與修補(bǔ)數(shù)據(jù)曲線對比及對應(yīng)殘差Fig.9 Curve comparison and corresponding residuals between original data and patched data

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

3 結(jié) 語

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

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