地球运动的长期演化范文

《天文学进展杂志》2016年第二期

摘要：

Milankovitch理论认为，地球表面日照量受地球运动变化的影响，这些改变又会导致地球气候的变化。因此，古气候学的研究需要相关地球运动变化的精确信息。目前，关于地球长期平动和转动的最新理论是La2004,La2010和La2011模型，它们主要通过数值积分行星轨道和地球自转运动方程得到；这些结果已被用于过去50Ma的地质记录研究中。而超过这一时间范围，解的准确性会由于混沌效应迅速降低。主要介绍这些天文理论的相关内容以及有关地球运动演化方面的最新进展。第2章是关于数值积分太阳系长期运动，第3章介绍长期岁差方面的内容，第4章给出有关地球运动的长周期分析，第5章则讨论数值解的稳定性和太阳系稳定性方面的情况，最后是对目前地球运动长期演化领域的总结以及对未来的展望。

关键词：

地球轨道；历表；岁差

1引言

在太阳、月亮和其他行星引力作用下，地球轨道及其自转轴在空间中的指向都随时间发生缓慢周期性变化。根据Milankovitch理论，这些变化导致地表日照辐射量的改变，从而引起地球气候的周期性变化[1]，而这种周期信息随后被保存在对气候变化敏感的沉积地层中。20世纪70年代，科学家通过分析深海沉积物的氧同位素验证了Milankovitch提出的地球轨道周期与冰期旋回关联的正确性[2]。随后，这种以地球轨道参数和自转变化规律来研究地球气候演变和地质年代的方法，逐渐成为地质学的重要分支学科——旋回地层学[3,4]。目前，旋回地层学在新生代(650Ma)研究中已取得重要进展，但对中生代及更老的地质时期，则需要更高精度、更长时间尺度的地球轨道和自转演化数据作为支撑。早期主要通过解析摄动法计算地球轨道。早在1781年，Lagrange利用分析方法研究了包括水星、金星、地球、火星、木星和土星在内的行星系统运动规律[5,6]。随后在1856年，LeVerrier发表了更高精度的地球轨道结果，他考虑了除海王星外所有大行星的作用[7]。1873年，Stockwell则在LeVerrier的工作基础上增加对海王星的处理[8]。而此前所有的工作都是针对一阶线性动力学方程，Hill(1897)[9]、Bretagnon(1974)[10]和Duriez(1977)[11]则开始考虑更高阶的作用项以提高计算结果的精度。然而所有这些方法都难以在相当长的时间跨度内给出太阳系运动较好的近似结果。随着计算机的出现，数值积分逐渐成为解决此类问题的一种更为行之有效的方法。最早直接数值积分太阳系真实模型(无平均化处理)得到地球平动及自转结果的工作由Quinn等人(QTD)(1991)完成，其积分时间为3Ma，并考虑到广义相对论和月球的影响[12]。随后Wisdom(1991)[13]和Sussman(1992)[14]通过含混和变量的辛积分器计算得到100Ma太阳系运动情况，并验证了太阳系的混沌性质。随着计算机技术的发展，太阳系演化计算的时间跨度不断增加，计算精度也越来越高。Ito和Tanikawa(2002)[15]针对不含月球影响的太阳系牛顿模型，积分得到几十亿年的最长演化结果。Varadi等人(2003)[16]给出范围50Ma的较高精度大行星历表。除此之外，法国巴黎天文台的天文学家Laksar也在相关方面取得巨大进展，其给出的太阳系长期演化模型是目前最权威的数值结果[17–19]。本文作者主要介绍他在地球轨道运动及自转方面的研究工作。1985年，Laskar给出针对太阳系自控系统的一种新的分析研究方法，以计算相应Lagrange解对应频率[20]。1986年则进一步考虑相对论效应和月球摄动作用，用数值方法计算10ka内行星轨道根数和地球自转运动，并给出其多项式表达式[21]。1987年和1988年，Laskar将之前的太阳系演化计算推广到30Ma，并利用傅里叶分析方法给出轨道的拟周期表达式[22,23]。随后，Laskar完善频谱分析理论，提出基本频率数值分析法(NAFF)[24,25]。1993年，Laskar基于之前的工作，用Kinoshita刚体自转理论[26]计算得到长时间地球自转参数，同时根据QTD[12]解里的潮汐耗散模型对数据进行修正，得到地球自转解La93[27]，并通过数值模拟的方法研究地月系统稳定性[28]。1997年，N´erondeSurgy和Laskar利用Hamilton的方法研究刚体地球自转，并考虑潮汐耗散、地核地幔摩擦、大气潮等作用，形成新的地球自转理论[29]。1999年，Laskar研究地球轨道数值计算的局限性[30]。考虑到辛积分器更适合于太阳系长期演化的计算，Laskar等人针对此类摄动哈密顿系统给出SABAC4型辛积分器[31]。2004年，利用此积分器计算太阳系演化，同时在1997年相关工作基础上计算得到地球自转数据La2004[17]。2010年则改进数值模型，给出高精度长时间尺度下地球轨道运动模型La2010[18]。随后La2011模型进一步优化数值积分算法，并以此研究小行星密近交会的现象[19]。与此同时，Laskar利用数值模拟的方法研究太阳系稳定性和混沌现象[32–39]。可以看出，Laskar通过不断提高积分器性能，精化频谱分析方法，改善太阳系动力学模型，从而不断提高长时间尺度下太阳系内大行星轨道与行星自转计算的精度，并研究其运动规律及稳定性。本文将介绍Laskar等科学家在太阳系长期演化及地球自转领域的最新研究成果。第2章介绍利用数值方法计算太阳系长期演化，包括La2010的动力学模型、积分器和所得数值计算结果及分析；第3章介绍有关地球长期岁差方面的内容，并介绍La2004所使用的地球刚体自转模型以及对潮汐耗散影响的处理等；第4章主要介绍地球轨道运动及自转变化的长周期分析，包括频谱分析方法NAFF，以及通过此分析方法得到的地球长期平动和转动模型；第5章则讨论数值积分的局限性，以及太阳系运动的稳定性；最后，我们就以上研究工作进行总结和展望。

2太阳系演化

太阳系是以太阳为中心，和所有受到太阳引力约束天体的集合，主要包括“八大行星”，及数以亿计的太阳系小天体。与太阳系运动相关的N体问题一直为世人关注而未被真正解决。针对N体问题常用的解法包括：分析方法、数值方法以及利用定性方法辅助分析。与其他方法相比，数值计算不失为一种简单高效的处理此类复杂问题的途径。目前随着计算机科技的高速发展，我们可以直接通过数值积分得到上亿年太阳系运动变化情况。目前有很多通用型积分器，例如HNBody[40]、Mercury[41]可以直接用来积分计算太阳系运动演化，而对于特别问题的研究则需要自行设计专用积分器。当前最为权威的太阳系长期演化结果是由Laskar提供的La2010d[18]和La2011[19]。其中La2011为Laskar等人利用历表INPOP10a[42]，并改进数值积分方法后得到的最新太阳系长期演化结果。其与La2010使用的数值模型及算法基本一致，这里我们主要介绍La2010的模型及结果。与历表计算尽可能考虑所有的摄动项不同，在研究天体轨道长期演化时，则一般采用长期理论模型或称为平均模型，此模型在运动方程中去掉全部短周期项[43]。La2010所使用的动力学模型为接近真实太阳系模型，主要包含所有太阳系八大行星，以及冥王星；为了获得精确的地月系统演化结果，月球被当作单独天体加以计算，并考虑地球和月球的主要引力场参数，以及地月系统的潮汐耗散作用；同时考虑最主要的相对论效应——太阳的一阶后牛顿效应[44]；此外还计算5颗小行星Ceres、Pallas、Vesta、Iris、Bamberga所带来的影响，以及针对剩余小行星和其他未模型化参数进行一定的修正。在天文动力学中所遇到的力学系统，很多都是哈密顿系统，或者带有少量耗散的拟哈密顿系统。近些年发展起来的以辛几何为理论框架的辛算法，已成为目前研究太阳系天体轨道长期演化的最佳数值方法。它能保持哈密顿系统的主要性质而不会引入人为的耗散，因而可以用大步长进行长时期的数值积分[45]。相对而言，其他诸多数值方法由于局部截断误差造成人工耗散，而不能保持能量守恒，导致最后计算结构失真[46]。La2010选取SABAC4型辛积分器[31]，此积分器将太阳系整体运动分为两部分：可积部分(行星绕日和月球绕地的开普勒运动)和微小扰动部分(各个行星之间相互作用)，并增加对月球轨道修正处理过程。积分初始条件来自特殊的INPOP历表[47](此历表将INPOP历表结果推广计算至1Ma)。根据上面给出的动力学模型及积分器构造，就可以计算任何时刻的天体运动状态。而数值积分的准确性，则是我们最为关心的问题。La2010针对不同积分模型及积分初值给出四个版本的积分结果La2010a、La2010b、La2010c、La2010d，以方便研究其数值积分的稳定性，表1给出这四个版本的差别。为估计积分结果的误差，采用对相同模型积分但积分步长略微不同的两次积分结果进行比较。图1为La2010x和相应积分La2010x*，x=a、b、c(详见表1)，所得地球轨道偏心率在100Ma内的比较。可以看到在大约50Ma内两者差别基本为零，但由于混沌效应，之外结果差别迅速增大。这也说明由于数值积分方法引入的误差，会导致La2010b,c的结果在55Ma前的结果不可采信，对于La2010a则为60Ma。此外，数值计算天体轨道长期演化通常采用太阳系简化模型，忽略真实模型中的短周期效应。但这种近似是否能真正代表太阳系天体运动的真实情况，特别是所得长期解是否接近天体长时期的真实轨道？我们不得而知。不过，通过改变动力学模型等积分条件，可以利用数值方法研究动力学模型对积分结果的影响，具体内容详见第5章。Laskar发现地球轨道偏心率的变化规律中，一项与木星和土星相关的405ka周期项(详见第4章)在整个250Ma的时间范围内保持稳定，受动力学模型影响较小，而被认为可能存在于地球轨道长期演化规律中。这一周期项也被地质学家用于古地质的研究中。Westerhold等人利用地层学方法得到：在早期古近纪(4754Ma)时期的地质数据与La2011和La2010d所给出偏心率变化周期相符[48]；也说明La2010d和La2011中偏心率参数的适用范围可达54Ma。虽然数值积分方法极大程度地推动了太阳系长期演化领域的研究，但如何得到更精确、时间跨度更长的演化结果仍然是一个十分困难的问题。通过上面的分析可以知道，辛积分器虽然很好地保持了哈密顿系统的特性，但积分精度较低，提高结果的数值精度需要进一步改进数值积分方法；同时为了提高结果的准确性，需要改进太阳系动力学模型，模型不宜过于复杂导致积分时间过长，但又能真实地反映太阳系演化的长周期规律。这样，通过太阳系长期演化的结果，就可以计算地球轨道长期变化规律——地球长期岁差。然后利用频谱分析的方法，进而得到长期岁差模型。

3地球长期岁差

地球岁差是描述在太阳、月球和其他行星的引力作用下，地球自转轴在空间中的长周期进动现象。国际天文联合会(IAU)2006决议B1规定，从2009年1月1日开始，采用P03岁差理论[49,50]取代之前IAU2000岁差章动模型中的岁差部分，作为最新IAU岁差章动模型使用[51]。P03岁差表达式为五次多项式，虽然精度很高，但只适用于J2000.0附近几千年的范围，其误差会在更长的时间范围内迅速增加。此外岁差表现的是一种复杂的长周期进动现象，其主要周期都在上万年以上，这点很难用多项式体现出来。于是长期岁差(longtermprecession)的概念应运而生。长期岁差主要指较长时间范围的地球岁差，其时间范围一般在几十万年以上。长期岁差所包含的周期一般都在一万年以上，较短周期则归为章动部分[52]。从这点来看，长期岁差能更好地体现出岁差的长周期特性。2011年，Vondra´k等人给出与IAU2006岁差模型最为接近的长期岁差理论[52]。根据太阳系天体对地球长期作用所产生两个不同表现，将岁差分为黄道岁差和赤道岁差。其中黄道岁差表示由于行星质点引力对地球轨道面产生的摄动，使得黄道面的空间取向发生变化，引起春分点沿赤道运动。另一方面，赤道岁差表示由于日月和行星对地球赤道突起部分的较差力矩而引起的赤道面(极轴)的长期进动[53]。过去的20年里，很多人尝试用不同方法计算地球赤道岁差。Kinoshita[26,54]、Kinoshita和Souchay[55,56]，以及N´erondeSurgy和Laskar[29]用哈密顿方法计算地球赤道岁差。Williams等人[57,58]利用月亮、太阳及其他行星轨道分析解，计算作用在椭圆型地球上的力矩。Hartmann和Soel[59]则利用数值积分得到的历表计算力矩。Bretagnon等人[60,61]根据最新观测数据(例如VLBI)，基于太阳系历表DE403和月球运行理论ELP2000，用解析的方法给出高精度刚体地球自转理论。最近，Klioner等人建立了相对论自转理论[62]。他们根据Klioner[63]给出的后牛顿自转方程，并考虑相对论力矩、相对论岁差/章动、多参考系及相应时间系统、相对论尺度变换等因素，给出建立在较为严格的相对论框架下的自转理论。但真实的地球并非刚体。很早人们就通过古地质时钟发现地球角动量的损失[64]，现如今通过激光测月技术精确测得月球潮汐加速度[65]。造成这些现象是多方面的，其中一个很重要的原因来自于潮汐耗散现象：月球等其他天体会造成地球的潮汐隆起反应，但地幔滞弹性和海洋能量耗散，使得潮汐变形反应有滞后性[66]，从而产生一个使地球自转角速度减小的力矩和反作用于月球上的力。根据Mignard(1979)[67]，Touma等人(1994)[68]的潮汐耗散模型，认为此力矩与地球形变恢复到平衡状态的延迟时间t成正比。N´erondeSurgy等人利用数学方法将潮汐耗散作用表示为对地球赤道岁差参数的影响，以及此耗散对月球轨道的反作用，具体潮汐耗散公式参考文献[17,29]。此外，地球长期岁差的计算需要考虑更为精细的地月系统模型，其他摩擦作用，例如地核与地幔间的摩擦作用、大气潮汐、地幔对流、气候摩擦(climatefriction)等，也会对地球自转产生影响。在La2004的动力学模型中，Laskar等人利用等效粘滞系数来处理地核与地幔间的摩擦作用[69]，其他效应则未被考虑[17]。于是，根据地月系统轨道参数就可以通过公式(2)和(3)计算得到地球赤道岁差参数ψ和。另外通过数值积分，可以知道此潮汐耗散模型对月球半长轴的改变大约是3.98cm/a，造成现阶段日长的变化为2.68ms/100a[17]。另一方面，现阶段的计算仅使用当前时期的动力学模型，并不能反映整个地质时期的真实情况。例如大陆板块漂移、地球两极冰块覆盖情况变化等因素都会对地球运动长期演化计算产生影响，而这些影响如何模型化也是非常困难的问题。基于相对论地球自转理论[62]，Tang等人(2015)通过数值方法在一阶后牛顿精度下计算地球姿态的长期演化，得到与相对论原理相符的长期岁差模型。通过分析相对论效应对长期岁差的影响，可以知道太阳的一阶后牛顿作用和测地岁差为影响地球长期岁差的主要相对论效应。而其他相对论效应的影响在1Ma这个尺度上很小，可以在大多数情况下不考虑[70]。

4周期分析

用数值方法积分太阳系运动或者地球自转，最终目的都是为了研究天体的运动规律，特别是地球轨道在其他行星摄动作用下的变化规律。对于轨道半长轴，带内行星的变化非常小，而带外行星则由于平均运动共振现象变化较为明显[71,72]。在地质历史时期，至少过去的250Ma内，地日平均距离没有明显变化[17]。虽然Poincar´e曾预言，由于混沌现象，太阳系中的行星轨道无法外推至无穷远时刻[76]。不过在有限时间内(例如几百万年)，我们仍然可以利用数值积分得到行星轨道数值解，并通过频谱分析方法给出轨道的拟周期近似解，而所得的结果也被广泛应用于古气候学的研究当中。在2004年和2010年，Laskar等人通过数值方法求得轨道频率gi和si，详见表2(第4列和第5列)[17,18]。其中4颗带内行星相关频率是利用La2004和La2010中20Ma的数据计算得到的，5颗带外行星则利用50Ma的数据求得。频率分析方法主要用来从动力学系统的数值解中找到主要的基本频率，然后得到拟周期函数以用来研究系统长期演化规律。Laskar等人改进相关频谱分析方法，针对研究太阳系混沌现象而给出NAFF方法[24,25]。整个分析方法过程大致如下：通过傅里叶变换公式找到具有最大作用的频率项，然后将此频率项的贡献从数据中去除，再重复上面的步骤从剩余部分中找到其他主要频率项。同时利用增加Hanning窗函数的方法提高寻找频率的准确度。此方法能精确找到所需的基本频率，以便更好地分析动力学系统的运动规律。

5稳定性

将地球轨道数值演化结果运用于地质研究中，目前已在新近纪时期(23.030Ma)研究中取得巨大成功[85]。下一步正努力提高天文轨道计算结果的精度，进而用于古近纪时期(6523Ma)。对于更遥远的地质年代，则很难获得精确的地球轨道运动解[17]。对于这种长时间尺度下的数值积分，积分结果的稳定性是我们最关心的内容。目前没有太多行之有效的手段来估算长时间积分的精度。除利用位置、能量或沿迹误差，还可以通过略微改变积分初始条件或者积分步长，然后对两次积分结果的差异进行比较评估。对于太阳系长期演化问题，积分误差来源有很多，包括数值积分方法的限制、舍入误差的影响、积分初始条件的误差、积分模型与真实模型的差别等。同时由于系统的混沌效应，这些误差会在一定时间Tv后以指数的形式迅速被放大。Laskar对不同情况下时间Tv的大小进行估算。对于由于积分方法和截断误差对太阳系长期演化的限制大约是56Ma。这里针对相同的动力学模型，通过改变积分步长和初始条件分别积分，以检验结果的稳定性。而提高精度的办法就是改进数值积分算法，或者是提高计算机计算精度[17]。另一项测试则针对动力学模型的不确定性，相应结果列在表5。当时认为最大的不确定性源于太阳的扁率，不过由于科学技术的发展，通过太阳和太阳风层探测器(SOHO)以及太阳全球振荡监测网(NG)的测量可以得到JS2=(2.18±0.06)×107[86]，也将这一限制推后至大约40Ma[17]。对于地球自转的计算，不确定性则源于潮汐耗散模型，以及由于冰川覆盖和地幔对流对地球力学椭率的影响等[17,27]。数值积分方法让我们从一个全新的角度去研究太阳系长期演化问题。自牛顿万有引力定律诞生以来，人们一直希望解决N体问题的计算难题。随后Poincar´e发现大多数情况下这个问题不可解[76]，即使对于最简单的三体问题，仍然不可能找到一个解析式将行星轨道外推至无穷远。随后出现的KAM定理则进一步研究了在微小扰动下的非退化哈密顿系统存在拟周期运动[87–89]，并扩展运用于行星科学研究其稳定性[90–92]以及轨道发散变化性[93–96]。这些理论虽不能直接用于真实的太阳系，不过目前普遍认为太阳系在其有限的生命中是稳定的[39]。随着计算机科技的发展，太阳系稳定性的研究也取得巨大进展。在计算机的帮助下，Laskar在1994年对数十亿年太阳系的多种可能演化结果进行模拟，发现带外行星运动较为稳定，而带内行星出现明显混沌和轨道漂移；并且水星和金星可能发生碰撞[38]。随后改进数值计算方法和动力学模型，Laskar与合作者Gastineau于2009年公布对太阳系各种可能演化情况的庞大模拟结果[38]。在2501种可能性中(见图2)，大多数结果表明在5Ga内的太阳系仍然像现在一样运转着，没有发生巨大变化。然而有1%的结果预测水星轨道偏心率增加，导致可能与金星或者太阳发生碰撞，同时其中一个结果造成火星轨道偏心率增加，进而引起带内行星运动的紊乱。太阳系演化的数值模拟结果一定程度上给出了太阳系稳定性的概率解。发生灾难性变化可能性为1%，这一可能性比之前预想的要高。更为令人惊讶的是，如果只考虑纯牛顿效应，这种混乱可能性会提高到60%，或许这也间接说明相对论的正确性[39]。

6总结与展望

近20年，随着计算机科学的发展，数值计算和数值分析的方法极大程度地提高了我们对于整个太阳系演化的认识和了解，从地质时期地球轨道的计算到未来太阳系稳定性的研究都取得巨大突破性进展。通过数值积分，我们可以得到数百万年内较为精准的太阳系演化结果，特别是地球轨道运动和自转结果。借助频谱分析的方法，也使我们对行星运动规律有了更深刻的了解。对太阳系演化可能性的模拟，也在一定程度上回答了太阳系稳定性的疑问。本文主要介绍了太阳系演化及地球自转方面的发展状况和趋势，包括Laskar等科学家在行星轨道计算及分析，以及太阳系稳定性方面的工作情况。目前已公布的最新研究成果是地球轨道运动解La2010和地球自转运动解La2004，相关数据可以在网页上下载。这些数据主要利用辛积分器直接积分太阳系真实模型得到。其中地球偏心率的结果适用于过去的54Ma。而超过这一时间范围，解的准确性会因混沌效应迅速降低。如何突破这一限制，提高地球长期演化解的精度和扩大适用时间范围都将成为未来此领域研究重点。一方面需要提高数值积分算法的速度和精度；另一方面，目前所使用的动力学模型为简化的太阳系模型，如何完善这一模型，考虑更微小的作用；同时加强与地质学家的合作，从地质记录中获得地质历史时期的精确地球轨道参数周期，进而对模型进行修正，也是今后工作的方向。我国于2013年启动“地时–中国(EARTHTIME-CN)”计划[97]，主要通过建立与国际相关组织的联系与合作，促使国内地质年代学、古地磁学和旋回地层学等实验室之间的联合研究，进一步提高我国高精度地质年代学和定量年代地层学研究的水平和国际影响[3]。该计划将促进中国的地质学家和天文学家在未来更紧密地合作，进而为地球科学和太阳系演化研究做出更大贡献。

作者：唐凯 唐正宏 陶金河 单位：中国科学院 上海天文台