BIC定阶球谐函数-武文俊

（3）
[2]
通过模型参数的逆函数估计法可以确定出 ARMA（p，q）模型的自回归系数和滑动平均系数 。下 面确定自回归阶数 p 和滑动平均阶数 q。 自回归 AR ( p ) 模型的参数又满足 Yule-Walker 方程 ，即有：
∗
[3]
⎡ φ1 ⎤ ⎡ γ 0 ⎢φ ⎥ ⎢ γ ⎢ 2⎥ ⎢ 1 ⎢ ... ⎥ = ⎢ ... ⎢φ ⎥ ⎢ ⎣ p ⎦ ⎣γ p −1
[2]
2期
武文俊等：利用时间序列模型预报电离层 TEC
143
3
数据的比较与分析
IGS 每天根据分布在全球、全天候观测得到的双频 GPS 资料归算出两种类型的 TEC 结果。一是每
天的快速解（推迟一天） ，公布在文件名为 igrddd0.yyi 的文件中；另一个是推迟两个星期公布的每天的 综合解，即综合了若干个 TEC 归算中心的结果而得到，公布在 igsgddd0.yyi 文件中（ddd 表示年积日， yy 表示简化年数，i 表示电离层） 。尽管有通过火箭等其他技术手段直接获得的电离层资料，但其时空覆 盖率都不能和 GPS 观测资料相比。因此 IGS 的 TEC 的结果是检验各种 TEC 预报模型预报效果的较好的 客观标准。 取 2007 年 1 月 1 日至 2007 年 5 月 31 日综合电离层资料为时间序列样本（实际工作中可以取快速 解中的电离层资料作为样本，它更新更快，更具有工程应用价值） ，用 ARMA（p，q）模型中 BIC 定阶 向后预报 15 天（即 2007 年 6 月 1 日至 2007 年 6 月 15 日）的 TEC 值。同时，利用 IRI2007 计算相同 时间、相同经纬度的电离层 TEC。IRI2007 计算电子密度剖面时把电离层分成 6 个区域，即顶部、 F2 层、
摘要：以 IGS（international GPS service）发布的电离层 TEC（total electron content） 资料为样本， 用时间序列模型对全球的电离层总电子含量进行了预报。 在时间序列预报模型中， 不同的定阶方法导致不同的预报结果；实践证明本文使用的 BIC 定阶准则较好地实现了电离层 总电子含量的预报。结果表明：对 10 d 左右的预报时间段，时间序列模型的 TEC 计算结果相对 精度高，预报相对精度优于 60%的网格点数在总网格点数中所占百分比可达 90%以上。 关 键 词： 时间序列；电离层总电子含量；预报; IGS（international GPS service） 文章编号：1674-0637(2008)02-0141-06
P = I y − II / II
（8）
其中：I y 为 ARMA 或 IRI 预报的电离层 TEC，I I 为 IGS 发布的 TEC。 将它们每天相对精度分别大于 90%， 大于 80%而小于 90%，大于 70%而小于 80%，大于 60%而小于 70%以及大于等于 60%的网格点数在总 网格点数中所占比例列于表 1（其中 100%~90%表示它们相对精度大于 90%而小于 100%的网格点数在 总网格点数中的百分比；余类推） 。
中图分类号： P352 文献标识码： A
1
引
言
在卫星单向授时、GPS 共视时间比对等技术中，信号会受到电离层的影响。信号传播过程中，其时 延取决于电离层总电子含量（TEC） 。因此，提高电离层 TEC 的预报精度，有利于提高单向授时和实时 时间比对的精度。目前国内外预报电离层 TEC 的方法主要有两类：一类是根据经验模式给出的电离层 电子浓度剖面，沿高度积分得到所需要的 TEC。例如，国际上著名的 IRI 模型是被广泛采用的经验模型， 该模型描述了无极光电离层在地磁宁静条件下的特定时间、特定地点上空 50 km 以上范围的电子密度、 电子温度、离子温度、离子成分、电子含量等月平均值。另一类方法是利用 TEC 观测数据直接建立 TEC 经验模式。例如，有些学者利用新乡站 1982—1989 年间测得的 TEC 数据，应用 TEC 与太阳黑子数滑动 均值的线性关系，建立了 TEC 的经验模式；有的研究者利用经验正交函数构造了武汉地区电子浓度总含 量的经验模式。这些经验模式主要是基于单站的 TEC 资料，利用 TEC 与太阳活动的相关性进行模拟的 。 时间序列分析是应用概率统计方法对一串随时间变化的随机数字序列进行分析的数学方法。 它主要 分时域分析和频域分析两方面，通过时域分析可以对时间序列进行预报，而通过频域分析可以掌握时间 序列的频域特性。本文使用时间序列的时域分析部分，利用 IGS（international GPS service）发布的全球 电离层 TEC 结果实现了全球未来十几天时间段内的 TEC 预报。
总第 31 卷 第 2 期 2008 年 12 月
时间频率学报
Journal of Time and Frequency
Vol.31 No.2 Dec., 2008
利用时间序列模型预报电离层 TEC1
武文俊 ，李志刚 ，杨旭海 ，程宗颐 ，王晓晗
1,3 1 1 1,2 4
（1. 中国科学院国家授时中心，陕西 西安 710600；2. 中国科学院上海天文台，上海 200030； 3. 中国科学院研究生院，北京 100039；4. 中国人民解放军 61541 部队，北京 100094 ）
*
γ1 γ0
...
... ... ... ...
γ p−2
1 N
N −k t =1
γ p−2 ⎥ ⎥ ... ⎥ ⎥ γ0 ⎦
γ p −1 ⎤
−1
⎡γ1 ⎤ ⎢γ ⎥ ⎢ 2⎥ ⎢ ... ⎥ ⎢γ ⎥ ⎣ p⎦
（4）
其中理论自协方差函数 γ k 为：
γk =
∑ xt + k xt
∗
p
（5）
（5）式中 N 为样本总数，通过上述关系可以求得 AR ( p ) 模型的白噪声方差估计
σ a = γ 0 − ∑ϕ jγ j
j =1
∧ 2
（6）
假定已经有阶数 p 的上界 p，假设 AP 模型的阶数是 k 时，引入 BIC 函数 BIC (k ) = ln σ a +
∗
∗
∧ 2
∗
k ln N N
∗ ∗
（7）
将 BIC ( k ) 在（0，1，2，…，p）中的第一个最小值点称为 AR ( p ) 的阶数。这样就确定了 AR ( p ) 模型的 阶数 p 。让 p +1 作为 ARMA（p，q）模型中的 p，1 作为 ARMA（p，q）模型中的阶数 q。 定阶以后采用平稳线性最小方差预报中的 ARMA（p，q）模型预报公式进行预报 。
142
时间频率学报
o o Байду номын сангаас o
总 31 卷
o
预报。 IGS 提供的电离层 TEC 资料是按全球每 2 h、 经度方向 （0 ～360 ） 每5 、 纬度方向 （87.5 ～-87.5 ） 每 2.5 的网格点给出的（即每 2 h 一幅全球电离层图，每图 71×73 = 5 183 个网格点） 。进行 ARMA（p， q）模型预报时，首先对电离层球壳上 5 183 个 TEC 值用下面的球谐函数拟合：
[1]
2
用 ARMA(p,q)模型预报电离层球面系数
利用 IGS 提供的电离层 TEC 资料，使用时间序列 ARMA（p，q）模型对未来时空点上的 TEC 进行
收稿日期：2007-12-26；修回日期：2008-06-13 基金项目：863 课题资助项目（2006AA12Z322） ；国家 973 计划资助项目（2007CB815503） 作者简介：武文俊，男，硕士，主要从事卫星轨道的定轨等方面的研究。
TEC (α , β ) =
max_ n n
o
n =0 m =0
∑ ∑ Pnm (sin α )(anm cos(mβ ) + bnm sin(mβ ))
（1）
其中， TEC (α , β ) 是对纬度 α 、经度 β 而言的总电子含量，max_n 是拟合时采用的最大阶数， Pnm (α ) 是归一化连带勒让德多项式， anm 及 bnm 则是对 5 183 个 TEC (α , β ) 值用最小二乘法按（1）式拟合得到的 球谐系数。不同时刻有不同的 anm 及 bnm 值，这样就得到了 anm 及 bnm 的时间序列。反之，如果已知某时 刻的 anm 及 bnm ，根据（1）式可计算该时刻对纬度 α、经度 β 而言的总电子含量。根据已有的 anm 及 bnm ， 用时间序列预报未来时刻的 anm 及 bnm ，算出相应时刻相应经、纬度的 TEC，实现电离层 TEC 的预报。 根据所用时间序列物理背景，下面选用 ARMA（p，q）模型。 设 {xi } 为上述球谐系数所构成的时间序列，且满足如下 ARMA（p，q）模型方程：
xt = ϕ1 xt −1 + ϕ2 xt − 2 + Lϕ p xt − p + ε t
（2）式中， ϕ 1 ， ϕ 2 ， L ， ϕ p 是 ARMA（p，q）模型的自回归系数， ε t 表达为：
（2）
ε t = at − θ1at −1 − L − θ q at − q
（3）式中， {at } 为白噪声序列， θ1 ， θ 2 ， L ， θ q 为 ARMA（p，q）模型的滑动平均系数。
表1 年 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ARMA 模型与 IRI 模型 TEC 预报各种相对精度的网格点数在总网格点数中的百分比 月 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 日 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 100%~90% ARMA 0.60 0.56 0.50 0.47 0.48 0.51 0.44 0.28 0.28 0.42 0.50 0.49 0.46 0.30 0.25 IRI 0.17 0.17 0.17 0.19 0.18 0.18 0.17 0.16 0.16 0.18 0.19 0.18 0.18 0.20 0.28 90%~80% ARMA 0.26 0.27 0.32 0.32 0.32 0.30 0.34 0.35 0.28 0.26 0.25 0.25 0.24 0.25 0.28 IRI 0.19 0.19 0.20 0.19 0.18 0.17 0.16 0.15 0.15 0.18 0.19 0.20 0.19 0.23 0.25 80%~70% ARMA 0.08 0.08 0.11 0.12 0.12 0.11 0.13 0.21 0.23 0.16 0.11 0.12 0.12 0.18 0.21 IRI 0.24 0.22 0.23 0.24 0.23 0.21 0.19 0.17 0.22 0.25 0.26 0.25 0.26 0.22 0.18 70%~60% ARMA 0.03 0.04 0.04 0.05 0.04 0.04 0.04 0.06 0.11 0.07 0.06 0.06 0.06 0.10 0.10 IRI 0.21 0.20 0.20 0.20 0.20 0.21 0.21 0.23 0.22 0.19 0.18 0.17 0.17 0.13 0.10 100%~60% ARMA 0.97 0.95 0.97 0.96 0.96 0.96 0.95 0.91 0.89 0.90 0.92 0.92 0.89 0.84 0.84 IRI 0.81 0.80 0.80 0.82 0.79 0.77 0.73 0.71 0.75 0.80 0.82 0.80 0.81 0.78 0.82