地球物理学进展  2017, Vol. 32 Issue (5): 1880-1885   PDF    
利用卡尔曼滤波技术提取地磁规则日变化
陈果1,2,3, 杜爱民1,2,3, 张莹1,2,3, 赵旭东4     
1. 中国科学院地质与地球物理研究所, 地球与行星物理重点实验室, 北京 100029
2. 中国科学院地球科学研究院, 北京 100029
3. 中国科学院大学, 北京 100049
4. 中国地震局地球物理研究所, 北京 100081
摘要:地磁规则日变化的有效提取能够提高地磁活动指数现报的准确性.本文基于1998—2013年北京十三陵(BMT)台站的地磁数据,利用卡尔曼滤波(KF)方法提取了地磁场规则日变化(SR)曲线.该方法提取的SR曲线具有明显的逐日变化特征,能够准确反映SR曲线的季节变化和地方时效应.根据KF方法得到的SR日变幅度逐日变化的结果,改进了K指数的量算方法,重新计算了十三陵台站的K指数以及Ak指数,并分析了SR日变幅度的逐日变化对K指数量算的改进程度.分析结果表明,当K≤3时,KF方法可以较好地改善SR曲线扰动幅度的计算,提高了地磁规则日变化提取的有效性,从而改进相应的地磁活动指数的量算.
关键词卡尔曼滤波    地磁规则日变化    地磁活动指数    
Kalman filter technique for defining solar regular geomagnetic variations
CHEN Guo1,2,3 , DU Ai-min1,2,3 , ZHANG Ying1,2,3 , ZHAO Xu-dong4     
1. Key Laboratory of Earth and Planetary Physics, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
2. Institutions of Earth Science, Chinese Academy of Sciences, Beijing 100029, China
3. University of Chinese Academy of Sciences, Beijing 100049, China
4. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
Abstract: Defining Solar Regular(SR) geomagnetic variance can improve the accuracy of geomagnetic indices estimation. Based on geomagnetic data from 1998 to 2013 at Beijing Ming Tombs Observatory (BMT) of China, Kalman filter technique is developed to define solar regular geomagnetic variation. The new algorithm can distinguish day-to-day (DTD) variations of geomagnetic field's regular variation and reflect both the diurnal and seasonal variations of magnetic disturbance. Based on DTD variations in the amplitude of SR curve obtained by Kalman filter, we improved the evaluation of K index, recalculated the K and Ak index at BMT, and evaluated the improvements of K index estimation after taking account of DTD variations in SR curve. The results showed that Kalman filter technique can improve the amplitude estimation of SR when K ≤ 3, then increase the accuracy of SR definition, which consequently improve the estimation of geomagnetic activity indices.
Key words: Kalman filter     solar regular geomagnetic variations     geomagnetic activity indices    
0 引言

地磁场可分为内源场和外源场两大部分,内源场起源于地表以下的磁性物质和电流,进一步分为主磁场、岩石圈磁场和感应场三部分,其中主磁场和岩石圈磁场变化较为缓慢,称为稳定磁场.外源场起源于地表以上的空间电流体系,这些电流体系随时间变化较快,所以外源磁场通常称为变化(瞬变)磁场(徐文耀等,2014).变化磁场又可分为规则变化和不规则变化.地磁场的规则变化主要包括太阳静日变化(Sq),太阴静日变化(L)和磁扰过后恢复期环电流等磁层电流缓慢变化(Dma)等成分(Mayaud, 1980).地磁场的不规则变化俗称扰动变化,简称磁扰.

地磁指数是表征地磁活动水平的重要指标,也是研究日地耦合过程的关键参数.地磁活动性关心的是不规则变化的磁扰活动特点.在计算地磁活动指数时,首先要消除的是地磁场的规则日变化(Solar Regular (SR) variations),残差部分当作地磁扰动的变化,然后按照一定的时间间隔和事先由统计对比分析得到的数值标尺换算成相应的地磁指数.因此,地磁规则日变化曲线的提取是计算地磁活动指数的最重要的一步.

常规的提取地磁场SR曲线的方法很简单:采用一个月内5个国际磁静日(IQD)的平均值作为该月每一天的标准静日变化(Chapman et al., 1940).然而,SR具有逐日变化特性,主要归因于太阳静日变化Sq的逐日变化,它是由电离层和磁层的几个电流体系引起的(Crooker and Siscoe, 1971; Kawasaki and Akasufu, 1971; Blanc and Richmond, 1980; Olson, 1984; Iyemori, 1990; McPherron, 1991; Iyemori and Rao, 1996).在地磁场平静日期间,Sq曲线的日变幅度和形态都存在不可忽视的逐日变化特性,有时,连续几个磁静日的日变幅可相差一倍以上(Xu, 2008).如果要真正认识地磁场的扰动变化,首先必须准确计算地磁场规则日变化,而规则变化的逐日变易性所带入的误差,在中小扰动时尤为明显.

目前,识别SR曲线逐日变化的方法可分为统计回归法、频率滤波法、正交基分解法等(Menvielle et al., 1995; Xu and Kamide, 2004),应用较为广泛的主要有:芬兰气象研究所研制的Finnish Meteorological Institute(FMI)方法(Pirjola et al., 1990)和Takahashi方法(Takahashi et al.,2001).Menvielle等(1995) 评估了多种计算机识别SR曲线的方法,结果表明FMI方法与手工标定的效果最为接近.FMI方法是一种变窗口回归方法,在计算当日的SR曲线时,需要衡量这段时间前后的磁场变化.该方法由每小时前后(m+n)分钟的数据计算该小时的平均值,其中,m与地方时有关,n与地磁活动性有关(Pirjola et al., 1990).由于考虑了每日地磁场实际变化的特点,FMI方法可以很好的识别SR曲线的逐日变化.但该方法需要使用每日前后时间段的数据进行均值计算,所以无法用于地磁活动性指数的现报.Takahashi方法采用了一种时序叠加的方法提取SR曲线:计算当日前10天地磁实测数据各小时的中值,再对该日二十四个数据点进行傅立叶变换,滤去5次以上谐波,采用样条函数插值成分钟值,即可得到该日的SR曲线.该方法相比Chapman的常规方法的优点是不需要人工标定磁静日,可用于地磁活动性指数的现报,但缺点是无法较好提取具有逐日变化的SR曲线(Wang et al., 2016).

本文基于1998—2013年北京十三陵(BMT)台站的地磁数据,利用卡尔曼滤波(KF)方法提取了地磁场规则日变化(SR)曲线.通过与FMI方法,Takahashi方法以及常规IQD方法结果的对比,我们检验了KF方法提取SR变化曲线的效果.该方法可有效识别地磁场规则变化的逐日变化特性,反映SR变化曲线的季节变化和地方时效应.根据KF方法得到的SR日变幅度逐日变化的结果,改进了K指数的量算方法,重新计算了K指数以及Ak指数,分析了SR日变幅度的逐日变化对K指数量算的改进程度.

1 卡尔曼滤波方法 1.1 卡尔曼滤波基本原理

卡尔曼滤波是一种递推线性最小方差估计(Kalman, 1960).它的主要思想是:首先给定预报模式,然后在模式的动力约束条件下根据当前k时刻模型的状态量预测k+1时刻的状态量.引入观测数据对k+1时刻的模式状态进行重新分析,在使观测和模式结果误差方差达到最小的条件下,得到k+1时刻状态量的最优估计.随着模式状态预报的持续进行和不断地将新的观测资料同化到动态模式系统中,使得预报过程可以不断向前推进(马建文等, 2013).

在卡尔曼滤波算法中,用状态方程和量测方程来描述系统.一般的线性离散系统可表示为

(1)
(2)

其中,式(1) 为状态方程,式(2) 为量测方程.XkRn×1k时刻的状态向量,FkRn×nk时刻到(k+1) 时刻的状态转移矩阵,WkRn×1k时刻的系统噪声向量;YkRm×1k时刻的观测向量,HkRm×nk时刻的预测输出转移矩阵,VkRm×1k时刻的观测噪声向量.WkVk是互不相关,均值为零,方差分别为QkRn×nVkRm×m的白噪声.

假定现有系统状态为Xk,没有使用该时刻观测Yk得到的状态估计值称为先验估计值,记为 (-表示先验,^表示估计);使用了该时刻观测Yk得到的状态估计值称为后验估计值,记为;先验估计和后验估计的误差协方差矩阵分别表示为Pk-Pk.

如果已知(k-1) 时刻的后验估计值及其误差协方差矩阵Pk-1,可以得到k时刻的先验估计值及其误差协方差矩阵Pk-,它们会满足的等式为

(3)
(4)

式(3)和式(4) 称为时间更新方程.

当得到k时刻的观测Yk后,通过使k时刻后验估计的误差协方差Pk最小,可以得到的公式为

(5)
(6)

其中Kk叫卡尔曼增益矩阵,它满足的表达式为

(7)

式(5)~(7) 称为量测更新方程.这样就得到了对随机线性离散系统k时刻状态的最优估计值及其误差协方差矩阵(Brown et al., 1992).

1.2 地磁规则日变化的提取

下面我们将利用卡尔曼滤波技术给出地磁规则日变化SR曲线的最优估计.这里采用的地磁数据为十三陵台站的分钟值数据,首先对每小时的分钟值数据计算得到该小时的平均值.

我们用状态向量Xk(ti)表示第k天的地磁场日变化,其中,ti表示小时,i=0, 1, 2, …23,观测向量Yk(ti)为地磁H分量的实际测量值.由于SR曲线为地磁场的规则变化,那么相邻两天的值差别应该不大,因此我们将(1) 式中的状态转移矩阵Fk-1设为单位矩阵I.由于H分量是地磁场日变化的实际测量值,我们将(2) 式中的预测输出转移矩阵Hk也设为单位矩阵I.地磁场规则日变化的状态演化方程和量测方程分别表示成式(8)和式(9) 为

(8)
(9)

其中,系统噪声向量W和观测噪声向量Vk满足高斯正态分布,公式为

(10)
(11)

这里σ经验地取值为10 nT,该值是对磁静日期间地磁规则日变曲线模型值的方差进行统计分析后给定的;Σk是第k天观测向量的误差协方差矩阵.

首先给定模式开始时刻的地磁规则日变曲线的初始估计值X0X0满足的条件为

(12)

其中,X0是高斯随机向量,期望值E(X0)为开始时刻的实际测量值Y0,协方差矩阵C(X0)为单位矩阵I.在完成初始条件的设定后,即可采用卡尔曼滤波算法使模式向前运行,具体步骤如下:

1) 利用地磁场规则日变化的时间演化模型,根据(3)(4) 式计算得到状态量Xk的先验估计值及其误差协方差矩阵Pk-.

2) 逐点计算观测数据Yk和先验估计值的差值,公式为

(13)

为了消除地磁场不规则扰动的影响,我们设定经验阈值为30 nT,舍弃超过阈值的数据点,且被舍弃的点不参与观测噪声向量的误差协方差矩阵Σk的计算.阈值的确定不仅与台站数据的逐日变化有关,而且还与所用台站的纬度有关.对于不同的台站,需经验地选取不同的阈值大小.

3) 计算第k天观测噪声向量的误差协方差矩阵Σk:当k < 30天时,Σk=measvar*I;其中,measvar为观测向量的误差协方差,这里经验地取为100.当k≥30天时,首先从地磁场H分量的分钟值数据估计出的初始K指数值,即由每3个小时中的最大值和最小值的差值确定K指数;如果当天的最大K指数≤3,那么Σk=measvar*I,其中measvar取为10,该值是统计分析磁静日期间观测值的方差后给定的;如果当天的最大K指数>3,那么Σk通过使用前30天观测数据小时均值的样本方差计算得到.

4) 利用地磁场观测数据Yk,根据(5)~(7) 式计算状态变量的后验估计值及其误差协方差矩阵Pk.

5) 对做傅立叶分析,滤去5次以上谐波,然后再使用样条函数插值到分钟精度,这样就得到第k天的地磁规则日变化SR曲线.

6) 将后验估计值作为第k+1天的地磁场规则日变化时间演化模型运行的初始值.重复步骤1)~5),直到完成所有需要处理的磁场数据.

2 计算结果与分析

为了检验KF方法提取地磁场规则日变化SR曲线的准确性,我们挑选了2008年2月24日和25日连续两天磁静日来测试KF方法的效果.图 1a给出了连续两个磁静日十三陵台站地磁H分量随地方时的变化,可以看出每日正午附近的地磁场变化均非常平缓,但是日变化幅度和相位都有明显差异.图 1b给出了分别用FMI方法(红线)、Takahashi方法(灰色实线)和KF方法(蓝线)提取的SR曲线的情况.由于我们选择的是磁静日数据,那么地磁场规则日变化SR曲线应接近于观测值(黑色点线).从图中可以看出,FMI方法提取的SR曲线最接近于实测值,然而我们发现对连续两天SR时均值数据作五次谐波拟合时的所用到的外推方法,会在两天数据的连接处产生“突跳”.相比FMI方法,Takahashi方法无法识别SR的逐日变化.我们采用KF方法提取的SR曲线比Takahashi方法更接近于实际观测值,并且具有明显的逐日变化,在两天数据的连接处也避免了出现类似FMI方法的数据跃变的情况.

图 1 (a)2008年2月24日和25日北京十三陵台站地磁H分量的实测值; (b)FMI方法(红线)、Takahashi方法(灰色实线)和KF方法(蓝线)得到的SR曲线计算值随时间的变化;黑色点线为地磁H分量的实测值 Figure 1 (a)Observed H component data at Beijing Ming Tombs observatory (BMT) on 24 and 25 February, 2008; (b) SR curves obtained by FMI method (red line), Takahashi method (grey solid line) and KF method (blue line); observed magnetic H component data (black dotted line)

我们对1998—2013年共16年的十三陵台站地磁数据进行了分析,对比了KF方法与常规IQD方法对SR曲线的提取情况.图 2给出了1998—2013年KF方法(a)和IQD方法(b)提取的SR曲线的季节平均值随地方时的分布情况.从图中可以看出,不管是在日变幅度还是在形态上KF方法得到的结果都与IQD方法的基本一致,而且均具有明显的季节变化和地方时效应.

图 2 KF方法(a)和IQD方法(b)提取的SR曲线随地方时的变化 Figure 2 SR curves obtained by KF method(a) and IQD method(b) vary with local time

不同方法识别SR曲线逐日变化的效果不同,Menvielle等(1995) 指出FMI方法得到的SR曲线与实测值曲线贴合太紧,在磁扰时会混入规则变化以外的成分,这样就会过高估计地磁规则日变化,导致地磁活动指数被低估.常规IQD方法以及Takahashi方法由于不能体现明显的逐日变化,在计算地磁活动指数时,不能完全消除地磁规则日变化.在大的磁暴期间,这些残留的地磁规则日变化远小于扰动变化,因而没有大的影响,但是在中小扰动,这些残留的地磁规则日变化就不可忽略.

为了分析SR日变幅度的逐日变化对K指数量算的改进程度,我们根据由KF方法得到SR日变幅度逐日变化的结果,对十三陵地磁台的K指数重新进行了计算(记作KF-K指数),并与常规IQD方法得到的K指数值(记作IQD-K指数)进行了对比.图 3给出了两种方法结果的对比,由上至下依次为:十三陵台站2001年7月H分量的分钟值曲线;KF方法得到的SR变化;常规IQD方法得到的SR变化;KF-K指数值;IQD-K指数值;两种方法得到的K指数之差.对2001年7月的统计结果表明,对于这个月248个地磁活动K指数,两种方法得到的K指数差值为0的情况有153个,占61.7%,±1的情况有90个,占36.3%,±2的情况有5个,占2%.两种方法计算得到的K指数约有1/3存在差异,因此考虑SR曲线的逐日变化是有必要的.

图 3 (a)北京十三陵台站2001年7月H分量的分钟值曲线;(b)KF方法得到的SR变化;(c)常规IQD方法得到的SR变化;(d)KF-K指数值;(e)IQD-K指数值;(f)两种方法得到的K指数之差 Figure 3 (a)Observed magnetic H component data at BMT station in July, 2001; (b)SR curve obtained by KF method; (c)SR curve obtained by regular IQD method; (d)K-index values obtained by KF method; (e)K-index values obtained by IQD method; (f)Differences between K-index values obtained by the two methods above

另外,我们将2003年十三陵地磁台站的KF-K值,与常规方法计算的IQD-K值的频次分布进行比较.从图 4a给出的十三陵台站K值的频次分布对比可看出二者具有相似的分布,为能更清晰地分辨出二者的差异,图 4b给出了两种方法计算得到的K值之差的频次分布.可以看出,|ΔK|的发生频次呈正态分布,且|ΔK|≤1的情况几乎达到98.79%.这说明了根据KF方法得到SR曲线的结果改进K指数的量算方法可以相当准确地反映K指数量值.

图 4 (a)十三陵台站2003年IQD-K值(灰色柱状)与KF-K值的频次分布(白色柱状);(b)基于两种方法计算的K值之差的频次分布 Figure 4 (a) Occurrence of K index of BMT for year 2003 by IQD method (grey bar) and KF method (white bar); (b) Percentage distribution of ΔK of year 2003

K指数与对应3小时时段最小变幅近似为半对数关系,使得对K指数求和以及求日均值在数学上变得毫无意义,不便于统计分析和定量计算.因此,为了进行统计分析,这里我们将K指数值转换为线性振幅值,计算了对应的三小时等效幅度ak指数,然后再计算一天8个ak指数的平均值得到每日等效幅度Ak指数.IQD-K与KF-K对应的Ak指数分别记作Ak_IQD与Ak_KF.

图 5ad分别给出了太阳活动上升期1999年,太阳活动较高年份2003年,太阳活动低年2008年,以及1998—2013年共16年的Ak_IQD与Ak_KF指数的散点分布图.可以看出,Ak_IQD与Ak_KF指数的相关性系数均在0.95以上.在太阳活动低年,改进的K指数量算方法对K指数的重新计算影响较大,从图 5c可以看出,重新计算的K指数在整体上高于用常规方法计算的K指数,对应的Ak指数的范围≤25, 相应的K指数的范围为≤3.因为对于常规方法,计算K指数的时候往往会残留一些地磁规则日变化的成分,尤其是对于K≤3的扰动,这些残留的成分会影响K指数的大小.这同时也说明了,改进的K指数的量算方法由于考虑了SR的逐日变化,在一定程度上可以消除常规方法难以去除的地磁规则日变化的残留成分,从而改善SR日变曲线扰动幅度的计算,提高了K指数的量算精度.

图 5 1999年(a),2003年(b),2008年(c),以及1998—2013年(d)的Ak_KF与Ak_IQD指数的散点分布 Figure 5 Scatter plots of Ak_KF index and Ak_IQD index. (a) year 1999; (b) year 2003; (c) year 2008; (d) years from 1998 to 2013
3 总结与讨论 3.1

地磁规则日变化的提取是目前地磁活动指数现报的误差来源之一,为了降低现报误差,提高地磁活动指数现报的准确性,我们采用KF方法提取了地磁场规则日变化.通过与FMI方法、Takahashi方法结果的对比,检验了KF方法的提取效果,结果表明KF方法提取的SR曲线具有明显的逐日变化特征,能够准确反映SR曲线的季节变化和地方时效应.

3.2

基于1998—2013年北京十三陵台站的数据,我们根据KF方法得到的SR日变幅度逐日变化的结果,改进了K指数的量算方法,重新计算了北京十三陵台站的K指数,并分析了SR日变幅度的逐日变化对K指数的量算的改进程度.在大的磁暴期间,地磁规则日变化远小于扰动变化,因而没有大的影响,但是在中小扰动,对地磁规则日变化的过高或过低估计就不可忽略了.通过分析KF-K值与IQD-K值之差的频次分布,可以看出改进的K指数的量算方法可以相当准确地反映K指数量值.

3.3

为了便于统计分析,我们将K指数值转换为线性振幅值,计算了相应的Ak指数,并对Ak_IQD与Ak_KF指数进行了相关性分析.分析结果表明,当K≤3时,KF方法可以较好地改善SR曲线扰动幅度的计算,提高了地磁规则日变化提取的有效性,从而改进相应的地磁活动指数的量算.本文提出的根据KF方法得到的SR日变幅度逐日变化的结果改进K指数的量算方法,也可以用来确定其他地磁活动指数,如Dst、AE等指数.

致谢 诚挚地感谢中国科学院地质与地球物理研究所区家明博士后和编辑部的大力支持.
参考文献
[] Blanc M, Richmond A D. 1980. The ionospheric disturbance dynamo[J]. J. Geophys. Res., 85(A4): 1669–1686. DOI:10.1029/JA085iA04p01669
[] Brown R G, Hwang P Y C. 1992. Introduction to random signals and applied Kalman filtering with Matlab[M]. second edition. John Wiley & Sons, Inc.
[] Chapman S, Bartels J. 1940. Geomagnetism[M]. New York: Oxford University Press.
[] Crooker N U, Siscoe G L. 1971. A study of the geomagnetic disturbance field asymmetry[J]. Radio Science, 6(4): 495–501. DOI:10.1029/RS006i004p00495
[] Iyemori T. 1990. Storm-time magnetospheric currents inferred from mid-latitude geomagnetic field variations[J]. J. Geomagn. Geoelectr., 42(11): 1249–1265. DOI:10.5636/jgg.42.1249
[] Iyemori T, Rao D R K. 1996. Decay of the Dst field of geomagnetic disturbance after substorm onset and its implication to storm-substorm relation[J]. Ann. Geophys., 14(6): 608–618. DOI:10.1007/s00585-996-0608-3
[] Kalman R E. 1960. A new approach to linear filtering and prediction problems[J]. J.Basic Eng., 82(1): 35–45. DOI:10.1115/1.3662552
[] Kawasaki K, Akasufu S I. 1971. Low-latitude DS component of geomagnetic storm field[J]. J.Geophys.Res., 76(10): 2396–2405. DOI:10.1029/JA076i010p02396
[] Ma J W. 2013. Research and Development of Data Assimilation Algorithm (in Chinese)[M]. Beijing: Science Press.
[] Mayaud P N. 1980. Derivation, Meaning, and Use of Geomagnetic Indices[M]. Geophysical Monograph 22. Washington, D.C:American Geophysical Union, 85-96.
[] McPherron R L. 1991. Physical processes producing magnetospheric substorms and magnetic storms[J]. Geomagnetism, 4: 593–739.
[] Menvielle M, Papitashvili N, Häkkinen L, et al. 1995. Computer production of K indices:Review and comparison of methods[J]. Geophys. J. Int., 123(3): 866–886. DOI:10.1111/gji.1995.123.issue-3
[] Olson W P. 1984. Introduction to the topology of magnetospheric current systems[A].//Potemra T A ed. Magnetospheric Currents[M]. Washington, D.C:American Geophysical Union, 28:49-66.
[] Pirjola R, Ryno J, Sucksdorff C. 1990. Computer production of K-indices by a simple method based on linear elimination[C]//Proceedings of the International Workshop on Geomagnetic Observatory Data Acquisition and Processing. Proceedings of the International Workshop on Geomagnetic Observatory Data Acquisition and Processing.
[] Takahashi K, Toth B A, Olson J V. 2001. An automated procedure for near-real-time Kp estimates[J]. J. Geophys. Res., 106(A10): 21017–21032. DOI:10.1029/2000JA000218
[] Wang G, Luo B X, Liu S Q, et al. 2016. An improved algorithm for nowcast of Kp index[J]. Chin. J. Space Sci. (in Chinese), 36(2): 153–166. DOI:10.11728/cjss2016.02.153
[] Xu W Y. 2008. Uncertainty in magnetic activity indices[J]. Sci. China Ser. E-Tech. Sci., 51(10): 1659–1664. DOI:10.1007/s11431-008-0261-z
[] Xu W Y. 2014. Introduction to Geomagnetic Activity (in Chinese)[M]. Beijing: Science Press.
[] Xu W Y, Kamide Y. 2004. Decomposition of daily geomagnetic variations by using method of natural orthogonal component[J]. J. Geophys. Res., 109(A5): A05218. DOI:10.1029/2003JA010216
[] 马建文. 2013. 数据同化算法研发与实验[M]. 北京: 科学出版社.
[] 王庚, 罗冰显, 刘四清, 等. 2016. 一种改进的Kp指数现报模式[J]. 空间科学学报, 36(2): 153–166. DOI:10.11728/cjss2016.02.153
[] 徐文耀. 2014. 地磁活动性概论[M]. 北京: 科学出版社.