算法系列之十八：用天文方法计算二十四节气（上）

二十四节气在中国古代历法中扮演着非常重要的角色，本文将介绍二十四节气的基本知识，以及如何使用VSOP82/87行星运行理论计算二十四节气发生的准确时间。

中国古代历法都是以月亮运行规律为主，严格按照朔望月长度定义月，但是由于朔望月长度和地球回归年长度无法协调，会导致农历季节和天气的实际冷暖无法对应，因此聪明的古人将月亮运行规律和太阳运行规律相结合制定了中国农历的历法规则。在这种特殊的阴阳结合的历法规则中，二十四节气就扮演着非常重要的作用，它是联系月亮运行规律和太阳运行规律的纽带。正是由于二十四节气结合置闰规则，使得农历的春夏秋冬四季和地球绕太阳运动引起的天气冷暖变化相一致，成为中国几千年来生产、生活的依据。

二十四节气起源于中国黄河流域。远在春秋时代，古人就开始使用仲春、仲夏、仲秋和仲冬四个节气指导农耕种植。后来经过不断地改进与完善，到秦汉年间，二十四节气已经基本确立。公元前 104年，汉武帝颁布由邓平等人制定的《太初历》，正式把二十四节气订于历法，明确了二十四节气的天文位置。二十四节气天文位置的定义，就是从太阳黄经零度开始，沿黄经每运行15度所经历的时日称为“一个节气”。太阳一个回归年运行360度，共经历24个节气，每个公历月对应2个节气。其中，每月第一个节气为“节令”，即：立春、惊蛰、清明、立夏、芒种、小暑、立秋、白露、寒露、立冬、大雪和小寒等12个节令；每月的第二个节气为“中气”，即：雨水、春分、谷雨、小满、夏至、大暑、处暑、秋分、霜降、小雪、冬至和大寒等12个中气。“节令”和“中气”交替出现，各历时15天，人们习惯上把“节令”和“中气”统称为“节气”。

为了更好地理解二十四节气的天文位置，首先要解释几个天文学概念。“天球”是人们为了研究天体的位置和运动规律而引入的一个假象的球体，根据观察点（也就是球心）的位置不同，可分为“日心天球”、“地心天球”等等。图（1）就是天球概念的一个简单示意图：

图（1）天球概念示意图

天文学中常用的一个坐标体系就是“地心天球”，它与地球同心且有相同的自传轴，理论上具有无限大的半径。地球的赤道和南北极点延伸到天球上，对应着天赤道和南北天极点。和地球上用经纬度定为位置一样，天球也划分了经纬度，分别命名为“赤经”和“赤纬”，地球上的经度用的是度（分秒）为单位，赤经以时（分秒）为单位。天空中的所有天体都可以投射到天球上，用赤经和赤纬定为天体在天球上的位置。“黄道（Ecliptic）”是地球绕太阳公转轨道的轨道平面与天球（地心天球）相交的大圆，由于地球公转受月球和其它行星的摄动，地球的公转轨道并不是严格的平面，因此黄道的严格定义是：地月系质心绕太阳公转的瞬时平均轨道平面与天球相交的大圆。黄道和天赤道所在的两个平面并不是重叠的，它们之间存在一个23度26分的交角，称为“黄赤交角”。由于黄赤交角的存在，黄道和天赤道就在天球上有两个交点，这两个交点就是春分点和秋分点。在天球上以黄道为基圈可以形成黄道坐标系，在黄道坐标系中，也使用了经纬度的概念，分别称为“黄经”和“黄纬”。天体的黄经从春分点起沿黄道向东计量，春分点是黄经0度，沿黄道一周是360度，使用的单位是度、分和秒。黄纬以黄道测量平面为准，向北记为0度到90度，向南记为0度到-90度。

黄道平面可以近似理解为地球绕太阳公转的平面，以黄道为基圈的黄道坐标系根据观测中心是太阳还是地球还可以区分为日心坐标系和地心坐标系，对应天体的黄道坐标分别被称为“日心黄经、日心黄纬”和“地心黄经、地心黄纬”。日心黄经和日心黄纬比较容易理解，因为太阳系的行星都是绕太阳公转，以太阳为中心将这些行星向天球上投影是最简单的确定行星位置关系的做法。但是人类自古观察太阳的周年运动，都是以地球为参照，以太阳的周年视运动位置来计算太阳的运行轨迹，使用的其实都是地心黄经和地心黄纬，要了解古代历法，理解这一点非常重要。图（2）就解释了造成这种视觉错觉的原因：

图（2）太阳黄道视觉位置原理图

古人定义二十四节气的位置，是太阳沿着黄道运行时的视觉位置，每个节气对应的黄道经度其实是地心黄经。从图（2）可以看出日心黄经和地心黄经存在180度的转换关系，同样可以理解，日心黄纬和地心黄纬在方向上是反的，因此可以很方便地将两类坐标相互转换，转换公式是：

太阳地心黄经 = 地球日心黄经 + 180° （3.1式）

太阳地心黄纬 = -地球日心黄纬 （3.2式）

了解了以上的天文学基础之后，就可以着手对二十四节气的发生时间进行计算。我们常说的节气发生时间，其实就是在太阳沿着黄道做视觉运动过程中，当太阳地心黄经等于某个节气黄经度数时的那个瞬间的时间。所谓的用天文算法计算二十四节气时间，就是根据牛顿力学原理或开普勒三大行星定律，计算出与历法密切相关的地球、太阳和月亮三个天体的运行轨道和时间参数，以此得出当这些天体位于某个位置时的时间。这样的天文计算需要计算者有扎实的微积分学、几何学和球面三角学知识，令广大天文爱好者望而却步。但是随着VSOP-82/87行星理论以及ELP-2000/82月球理论的出现，使得天文计算变得简单易行，本文就是以VSOP-82/87行星理论为计算依据，计算二十四节气的准确时间。

古代天文学家在对包括地球和月亮在内的行星运行轨道精确计算后发现，天体的运行因为受相近天体的影响，并不严格遵循理论方法计算出来的轨道，而是在理论轨道附近波动。这种影响在天文学上被称为摄动，摄动很难被精确计算，只能根据经验估算。但是经过长期的观测和计算，天文学家发现行星轨道因为摄动影响而产生的波动其实也是有规律的，即在相当长的时间内呈现出周期变化的趋势。于是天文学家开始研究这种周期变化，希望通过一种类似曲线拟合的方法，对一些周期计算项按照某种计算式迭代求和计算代替积分计算来模拟行星运行轨迹。这种计算式可以描述为：a + bt + ct² + … xcos(p + qt + rt² + …)，其中t是时间参数，这样的理论通常被称为半解析（semi-analytic）理论。其实早在十八世纪，欧洲学者Joseph Louis Lagrange就开始尝试用这种周期项计算的方法修正行星轨道，但是他采用的周期项计算式是线性方程，精度不高。

1982年，P.Bretagnon公开发表了VSOP行星理论（这个理论的英文名称是：Secular Variations of the Planetary Orbits，VSOP的缩写其实是源于法文名称：Variations Séculaires des Orbites Planétaires），VSOP理论是一个描述太阳系行星轨道在相当长时间范围内周期变化的半分析（semi-analytic）理论。VSOP82理论是VSOP理论的第一个版本，提供了对太阳系几大行星位置计算的周期序列，通过对周期序列进行正弦或余弦项累加求和，就可以得到这个行星在给定时间的轨道参数。不过VSOP82由于每次都会计算出全部超高精度的轨道参数，这些轨道参数对于历法计算这样的民用场合很不适用。1987年，Bretagnon 和 Francou 创建了VSOP87行星理论，VSOP87行星理论不仅能计算各种精密的轨道参数，还可以直接计算出行星的位置，行星位置可以是各种坐标系，包括黄道坐标系。VSOP87行星理论由6张周期项系数表组成，分别是VSOP87、VSOP87A、VSOP87B、VSOP87C、VSOP87D和VSOP87E，其中VSOP87D表可以直接计算行星日心黄经（L）、日心黄纬（B）和到太阳的距离（R），此表计算出的结果适用于节气位置判断。

VSOP87D表包含了三部分数据，分别是计算行星日心黄经的周期项系数表（L表）、计算行星日心黄纬的周期项系数表（B表）和计算行星和太阳距离的周期项系数表（R表）。VSOP87D表有太阳系8大行星的数据，本文的计算只关心与地球相关的数据。L表由L0－L5六部分组成，每一部分都包含若干个周期项系数条目，比如L0表有559个周期项系数条目，L1表有341个条目等等。L表的每个周期项系数条目包含若干个参数，用于计算各种轨道参数和位置参数，计算地球的日心黄经只需要用到其中三个系数。计算所有的周期项系数并不是必须的，有时候减少一些系数比较小的周期项可以减少计算所花费的时间，当然，这会牺牲一点精度。假设计算地球日心黄经的三个系数是A、B和C，则每个周期项的计算表达式是：

A * cos(B + Cτ) （3.3式）

其中τ是儒略千年数，τ的计算公式如下：

τ = (JDE - 2451545.0) / 365250 （3.4式）

JDE是计算轨道参数的时间，单位是儒略日，2451545.0是公元2000年1月1日 12时的儒略日数，关于儒略日的概念，请参考“日历生成算法”的第一篇《中国公历（格里历）》中的说明以及计算方法。以L0表的第二个周期项为例，这个周期项数据中与日心黄经计算有关的三个系数分别是A= 3341656.456，B=4.66925680417，C=6283.07584999140，则第二个周期项的计算方法是：3341656.456 * cos(4.66925680417 + 6283.0758499914 * τ)。对L0表的各项分别计算后求和可得到L0表周期项总和L0，对L表的其它几个部分使用相同的方法计算周期项和，可以得到L1、L2、L3、L4和L5，然后用用3.5式计算出最终的地球日心黄经，单位是弧度：

L = (L0 + L1 * τ+ L2 * τ² + L3 * τ³ + L4 * τ⁴ +L5 * τ⁵) / 10⁸（3.5式）

用同样的方法对地球日心黄纬的周期项系数表和计算行星和太阳距离的周期项系数表计算求和，可以得到地球日心黄纬B和日地距离R，B的单位是弧度，R的单位是天文单位（AU）^[1]。由于3.5式的计算方法需要多次计算τ的乘方，浮点数的乘方计算的速度比较慢，实际计算时，通常对3.5式进行变换，用乘法和加法代替直接的乘方计算，这是一种常用的转换：

L = (((((L5 * τ + L4) * τ + L3) * τ + L2) * τ + L1) * τ + L0) / 10⁸（3.6式）

本文就是使用3.6式代替3.5式进行计算。

VSOP82/87行星理论中的周期项系数对不同的行星具有不同的精度，对地球来说，在1900-2100年之间的200年跨度期间，计算精度是0.005"。前文曾说过，对于不需要这么高精度的计算应用时，可以适当减少一些系数比较小的周期项，减少计算量，提高计算速度。Jean Meeus在他的《天文算法》一书中就给出了一套精简后的VSOP87D表的周期项，将计算地球黄经的L0表由原来的559项精简到64项，计算地球黄纬的B0表甚至被精简到只有5项，从实际效果看，计算精度下降并不多，但是极大地减少了计算量。

使用VSOP87D周期项系数表计算得到的是J2000.0平黄道和平春分点（mean dynamic ecliptic and equinox）为基准的日心黄经（L）和日心黄纬（B），其值与标准FK5系统略有差别，如果对精度要求很高可以采用下面的方法将计算得到的日心黄经（L）和日心黄纬（B）转到FK5系统^[2]：

首先然后L'，单位是度：

L' = L - 1.397 * T - 0.00031*T² （3.7式）

3.7式中的T是儒略世纪数，它与儒略千年数τ的关系是：T = 10 *τ。然后使用L'计算L和B的修正值ΔL和ΔB：

ΔL = -0.09033 + 0.03916 * ( cos(L') + sin(L') ) * tan(B) （3.8式）

ΔB = +0.03916 * ( cos(L') - sin(L') ) （3.9式）

ΔL和ΔB的单位都是"，是度分秒角度单位体系，需要将3.6式计算出得L和B转换成以度（°）为单位的值后再进行修正。

CalcSunEclipticLongitudeEC()函数就是使用VSOP87行星理论计算行星日心黄经的代码实现，整个计算过程和前文描述一样，首先根据VSOP87D表的数据计算出L0-L5，然后用3.6式计算出地球的日心黄经，3.6式计算出来的单位是弧度，因此转换成度分秒单位，最后使用3.1式将结果转换成太阳的地心黄经：

double CalcSunEclipticLongitudeEC(double dt)

{

double L0 = CalcPeriodicTerm(Earth_L0, sizeof(Earth_L0) / sizeof(VSOP87_COEFFICIENT), dt);

double L1 = CalcPeriodicTerm(Earth_L1, sizeof(Earth_L1) / sizeof(VSOP87_COEFFICIENT), dt);

double L2 = CalcPeriodicTerm(Earth_L2, sizeof(Earth_L2) / sizeof(VSOP87_COEFFICIENT), dt);

double L3 = CalcPeriodicTerm(Earth_L3, sizeof(Earth_L3) / sizeof(VSOP87_COEFFICIENT), dt);

double L4 = CalcPeriodicTerm(Earth_L4, sizeof(Earth_L4) / sizeof(VSOP87_COEFFICIENT), dt);

double L5 = CalcPeriodicTerm(Earth_L5, sizeof(Earth_L5) / sizeof(VSOP87_COEFFICIENT), dt);

double L = (((((L5 * dt + L4) * dt + L3) * dt + L2) * dt + L1) * dt + L0) / 100000000.0;

/*地心黄经 = 日心黄经 + 180度*/

return (Mod360Degree(Mod360Degree(L / RADIAN_PER_ANGLE) + 180.0));

}

Mod360Degree()函数将大于360°或小于0°的值调整到0－360°之间，便于转换显示。CalcPeriodicTerm()函数使用3.3式对一个周期项系数表进行求和计算，可以指定需要计算的周期项数：

double CalcPeriodicTerm(const VSOP87_COEFFICIENT *coff, int count, double dt)

{

double val = 0.0;

for(int i = 0; i < count; i++)

val += (coff[i].A * cos((coff[i].B + coff[i].C * dt)));

return val;

}

同样的方法计算太阳的地心黄纬：

double CalcSunEclipticLatitudeEC(double dt)

{

double B0 = CalcPeriodicTerm(Earth_B0, sizeof(Earth_B0) / sizeof(VSOP87_COEFFICIENT), dt);

double B1 = CalcPeriodicTerm(Earth_B1, sizeof(Earth_B1) / sizeof(VSOP87_COEFFICIENT), dt);

double B2 = CalcPeriodicTerm(Earth_B2, sizeof(Earth_B2) / sizeof(VSOP87_COEFFICIENT), dt);

double B3 = CalcPeriodicTerm(Earth_B3, sizeof(Earth_B3) / sizeof(VSOP87_COEFFICIENT), dt);

double B4 = CalcPeriodicTerm(Earth_B4, sizeof(Earth_B4) / sizeof(VSOP87_COEFFICIENT), dt);

double B = (((((B4 * dt) + B3) * dt + B2) * dt + B1) * dt + B0) / 100000000.0;

/*地心黄纬 = －日心黄纬*/

return -(B / RADIAN_PER_ANGLE);

}

AdjustSunEclipticLongitudeEC()函数根据3.8式计算黄经的修正量，longitude和latitude参数是由VSOP87理论计算出的太阳地心黄经和地心黄纬，单位是度，dt是儒略千年数，返回值单位是度：

double AdjustSunEclipticLongitudeEC(double dt, double longitude, double latitude)

{

double T = dt * 10; //T是儒略世纪数

double dbLdash = longitude - 1.397 * T - 0.00031 * T * T;

// 转换为弧度

dbLdash *= RADIAN_PER_ANGLE;

return (-0.09033 + 0.03916 * (cos(dbLdash) + sin(dbLdash)) * tan(latitude * RADIAN_PER_ANGLE)) / 3600.0;

}

经过上述计算转换得到坐标值是理论值，或者说是天体的几何位置，但是FK5系统是一个目视系统，也就是说体现的是人眼睛观察效果（光学位置），这就需要根据地球的物理环境、大气环境等信息做进一步的修正，使其和人类从地球上观察星体的观测结果一致。

【下篇将介绍修正理论和修正算法，请继续关注】

小知识1：天文单位

天文单位（英文：Astronomical Unit，简写AU）是一个长度的单位，约等于地球跟太阳的平均距离。天文单位是天文常数之一，是天文学中测量距离，特别是测量太阳系内天体之间的距离的基本单位。地球到太阳的平均距离大约为一个天文单位，约等于1.496亿千米。 1976年，国际天文学联会把一天文单位定义为一颗质量可忽略、公转轨道不受干扰而且公转周期为365.2568983日（即一高斯年）的粒子与一个质量相等约一个太阳的物体的距离。当前普遍被接受并使用的天文单位的值是149,597,870,691±30米（约一亿五千万公里）。

小知识2：FK5系统

FK5常用的目视星表系统，又称第五基本星表，是在FK4（第四基本星表）的基础上发展出来的，对FK4星表进行了修正，于1984年正式启用。它定义了一个以太阳质心为中心，J2000.0平赤道和春分点为基准的天球平赤道坐标系。近年来国际上又编制了FK6星表（第六基本星表），但是还没有被正式启用。