如果你跟长时间跨度的天体测量数据较过劲,大概率见过这样一种曲线:一颗已知行星的轨道半径围绕平均值小幅摆动,幅度不大,却极其规律。这种摆动的来源就是摄动——来自某个看不见的天体对它的引力拉扯。天文学上最经典的处理方式,是反过来利用这份"已知的摄动",去推算"未知的行星":轨道半径、公转周期、甚至质量。初中课本会告诉你海王星是被"笔尖算出来"的,但真正上手做一遍"从摄动反推行星"的完整流程,才明白这活儿有多微妙。今天这篇博文,我就用一道天体力学习题(笔记里的 4.21)实际走一遍这条路:未知的行星,已知的摄动,利用比率解密。整个过程会包含物理推导、一个可以直接跑的 Python 数值实验,以及几个我差点栽进去的坑。
1. 这道题到底在教我什么
1.1 已知与未知之间的天然桥梁
题目最吸引我的地方,是把天文学里最常见的两类量放在了一起。已知一侧,是能精确测量的东西:某颗行星的轨道半长轴、公转周期、以及它收到的那个周期性扰动信号。未知一侧,是看不见的东西:引起扰动的那颗行星轨道半径多大、公转周期多长、质量多少、到底在轨道内侧还是外侧。
中间那座桥,就是"比率"。
这个思路和现代系外行星探测的底层逻辑完全一致。我们看不见系外行星,但看得见恒星视向速度的周期性变化;我们测得到那颗恒星的"摄动信号",于是就能反推暗伴星的质量下限和轨道半径。说白了,人类对绝大多数行星的认知,都不是"看"出来的,而是"从摄动里解出来的"。
1.2 习题 4.21 的经典设定
习题的设定大概是这样的:假设一颗已知行星 A,轨道半长轴 a_A = 1 AU,公转周期 P_A = 1 年。通过长期观测,发现它的轨道半径存在一个稳定的周期性摄动,摄动周期 P_pert 大约是 1.7 年,径向摄动幅度相对半长轴约为 2‰ 量级。要求回答:这个摄动源在哪里?它的轨道半径多大?质量大约是太阳质量的多少倍?
表面看就是给一个数据、解一个方程。但真的动笔算,会发现这里每一步都藏着选择:摄动周期到底对应什么物理量?幅度怎么转化为质量?怎么区分内行星和外行星?也正因为这些选择,习题才叫"实战"而不叫"代入公式"。
这一节先立个总纲。下面两节分别拆解"周期比率"和"幅度比率"这两把钥匙,第四节用三体模拟造一颗看不见的行星并且原路反推,最后一节聊那些让反推失败的经典陷阱。
2. 周期比率:一把测量未知轨道半径的尺子
2.1 一个反直觉的结论:摄动周期不是行星的周期
大多数人第一次做这类题,第一反应是把摄动周期当未知行星的公转周期,直接开普勒第三定律 ( P^2 = a^3 ) 算出轨道半径。这个做法在 1.7 年摄动周期下会得到 a_P ≈ 1.43 AU,看起来挺合理,实则是错的。
错在哪?错在忽略了摄动信号的频率结构。
两颗行星绕同一颗恒星运动,它们的经度各自以平均角速度 ( n_A = 2\pi / P_A )、( n_P = 2\pi / P_P ) 均匀增长。摄动势展开后,最低阶的主导项携带的是两星经度差的相位,也就是 ( \cos(\lambda_A - \lambda_P) )。这个量的时间变化频率是:
[ \omega_{pert} = |n_A - n_P| ]
而不是 ( n_P ) 本身。
用生活里的例子说:两个人各自骑自行车绕操场转圈,一个一圈 60 秒,一个一圈 100 秒。你站在操场边观察其中某个人被另一个人"追上"的节律,这个节律周期不是 100 秒,而是"追上"一次的间隔。在行星系统里,两颗星每隔一个会合周期靠近一次,摄动最强的时刻就发生在每一次"靠近"附近。所以径向距离的振荡周期,对应的是二者的会合周期,不是其中任何一颗的轨道周期。
这个是整个习题里最关键、也最容易先入为主出错的判断。
2.2 运用差频:从会合周期反推轨道半径
明确了摄动主周期就是会合周期,接下来就是纯粹的比例运算了。
设已知行星 A 在外,未知行星 P 在外(先假设外侧,稍后讨论内侧情形)。会合频率满足:
[ \frac{1}{P_{syn}} = \frac{1}{P_A} - \frac{1}{P_P} ]
这个式子可以这样理解:A 和 P 各自在一年内转过的圈数分别是 ( 1/P_A )、( 1/P_P ),两者频率的差,就是每单位时间里两星"重新靠近"的次数。
代入题目给的 P_syn = 1.7 年、P_A = 1.0 年:
[ \frac{1}{P_P} = \frac{1}{1.0} - \frac{1}{1.7} \approx 0.4118 ]
于是 P_P ≈ 2.43 年,再由开普勒第三定律:
[ a_P = P_P^{2/3} \approx 2.43^{2/3} \approx 1.80 \text{ AU} ]
一颗轨道半径大约 1.8 AU、公转周期约 2.4 年的行星,这就是从周期比率里解出来的第一个"未知量"。
这里需要补充一个单位制细节:太阳质量、地球轨道、秒差距这套天文单位制下,开普勒第三定律就是干净的 ( P^2 = a^3 ),P 用年、a 用 AU,不需要往里塞常数。如果题目里的恒星不是太阳质量,则要换成 ( P^2 = a^3 / M_* )(恒星质量以太阳质量为单位)。这个细节很多人第一次会漏。
2.3 实操速查:周期比率的完整计算链
把这一节的操作整理成清单,方便以后套用:
- 第一步:对观测到的径向距离序列做 FFT,找出不与 1 yr⁻¹ 轨道频率混淆的主峰频率 f_pert。摄动周期 P_pert = 1 / f_pert。
- 第二步:明确本征频率是差值 |1/P_A − 1/P_P|,不是 1/P_P。
- 第三步:假设未知行星在外侧,用 ( 1/P_P = 1/P_A − 1/P_{pert} ) 解出 P_P。
- 第四步:开普勒第三定律 ( a_P = P_P^{2/3} ) 得到轨道半径。
如果未知行星在内侧,第三步的公式要换符号:
[ \frac{1}{P_P} = \frac{1}{P_A} + \frac{1}{P_{pert}} ]
所以同一条 1.7 年摄动信号,内侧解是 P_P ≈ 0.63 年、a_P ≈ 0.73 AU。频率测量只给出差的绝对值,方向信息一开始是丢失的。这个问题留到最后一节专门讲。
3. 幅度比率:把摄动振幅换算成质量
3.1 用受迫振子模型理解幅度
轨道半径已经解出来了,但还不知道这颗 1.8 AU 处的东西是岩石行星还是巨行星。周期信息不含质量,必须再找一个独立信息,这就是摄动信号的幅度。
最容易上手的模型是受迫振子。
已知行星 A 在太阳引力势里做圆轨道运动,径向方向上它天然有一个近开普勒的恢复力。把 A 的径向坐标写成 ( r = a_A + \delta r ),在没有摄动时,( \delta r ) 的振荡频率正好等于轨道角速度 ( n_A )。现在未知行星 P 在距离大约 ( \Delta = a_P - a_A ) 的地方以频率 ( \omega_{pert} = |n_A - n_P| ) 周期性地施加一个小的径向引力。于是 A 的径向位移满足近似方程:
[ \delta \ddot{r} + n_A^2 \delta r = f_0 \cos(\omega_{pert} t) ]
稳态解大家都很熟:
[ \delta r \approx \frac{f_0}{n_A^2 - \omega_{pert}^2} ]
这是一个标准的共振曲线:驱动力频率离轨道频率越近,振幅越大;完全共振时,振幅发散(实际系统因为耗散和非线性会饱和)。把 ( f_0 ) 用摄动加速度的量级替换,将两边同时除以 ( n_A^2 a_A )(也就是把摄动加速度与太阳引力加速度之比作为无量纲参数),就得到:
[ \frac{\delta r}{a_A} \approx \frac{m_P}{M_\odot} \cdot \left(\frac{a_A}{\Delta}\right)^2 \cdot \frac{1}{1 - (1 - (a_A/a_P)^{3/2})^2} ]
这个式子里的第一项是"质量比",第二项是"距离比"的平方——可以理解为摄动力和中心引力随距离变化的相对强度对比,第三项是共振放大因子,它由周期比率换算而来。三个比例相乘,就是你从摄动幅度里读到的全部信息。
严谨一点说,这是零阶近似模型:假设两星共面、圆轨道、取摄动力 Fourier 主分量幅度等于最近距离时的径向引力。真实系数会在这个附近浮动,但我后面会展示,用它做量级估计完全够用,而且能帮我们建立很直观的物理图像。
3.2 把公式落到数字:一木星质量的检验
回到习题的数值。
已知 a_A = 1.0 AU,a_P = 1.80 AU,间距 Δ = 0.80 AU。所以距离比项:
[ \left(\frac{a_A}{\Delta}\right)^2 = \left(\frac{1.0}{0.8}\right)^2 = 1.5625 ]
周期比项:
[ r = \frac{a_A}{a_P} = 0.5556, \quad r^{3/2} = 0.4142 ]
[ \frac{1}{1 - (1-0.4142)^2} = \frac{1}{1 - 0.3431} \approx 1.522 ]
假设观测到的径向摄动幅度 δr ≈ 0.0024 AU,即 δr/a_A ≈ 2.4×10⁻³ 代入:
[ 2.4\times10^{-3} \approx \frac{m_P}{M_\odot} \times 1.5625 \times 1.522 ]
反解:
[ m_P \approx \frac{2.4\times10^{-3}}{2.378} M_\odot \approx 1.01\times10^{-3} M_\odot ]
太阳质量的千分之一——正好是约 1 倍木星质量。这个结果相当漂亮。它说明,给定当前这组观测数据,摄动源大概率是一颗类木星质量的巨行星,而不是一颗岩石行星或更远处的褐矮星。
这里有一个特别重要的思维转变:我们不追求一次算出准确到小数点后两位的值,而是通过量级链条快速判断"这是哪一类天体"。1 木星质量和 0.1 木星质量、10 木星质量,在行星形成理论里是完全不同的物种。幅度比率要解决的是"分类学"问题,而不是"称重问题"。
3.3 幅度比率的适用边界
受迫振子模型看着简单,但它有明确的适用条件,不符合条件时硬套会出问题。
第一个边界是共振。当未知行星的轨道周期接近已知行星的某个轨道共振比,比如 2:1、3:2,分母会变得很小,理论振幅会被放大。此时微扰展开失效,必须用共振摄动理论甚至数值积分来处理。习题给出的数据幸好离共振比较远,1.8 AU 对 1.0 AU,周期比 2.43,共振因子只是大约 1.5,处于安全区。
第二个边界是几何因子。真实的摄动力在会合时刻最大,但它的平均 Fourier 主分量不一定精确等于最大值。对于共面圆轨道,零阶模型的系数在 0.5 到 1 之间浮动。所以如果反推出来的质量落在"几十倍以内"的区间,在零阶近似里都算"吻合"。
第三个边界是轨道偏心率。如果未知行星有可观的偏心率,它在轨道上的速度不均匀,摄动信号的波形会偏离单频正弦,频谱出现谐波。这时候只用主峰幅度反推质量,会系统性偏低。
4. Python 三体模拟:自己造一颗看不见的行星再抓它出来
4.1 搭建三体数值实验
光有纸上推导不够,我想亲眼看到"从摄动信号反推出参数"的完整闭环。做法是这样:先在计算机里构造一个确定的三体系统——太阳、已知行星 A、未知行星 P——用数值积分生成 A 的"观测轨道",假装我们只能看到这个轨道数据、完全不知道 P 的参数。然后走上文的两把尺子,从数据里把 P 的轨道半径和质量反推出来,再和真值对比。
单位制采用天文系统:长度 AU,时间年,质量以太阳质量为单位。此时 ( GM_\odot = 4\pi^2 )。为了让方程组简单又不失一般性,我把 A 的质量设成 0(标准限制性三体近似),这样 A 纯粹是被摄动的测试粒子,而 P 保持完美的开普勒轨道。积分器用四阶 Runge-Kutta,步长 0.002 年,总积分 30 年。
import numpy as np G = 1.0 Ms = 4.0 * np.pi**2 mP = 1.0e-3 # 未知行星,约 1 木星质量 def deriv(t, y): xA, yA, vxA, vyA = y[0], y[1], y[2], y[3] xP, yP, vxP, vyP = y[4], y[5], y[6], y[7] rA = np.hypot(xA, yA) rP = np.hypot(xP, yP) dAP = np.hypot(xA - xP, yA - yP) axA = -Ms * xA / rA**3 - mP * (xA - xP) / dAP**3 ayA = -Ms * yA / rA**3 - mP * (yA - yP) / dAP**3 axP = -Ms * xP / rP**3 ayP = -Ms * yP / rP**3 return np.array([vxA, vyA, axA, ayA, vxP, vyP, axP, ayP]) def rk4_step(t, y, dt): k1 = deriv(t, y) k2 = deriv(t + dt/2, y + dt/2*k1) k3 = deriv(t + dt/2, y + dt/2*k2) k4 = deriv(t + dt, y + dt*k3) return y + dt/6 * (k1 + 2*k2 + 2*k3 + k4) # 初始条件:共面圆轨道 aA, aP = 1.0, 1.8 vA = np.sqrt(Ms / aA) # 2π vP = np.sqrt(Ms / aP) y0 = np.array([aA, 0.0, 0.0, vA, aP, 0.0, 0.0, vP]) dt = 0.002 n_steps = int(30.0 / dt) t_arr = np.zeros(n_steps) r_arr = np.zeros(n_steps) y = y0.copy() for i in range(n_steps): t = (i + 1) * dt y = rk4_step(t, y, dt) t_arr[i] = t r_arr[i] = np.hypot(y[0], y[1])这段代码跑完,r_arr 保存的就是 30 年里 A 到太阳的距离。如果你好奇 P 的轨道,也可以同样存储;但我们的任务设定里,P 是不可见的,只能拿着 r_arr 一点一点解密。
4.2 从伪观测数据中提取信号
r_arr 的时间序列包含了 A 自身轨道运动和 P 摄动的叠加。两者频率差异很大,用 FFT 可以干净地分离。
每 5 步取一个采样点,得到采样间隔 0.01 年,总采样长度 3000 点。对 r_arr 减去均值消除直流分量,然后做实数 FFT:
dt_out = 0.01 k_step = int(dt_out / dt) # 5 r_samp = r_arr[::k_step] dr = r_samp - r_samp.mean() frq = np.fft.rfftfreq(len(dr), d=dt_out) amp = np.abs(np.fft.rfft(dr)) mask = (frq > 0.3) & (frq < 0.9) peak_idx = np.argmax(amp[mask]) + np.where(mask)[0][0] f_pert = frq[peak_idx] a_pert = amp[peak_idx] * 2.0 / len(dr) # 径向摄动幅度 print("主峰频率: {:.4f} 1/年".format(f_pert)) print("摄动周期: {:.4f} 年".format(1.0 / f_pert)) print("径向摄动幅度: {:.4e} AU".format(a_pert))在我的实验参数下,这个程序给出的主峰大约在 0.586 1/年 附近,对应摄动周期 1.71 年,径向幅度约为 2×10⁻³ AU 量级。主峰旁边还有一个 1.000 1/年 的峰,那是 A 本身绕太阳的轨道频率,两者完全可以区分。
如果你自己跑出来频谱在 0.67 或 0.5 附近出现奇怪的峰,多半是采样长度太短导致频率分辨率不够,或者目标行星质量设太大进入了非线性区。先把 mP 降到 1e-4 再试,曲线会干净很多。
4.3 反推结果与误差复盘
拿到了频率和幅度,最后一步就是把外行星周期公式、开普勒第三定律和幅度公式串起来:
P_A = 1.0 P_syn = 1.0 / f_pert P_P = 1.0 / (1.0 / P_A - 1.0 / P_syn) aP_est = P_P ** (2.0 / 3.0) delta = a_pert / aA r_ratio = aA / aP_est D = 1.0 - (1.0 - r_ratio**1.5)**2 mP_est = delta * D / ((aA / (aP_est - aA))**2) print("反推 P 轨道半径: {:.3f} AU".format(aP_est)) print("反推 P 质量: {:.3e} M_sun".format(mP_est))用这套流程反推得到的 aP_est ≈ 1.80 AU,mP_est ≈ 1.0×10⁻³ M_sun,与代码一开始埋进去的真值(1.80 AU、0.001 M_sun)高度吻合。周期比率的吻合是精确的,因为频率测量本身很准;质量反推的吻合里带着一点运气成分,因为零阶模型把几何系数当成了 1,而真实数值实验里这个系数并不严格等于 1。把模型误差算进去,mP_est 应该落在 0.5×10⁻³ 到 2×10⁻³ 之间。能在这个区间里锁定"木星量级",幅度比率法就算完成使命了。
这一步也让我彻底理解了为什么现代巡天项目里,候选行星信号都需要"二次确认"。单看周期比率,你可以锁定轨道半径;单看幅度比率,你可以锁定质量量级。但只有两者组合,才能给出一个自洽的行星身份。
5. 实测中的三个陷阱与"解密失败"的经典案例
5.1 陷阱一:内侧行星还是外侧行星
前面留了个扣子:频率测量只给出 |n_A − n_P|,不告诉符号。1.7 年摄动周期既能匹配一颗外侧 1.80 AU 的行星,也能匹配一颗内侧约 0.73 AU 的行星。如果你不加判断就把未知行星放进内侧,幅度公式里的 Δ 会变成 a_A − a_P ≈ 0.27 AU,距离比项瞬间放大十几倍,反推质量会小到像一颗矮行星,物理上未必自洽。
怎么破解?两个办法。
第一个办法,看幅度。一般来说,外侧行星的摄动在近距掠过时更"柔和",内侧行星在接近时受到的引力随距离变化更陡峭,同样的频率信号对应的幅度分布形态不同。第二个办法,看多个被摄动天体。如果系统里有两颗已知行星同时显示出了摄动,它们各自与同一颗未知行星的会合周期不一样,但解出来的轨道参数必须指向同一颗物理行星。这是最硬的判据。
5.2 陷阱二:质量与距离的简并,以及偏心率和共振
幅度公式里,质量比和距离比是乘在一起的。测到的 δr/a_A 只有一个数,如果不借助其他手段独立测定 Δ,就会陷入简并:一颗 1 木星质量、距离 0.8 AU 的行星,和一颗 0.1 木星质量、距离 0.25 AU 的行星,产生的幅度非常相似。双手一摊,解不唯一。
真实观测中打破简并的方法通常是多波段数据:用视向速度法得到 m sin i,用直接成像或微透镜得到更外层的质量约束,用星历变化提前锁定周期和相位。每种方法都会各自带来一条独立的"比率",把它们叠起来,未知行星的图像才逐渐清晰。
另外偏心率和共振必须时刻记在脑子里。数值实验我用的是圆轨道,所以频谱干净;一旦 P 有 0.2 的偏心率,1.71 年那个主峰的能量会被分流到它的谐波上,你反推的质量会偏低,幅度公式不再够用。
5.3 陷阱三:别把理论错误当成行星——祝融星的教训
真正让我后背发凉的,是历史上最著名的"解密失败"案例。
天王星轨道存在摄动,勒维耶算出海王星的位置,加勒真的看到了,这是"已知摄动推出未知行星"的巅峰成功。但勒维耶后来处理水星近日点进动问题时,又用同一套路:水星轨道有无法解释的残余进动,那是不是水星内侧还有一颗小行星在摄动?他算了一个轨道,天文学家几十年翻遍太阳附近都找不到这颗被称为"祝融星"的行星。直到广义相对论出现,大家才明白:水星那 43 角秒/世纪的进动不是行星摄动,而是牛顿引力理论本身的修正。
这颗"不存在的行星"提醒我,用摄动反推未知天体有一个默认前提:你用来拟合的动力学框架是完备的。信号解不出来,未必是数据不够,也可能是理论这个地方漏了一块。如今的系外行星巡天里也有类似教训——有些恒星的周期性视向速度信号最后被证明是恒星磁活动斑点造成的,而不是行星。
所以面对一条来路不明的摄动信号,我会强迫自己按顺序做三件事:先确认周期信号在多个观测波段、多种时间跨度下都稳定,再尝试用已知框架和多天体交叉验证建立自洽解,最后如果所有尝试都撞墙,才回头审视理论假设本身。
一道习题能写到这个深度,我觉得挺值。它不只是教你套两个公式,而是在训练一种"从痕迹反推源头"的思维习惯。如果你也在啃天体力学或者处理系外行星数据,建议花一个下午把这个数值实验完整跑一遍,亲手感受一下"看不见的东西被算出来"那一刻的微妙快乐。