鄧方進， 王緒本， 李德偉

(成都理工大學 地球物理學院,成都 610059)

ARIMA模型在LMT數(shù)據(jù)處理中的應用

鄧方進， 王緒本， 李德偉

(成都理工大學 地球物理學院,成都 610059)

為了解決LMT時間序列出現(xiàn)的缺失和強干擾現(xiàn)象，根據(jù)實測資料數(shù)據(jù)量大、非線性、非平穩(wěn)性等特點，首次采用ARIMA模型進行預測和填補，基于平穩(wěn)性檢驗和貝葉斯信息準則確定模型階數(shù)，采用最小二乘原理確定模型參數(shù)，建立雙向預測模型和線性合并方法進行預測，并對比ARIMA模型和AR模型預測數(shù)據(jù)的準確度。實例表明，ARIMA模型預測結果準確，精度比AR模型高，且誤差不會累積，解決了原始資料的不連續(xù)性和強干擾的問題。

ARIMA模型; 預測; 大地電磁測深; 資料處理; 缺失數(shù)據(jù)

0 引言

在大地電磁測深(Magnetotelluric，MT)資料處理過程中，為了得到研究區(qū)域合理的地電屬性特征，就必須保證有較高質量的野外原始觀測數(shù)據(jù)。方勝[1]總結了在大地電磁野外數(shù)據(jù)采集過程中，如何保證數(shù)據(jù)質量的操作經驗。但是實際的長周期大地電磁測深(Long period Magnetotelluric，LMT)的數(shù)據(jù)測量位置選在遠離城區(qū)電磁干擾較小的平坦地區(qū)，測量周期長達半個月甚至數(shù)個月，在LMT儀器的運行過程中，由于儀器的響應問題，采集到的電磁場時間序列偶爾會出現(xiàn)部分缺失和隨機強干擾現(xiàn)象。數(shù)據(jù)的缺失對長周期資料處理影響非常大，長周期數(shù)據(jù)地處理需要半個月甚至幾個月的數(shù)據(jù)進行分析和處理，并且時間序列是需要連續(xù)不間斷的。為了保證大地電磁測深資料的準確度，就需要選擇一個適當?shù)姆椒▽θ笔?shù)據(jù)進行填補。

Box&Jenkins[2]提出著名的差分自回歸滑動平均模型(Autoregressive Integrated Moving Average Model, ARIMA)時間序列預測方法，能夠準確地預測非平穩(wěn)時間序列，只需要給出時間序列的前N個連續(xù)數(shù)據(jù)，就可以采用該預測模型進行數(shù)據(jù)的預測，并針對澳大利亞1972年到1991年的居民消費價格指數(shù)(Consumer Price Index，CPI)，根據(jù)自相關函數(shù)和偏自相關函數(shù)的性質確定模型階數(shù)建立ARIMA模型，預測未來四年的CPI數(shù)值；許立平等[3]根據(jù)1973年1月至2010年11月倫敦黃金月度價格建立ARIMA模型，預測2011上半年的黃金價格走勢，并得出短期上漲的結論；張勃等[4]將ARIMA模型用于生態(tài)足跡動態(tài)模擬和預測。關于ARIMA模型在地球物理領域的應用，唐恒專[5]應用AR模型對地震信號P波初至時段及其后一時段進行建模，求出預測誤差并將其作為特征判據(jù)對地震信號進行分類和識別。

1 MT數(shù)據(jù)缺失和強干擾現(xiàn)象

MT資料處理的目的是從采集到的大地電磁時間序列中，提取出能夠反映研究區(qū)域的地電特性的各種頻率成分，通過計算得到準確的阻抗張量要素、視電阻率、相位以及其他相關的大地電磁響應函數(shù)，從而獲得研究區(qū)域的電性結構，為資料解釋提供依據(jù)[6]。目前，大地電磁測深數(shù)據(jù)處理方法有遠參考法[7]、Robust法[8]、Hilbert-Huang變換法[9]、S變換[10]等。常使用的方法是采用傅里葉變換將大地電磁場的時間序列向頻率域轉換，這些處理方法都是針對原始資料的時間序列是連續(xù)的。

LMT數(shù)據(jù)采集儀器將采集到的數(shù)據(jù)保存到TF存儲卡里，存儲卡中的數(shù)據(jù)包含了表征地電結構的原始數(shù)據(jù)和其他干擾數(shù)據(jù)，其中干擾數(shù)據(jù)來源于外界影響和內部影響，其中外界影響包含了人為電磁場干擾和自然電磁場干擾，如高壓輸電線、工廠、磁暴和雷電等；內部影響由于儀器內部元件的影響，如電極、磁棒和儀器主機的故障。在資料整理過程中，需要盡量剔除噪聲干擾數(shù)據(jù)，常用平滑濾波的方法篩選出質量較高的數(shù)據(jù)，將挑選出來的數(shù)據(jù)再進行下一步處理。但是實際的野外工作中，隨機干擾數(shù)據(jù)可能會導致原始資料在某一段時間內出現(xiàn)電磁場數(shù)值特別大或者特別小的情況。對于這種情況，采用濾波方法進行消除就會對資料的有效信號產生影響，所以可以采用一種合理的方法將隨機干擾數(shù)據(jù)進行置換。

為了得到低頻段的資料，要盡可能減小截斷效應，增大記錄長度，保證低頻段也有足夠多的疊加次數(shù)，由于探測深度的不同，長周期大地電磁信號對應的周期范圍在10 s～20 000 s，通常要保證記錄長度不低于最低頻率所對應周期的6倍，則需要在測點上采集足夠長的時間。經過長期的野外測量工作，特別是針對采集時間長達半個月甚至數(shù)個月的LMT探測會出現(xiàn)這樣的問題：在野外的數(shù)據(jù)采集過程中，LMT測量儀器會偶爾出現(xiàn)故障，儀器工作不穩(wěn)定的狀況，這種偶然情況就會導致記錄到的電磁場數(shù)據(jù)出現(xiàn)缺失，缺失數(shù)據(jù)1個～30個不等。在進行資料處理時，對于這種原始資料缺失情況導致時間序列不連續(xù)，可能會造成在記錄狀態(tài)下的高質量數(shù)據(jù)無法利用。因此，找到一種合理的預測模型對有缺失的LMT數(shù)據(jù)進行填補，彌補由于儀器硬件原因造成的不足，可以提高LMT資料的利用率。

對LMT數(shù)據(jù)的分析及處理，張偉[11]在信號去噪和時頻分析等方面做了很多研究，張春虹等[12]將無激勵的AR預測模型用于預測大地電磁測深缺失數(shù)據(jù)，但是存在預測精度和預測個數(shù)的不足。

2 大地電磁場信號特征分析

2.1 信號的分類

在實際信號的過程中，盡管信號的處理方法和基本原理都相同，但是針對不同學科領域的不同信號或者同一學科的不同信號，處理過程也會存在很多差異。所以在做信號處理之前，首先要明確信號的類別，現(xiàn)代信號處理中對信號的特征將信號分為以下幾類[13]：

1)線性信號和非線性信號。線性信號就是基于線性系統(tǒng)的信號，該系統(tǒng)具有齊次性和可加性兩個性質；而非線性信號不滿足上述線性系統(tǒng)的性質。

2)平穩(wěn)信號和非平穩(wěn)信號。平穩(wěn)信號滿足平穩(wěn)性特征，最大特點是它的統(tǒng)計特征與時間分布無關。其特征為：①序列均值為常數(shù); ②序列二階矩有界; ③序列協(xié)方差函數(shù)與時間無關，只與時間間隔有關; ④序列聯(lián)合概率分布函數(shù)與時間無關，只與時間間隔有關。不滿足以上特征的信號為非平穩(wěn)信號。

3)高斯信號和非高斯信號。高斯信號定義為概率密度分布為正態(tài)分布的隨機信號。在隨機過程中使用偏斜度和峭度兩個常量來描述該隨機過程，高斯信號的偏斜度和峭度都為零，而不滿足該類條件的隨機信號都為非高斯信號。

通過對信號進行分類，確定信號的類型，在處理數(shù)據(jù)時，根據(jù)信號特征選擇合適處理的方法，才能“對癥下藥”，準確地對信號進行處理，得到合理的處理結果。

2.2 實測大地電磁信號特征分析

大地電磁場是指地球天然電磁場中隨時間變化的電場和磁場，頻率范圍寬。它們是由太陽不斷射出的粒子流與電磁輻射在地球周圍空間所引起的電磁效應所形成的。通過采集天然變化的電磁場分量，將電磁場信號轉換成視電阻率曲線和相位曲線，反演求得各地層的電阻率和厚度值。

20世紀中期，進行大地電磁數(shù)據(jù)處理前，通常假設天然場的大地電磁信號是具有線性、因果性、平穩(wěn)性和最小相位性等特征的隨機信號，并且統(tǒng)計特征滿足高斯分布。21世紀初期，隨著現(xiàn)代信號處理方法的發(fā)展，對信號特征的分析研究取得了新的進展，大地電磁信號的研究思路也發(fā)生了重大轉變。王書明等[14]根據(jù)大量實測資料進行分析研究，其結果表明大地電磁信號的零極點圖不滿足最小相位系統(tǒng)的要求；分析實測資料的高階統(tǒng)計量，證明了大地電磁信號是非高斯性，通過相關系數(shù)證明信號是非平穩(wěn)的，并且不滿足線性要求。

LMT野外資料采集時間長，本次試驗采樣頻率為1 Hz，以200 h的觀測時長進行計算，那么采集的電磁場一個分量數(shù)據(jù)達到了720 000個，電磁場包含了電磁的5個分量，那么五道數(shù)據(jù)就共計3 600 000個，數(shù)據(jù)文件大小近百兆。

圖1為陜西某地LMT采集資料中某一天的記錄數(shù)據(jù)，使用LEMI-417儀器，采樣率為1 Hz，采集大地電磁場五分量數(shù)據(jù)，包含兩個電場分量Ex和Ey，三個磁場分量Hx、Hy和Hz。

圖1 陜西某地實測一天的LMT數(shù)據(jù)曲線圖Fig.1 Graph of LMT data in Shanxi province

由圖1可以看出，每一條曲線圖沒有一個明確的趨勢，很難用一個確切的函數(shù)來表達。整個時間序列中，每一道數(shù)據(jù)的振幅的幅值極值相差不大，變化較小，相對來說，磁場的變化比電場的變化大，電場的曲線趨勢相對平穩(wěn)。但是在其中50 000～55 000數(shù)據(jù)區(qū)間中，見紅色虛框內，電場兩個方向的都出現(xiàn)了一個明顯的跳動(這是偶然現(xiàn)象，出現(xiàn)次數(shù)不多)，出現(xiàn)跳動的時間點一致，證明這些時間點的數(shù)據(jù)存在信號干擾。將Hx道中數(shù)據(jù)中20 000～210 000中的1 000個數(shù)據(jù)(約17 min)提出來另成圖像，可以看出，在相對上升趨勢的曲線中，也存在諸多類似正余弦信號的疊加信號。

隨著時間序列非平穩(wěn)性的提出，單位根檢驗目前已經成為宏觀數(shù)據(jù)建模前首先要進行的工作。為此，Dickey等[15]提出了著名的單位根檢驗方法(Augmented Dickey-Fuller test， ADF檢驗)，用于檢驗序列的平穩(wěn)性。

ADF檢驗是通過下面三個模型完成的：

模型1：

(1)

模型2：

(2)

模型3：

(3)

式中：t為時間變量；α為常數(shù)項；βt為趨勢項；εt為殘差Xt項；Xt為原始序列；為ΔXt差分序列。原假設H0:δ=0，檢驗時從模型3開始，然后模型2，最后模型1，若檢驗拒絕H0:δ=0，則原始時間序列不存在單位根，即序列為平穩(wěn)序列，即可停止，若檢驗接受H0:δ=0，則序列為非平穩(wěn)序列。通過對大地電磁場的電磁數(shù)據(jù)進行平穩(wěn)性分析，檢驗結果為非平穩(wěn)。

長周期大地電磁儀所采集的五分量大地電磁場時域數(shù)據(jù)，屬于時間連續(xù)的隨機非周期信號，并且具有非線性、非高斯性、非最小相位性等特點。野外測量時間長，在時域表現(xiàn)出數(shù)據(jù)量龐大和非平穩(wěn)性的特征，這些特征為后續(xù)信號處理方法的選擇提供了依據(jù)。

3 ARIMA模型及其預測方法

3.1 ARIMA預測模型

ARIMA 模型的基本原理：把時間序列視為隨機過程，用一個數(shù)學模型來描述或模擬；可用該時間序列的過去值和現(xiàn)值來預測未來值。該模型考察了時間序列的動態(tài)、持續(xù)特征，揭示了時間序列過去與現(xiàn)在、將來與現(xiàn)在的相互關系。

若序列{X}能通過d次差分后變成平穩(wěn)序列{Y}，則：

yt=Δdxt

(4)

那么可以建立模型

(5)

(6)

3.2 ARIMA模型階數(shù)的確定及參數(shù)估計

對于非平穩(wěn)序列，相應的措施是進行差分處理使其平穩(wěn)。在連續(xù)變化的時間范圍內，變量X關于時間t的變化率用dX/dt來表示；對于離散的變量{X}，我們常取在規(guī)定的時間區(qū)間上的差商ΔX/Δt來表示變量的變化率，如果選擇Δt=1，則

ΔXt=Xt-Xt-1

(7)

其中：時間t取非負數(shù)；ΔXt則為差分序列。對差分序列再進行ADF檢驗，若檢驗結果還是非平穩(wěn)，則繼續(xù)進行差分處理，直到檢驗結果為平穩(wěn)。期間所進行的差分次數(shù)則為ARIMA模型的差分階數(shù)d。

確定了差分次數(shù)，接下來確定階數(shù)p(AR模型的階數(shù))和q(MA模型的階數(shù))。常用的定階方法：最終預測誤差[16]、阿凱克信息準則(Akaike information criterion，AIC)[17]、BIC[18]以及奇異值分解[19]等。

AIC準則法的原理是使得式(8)最小。

(8)

筆者采用貝葉斯信息準則法(Bayesian Information Criterion，BIC)，是一種改進的AIC準則法。原理是使得式(9)為最小，BIC的準則比AIC更嚴格。

(9)

根據(jù)經驗法[20]，當N=20～50時，AR模型最大階數(shù)定為N/2,當N=50～100時，最大階數(shù)選在N/3～N/2,當N=100～200時，最大階數(shù)定在2N/ln(2N)效果會相對比本文好，所以在進行模型定階時，可以依據(jù)經驗法縮小階數(shù)確定范圍。

模型參數(shù)的精細估計采用最小二乘準則來進行，并有近似極大似然估計等估計方法。

3.3 雙向預測和線性合并

通過對比大量的預測數(shù)據(jù)與真實數(shù)據(jù)，發(fā)現(xiàn)隨著時間地增長，其預測誤差會越來越大。結合LMT缺失數(shù)據(jù)的特性，數(shù)據(jù)的缺失部分不會出現(xiàn)在數(shù)據(jù)的開頭和結尾部分，那么就可以結合從前往后的預測結果和從后往前的預測結果綜合決定最后的結果。具體方法是首先從缺失數(shù)據(jù)前的N個數(shù)據(jù)預測k個數(shù)據(jù)，采用ARIMA模型的預測結果為ans1(k)，然后從缺失數(shù)據(jù)后的N個數(shù)據(jù)反向預測k個數(shù)據(jù)，預測結果為ans2(k)，設計線性比例a(k)和b(k)：

(10)

與

(11)

通過線性比例，將兩次前后分別進行的預測值進行合并，最后的合成結果為ans(k)：

ans(i)=a(i)*ans1(i)+b(i)*ans2(i)

i=1,2,…,k

(12)

將此序列作為最終的預測結果，因此隨著預測個數(shù)的增加，預測值和實測值誤差不會累積，理論上預測數(shù)據(jù)的開頭和結尾部分誤差相對較小。

3.4 ARIMA模型預測步驟

采用ARIMA預測模型對含有缺失數(shù)據(jù)的序列進行預測，對需要置換的數(shù)據(jù)進行修補，程序流程如圖2所示。

圖2 ARIMA模型數(shù)據(jù)預測流程圖Fig.2 The flow chart of predict data using ARIMA model

1)輸入需要進行預測的大地電磁序列，確定序列個數(shù)N，預測數(shù)據(jù)個數(shù)k，提取預測區(qū)間前N列和后N列，分離有用數(shù)據(jù)存入二維數(shù)組x[N][10]。

2)將預測區(qū)間后N列數(shù)據(jù)進行反向處理，然后對每列數(shù)據(jù)進行平穩(wěn)性檢驗，如果檢驗結果非平穩(wěn)，則進行差分處理，直到序列平穩(wěn)。

3)通過BIC準則確定模型階數(shù)，采用最小二乘方法進行參數(shù)估計，再進行參數(shù)估計及數(shù)據(jù)預測，并通過雙向預測進行數(shù)據(jù)合并。

4)進行預測數(shù)據(jù)與真實數(shù)據(jù)的誤差分析，并輸出結果。

4 ARIMA預測模型應用實例

4.1 大地電磁缺失數(shù)據(jù)的預測

以LEMI-417長周期大地電磁儀所采集的五分量大地電磁場數(shù)據(jù)為例，利用ARIMA模型分別對LMT資料的時間序列中磁場x分量(Hx)、磁場y分量(Hy)、磁場z分量(Hz)、電場x分量(Ex)和電場y分量(Ey)中缺失的時域數(shù)據(jù)進行預測，筆者以采樣率為1 s的LMT時域序列作為實驗樣本數(shù)據(jù)，隨機取430個約4 min的數(shù)據(jù)作為分析模型(圖3)。

圖3中包含了電磁場五個分量的曲線圖，每個分量的形態(tài)不一，有相對平穩(wěn)序列如Ex；非平穩(wěn)序列，如Hz存在下降趨勢；數(shù)據(jù)波動較大的曲線如Ey。圖3中紅色虛線內的數(shù)據(jù)(200-230)作為本次試驗的預測目標數(shù)據(jù)，ARIMA預測模型所使用到的數(shù)據(jù)包含預測區(qū)間前200個數(shù)據(jù)和預測區(qū)間后200個數(shù)據(jù)。

圖3 原始LMT時間序列曲線圖Fig.3 Time series graph of raw LMT data

4.2 預測結果分析

通過ARIMA模型對時間序列缺失的30個數(shù)據(jù)進行預測填補后，將預測結果與實測數(shù)據(jù)進行對比，為了直觀反映預測結果與實測數(shù)據(jù)進行對比，圖4給出了ARIMA模型預測后的各分量的缺失數(shù)據(jù)和實際樣本數(shù)據(jù)的散點圖比較。

圖4 預測數(shù)據(jù)和實測數(shù)據(jù)誤差對比圖Fig.4 The comparison of predicted data and observed data

由圖4可以看出，ARIMA預測的數(shù)據(jù)和實測數(shù)據(jù)基本重合，其曲線走向基本一致。AR模型預測的數(shù)據(jù)當曲線走勢較為平穩(wěn)時，其效果較好，但是當曲線有明顯的上升或者下降的趨勢時(如Ey)，其預測結果前5個的準確度較高，而其后面的數(shù)據(jù)誤差則越來越大，這是因為誤差累積造成的。所以ARIMA模型預測的結果準確度比采用AR模型預測的結果高，隨著預測個數(shù)的增加，預測數(shù)據(jù)值和真實值的誤差并沒有像常規(guī)預測方法那樣會出現(xiàn)誤差的累積而造成誤差逐漸增大。對于數(shù)據(jù)本身波動較大的情況，預測誤差相對較大。

表1為ARIMA模型預測數(shù)據(jù)與真實數(shù)據(jù)的誤差值，表1中ΔHx、ΔHy和ΔHz分別表示磁場的x分量、y分量和z分量預測值和真實值的誤差絕對值，ΔEx和ΔEy分別表示電場x分量和y分量預測值和真實值的誤差絕對值。

表1 預測數(shù)據(jù)與實測數(shù)據(jù)誤差表

從表1中可以看出各個分量的整體誤差值都比較小，通過計算平均誤差為0.06，可以滿足作為資料處理的要求。

為了直觀地表示ARIMA預測模型和AR預測模型預測數(shù)據(jù)與實測數(shù)據(jù)的誤差對比，對兩者的誤差數(shù)據(jù)繪制誤差散點圖，如圖5所示。

圖5 預測數(shù)據(jù)與實測數(shù)據(jù)誤差散點圖Fig.5 Error scatterplot of predicted data and observed data

由圖5可以看出，ARIMA模型預測數(shù)據(jù)的誤差絕對值都保持在0.2以下，并且誤差值不會隨著預測個數(shù)的增加而增加。而采用AR模型預測的結果的誤差總體偏大，最大可達到0.8，由于誤差累積，隨著預測個數(shù)的增加，整體呈現(xiàn)上升的趨勢。通過圖5的對比可知，采用ARIMA預測模型預測的結果誤差值比AR預測模型的結果誤差值小，且誤差相對穩(wěn)定。所以采用ARIMA預測模型對LMT缺失的30個數(shù)據(jù)進行預測填補的方法是可行，且結果較好。

為了進一步驗證ARIMA模型預測數(shù)據(jù)的準確定，反映經過ARIMA模型對缺失數(shù)據(jù)進行填補后的數(shù)據(jù)在大地電磁測深資料處理中的影響，分別對實測原始數(shù)據(jù)和經過ARIMA模型預測填補后的數(shù)據(jù)進行頻譜分析。原始數(shù)據(jù)前后平均選取共計2 048(211)個數(shù)據(jù)，預測數(shù)據(jù)根據(jù)原始數(shù)據(jù)選取的原則在相同位置選取相同個數(shù)的數(shù)據(jù)，進行傅里葉變換，將時域數(shù)據(jù)變換到頻率域，則得到預測填補后的數(shù)據(jù)與實測數(shù)據(jù)的頻譜對比(圖6)。

圖6 預測填補后的數(shù)據(jù)與實測數(shù)據(jù)頻譜對比Fig.6 Comparision of the spectrum of observed data and the predicted data

由圖6可以看出，紅色實線和藍色實線基本重合，表明經ARIMA模型預測填補后的數(shù)據(jù)頻譜與原始數(shù)據(jù)的頻譜差值很小，基本可以忽略，故經預測填補后的數(shù)據(jù)對整體數(shù)據(jù)頻譜基本不會產生影響。所以采用ARIMA模型與缺失數(shù)據(jù)進行預測和填補是可行的。對于強干擾數(shù)據(jù)(如圖1中的電場跳點)，就可以直接將干擾數(shù)據(jù)定義為缺失數(shù)據(jù)進行處理，采用預測數(shù)據(jù)進行置換，但是當野外采集到的數(shù)據(jù)出現(xiàn)長時間的數(shù)據(jù)缺失現(xiàn)象或者干擾數(shù)據(jù)較多時，就需要重新對數(shù)據(jù)進行采集，以便得到更加準確的原始數(shù)據(jù)。

5 結論

筆者將ARIMA預測模型應用到LMT缺失數(shù)據(jù)的填補中，可以較好地還原野外數(shù)據(jù)采集過程中缺失的數(shù)據(jù)。實驗結果表明，ARIMA預測模型預測結果精度比AR模型高，隨著預測個數(shù)的增加，預測誤差值不會逐漸累加而增大，而是保持在一定范圍之內。ARIMA預測模型對原始LMT數(shù)據(jù)進行預處理，可以得到完整的野外觀測資料，置換由于強干擾導致數(shù)據(jù)出現(xiàn)的跳點，滿足大地電磁資料處理過程中對數(shù)據(jù)完整性的要求，提高了野外采集數(shù)據(jù)的利用率。因ARIMA預測模型對于缺失數(shù)據(jù)個數(shù)較多的數(shù)據(jù)其預測結果誤差較大，對于出現(xiàn)長時間的數(shù)據(jù)缺失或者干擾數(shù)據(jù)較多時，建議對該測點進行從新測量，以獲得更準確的原始資料。

[1] 方勝.如何保證大地電磁測深原始資料的質量[J].地球科學,1990(S1):14+80.

FANG S. The way to guarantee the quality of magnetotelluric sounding raw data[J]. Earth Science,1990(S1):14+80.(In Chinese)

[2] BOX G E P,JENKINS G.Time series analysis forecasting and control[J].Journal of the Operational Research Society,1976,22(2):199-201.

[3] 許立平,羅明志.基于ARIMA模型的黃金價格短期分析預測[J].財經科學, 2011(1):26-34.

XU L P, LUO M Z. Short-term analysis and prediction of gold price based on ARIMA model[J].Finance&Economics, 2011(1):26-34. (In Chinese)

[4] 張勃,劉秀麗.基于ARIMA模型的生態(tài)足跡動態(tài)模擬和預測——以甘肅省為例[J].生態(tài)學報,2011,31(20):6251-6260.

ZHANG B, LIU X L. Dynamic ecological footprint simulation and prediction based on ARIMA model: a case study of Gansu province, China[J]. Acta Ecologica Sinica, 2011,31(20):6251-6260. (In Chinese)

[5] 唐恒專.AR模型預測誤差比在地震事件識別中的應用[J].數(shù)據(jù)采集與處理,2004,19(4):446-449.

TANG H Z. Predication error ratio of AR model and its application to seismic event discrimination[J]. Journal of Data Acquisition & Processing,2004,19(4):446-449. (In Chinese)

[6] TRAVASSOS J M,BEAMISH D.Magnetotelluric data processing—a case study[J].Geophysical Journal International,1988, 93(2):377-391.

[7] 田紹耕.遠參考法大地電磁測深[J].石油地球物理勘探,1983,18(3):232-235.

TIAN S G. Magnetotelluric using remote reference method[J]. Oil Geophysical Prospecting, 1983,18(3):232-235. (In Chinese)

[8] SUTARNO D, VOZOFF K. Robust M-estimation of magnetotelluric impedance tensors[J]. Exploration Geophysics, 1989, 20(20):383-398.

[9] 湯井田,化希瑞,曹哲民,等.Hilbert-Huang變換與大地電磁噪聲壓制[J].地球物理學報,2008,51(2):603-610.

TANG J T,HUA X R,GAO Z M,et al.Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data[J].Chinese J.Geophys,2008,51(2):603-610. (In Chinese)

[10] 景建恩,魏文博,陳海燕,等.基于廣義S變換的大地電磁測深數(shù)據(jù)處理[J].地球物理學報,2012,55(12):4015-4022.

JING J E, WEI W B, CHEN H Y, et al. Magnetotelluric sounding data processing based on generalized S transformation[J]. Chinese J. Geophys,2012,55(12):4015-4022. (In Chinese)

[11] 張偉.長周期大地電磁法信號處理關鍵技術研究與應用[D].成都:成都理工大學,2011.

ZHANG W. Research and application in key single processing technology of Long-period Magnetotelluric[D]. Chengdu:Chengdu university of technology,2011. (In Chinese)

[12] 張春虹,張剛.AR模型在大地電磁測深資料處理中的應用[J].地球物理學進展,2013,28(3):1227-1233.

ZHANG C H, ZHANG G. The application of AR model in magnetotelluric data processing[J]. Progress in Geophys., 2013,28(3):1227-1233. (In Chinese)

[13] 張賢達. 現(xiàn)代信號處理:第二版[M]. 北京:清華大學出版社, 2002.

ZHANG X D. Modern signal processing: second edition [M]. Bejing: Tsinghua University Press,2002. (In Chinese)

[14] 王書明, 王家映. 大地電磁信號統(tǒng)計特征分析[J]. 地震學報, 2004, 26(6):669-674.

WANG S M, WANG J Y. Analysis on statistic characteristics of magnetotelluric signal[J]. Acta Seismologica Sinica, 2004, 26(6):669-674. (In Chinese)

[15] DICKEY D A, FULLER W A. Distribution of the estimators for autoregressive time series with a unit root[J]. Journal of the American Statistical Association, 1979, 74(366a):427-431.

[16] AKAIKE H. A new look at the statistical model identification[J]. Automatic Control IEEE Transactions on, 1974, 19(6):716-723.

[17] AKAIKE H. Fitting autoregressive models for prediction[J]. Annals of the Institute of Statistical Mathematics, 1969, 21(1):243-247.

[18] SCHWARZ G. Estimating the dimension of a model[J]. Annals of Statistics, 1978, 6(2):15-18.

[19] KANJILAL P P, PALIT S, SAHA G. Fetal ECG extraction from single-channel maternal ECG using singular value decomposition[J]. IEEE transactions on bio-medical engineering, 1997, 44(1):51-9.

[20] 何正嘉,訾艷陽,張西寧.現(xiàn)代信號處理及工程應用[M].西安:西安交通大學出版社,2007.

HE Z J, ZI Y Y, ZHANG X N. Modern signal processing and application in engineering[M]. Xian: Xi'an Jiaotong University Press,2007. (In Chinese)

TheapplicationofARIMAmodelinLMTdataprocessing

DENG Fangjin, WANG Xuben, LI Dewei

(Chengdu University of Technology,Chengdu 610059， China)

To solve the deletion or skipping of LMT time series which is is large, nonlinear, non-stationary based on observed data, the ARIMA method is applied to predict the missing data and replace the jump-point. The model order number is obtained through the stationarity test and the bayesian information criteria(BIC), the model parameters is determined using the LS theory. Establishing the bi-directional prediction model and linear combination method to forecast the missing data, then compare the predicted data accuracy of ARMA and AR. The result shows that the ARMA method with higher accuracy and without error accumulation, can solve the problem that the observed data is discontinuity and strong interference.

ARIMA model; predict; magnetotelluric; data processing; missing data

P 631.2

A

10.3969/j.issn.1001-1749.2017.05.05

2017-02-17 改回日期： 2017-03-31

“深地資源勘查開采”重點專項(2106YEC0600302);國家高技術研究發(fā)計劃(2014AA06A612)

鄧方進(1992-)，男，碩士，研究方向為固體地球物理學，E-mail：fangjindeng@126.com 。

1001-1749(2017)05-0612-09