于 畅，张 羽

（1.吉林大学 地球科学学院，吉林 长春 130061；2.长春合隆地震台，吉林 长春 130200；3.吉林省地震局，吉林 长春 130117）

0 引言

在工程实际中常出现随时间变化又相互关联的数据，一般称之为时间序列。在地震分析预报探索中，将地磁观测数据作为时间序列进行数学模型拟合，可以实现数据的建模、预报和系统控制[1]。研究地球磁场及其长期变化，有助于了解地球内部运动状况及地球的演化历史。利用局部地区地磁场随时间的变化异常可进行地震与火山等自然灾害的预测分析[2]。

长春地磁台位于吉林省长春市三岗县光荣村，处于长春市西北方向，距离市中心约60 公里。台基为白垩纪松花江群巨厚沉积地层，由无磁性和弱磁性泥岩、细粉砂岩和砂岩组成。下伏为侏罗纪火山岩类，构成盆地的磁性基地第四纪盖层厚约20 米。地面磁场水平梯度小于1nT/m。长春地磁台为吉林省地磁基准台，其FHDZ-M15 地磁相对记录仪记录的数据连续、准确地追踪地磁场的变化，为吉林省乃至全国地磁学科建设提供基础性服务。

MATLAB 是由美国MathWorks 公司出品的商业数学软件，在数据可视化、数据建模及数据分析等方面提供了全面的解决方案和完备的运行环境。它在数值计算领域首屈一指，在计算矩阵、绘制函数和实现算法等方面具有极快的运行效率。因而选择MATLAB 作为数学建模的运行平台能够更好地提高效率，实现模型可视化。

地磁数据F 分量记录的是地磁场强度，在无外界磁干扰的情况下地磁场应无大幅度变化，若F 分量数据在一定时间范围内出现陡升、骤降或其他明显变化，可对应地磁分析预报中的多种数学算法进行计算，试探引起F 值变化的原因，并探索可能产生的结果。

本文即以长春地磁台2007—2018 年地磁相对记录F 分量为样本，采用ARMA 模型为拟合手段，对F 分量进行处理与分析，旨在寻求能真实反映样本变化趋势的回归函数，为研究历史数据的相关性及预测未知数据的变化趋势提供铺垫，从而达到对异常数据范围的控制。

1 ARMA 模型的拟合原理

在应用数学中谱密度是用于表征信号处理的通用概念，它表示每赫兹的功率。计算数学模型的谱密度是研究平稳随机过程的典型方法。由于ARMA 过程的谱密度是有理函数，它能逼近任何连续谱密度，同时ARMA 过程具有线性结构，可用求一阶导数求得序列中某一时间段内数据的变化趋势，可用求二阶导数求得数据变化量的极值、最值及取值范围。这将会对地磁异常预报起到辅助作用，因此运用这一模型对地磁相对记录数据进行拟合分析具有可行性及重要意义。

ARMA 模型拟合分为以下几个步骤：

（1）数据的平稳性检验。二维图形中平稳的时间序列通常表现为数据值在均值附近上下浮动，而非持续上升或持续下降。如果时间序列呈现出明显的持续上升或持续下降的变化趋势，则需要使用一阶差分或n 阶差分的手段，直至变化趋势消失。一般地，一阶差分或二阶差分即可消除变化趋势。

（2）模型初步识别。由平稳时间序列的样本函数计算其自相关函数和偏自相关函数，从函数图像中可分别判断出自相关函数与偏自相关函数的拖尾性和截尾性，由函数的拖尾性和截尾性可以进行模型的初步估计和阶的初估计。

（3）模型的参数估计。确定模型类别后能够对应到1.2 中列出的三项公式中，运用MATLAB 软件中的AIC 函数准则法定阶并确定参数φ，θ，σ2的值。

（4）模型拟合优度检验。对初步建立的模型运用Q 值检验法进行合理性检验。Q 值检验法的通过条件是数据的Q 值小于置信表内置信

区间的值，若不成立，则需要从获得平稳性时间序列开始，重复上述三项拟合步骤，选取新的拟合模型，再次确定模型的阶数和参数，直至所得模型通过Q 值检验法[1]。

1.1 平稳性检验

数据的平稳性检验需要经过两个步骤，一是其数学期望和方差是否为常数；二是其自协方差函数是否只与时间间隔有关，与间隔端点的位置无关。

对于长度为N 的样本序列X1，X2，…，XN，可构造以下统计量：

1.2 模型初步识别与参数估计

对于平稳时间序列，可用样本自相关函数和偏自相关函数的统计特性识别模型。ARMA模型自相关和偏自相关函数的统计特性如表1所示[3]。

表1 自相关函数与偏自相关函数的统计特性

从自相关和偏自相关函数图像中可分别判断函数的截尾性和拖尾性，按照统计特性可以初步确定样本模型的类别。

ARMA 模型的三种表现形式如下：

（1）AR（p）模型

其中c 为常数，φt为自回归系数，随机误差εt是均值为0，方差为σ2的白噪声序列。

（2）MA（q）模型

其中c 为常数，θi为滑动平均系数，随机误差εt 是均值为0，方差为σ2的白噪声序列。

（3）ARMA（p，q）模型

其中c 为常数，φt和θi为自回归系数和滑动平均系数，随机误差εt是均值为0，方差为σ2的白噪声序列。

依托MATLAB 软件并运用AIC 准则法对模型定阶，探求使得AIC 值最小的参数p 和q。Box-Jenkins 定阶方法中提到实际应用中混合模型的阶数是很低的，所以在确定模型参数时，从低阶到高阶进行尝试是可以找到最佳参数的。结合MATLAB 中AIC 函数的使用方法，得到使得AIC 值最小的模型参数并不困难。

1.3 拟合优度检验

ARMA 模型的识别与参数估计需要反复设置、检验，通常有假设检验、Q 值检验、T 值检验等方法验证模型的拟合优度。本文运用Q值检验法检验模型的拟合优度。如果对数据求得的Q 值小于置信表中置信区间对应的值，则可验证拟合模型为有效模型。若大于置信表中对应的值，表明模型有效程度不高，需要重新判断模型参数，提高拟合优度。

2 应用实例

以长春地磁台2007—2018 年地磁相对记录F 分量分数据作为原始样本序列，绘制出原始序列曲线，如图1 所示。

可见，F 分量呈逐年上升趋势，为消除年变递增趋势，对原始序列进行一阶差分运算，即运用y（k）=x（k+1）-x（k）方法得出新序列。绘制一阶差分后的序列，由图2 可见一阶差分后的数据上下起伏均匀，已无明显年变趋势。

接下来需要按照1.1 中提出的平稳性检验方法验证一阶差分序列的平稳性。

图1 2007—2018 年长春地磁台地磁相对记录F 分量原始序列曲线Fig.1 Curve for element F in geomagnetic data of Changchun from 2007 to 2018

图2 2007—2018 年长春地磁台地磁相对记录F 分量一阶差分曲线Fig.2 Curve for 1st difference of F in geomagnetic data of Changchun from 2007 to 2018

其中x（t）为常数，x（t-τ）为t 的一元函数，令x（t-τ）=a（t-τ）+b，代入积分中可得R（τ）为τ 的一元函数。至此满足平稳性时间序列的所有条件，可表明此时的一阶差分序列为平稳性时间序列。

由图3 所示，一阶差分序列的自相关函数变化值逐渐接近0 又不等于0，故其具有拖尾性。从图4 中可以看出，偏自相关函数的变化趋势从平缓到最终全部为0，可以得出具有截尾性。

根据统计特性自相关函数具有拖尾性，偏自相关函数具有截尾性可初步判断序列为AR（p）模型。由序列的自相关和偏自相关函数的拖尾性和截尾性也可反证该序列为平稳性时间序列。

对AR（p）模型的序列使用AIC 准则法定阶，若将阶数限定在5 阶以内，则使得AIC 值最小的阶数为5 阶。

求解的主程序如下：

此时对模型的初步估计为AR（5）。运用最小二乘估计得到AR（5）模型的表达形式如下，Xt=-0.4169Xt-1-0.0243Xt-2+0.0112Xt-3+0.2032Xt-4+0.2351Xt-5+εt，t=6，…，n

主程序如下：

为验证模型的可信度，对残差序列进行检验是必要的。运用Q 值检验法，将一阶差分序列按递增顺序排列，取出最大值和最小值并求得极差可得m=60.3，选出序列中的可疑值并求得其与相邻数据差的绝对值可得n=4.1，最终得出Q=n/m=0.0667，规定所求函数的置信水平在95%以上，需满足所求Q 值小于Q 表中对应的置信区间的值，比较后得出所求Q 值远远小于Q 表中对应的值，故可认为在显著情况下残差序列为白噪声序列，所得模型符合拟合优度要求。可运用模型表达式进行数据的初步预测。

图3 F 分量一阶差分序列自相关函数图像Fig.3 Autocorrelation of first order difference for element F

图4 F 分量一阶差分序列偏自相关函数图像Fig.4 Partial autocorrelation of first order difference for element F

3 结论

本文以地磁相对记录F 分量预处理分数据为例建立ARMA 数学模型拟合方程。根据数据的统计特性计算预处理数据的拖尾性和截尾性并判断模型的类别，运用AIC 准则法确定模型阶数及模型的拟合参数，对方程进行Q 统计量检验，最终得出数据的拟合模型。

在模型定阶与参数估计的过程中需要对阶数由小到大尝试多次，若依赖计算机程序进行足够次数的迭代与优化，更方便得到精度较高的拟合方程。对高精度的拟合方程求一阶导数可以得到某个时间段内F 值的变化趋势，对骤变时间段内的F 值进行地磁异常分析，这是对原有地磁异常判定方法的补充。对二阶导数置零可以求得F 值在某一时间段内的极值，极值点对应的F 值为对应时间段内F 的最大或最小值，求得某一时间段内的最大值和最小值也同样是对地磁异常核实的补充，因而运用ARMA模型拟合的方法获得表征数据的拟合方程，这种方法为地震异常的预测和分析提供一种新思路。

在实际应用前对原始数据进行有效预处理是必要的，保证数据的连续性和准确性，去除原始数据中明显的线性趋势，这对得到平稳时间序列、算取更为精确和优化的拟合模型具有重要作用。