天文观测中通常采用地球定向参数(Earth orientation parameters, EOP)来描述地球自转运动, 包括三个部分, 一是地球自转轴相对于地球本体的位置发生的周期性变化, 称为地极移动; 二是地球自转轴相对于空间的变化, 周期性变化部分称为章动, 长期性变化称为岁差; 三是地球自转轴的旋转, 直接反映的观测量是地球自转速率(日长)变化. 地球自转运动反映了固体地球与大气、海洋、地幔与地核在各种时空尺度上的相互作用过程^[1–4], 研究地球自转运动的整体行为演化规律, 有助于加深我们对地球自转运动过程的理解, 对于构建一个能定量的、与观测相吻合的地球自转理论模型具有十分重要的理论意义^[5,6]. 地球自转运动的非线性动力学特征研究是天文地球动力学领域中一项十分重要的研究内容, 一方面可以加强我们对地球自转运动产生和演化过程的理解, 有助于深入研究地球自转变化规律, 为建立地球自转变化理论模型提供科学依据, 从而为研究地球动力学过程提供重要的测量约束; 另一方面, 采用新颖的混沌和分形理论来研究地球自转运动的非线性特征, 对于地球自转运动在不同时间尺度上的预测与重构是十分重要的. 本文以反映地球自转速率变化的日长变化(length of day changes, ΔLOD)观测量作为研究对象, 分析地球自转速率变化的周期、混沌和分形非线性特性.

国际地球自转与参考系统服务组织(International Earth Rotation and Reference System Service, IERS)发布包含ΔLOD的长期EOP观测资料^[7], 大部分学者在分析ΔLOD观测资料时, 最关注的是ΔLOD的周期性变化. 研究发现, ΔLOD观测序列中存在多个显著的周期性信号, 如准两年、周年、半周年、周月和半月信号, 具有多时间尺度特性^[8–11]. 除了上述比较显著的周期性信号, ΔLOD序列中还包含一些频率较低的年际和年代际长周期性信号, Yu等^[12]利用小波分析方法从1962—2012年的ΔLOD序列检测到周期为6 a, 13 a和18.6 a的低频信号; Chao 等^[13]运用集合经验模态分解方法从1962—2015年的ΔLOD序列中得到周期为6 a 和13.69 a的信号分量; Shen 和Peng^[14]采用Morlet小波变换方法从1760—2018年的ΔLOD序列中检测到8.5 a, 11 a, 13.5 a, 18.6 a, 22.3 a, 33 a, 68 a和149 a 8个长周期性信号. 近年来, Ding^[15], Duan和Huang^[16]分别运用标准小波变换和AR-z谱分析方法对时间跨度长短不同的ΔLOD观测资料进行分析, 发现ΔLOD序列中除存在约6 a的年际变化信号, 还呈现8.6 a的亚十年尺度波动. 虽然ΔLOD的周期性特征取得了显著性成果, 但周期信号分析方法只能获得信号的时域或频域方面的信息.

20世纪60年代, 混沌与分形非线性科学迅速发展, 混沌是外在表现为无序的随机现象而内部存在确定性的规律, 是有别于周期信号与随机信号的第三类系统^[17,18]. 分形理论打破了传统几何学所不能描述的整体和局部自相似的几何图形. 因此, 混沌和分形被认为是继相对论与量子力学之后的第三次物理学革命. 混沌与分形理论被广泛应用于气地球科学、空间科学和天文学等科学领域^[17–20], 拓展了这些科学领域研究的广度与深度. 地球自转运动具有复杂的非线性特征, 目前这些研究仅单一地侧重于多时间尺度一方面来分析, 对于地球自转运动的混沌和分形特征的研究却很少, 无法全面地揭示地球自转运动的本质规律. 目前混沌识别的主要方法使用K熵和最大Lyapunov指数, 但两个物理量的计算结果易受数据长度、观测噪声及参数设置的影响, 导致计算结果不可靠. 定量递归分析(recurrence quantification analysis, RQA)法克服了传统方法的缺点, 无需考虑数据长度、观测噪声及数据分布等信息, 计算结果具有较高的可靠性, 是定性与定量研究非线性动力系统混沌特征的有效方法^[21], 被广泛应用于机械故障识别、大地电磁信噪辨识太阳活动动力学特征研究等领域. 分形维是一种用于描述系统动态复杂性的量度, 它与动力系统的混沌理论交叉结合、相辅相成, 为描述系统动态行为的复杂度提供了新的视角和方法^[22]. 分形维有多种定义, 如盒子维^[23]、豪斯道夫(Hausdorff)维^[24]和关联维^[25], 其中关联维是从时间序列数据中直接计算得出的, 使其成为刻画动力系统是否具有分形特征的重要定量指标, 本文采用Grassberger和Procaccia提出的GP算法计算关联维.

为从多角度视角探讨ΔLOD的演化规律, 本文将周期分析方法和非线性分析方法相结合, 从周期性、混沌动力学和分形特性对ΔLOD的非线性特征进行了全面的分析, 首先采用自适应噪声完备经验模态分解(complete ensemble empirical mode decomposition with adaptive noise, CEEMDAN)算法分析ΔLOD的多时间尺度特性, 它能将ΔLOD序列分解为趋势项、周期项和残差项, 如果ΔLOD具有混沌特性, 混沌成分应该包含在扣除趋势项和周期项的残差项中; 然后将剩余项当作新的时间序列, 通过RQA 和GP算法对剩余序列的混沌与分形特性进行分析, 并对比分析原ΔLOD序列的分析结果. 通过本文对ΔLOD的非线性动力学特性的分析, 为全方位、多角度理解地球自转运动提供一个新视角. 同时若ΔLOD具有混沌与分形特性, 就可以借助非线性动力学模型来建立地球自转速率变化理论模型, 为地球自转速率变化建模和预测提供理论支持.

IERS提供长期EOP观测数据, 这些数据早期由光学天体测量方法获得, 随着上世纪70年代甚长干涉基线、全球定位系统和人造卫星激光测距等现代空间测地资料逐渐加入到EOP测量中, EOP测量精度取得了很大提高, 目前ΔLOD的观测精度可以达到μs水平^[26], 能够用来分析ΔLOD的非线性动力学特性. 本文选取IERS EOP 14C04产品中1962-01-01至2023-12-31的ΔLOD观测序列进行分析^[26], 数据来源于https://www.iers.org/IERS/EN/DataProducts/EarthOrientationData/eop.html, 采样间隔为一天, 观测序列如图1所示.

CEEMDAN算法是在经验模态分解算法的基础上发展的, 它通过在原始时间序列中自适应地 加入相反的高斯白噪声, 以解决经验模态分解和集合经验模态分解的模态混叠和重构误差大的问题, 该算法的具体步骤如下所述^[8], 其具体分解流程如图2所示.

步骤1　计算第一阶本征模态函数(intrinsic mode function, IMF)分量IMF₁. 向给定时间序列$ x\left( t \right) $($ t=1, 2, \cdots , n) $中添加高斯白噪声得到新序列$ {x_I}\left( t \right) = x\left( t \right) + \beta {\xi _I}\left( t \right) $ , 其中, β为信噪比, $ {\xi _I}\left( t \right) $为第I次添加的服从标准正态分布的高斯白噪声, $ I = 1, 2, \cdots , L $. 对产生的L个新序列$ {x_I}\left( t \right) $分别进行经验模态分解, 获得L个$ {\text{IM}}{{\text{F}}_{I, 1}} $分量, 进行总体平均后得到第一阶IMF分量IMF₁, 则$ {\text{IM}}{{\text{F}}_1}\left( t \right) = \dfrac{1}{L}\displaystyle\sum\nolimits_{I = 1}^L {{\text{IM}}{{\text{F}}_{I, 1}}} $, 余项$ {r_1}\left( t \right) = x\left( t \right) - {\text{IM}}{{\text{F}}_1}\left( t \right) $.

步骤2　计算第二阶IMF分量IMF₂. 向余项$ {r_1}\left( t \right) $中添加L次高斯白声后产生L个新序列$ {r_{1, I}}\left( t \right) = {r_1}\left( t \right) + \beta {\xi _I}\left( t \right) $, 对其进行经验模态分解得到相应新序列的第一阶IMF分量$ {\text{IM}}{{\text{F}}_{I, 1}} $, 则$ {\text{IM}}{{\text{F}}_2}\left( t \right) = \dfrac{1}{L}\displaystyle\sum\nolimits_{I{\text{ = 1}}}^L {{\text{IM}}{{\text{F}}_{I, 1}}} $, 余项$ {r_2}\left( t \right) = {r_1}\left( t \right) - {\text{IM}}{{\text{F}}_2}\left( t \right) $.

步骤3　重复以上处理过程, 得到满足要求的IMF分量与相应的余项, 直至余项为单调函数不满足经验模态分解条件, 终止计算, 则原时间序列最终被分解为K个IMF分量和一个余项之和, 即

频谱分析通过对时间序列的幅值进行分析, 找出序列中存在的显著周期项, 这些周期项在频谱图上表现为具有较大的幅值特征. 具体而言, 频谱分析通过傅里叶级数展开式来计算时间序列的频谱值, 对于离散傅里叶级数来说^[27]:

式中, $ f\left( k \right) $表示频谱值.

根据(2)式可计算序列各数据点的频谱值, 进而画出序列的频谱图. 为降低计算复杂度, 可通过选择合适的时间序列长度n, 使其满足n = 2^c, $ c = 0, 1, 2, \cdots $, 对于长度不满足2的幂次要求的, 可通过向时间序列中增加零值, 使序列的长度满足2的幂次要求, 再利用快速傅里叶变换(fast Fourier transform, FFT)算法计算频谱值. FFT算法的基本流程^[28]: 1) 将n点的离散傅里叶变换分解为两个$ {n {/ } 2} $离散傅里叶变换和一个旋转因子的计算; 2) 递归地应用步骤1方法, 直至所有数据点的离散 傅里叶变换都计算完毕; 3) 合并所有的离散傅里叶变换结果, 获得完整的频域表示. 限于篇幅, FFT算法的详细过程可参考Yeh等^[27]报道.

RQA方法通过在高维相空间中构建递归图及提取递归图中递归点的分布统计规律对非线性动力系统的行为进行定性和定量描述. 递归图是由Eckmann等^[29]于1987年提出以相空间重构技术为基础, 通过吸引子结构特性, 运用二维图形来揭示时间序列的内部结构. RQA方法的具体步骤如下所述.

1)对给定时间序列$ \left\{ {\left. {x\left( 1 \right), x\left( 2 \right), \cdots , x(n)} \right\}} \right. $, 重构m维相空间中的相点, 即$ {\boldsymbol{X}}\left(i\right)=\left[x\left(i\right), x\left(i+\tau \right), \cdots , x\left(i+\left(m-1\right)\tau \right)\right] $, $ i=1, \text{ }2, \cdots , N $, 其中m为嵌入维数, τ为延迟时间, N为重构相空间的相点个数, 且$ N = n - (m - 1)\tau $.

2)计算相空间中任意两点的

式中, $ \left\| {} \right. \cdot \left. {} \right\| $为距离范数; λ为截止距离; $ \theta \left( x \right) $为Heaviside阶跃函数, 当$ x \geqslant 0 $时, $ \theta \left( x \right) = 1 $, 当$ x < 0 $时, $ \theta \left( x \right) = 0 $. 若$ {R_{i, j}} = 1 $, 则在$ N \times N $的坐标平面的相应位置$ \left( {i, j} \right) $处描点, 这样就能够得到一副递归图. 不同系统产生的时间序列使得相空间中吸引子的表现不一样, 递归图展示的结构也会不同, 若递归图中除主对角线外表现为无规则的离散点, 则时间序列为随机信号, 如图3(a)所示; 若图中表现为一些平行主对角线的线段, 则时间序列为混沌信号, 如图3(b)所示; 若图中表现为棋盘或网格结构, 则时间序列为周期信号, 如图3(c)所示.

虽然递归图可以对系统动力学特性进行可视化分析, 但缺少相对客观的定量分析. Zbilut和Webber在递归图的基础上, 根据递归点的特点, 提出了一些定量分析指标, 鉴于本文目标, 选取确定性(DET)作为定量分析指标, 它描述平行主对角线线段的递归点数与总递归点数的比值:

式中, $ P\left( l \right) $是长度为l的对角线分布概率, $ {l_{\min }} $为最小对角线长度, 一般取2. 若DET越接近0, 表示信号随机性成分越多, 确定性成分越少; 若DET越接近1, 表示信号确定性成分越多, 随机性成分越少.

Grassberger和Procaccia^[30]在相空间重构理论基础上提出从时间序列中计算关联维的GP算法. 首先对给定时间序列$ \{ { {x( 1 ), x( 2 ){, } \cdots {, }x( n ){\text{ }}} \}} $, 重构m维相空间中的相点, 即$ {\boldsymbol{X}} (i) = [ x (i), x(i + \tau ), \cdots , x (i + (m - 1)\tau ) ] $; 然后得到关联积分:

式中, $ C\left( {m{, }\lambda } \right) $是一个累积分布函数, 表示相空间中的两个相点之间距离小于截止距离λ的概率. 关联维D与关联积分存在$ C\left( {m{, }\lambda } \right) \propto {\lambda ^{D\left( m \right)}} $, 关联维数D会随着嵌入维数m的增大而逐渐收敛, 最终达到某一收敛值, 收敛时的维数即为饱和关联维, 即

实际计算中, 利用最小二乘法拟合数据点对($ {\text{In}}C\left( {m, \lambda } \right) $, $ {\text{In}}\lambda $), 拟合直线的斜率即为关联维D. GP算法估计关联维的流程如图4所示.

对1962-01-01至2023-12-31的ΔLOD序列进行CEEMDAN分解, 获得9个不同时间尺度的IMF分量及一个表示ΔLOD总体随时间变化的趋势分量(图5), 这些IMF分量按频率从高频到低频依次排列, 它们隐含了ΔLOD不同时间尺度的局部特征信息, 其中, IMF₁和IMF₂分量波动较为剧烈, 属于ΔLOD序列中存在的几天至几十天准周期波动的高频尺度分量或亚季节性分量, IMF₃分量呈现貌似随机特性, 可能为随机时间序列, 亦有可能是混沌时间序列, IMF₄和IMF₅分量为周年和半周年等季节性变化信号, IMF₆分量为准两年振荡周期性信号, IMF₇分量为尺度约5 a—10 a的年际变化信号, IMF₈和IMF₉分量为尺度数十年的年代际变化信号. 各模态分量和趋势分量在ΔLOD序列中所占的比例并不相同, 如果采用方差贡献率来衡量不同分量在ΔLOD序列中所占的比例(表1), 趋势分量对ΔLOD的贡献率最大, 达到26.91%, 其次为年代际低频变化信号(IMF₈和IMF₉分量), 两者的方差贡献率之和为24.45%; 高频分量或亚季节性分量IMF₁和IMF₂的贡献率分别为10.43%和6.03%; 周年和半周年等季节性变化信号(IMF₄和IMF₅分量)对ΔLOD的贡献率之和为19.59%; 准两年振荡信号IMF₆分量和亚十年变化信号IMF₇分量的贡献率相当, 分别为4.45%和4.29%; IMF₃分量(可能是随机信号或混沌成分)的贡献率最低, 为3.85%.

通过CEEMDAN方法将ΔLOD分解成趋势信号、周期信号和随机信号(可能是混沌信号), 现有的文献在对ΔLOD进行周期分析时, 大部分并没有借助类似经验模态分解方法先分解出周期成分, 再分析其周期性. 本文利用FFT算法对IMF₁—IMF₉分量进行频谱分析, 获得各分量中的平均周期如图6所示. 除IMF₃分量外, 其余IMF分量均呈现周期性波动, 表明ΔLOD具有多尺度特征, 其中, IMF₁分量包含平均周期6—15天的高频振荡信号, IMF₂分量包含2个周期30天左右的亚季节性信号, 这些高频振荡和亚季节性信号主要受大气高频振荡激发^[30]; IMF₄分量包含1/4, 1/3和1/2周年季节性信号, 平均周期分别为91.35天, 121.8 天和182.69天, 主要激发因素是大气、海洋等地表流体与固体地球之间角动量的交换^[31], 由于IMF₁, IMF₂和IMF₄分量中存在的周期信号的频率比较接近, CEEMDAN方法未能将这些频率接近的周期信号有效地分离, 仍然存在EMD方法的模态混叠缺陷^[10]; IMF₅—IMF₇, IMF₉分量为单一信号成分, 平均周期分别为365.39 d, 2.4 a, 6 a和31 a, IMF₈分量存在平均周期为13.4 a和18.6 a的两个振荡信号, 其中, 365.39天的周年信号与季节性信号的激发机制相同, 均源于大气、海洋等地表流体与固体地球之间角动量的交换^[31], 准两年信号主要与大气平流层的准两年振荡有关^[32], 周期6 a的亚十年尺度变化和周期13.4 a的年代际变化可能与地核和固体地球之间角动量的交换有关, 但二者之间的耦合机制尚不明确^[11]; 18.6 a和31 a的低频信号分别由日、月引潮力产生的固体地球周期性形变和冰期后地壳反弹引起^[13,33]. IMF₄和IMF₅分量包含的季节性信号的频谱幅度较大, 二者的方差贡献率之和为19.59%, 大于其他IMF分量, 说明ΔLOD主要以大气和海洋主要激发源引起的季节性振荡变化为主^[1].

周期信号混杂在趋势信号与随机信号(可能为混沌信号), 是否会影响ΔLOD的周期分析结果? 为得到结论, 直接用FFT算法对扣除可能为随机信号或混沌成分的IMF₃分量的ΔLOD序列进行频谱分析, 周期分析结果如图7所示. 分析结果发现, 扣除混沌成分(可能为随机信号)的ΔLOD序列与经过CEEMDAN分解获得的IMF分量得到相同的周期分析结果, 说明随机信号(可能为混沌信号)虽然混杂在周期信号中, 但对周期信号分析结果没有影响.

需要指出的是, ΔLOD观测序列中除广泛熟知的高频信号、亚季节性和季节性信号以及18.6 a信号外, 6 a, 13.4 a和31 a周期信号均为近年来在ΔLOD观测序列中新发现的信号, 如Shi等^[12]利用小波分析方法从1962—2012年的ΔLOD序列中检测到6 a和13 a周期信号, Duan等^[5]利用标准时频变换方法也从同样时段的ΔLOD序列中检测到6 a周期信号. Shen和Peng^[14]利用集合经验模态分解方法从1962—2015年的ΔLOD序列中提取到6 a和13.69 a周期信号. Ding等^[15,17,34]先后利用Morlet小波变换、AR-z谱分析和最优序列估计方法从长期ΔLOD序列中探测到5.9 a, 13.5 a和33 a周期信号. 由于采用ΔLOD观测序列的时间跨度不同以及信号分析方法不同, 周期信号的分析结果也有差异. 因此, 本文6 a, 13.4 a和31 a周期分析结果与已有研究结果之间存在的略微差异, 可能是由信号分析方法的频率分辨率和数据分析时段不同导致的.

为了得到准确可靠的结论, 在递归分析中, 选取4组不同的参数λ, 得到4组递归图, 其中每列为一组, 如图8所示, 从上至下每列依次表示原ΔLOD序列的递归图和扣除趋势项和周期项后的残差项(IMF₃分量)的递归图. 从这8幅图中可以看到, 有一些线段与主对角线平行, 具有一定的规律性, 通过递归图的定义可以确定原ΔLOD序列及其扣除趋势项和周期项后的残差项均具有混沌特性; 在每组2幅图中, 相对残差项, 原ΔLOD序列包含相对较多的平行主对角线结构, 说明其确定性成分较多, 随机性成分较少, 可预测性相对强, 又由于其具有混沌特性, 进而其混沌程度相对弱; 而残差项包含相对较少的平行主对角线结构, 说明其随机性较强, 确定成分较少, 可预测性相对弱, 具有相对更强的混沌程度.

由于递归图只能对系统进行定性的动力学特性分析, 使用RQA方法中的DET指标来定量刻画递归图中包含相对较多平行主对角线结构的特征, 这为RQA方法研究ΔLOD的非线性特性提取提供了有利条件. 计算了4组DET的平均值, 原ΔLOD序列及其扣除趋势项和周期项后的残差项的DET的平均值分别为0.987±0.007和0.979±0.004, 所有DET都接近于1, 说明DET指数所表征的地球自转速率变化含有较高的确定性成分, 不是随机信号. 为定量对比原ΔLOD序列及其扣除趋势项和周期项后的残差项的混沌特性, 将两者的DET数值进行相对相差分析, 相对相差约为0.82%, 说明原ΔLOD序列及其扣除趋势项和周期项后的残差项的混沌特性差别不是很明显. 接下来从统计学的角度来客观判断它们是否具有显著性差异, 因为DET为不同λ时($ \lambda =0.4, \text{ }0.5, \text{ }0.6, \text{ }0.7 $), 4组数据的平均值加减标准差, 所以可以用统计学的T检验进行分析, 结果发现原ΔLOD序列及其扣除趋势项和周期项后的残差项的DET的统计概率为0.03, 小于显著性水平0.05, T检验结果说明, 原ΔLOD序列及其扣除趋势项和周期项后的残差项的DET无显著性差异. 由于它们都是混沌信号, 而且包含不同的确定性成分, 则它们的可预测性强弱会有所不同.

图9给出了ΔLOD序列分形关联维D与嵌入维数m之间的关系, 首先在图9(a)不同嵌入维数m下对应的颜色曲线中找出存在直线的部分, 然后求出每条曲线中直线区域的斜率, 所得直线的斜率就是每个嵌入维数m所对应的分形关联维, 最后绘制由此得到的分形关联维与嵌入维数m之间的关系(图9(b)). 从ΔLOD序列的关联维随嵌入维数m的变化曲线图中可以看出, 起始关联维随嵌入维数 的增大而增大, 当嵌入维数$ m{\text{ > 14}} $时, 关联维的变化趋于缓和, 在一定误差范围内基本不变, 趋向一个饱和值, 该值即为饱和关联维, 即D = 4.678 ± 0.139.

ΔLOD扣除趋势项和周期项的残差项分形 关联维D与嵌入维数m之间的关系如图10所示, 同理, 用上述同样的方法, 可以获得ΔLOD扣 除趋势项和周期项的残差项的饱和关联维为D = 4.511 ± 0.172, 最小嵌入维数为12, 因此最少需要5个独立变量才能描述动力系统的变化规律, 并说明地球自转速率变化的混沌吸引子具有分形结构. 根据混沌动力学理论, 关联维若不存在饱和关联维数, 则信号为随机信号, 反之若存在饱和关联维, 则信号为混沌信号, 可以进一步确认原ΔLOD序列及其扣除趋势项和周期项的残差项均含有混沌成分, 具有混沌特性. 原ΔLOD序列比其扣除趋势项和周期项后的残差项的饱和关联维大, 关联维相对相差约为3.7%, 表明原ΔLOD序列的复杂度强于其扣除趋势项和周期项后的残差项, 这是由于趋势项和周期项与混沌信号叠加在一起, 分形结构的复杂性相对更强. 由于关联维D是经过回归分析得到的, 故采用邹氏检验来验证这些分形维的差异性, 结果发现, 原ΔLOD序列和其扣除趋势项和周期项后的残差项的分形维的F值大于F_0.05(2, 23)的临界值3.422, 说明原ΔLOD序列和其扣除趋势项和周期项后的残差项的分形维具有显著性差异.

地球自转变化非常复杂, 已有研究仅从多时间尺度和复杂性其中某一方面来分析, 忽略了其他部分的有效信息, 且传统的非线性分析方法由于受噪声、非平稳性、数据长度以及参数的影响, 无法准确描述与刻画其本质特征, 导致虚假结果, 本文选择反映地球自转速率变化的ΔLOD观测量作为研究对象, 从多时间尺度、混沌和分形特性方面对其非线性特性进行全面分析. 首先利用数据驱动的自适应CEEMDAN算法对1962-01-01至2023-12-31的ΔLOD观测序列分析其多时间尺度特性, 然后选用新颖的RQA方法来定性和定量研究ΔLOD的混沌特性, 该算法简单, 计算速度快, 具有很强的抗干扰能力, 再结合快速计算关联维的GP算法进行分形特性分析. 结果发现, ΔLOD序列由趋势成分、周期成分及混沌成分构成, 同时包含线性和非线性信号, 是极其复杂的系统, 该系统具有确定性和分维性等非线性动力学特性, 不能较好地适用于传统如自回归时间序列分析和自适应滤波器等基于线性系统的建模与预测方法, 这与已有研究的结论并不矛盾, 而是从不同角度印证了地球自转变化的复杂性. 在进行多时间尺度分析和混沌、分形特性识别时, 由于趋势成分、周期成分及混沌成分混杂在一起, 可能影响最终的分析结论, 因此又对比分析了原ΔLOD序列和其扣除周期成分或混沌成分后的时间序列分析结果. 得到如下结论.

1)扣除混沌成分后仅保留趋势项和周期项的数据与原ΔLOD序列的周期完全相同, 这种结果可能与混沌成分本身就无周期性有关, 本文求得的平均周期与已有研究结果相比较一致, 支持Ding等^[15,17,34]利用Morlet小波变换、AR-z谱分析和最优序列估计方法分析周期的结论.

2) ΔLOD序列貌似具有随机属性, 但根据RAQ方法可知其实是混沌时间序列, 存在内在的、有序的混沌动力学机制, 是混沌动力学系统演化的结果, 具有可预测性, 扣除趋势成分、周期成分仅保留混沌成分的数据与原ΔLOD序列的确定率DET无显著性差异.

3)根据分形关联维可知, 至少需要5个以上的变量才能有效描述地球自转速率变化的动力学特征, 从混沌时间序列相空间重构的角度可初步认为模拟地球自转变化动力学系统所需的独立变量至少需要5个, 这从侧面印证地球自转变化的激发因素是多源的.

4) 原ΔLOD序列的关联维大于其扣除趋势项和周期项后的时间序列的关联维, 表明前者具有更复杂的分形结构, 因此趋势成分和周期成分叠加在混沌成分中, 增加了动力学系统的复杂性, 可预测的时间更短, 预测难度更大.

本文从周期、混沌与分形多角度对ΔLOD观测序列的非线性特性进行了研究, 获得了定量的分析结果, 对于深入了解地球自转变化的本质规律具有重要理论意义. ΔLOD若存在混沌与分形特性, 就可以借助非线性动力学模型来建立地球自转速率变化理论模型, 如基于混沌理论通过相空间重构技术将复杂的经验参数建模转化为以ΔLOD固有动力学特性实现建模, 或基于奇异吸引子的分形结构构建迭代函数跟踪混沌运动轨迹. 本文结果为建立地球自转速率变化动力学数学模型与开展地球自转速率变化重构及预测提供科学依据, 数据去除周期项或混沌项前后非线性动力学特性的异同, 也可以为其他领域的研究者提供一个参考.