太阳黑子数时间序列的奇异谱分析和小波分析

第32卷第6期2007年11月测绘科学Science of Surveying and M app ingVol 132No 16Nov 1作者简介:徐克红(19822),女,山东泰安人,辽宁工程技术大学与中国测绘科学研究院联合培养硕士研究生,主要研究方向为卫星轨道确定。

E 2mail:xukehong0719@1631com 收稿日期:2007206228太阳黑子数时间序列的奇异谱分析和小波分析徐克红①②,程鹏飞①,文汉江①(①中国测绘科学研究院,北京　100039;②辽宁工程技术大学,辽宁阜新　123000)【摘　要】本文对小波变换和奇异谱分析方法进行了简要介绍,对离散小波的分解和重构、奇异谱分析的重构进行了详细阐述。

结合太阳黑子数1749年至2007年3月期间的月平均值时间序列进行了小波变换的分解和重构及SS A 方法的重构,提取了其主要的周期特性,并对两种分析方法进行了比较。

【关键词】小波分析;离散小波的分解与重构;奇异谱分析;太阳黑子数【中图分类号】P228 【文献标识码】A 【文章编号】100922307(2007)06200352041　引言太阳黑子是太阳光球上经常出现的阴暗斑点,是太阳活动的羁绊标志,是反映太阳辐射变化的重要指标,一般用太阳黑子数表示。

太阳黑子数反映了太阳活动强弱的变化,对地球的影响很大,诸如地磁变化、大气运动、气候异常、海洋变化等,都和太阳黑子数变化有着不同程度的关系。

因此研究太阳黑子数的变化有利于深入了解它对卫星轨道、定位等方面的影响。

对太阳黑子数变化的研究已有很多,韩延本,韩刚用小波分析的方法对太阳黑子数变化进行研究,验证了小波分析方法的可行性,并得到太阳黑子数变化包含多种周期分量的结论。

郝立生,李新,李月英利用Morlet 小波变换对太阳活动变化进行了研究,得到太阳活动存在141和106a 的变化周期。

小波变换的概念是1984年法国地球物理学家J 1Morlte 在分析处理地球物理勘探资料时提出来的。

其数学基础是19世纪的傅里叶变换,其后理论物理学家A 1Gr oss man 采用平移和伸缩不变性建立了小波变换的理论体系。

1989年S 1Mallat 提出了多分辨率分析概念,统一了在此之前的各种构造小波的方法,特别是提出了二进小波变换的快速算法,使得小波变换完全走向实用性[8]。

奇异谱分析(SS A )是对一维的时间序列进行分析的主成分分析方法。

该方法适用于从短噪声时间序列中提取信息。

SS A 在时空域中,通过将序列分解成元素行为模式的方法,将含在延迟坐标相空间的信息拆开,通过使用数据适应滤波器来帮助将时间序列分开为统计的独立成分,这些成分可以当作趋势、振动或噪声来进行分类。

本文选用太阳黑子数月平均值,采用小波变换和奇异谱分析的方法对该时间序列进行分析,同时对两种分析方法进行比较。

2　奇异谱分析主成分分析(PCA,Princi pal Component Analysis ),也称为经验正交函数(E OF,E mp irical O rthogonal Functi on ),可以由多维的时间序列中获取时间序列的主要成分,是常用的多元统计分析方法之一,主要将多个彼此相关的指标变换为少数几个彼此独立的综合指标即主成分,并要求主成分能反映原始数据的几乎全部信息,其中,常用于对一维的时间序列进行分析的方法称为奇异谱分析(SS A,Sin 2gular s pectru m analysis )。

奇异谱方法(SS A )是一种特别适合于研究周期振荡行为的分析方法,它是从时间序列的动力重构出发,并与经验正交函数相联系的一种统计技术,是E OF 分解的一特殊应用。

分解的空间结构与时间尺度密切相关,可以较好地从含噪声的有限尺度时间序列中提取信息,目前已应用于多种时间序列的分析中。

SS A 的具体操作过程是,将一个样本量为n 的时间序列按给定嵌套空间维数(即窗口长度)构造一资料矩阵。

当这一个资料矩阵计算出明显成对的特征值,且相应的E OF 几乎是周期性或正交时,通常就对应着信号中的振荡行为,可见SS A 在数学上相应于E OF 在延滞坐标上的表达。

对给定的X 1,X 2,…,X n 的时间序列,给定嵌套维数M ,M <N /2,建立时滞矩阵:X =x 1x 2…x N -M +1x 2x 3…x N -M +2…………x m x M +1…x N(1) (1)式的滞后自协方差是一个M ×M 的矩阵:S =s (0)s (1)…s (M -1)s (1)s (0)…s (M -2)…………s (M -1)s (M -2)…s (0)(2) S 为对称阵且主对角线为同一常数,称为Toep litz 矩阵,其特征值为:λ1≥λ2≥…≥λm ≥0,(3)λ1≥λ2≥…≥λm ≥0,(4) (4)式即为序列{x i }的奇异谱,称对X 的奇异值运算为奇异谱分析。

最大的特征值对应的特征向量看为第一阶模式,第二大的特征值对应的特征向量看为第二阶模式,依次类推。

第一阶模式代表了信号的最大变化趋势,第二阶模式代表了与第一阶模式无关的剩余信号量的最大变化趋势,依此类推。

在实际分析过程中,通常只选取前面的低阶模式进行分析。

计算出S 的特征值λk 和相应的特征向量,序列的SS A 展开为:x i+j =∑Mk =1aikE kj (5) 式中,i =1,2,…N 2M +1;j =1,2,…M 。

E kj 为时间E OF,记为T 2E OF,a ik 为时间主分量,记为T 2PC:a ik =∑Mj =1x i+j E kj ,(6) SS A 最重要的应用是重建成分RC (Reconstructi on )。

由第k 个的T 2E OF 和T 2PC 重建x i 的成分记为x k i ,公式为xk i=1i∑ij =1aij ,kE kj ,1≤i ≤M -1;1M∑Mj =1aij ,kE kj ,M ≤i ≤N -M +1;1N -i +1∑Mj =i-n +maij ,kE kj ,N -M +2≤i ≤N 1(7) 截取前k 个占比重大的成分近似表示原序列,x i =∑kk =1xik,i =1,2,…,N (8) 同时也可对各分量进行重建,用于对原信号的分析。

SS A 可以提取具有显著振荡行为的信号分量,并可选择若干有意义的分量进行序列重建。

其中低频信号的重建分量,显示了原始序列的主要演变特征。

3　小波分析311　小波分析小波分析是目前分析时间序列的有效工具,它可以获取时间序列的时间—频率特征,该分析方法是一种窗口大小(即窗口面积)固定但其形状可改变,时间窗和频率窗都可以改变的时频局域化分析方法,即在低频部分具有较高的频率分辨率和较低的时间分辨率,在高频部分具有较高的时间分辨率和较低的频率分辨率,所以被誉为数学显微镜。

正是这种特性,使小波变换具有对信号的自适应性。

将基本小波函数ψ(t )做位移τ后,再在不同尺度a 下与待分析信号x (t )做内积得到x (t )的小波变换:W T x (a,τ)=1a∫+∞-∞x (t )ψ3(t -τa)d t a >0(9) 式中,ψ3(t -τa )是ψ(t -τa)的逆变换。

小波变换将时间函数投影到二维的时间—尺度相平面上,提供了信号随时间变化的频谱特征,这对提取信号的某些本质特征创造了有利条件。

312　离散小波分解与重构离散小波变换定义(D iscrete W avelet Transf or m,简称DW T )为:W T f (a j 0,k τ0)=∫f (t )ψ3a j 0,k τ0(t )d t j =0,1,2,…,k ∈Z (10) 在多分辨率分析理论的基础上,1988年S 1Mallat 提出了离散小波分解与重构的快速算法,称为Mallat 算法[9]。

Mallat 分解算法使采样信号通过离散小波变换得到原始信号的低频系数和高频系数,然后保留高频部分,对低频系数进行再次的离散小波变换,由此一直下去,直至符合实际的要求。

重构算法正好与此过程相反。

见图1、图2。

其中,x 为原始信号,c N 为分解后的第N 层低频系数,d N为分解后的第N 层高频系数。

Lo_D 为小波分解低频滤波器,H i_D 为小波分解高频滤波器,Lo_R 为小波重构低频滤波器,H i_R 为小波重构高频滤波器。

4　太阳黑子数时间序列的奇异谱分析与小波分析人类对太阳的活动做了长期的观测,记录了太阳黑子数随时间变化的资料,这些资料记录了每天的、每月的和每年的太阳黑子数,而在实际应用中,它的月均值、月滑值和年均值才有较大的实用意义。

目前被认为较系统的太阳黑子相对数月值的时间序列起始于1749年,采用1749年22007年4月太阳黑子数时间序列,时间间隔为一个月,共3100个数据。

所有数据均来源于比利时皇家天文台(SI D C )[1]。

411　奇异谱分析(SSA)对太阳黑子数进行奇异谱分析,分别取窗口长度M 为100和150,单位为月。

首先,将M 取不同值时,所得的特征向量取前8个,将其变化曲线画出,见图3、图4。

由图可见,曲线呈现出十分明显的周期。

图3　成对特征向量(M =100)横坐标为:窗口长度(单位:月)　纵坐标为:特征向量图4　成对特征向量(M =150)图3和图4分别为窗口宽度取100和150时太阳黑子数SS A 的特征向量所表示的倾向或准周期信号的型态变化特征。

图3中a 图的1与2特征向量代表了倾向的变化特点,图3中b 、d 中各自两特征向量的形态变化均有些差别,其周期长度也有差别。

图4中的a 图同样代表了变化的倾向,而c 、d 中各自两特征向量的形态变化却很相似,周期长度也差别不大。

图3与图4相比,周期特性更明显,反映了1014a 的周期。

由于所取窗口长度度的限制,所以无法反映更长的周期特性。

然后,对前3阶的各分量进行了重构,并画出其对应的功率谱图,结果见图5和图6。

其中,M =100时,各模式所占比重为:第一阶模式,38166%;第二阶模式,63测绘科学 第32卷　31186%;第三阶模式,14134%;同理,取M =150时,所占比重为:第一阶模式,34186%;第二阶模式,28101%;第三阶模式,16176%。

图5　奇异谱分析(M =100)图6　奇异谱分析(M =150)由以上两图,各分量重构的结果所侧重细节差异很大,由所得功率谱图可见,所反映的主要周期是一致得,如两图中b 和d,均反应了1014a 的主周期,与准11a 的周期相一致。

两图中的f,与以上所反映的周期不同,为9216a 周期。

最后,分别对M =100,和M =150的前3阶模式进行重构,并与原始信号进行比较,并画出其功率普图,结果见图7和8。