如果你跟长时间跨度的天体测量数据较过劲大概率见过这样一种曲线一颗已知行星的轨道半径围绕平均值小幅摆动幅度不大却极其规律。这种摆动的来源就是摄动——来自某个看不见的天体对它的引力拉扯。天文学上最经典的处理方式是反过来利用这份已知的摄动去推算未知的行星轨道半径、公转周期、甚至质量。初中课本会告诉你海王星是被笔尖算出来的但真正上手做一遍从摄动反推行星的完整流程才明白这活儿有多微妙。今天这篇博文我就用一道天体力学习题笔记里的 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 AUa_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 减去均值消除直流分量然后做实数 FFTdt_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), ddt_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 AUmP_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 角秒/世纪的进动不是行星摄动而是牛顿引力理论本身的修正。这颗不存在的行星提醒我用摄动反推未知天体有一个默认前提你用来拟合的动力学框架是完备的。信号解不出来未必是数据不够也可能是理论这个地方漏了一块。如今的系外行星巡天里也有类似教训——有些恒星的周期性视向速度信号最后被证明是恒星磁活动斑点造成的而不是行星。所以面对一条来路不明的摄动信号我会强迫自己按顺序做三件事先确认周期信号在多个观测波段、多种时间跨度下都稳定再尝试用已知框架和多天体交叉验证建立自洽解最后如果所有尝试都撞墙才回头审视理论假设本身。一道习题能写到这个深度我觉得挺值。它不只是教你套两个公式而是在训练一种从痕迹反推源头的思维习惯。如果你也在啃天体力学或者处理系外行星数据建议花一个下午把这个数值实验完整跑一遍亲手感受一下看不见的东西被算出来那一刻的微妙快乐。
