董寶偉 錢秋亮 任亞飛,2 陶秋喆,2 邵建龍

1 昆明理工大學(xué)信息工程與自動化學(xué)院，昆明市景明南路727號，650504 2 云南省麻栗坡民族中學(xué)，云南省文山州文麻路193號，663600

地震地磁數(shù)據(jù)預(yù)處理是地磁數(shù)據(jù)挖掘和地震預(yù)報預(yù)測的首要工作[1]。地磁觀測數(shù)據(jù)屬于時間序列，其計(jì)算方法要求數(shù)據(jù)序列連續(xù)完整，因此對缺失數(shù)據(jù)進(jìn)行插值是前兆數(shù)據(jù)應(yīng)用前必須進(jìn)行的預(yù)處理工作[2]。國內(nèi)外對缺失值進(jìn)行插值處理的方法較多，如回歸法、最近鄰域替代法、人工填補(bǔ)法和期望值最大化法等[3-5]。武艷強(qiáng)等[6]利用三次樣條插值方法對時間序列中數(shù)據(jù)缺失較多等問題展開研究；楊登科[7]采用拉格朗日插值等6種方法對時間序列的插值效果進(jìn)行對比。

在實(shí)際地磁偏角觀測工作中，由于地磁觀測數(shù)據(jù)[8-10]種類眾多，往往會出現(xiàn)單點(diǎn)缺失和連續(xù)多點(diǎn)缺失的現(xiàn)象，而現(xiàn)有的數(shù)據(jù)插值方法在地磁觀測數(shù)據(jù)的插值應(yīng)用中存在諸多局限性。目前針對地磁偏角數(shù)據(jù)插值方法的研究較少，本文根據(jù)地磁數(shù)據(jù)的記錄特點(diǎn)，將時間序列模型[11-12]應(yīng)用于地磁數(shù)據(jù)缺失處理。將ARMA預(yù)測模型[13-14]用于地磁數(shù)據(jù)的插值處理，可為地磁觀測數(shù)據(jù)預(yù)處理提供一種可行的方法。

1 插值方法

插值法是數(shù)據(jù)預(yù)處理的重要方法，在數(shù)據(jù)分析和數(shù)據(jù)挖掘中發(fā)揮著重要作用。常見的插值法有均值插值、線性插值等，不同方法的適用性也不同，如均值插值適用于原始數(shù)據(jù)方差較小、序列趨于平穩(wěn)的數(shù)據(jù)，而本文所采用的地磁觀測數(shù)據(jù)屬于非平穩(wěn)數(shù)據(jù)；線性插值為直線擬合運(yùn)算，適用于趨勢性較小的數(shù)據(jù)。但在處理連續(xù)多點(diǎn)缺失的數(shù)據(jù)時，上述方法的插值效果均不理想。本文嘗試將ARMA模型預(yù)測方法應(yīng)用于地磁序列插值。

1.1 均值插值

均值插值就是計(jì)算數(shù)據(jù)的平均值并代替原始數(shù)據(jù)中的缺失部分，假設(shè)求解Ti和Ti+1之間任意一點(diǎn)T，則直接取T為Ti和Ti+1的平均值。對于數(shù)據(jù)集data=[45.20,26.21,42.10,NaN,NaN,NaN,52.20,32.52]，Python可調(diào)用均值mean()函數(shù)計(jì)算均值，利用自帶fillna()函數(shù)進(jìn)行mean插值計(jì)算，計(jì)算方法如下：

data = data.fillna(data.mean())

(1)

1.2 線性插值

線性插值是通過已知量來求解兩者之間的未知量，由于其為直線擬合，適合處理較為平緩的數(shù)據(jù)。在Python中可用interpolate()函數(shù)進(jìn)行線性運(yùn)算，利用fillna()函數(shù)進(jìn)行interpolate插值，計(jì)算方法如下：

data = data.fillna(data.interpolate())

(2)

或

data=data.interpolate(method=‘linear’)

(3)

1.3 自回歸移動平均(ARMA)模型

本文嘗試將ARMA模型作為插值方法引入到地磁觀測數(shù)據(jù)插值預(yù)處理中。將缺失數(shù)據(jù)作為時間序列預(yù)測數(shù)據(jù)，利用缺失序列之前的數(shù)據(jù)預(yù)測缺失數(shù)據(jù)，插值效果較好。

ARMA模型作為一種時間序列經(jīng)典模型，由自回歸模型(AR)和移動平均模型(MA)組合而成，其結(jié)構(gòu)為：

xt=φ0+φ1xt-1+φ2xt-2+…+φpxt-p+

εt-θ1εt-1-θ2εt-2-…-θqεt-q

(4)

式中，隨機(jī)變量Xt的取值xt不僅與前p期的序列值有關(guān)，還與前q期的隨機(jī)擾動有關(guān)，p為自回歸模型階數(shù)，q為移動平均模型階數(shù)。

ARMA(p,q)模型預(yù)測步驟如下。

1)非白噪聲檢驗(yàn)：使用Python自帶函數(shù)庫statsmodels中acorr_ljungbox函數(shù)進(jìn)行檢驗(yàn)，若p值小于顯著性水平α，則該序列為非白噪聲，具備建模分析價值。

2)平穩(wěn)性檢驗(yàn)：使用Python自帶函數(shù)庫statsmodels中adfuller函數(shù)進(jìn)行檢驗(yàn)，也可通過地磁序列的自相關(guān)函數(shù)(ACF)和偏自相關(guān)函數(shù)(PACF)來確定序列是否平穩(wěn)，若不平穩(wěn)可使用差分方法使數(shù)據(jù)序列達(dá)到平穩(wěn)，進(jìn)而確定模型階數(shù)(p、q值)。

3)估計(jì)：對模型參數(shù)進(jìn)行評估。

4)預(yù)測：利用擬合模型進(jìn)行預(yù)測分析。

為驗(yàn)證利用ARMA模型對地磁數(shù)據(jù)進(jìn)行插值的可行性和有效性，將常用的均值插值、線性插值的插值結(jié)果作為參照，與ARMA模型的插值結(jié)果進(jìn)行對比研究。

2 地磁數(shù)據(jù)插值實(shí)驗(yàn)設(shè)計(jì)

2.1 前期地磁數(shù)據(jù)處理

ARMA模型要求數(shù)據(jù)序列必須為平穩(wěn)數(shù)據(jù)，對于地磁觀測數(shù)據(jù)而言，數(shù)據(jù)序列由各種趨勢、周期和隨機(jī)干擾組合而成，因此在建立ARMA模型之前需要去除數(shù)據(jù)中的趨勢和周期成分。通過分析地磁偏角觀測數(shù)據(jù)模型，可將地磁數(shù)據(jù)的曲線形態(tài)分為三類：1)隨機(jī)干擾；2)地球背景場或太陽輻射周期；3)周期或趨勢成分。由于地震儀器所處的臺站和自然環(huán)境不同，其記錄的地磁數(shù)據(jù)均包含不同形態(tài)的趨勢和周期成分。

去除趨勢或周期、地球背景場或太陽輻射周期最簡單的方法是對地磁觀測序列進(jìn)行差分處理。對于多數(shù)觀測序列而言，一階差分即可去除大部分干擾成分，使數(shù)據(jù)序列達(dá)到平穩(wěn)，若一階差分?jǐn)?shù)據(jù)仍未達(dá)到平穩(wěn)，可進(jìn)行多次差分。對大量觀測數(shù)據(jù)序列進(jìn)行統(tǒng)計(jì)分析發(fā)現(xiàn)，約92.43%的地磁觀測數(shù)據(jù)序列為非平穩(wěn)狀態(tài)，經(jīng)過一階差分后，96.78%的非平穩(wěn)序列均可達(dá)到平穩(wěn)[1]。在對地磁偏角進(jìn)行處理時發(fā)現(xiàn)，為了使觀測數(shù)據(jù)趨于平穩(wěn)化，對觀測數(shù)據(jù)進(jìn)行二階差分，非平穩(wěn)數(shù)據(jù)基本可達(dá)到平穩(wěn)化。

2.2 地磁插值實(shí)驗(yàn)方法設(shè)計(jì)

2.2.1 非震異常數(shù)據(jù)數(shù)量≤10

為比較ARMA模型與其他插值方法的插值效果，選取完整的地磁觀測數(shù)據(jù)序列，人為構(gòu)造缺失數(shù)據(jù)。數(shù)據(jù)缺失點(diǎn)包括單點(diǎn)缺失和連續(xù)多點(diǎn)缺失，本文設(shè)計(jì)的連續(xù)多點(diǎn)缺失為2～10點(diǎn)缺失，利用3種不同的插值方法對缺失數(shù)據(jù)序列進(jìn)行插值處理。選擇缺失值之前的400個數(shù)據(jù)點(diǎn)作為ARMA模型的基礎(chǔ)，并確定參數(shù)最佳的ARMA模型，然后將預(yù)測值作為缺數(shù)填補(bǔ)結(jié)果。采用標(biāo)準(zhǔn)誤差作為評價指標(biāo)，標(biāo)準(zhǔn)誤差越小，插值效果越優(yōu)。標(biāo)準(zhǔn)誤差計(jì)算公式如下：

(5)

式中，n為缺值個數(shù)，xi為插值結(jié)果，yi為實(shí)際觀測值。

2.2.2 插值實(shí)驗(yàn)及分析

選取地磁偏角D數(shù)據(jù)，對地磁數(shù)據(jù)各種缺失情況進(jìn)行插值實(shí)驗(yàn)，步驟如下：

1)在地磁偏角數(shù)據(jù)中選取10個位置點(diǎn)作為缺失數(shù)據(jù)的起始點(diǎn)。

2)依據(jù)起始點(diǎn)連續(xù)構(gòu)造n(n取值1～10，分別對應(yīng)10種缺失情況)對缺失數(shù)據(jù)序列進(jìn)行差分處理。

3)采用3種不同方法進(jìn)行插值處理。

4)利用式(5)計(jì)算3種不同方法的標(biāo)準(zhǔn)誤差。

表1為實(shí)驗(yàn)結(jié)果，從表中可以看出，3種插值方法中ARMA預(yù)測插值法的插值效果較好，其次為線性插值法和均值插值法。圖1為3種插值方法的平均標(biāo)準(zhǔn)誤差曲線，從圖中可以看出，在不同缺失值情況下，3種方法的標(biāo)準(zhǔn)誤差曲線趨勢表明，隨著缺失值個數(shù)的增加，插值結(jié)果與實(shí)際值的標(biāo)準(zhǔn)誤差也在增加，且不同插值方法標(biāo)準(zhǔn)誤差的增長速度有所不同，相比于線性插值和ARMA預(yù)測插值，均值插值標(biāo)準(zhǔn)誤差的增長速度最快；ARMA預(yù)測插值的效果最好，其標(biāo)準(zhǔn)誤差最小，標(biāo)準(zhǔn)誤差曲線的變化幅度小于其他方法。

表1 三種不同插值方法在不同缺失條件下的平均標(biāo)準(zhǔn)誤差

圖1 不同缺失值個數(shù)情況下不同插值方法的平均標(biāo)準(zhǔn)誤差比較

圖2為不同單點(diǎn)缺失位置3種插值方法對應(yīng)的標(biāo)準(zhǔn)誤差變化曲線，從圖中可以看出，ARMA預(yù)測插值相比于其他2種方法的標(biāo)準(zhǔn)誤差更小。

圖2 不同單點(diǎn)缺失位置3種方法對應(yīng)的平均標(biāo)準(zhǔn)誤差變化曲線

為了能夠更加直觀地說明3種插值方法插值效果的差異，本文以2008-05-12汶川地震北京臺的秒采樣數(shù)據(jù)為例，給出磁偏角D在連續(xù)缺失10個數(shù)值時3種插值方法的插值效果(圖3)。

圖3 連續(xù)缺失10個數(shù)值的插值情況

實(shí)驗(yàn)表明，在均值插值、線性插值和ARMA預(yù)測插值中，ARMA模型的插值效果最優(yōu)。本文以2017-08-08九寨溝地震滿洲里地震臺連續(xù)缺失10個數(shù)值的地磁偏角秒數(shù)據(jù)為例進(jìn)行數(shù)據(jù)擬合，結(jié)果如圖4所示。

圖4 九寨溝地震滿洲里地震臺連續(xù)缺失10個數(shù)值的地磁偏角秒數(shù)據(jù)ARMA模型擬合

綜上分析可知，不同插值方法對不同缺失情況的插值效果各不相同，表現(xiàn)為曲線形態(tài)特征存在差異。ARMA預(yù)測插值法的標(biāo)準(zhǔn)誤差較小，且該模型中插值數(shù)據(jù)與實(shí)際觀測值擬合度高，而均值插值和線性插值會改變數(shù)據(jù)原始趨勢，因此ARMA插值對地磁偏角D數(shù)據(jù)的處理效果相對較好。

2.2.3 非震異常數(shù)據(jù)數(shù)量>10

在數(shù)據(jù)缺失較多時不應(yīng)采用插值方法，而應(yīng)直接刪除。但本文研究的地磁數(shù)據(jù)屬于時間序列數(shù)據(jù)，為了后期可檢測到數(shù)據(jù)異常點(diǎn)，本文采用統(tǒng)一合理值來代替大量缺失數(shù)據(jù)和非震異常數(shù)據(jù)。非震異常數(shù)據(jù)指地震儀器記錄錯誤或未記錄到的數(shù)據(jù)，其一般為多個連續(xù)數(shù)值88 888.0或99 999.0，或者為多個連續(xù)的小于 99 999.0且遠(yuǎn)大于正常地磁峰值的數(shù)值。由于地磁由多個變化緩慢的磁場組成，其數(shù)值為隨時間變化的時間序列，不可能出現(xiàn)多個連續(xù)相同且遠(yuǎn)大于正常峰值的數(shù)值。因此，該部分?jǐn)?shù)值與真實(shí)值無相關(guān)性，會對地磁的趨勢信息造成不可恢復(fù)的影響。但這部分?jǐn)?shù)據(jù)是時間序列的一部分，而時間序列中每個數(shù)值所在的時間點(diǎn)也是重要信息，因此應(yīng)用合理的數(shù)值代替該部分?jǐn)?shù)據(jù)，有利于保持時間序列的完整性。處理這些數(shù)據(jù)的算法如下：

1)遍歷1 d內(nèi)所有的地磁偏角D數(shù)據(jù)，統(tǒng)計(jì)數(shù)值為88 888.0、99 999.0及其以上數(shù)值的個數(shù)，若個數(shù)小于10，采用§2.2.1算法進(jìn)行插值；若個數(shù)大于10，則采用鄰域值進(jìn)行代替，直到出現(xiàn)正常值。

2)遍歷1 d內(nèi)所有的地磁偏角D數(shù)據(jù)，若數(shù)值出現(xiàn)0或“NaN”的個數(shù)小于10，采用§2.2.1算法進(jìn)行插值；若個數(shù)大于10，則采用鄰域值進(jìn)行代替，直到出現(xiàn)正常值。

3)遍歷1 d內(nèi)所有的地磁偏角D數(shù)據(jù)，若重復(fù)數(shù)值出現(xiàn)的個數(shù)小于10，采用§2.2.1算法進(jìn)行插值；若個數(shù)大于10，則采用鄰域值進(jìn)行代替，直到出現(xiàn)正常值。

以2008-05-12汶川地震滿洲里臺站作為測試臺站，其缺數(shù)中88 888.0、99 999.0、0和”NaN”個數(shù)大于10，數(shù)據(jù)為地磁偏角D的1 d秒采樣數(shù)據(jù)，結(jié)果如圖5所示。

圖5 汶川地震滿洲里臺站多種缺數(shù)大于10

由于包含遠(yuǎn)大于正常值的數(shù)值88 888.0和99 999.0，從圖5中無法看出數(shù)據(jù)走勢，該部分?jǐn)?shù)據(jù)對于后期數(shù)據(jù)異常分析并不重要，但對異常時間點(diǎn)分析必不可少，因此本文采用鄰域值代替該部分?jǐn)?shù)據(jù)(圖6)。鄰域代替屬于就近代替，其數(shù)值差距較小，因此處理后的數(shù)據(jù)和原始數(shù)據(jù)相差較小。該處理過程是為了保證地磁序列的時序連貫性，并不會對后續(xù)分析產(chǎn)生影響。

圖6 鄰域值代替缺失數(shù)據(jù)擬合圖

3 結(jié) 語

為了處理地震儀器記錄的非震異常數(shù)據(jù)，本文將ARMA模型用于數(shù)據(jù)插值處理，并與均值插值、線性插值的實(shí)驗(yàn)效果進(jìn)行對比分析，得到以下結(jié)論：

1)無論是單點(diǎn)插值還是連續(xù)多點(diǎn)插值(2～10個)，ARMA模型的插值效果均優(yōu)于均值插值和線性插值。

2)ARMA模型能夠較好地保持實(shí)際觀測序列的曲線形態(tài)，而均值插值和線性插值會改變其形態(tài)。

3)均值插值適用于方差小、穩(wěn)定性好的數(shù)據(jù)；線性插值適用于曲線形態(tài)較為平坦的序列；ARMA模型適用于規(guī)律性變化、干擾較小的序列。

4)ARMA模型屬于外插值法，采用該模型進(jìn)行預(yù)測時需要建模，因此與均值插值和線性插值相比耗時較多。由于存在較多的歷史數(shù)據(jù)作為參考，即使在多點(diǎn)缺失情況下，ARMA的預(yù)測效果仍高于其他兩種方法。

5)ARMA模型有望在電磁學(xué)領(lǐng)域數(shù)據(jù)預(yù)處理中發(fā)揮作用。