CN108983321B - 一种基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法 - Google Patents

Abstract

本发明涉及一种基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法，属天文技术和时频处理领域。本发明采用采用结合同步压缩小波变换的方法对太阳黑子数和地磁Ap指数的周期分量进行提取，较精确地对太阳黑子数和地磁Ap指数的周期分量进行了提取，较好的解决了之前几种时频分析方法在在提取太阳黑子数和地磁Ap指数的周期分量结果较差，不完整的问题。

Description

一种基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数 的周期分量的方法

技术领域

本发明涉及一种基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法，属天文技术和时频处理领域。

背景技术

地磁活动指数周期变化的研究，对认识地磁活动的规律，开展空间天气预报起着重要作用。地磁Ap指数是全球的全日地磁扰动强度的指数，是描述全球地磁活动水平的重要指数。太阳黑子数是描述太阳活动水平高低的重要指数，在分析太阳活动周期是采用最多的数据就是太阳黑子数。太阳活动对地球磁场变化有着重要作用。通过对地磁活动指数和太阳活动参数的统计分析，有助于建立地磁与太阳活动的相应关系。该方法有助于更好研究的地磁活动和太阳活动之间关系、预测地磁活动和气候系统的变化。

描述全球地磁活动水平的重要指数Ap和描述太阳活动水平高低的重要指数太阳黑子数的变化特征和相关关系一直受到特别关注。我们发现这些方法主要包含两个过程：(1)对采样数据进行小波分析或者是傅里叶变换；(2)紧接着比较太阳黑子数和地磁Ap指数的相同的周期特性明显的分量。通过相关文献对这几种方法的介绍，可以看到这些方法对地磁活动和太阳活动之间关系的研究存在着以下的问题：(1)以Fourier变换作为理论基础，受到Heisenberg/Gabor不确定原理的影响，其时频分辨率受限；(2)时频分析得到的结果谱污染较为严重；(3)将噪声较为敏感，造成误检。然而，太阳黑子数和地磁Ap指数作为天文数据，天文望远镜获取过程中难免受到噪声或一些其他数据污染，具有明显的非平稳、时变特点。基于这些因素，依靠传统的时域和频域分析方法难以有效提取太阳黑子数和地磁Ap指数的周期特征。

发明内容

本发明提供了一种基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法，以用于提取太阳黑子数和地磁Ap指数的周期分量。

本发明的技术方案是：一种基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法，所述方法的具体步骤如下：

所述方法的具体步骤如下：

步骤1：对太阳黑子数数据和地磁Ap指数数据分别进行连续小波变换，得到小波系数，通过对小波系数求偏导得到瞬时频率；

步骤2：对瞬时频率进行同步压缩，得到同步压缩后的时间周期图；

步骤3：对同步压缩后的时间周期图中除COI区域以外的区域采用启发式方法连续的提取最大能量曲线，得到周期分量；

步骤4：对步骤3的结果进行滤波处理，得到滤波处理后的周期分量和其能量曲线；

步骤5：通过步骤4中能量曲线对周期分量进行重建，得到该周期分量的时间幅度图；

步骤6：对步骤5得到的太阳黑子数和地磁Ap指数周期分量的时间幅度图中的相同的周期分量的幅值进行归一化处理，得到归一化的周期分量的时间幅度图；

步骤7：对步骤5或者步骤6的时间幅度图中相同的周期分量进行相关性分析，得到相关系数、平均滞后值以及响应关系。

所述步骤1中，小波变换是以Morlet小波作为小波母函数；小波变换的声音数nv＝32。

所述步骤3中，采用启发式方法连续的提取最大能量曲线具体为：采用启发式方法对一个能量曲线提取完成时，将它们的相关能量对应的频率在同步压缩后的时间周期图中被清零，然后对更新后的同步压缩后的时间周期图重复以上步骤进行曲线提取，直到提取到预定的能量曲线为止。

所述步骤3中，提取能量曲线后，在同步压缩后的时间周期图中被清零，清零的误差值clwin＝10。

所述步骤6中，对于周期分量进行信号重建时幅值进行归一化处理，使其值控制在[-1 1]。

本发明的有益效果是：本发明采用采用结合同步压缩小波变换的方法对太阳黑子数和地磁Ap指数的周期分量进行提取，较精确地对太阳黑子数和地磁Ap指数的周期分量进行了提取，较好的解决了之前几种时频分析方法在在提取太阳黑子数和地磁Ap指数的周期分量结果较差，不完整的问题。

附图说明

图1是本发明基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法的总体流程图；

图2是本发明中采用1932年1月到2017年10月的太阳黑子数和地磁Ap指数的月均值数据得到的图像；

图3是本发明中对图2采用归一化处理之后的图像；

图4是本发明中对图2中地磁Ap指数数据经过同步压缩变换之后的图像；

图5是本发明中对图2中太阳黑子数数据经过同步压缩变换之后的图像；

图6是本发明中对图4经过提取能量曲线得到的图像；

图7是本发明中对图5经过提取能量曲线得到的图像；

图8是本发明中对图6经过带通滤波得到的特定的周期分量的图像；

图9是本发明中对图7经过带通滤波得到的特定的周期分量的图像；

图10是本发明中对图8、图9经信号重建并归一化处理后得到的特定周期分量的时间幅度图。

具体实施方式

实施例1：如图1-10所示，一种基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法，所述方法的具体步骤如下：

步骤1：读取太阳黑子数和地磁Ap指数数据(如图2所示)，并进行归一化处理：Ap指数与太阳黑子数量纲不同，为了方便比较，对这两个参数的月均值进行归一化处理；在这里所使用的公式如下：

y＝2*(x-min)/(max-min)-1

其中min为x的最小值，max为x的最大值，输入向量为x，归一化后的输出向量为y。其结果如图3所示；

通过该步骤可以有效地证明本发明研究的太阳黑子数和地磁Ap指数两种数据相互存在规律。

步骤2：对太阳黑子数数据和地磁Ap指数数据分别进行连续小波变换，得到小波系数，通过对小波系数求偏导得到瞬时频率：首先对时间域信号f(t)进行连续小波变换，得到小波系数W_f(a,τ)：

式中：a为尺度因子，τ为时间平移因子，ψ^*(t)是小波母函数ψ(t)的共轭函数。如果

的负频率域趋于0，且主频ξ＝ω₀，则理论上小波系数W_f(a,τ)应该在时间尺度平面上的水平线a＝ω₀/ω周围的区域上展开。然而，实际上小波系数往往在尺度方向上发生扩散，聚焦效果不理想，使得时频图变得模糊。虽然小波系数在尺度方向上有扩散，但是其相位保持不变。这得出了对于W_f(a,τ)≠0的任何(a，τ)，信号f(t)可以用小波系数通过求偏导来计算瞬时频率，其公式表示为：

步骤3：瞬时频率是时频表示中的脊，为了减少淡化，我们必须压缩这些脊周围的频率，即对瞬时频率进行同步挤压：根据步骤2，每个点(a，τ)被转换为(ω_f(a，τ)，τ)，对(ω_f(a，τ)，τ)做同步挤压操作。由于尺度a和时间τ是离散值，频率变量ω和尺度变量a被“封装(binned)”在一起，即仅在离散值a_k下计算W_f(a,τ)，其中a_k-a_k-1＝(Δa)_k。当从时间-尺度平面映射到时间-频率平面(a，τ)→((ω_f(a，τ)，τ)时，同步压缩变换T_f(ω_f,τ)只有在频率范围[ω_l-Δω/2,ω_l+Δω/2]中心的ω_l才是有效的，其中ω_f是ω_l的集合，Δω＝ω_l-ωl-1：

得到的同步压缩后的时间周期图如图4、5所示。

步骤4：对同步压缩后的时间周期图中除COI区域以外的区域采用启发式方法连续的提取最大能量曲线，得到周期分量：假设时间域信号f(t)信号的能量有限，根据Parseval定理，信号时域和频域上的能量相等：

因此，我们称|T_f(ω_f,τm)|²为信号的能量谱密度，即单位Hz里的能量。但是在COI区域(Cone ofinfluence，小波影响锥，表示连续小波系数受到远处点的信号值影响较大的区域)内的周期分量不予提取，因为可信度较低，不予提取。从同步压缩后的时间周期图中连续提取最大能量曲线。为了提取离散曲线c*，我们最大化曲线能量的函数，然后使用一种启发式(贪心)方法，提取最大能量曲线。其公式为：

其中E_f(ω_l,τ_m)＝log(|T_f(ω_l,τ_m)|²)是T_f(ω_f,τ_m)的归一化能量，τ是τ_m的集合。当一个能量曲线提取完成时，它们的相关能量对应的频率在同步压缩后的时间周期图T_f(ω_f,τ)中被清零。然后再次对更新后的同步压缩后的时间周期图再一次进行曲线提取。一直重复以上步骤，直到所有想要得到的能量曲线被提取完成为止(如设置能量少于阈值0.2为终止条件)。所得到的能量曲线就是所要提取的周期分量，能量大小则是该周期分量的幅度大小；提取能量曲线得到的图像如图6、7所示(如图6所示，可知u＝1，u＝2，u＝3，u＝4，...，u＝9表示进行了9次曲线提取，如图7所示，进行了7次提取)。

步骤5：为得到特定的周期分量，对步骤4的结果进行滤波处理：步骤4中抽取的能量曲线很多，为了得到特定的周期分量，对步骤4的结果进行滤波处理，结果会更清楚；得到的图像如图8、9所示。

步骤6：为了研究太阳黑子数和地磁Ap指数之间的响应关系，可以抽取太阳黑子数和地磁Ap指数的相同的周期分量进行信号重建。对步骤5中的结果通过能量曲线对周期分量进行重建，得到该周期分量的时间幅度图：由步骤5提取的能量曲线进行周期分量的重建。同步压缩变换并不影响信号的重建。信号重建的公式为：

式中：ξ是角频率，

是ψ(t)的傅里叶变换。抽取需要仔细研究的周期分量，然后进行信号重建，得到该周期分量的时间幅度图。

步骤7：太阳黑子数和地磁Ap指数的相同的周期分量的幅值进行归一化处理：为了方便观察，对这些周期分量进行信号重建时幅值可能不在同一数量级的，进行归一化处理；得到的图像如图10所示。

步骤8：对步骤6或者步骤7的时间幅度图中太阳黑子数和地磁Ap指数的相同的周期分量进行相关性分析：在相同周期分量中，研究太阳黑子数和地磁Ap指数的相关性，计算相关系数、平均滞后值以及响应关系。其公式为：

同时也在每个太阳周中进行研究，以便得到规律，推断下一个太阳周的太阳黑子数和地磁Ap指数之间的响应关系。

上述实施例中，所述步骤2中，小波变换是以Morlet小波作为小波母函数；小波变换的声音数nv＝32。所述步骤4中，提取能量曲线后，在同步压缩后的时间周期图中被清零，清零的误差值clwin＝10。所述步骤1和步骤7中，对于周期分量进行信号重建时幅值进行归一化处理，使其值控制在[-1 1]。

上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

Claims (5)

1.一种基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法，其特征在于：所述方法的具体步骤如下：

步骤1：对太阳黑子数数据和地磁Ap指数数据分别进行连续小波变换，得到小波系数，通过对小波系数求偏导得到瞬时频率；

步骤2：对瞬时频率进行同步压缩，得到同步压缩后的时间周期图；

步骤3：对同步压缩后的时间周期图中除COI区域以外的区域采用启发式方法连续的提取最大能量曲线，得到周期分量；

步骤4：对步骤3的结果进行滤波处理，得到滤波处理后的周期分量和其能量曲线；

步骤5：通过步骤4中能量曲线对周期分量进行重建，得到该周期分量的时间幅度图；

步骤6：对步骤5得到的太阳黑子数和地磁Ap指数周期分量的时间幅度图中的相同的周期分量的幅值进行归一化处理，得到归一化的周期分量的时间幅度图；

步骤7：对步骤5或者步骤6的时间幅度图中相同的周期分量进行相关性分析，得到相关系数、平均滞后值以及响应关系。

2.根据权利要求1所述的基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法，其特征在于：所述步骤1中，小波变换是以Morlet小波作为小波母函数；小波变换的声音数nv＝32。

3.根据权利要求1所述的基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法，其特征在于：所述步骤3中，采用启发式方法连续的提取最大能量曲线具体为：采用启发式方法对一个能量曲线提取完成时，将它们的相关能量对应的频率在同步压缩后的时间周期图中清零，然后对更新后的同步压缩后的时间周期图重复以上步骤进行曲线提取，直到提取到预定的能量曲线为止。

4.根据权利要求3所述的基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法，其特征在于：所述步骤3中，提取能量曲线后，将它们的相关能量对应的频率在同步压缩后的时间周期图中清零，清零的误差值clwin＝10。

5.根据权利要求1所述的基于同步压缩小波变换的提取太阳黑子数和地磁Ap指数的周期分量的方法，其特征在于：所述步骤6中，对于周期分量进行信号重建时幅值进行归一化处理，使其值控制在[-1 1]。

Images