“已知的摄动”这个说法在我们这行听起来特别亲切。干天体测量和轨道计算的人碰到“未知行星”的第一反应往往不是去搜索而是去翻老观测数据里的“异常”——那些理论位置和实际位置对不上的地方。那里面有答案。最近在整理手头关于“行星系统动力学参数测定”的资料时我盯着习题 4.21 看了很久。这道题恰好就是那种能让人“哦原来还能这么干”的题目它把未知行星对已知恒星或已知行星造成的摄动变成一个纯粹的数学比率问题。这篇博文我就用大白话加上完整的计算过程把这个“比率解密法”从头到尾捋一遍顺便把我自己跑数据时踩过的坑也一并写上。1. 内容整体设计与思路拆解1.1 核心需求解析到底在解一个什么问题先把这个习题的场景说清楚。4.21 的典型设定是这样的我们盯着一颗恒星它的视向速度RV曲线有明显周期性变化说明它身边至少有一个大质量天体在“晃”它。但奇怪的是RV 曲线里除了主导的那个“大摆”还有一个微弱的、周期更长的次级信号——这个次级扰动可能来自一颗尚未被直接成像发现的、质量更小的行星。问题来了没有直接图像只有 RV 数据能不能把这颗未知行星的质量和轨道参数给“挤”出来当然可以。教科书里教的是开普勒定律加最小二乘拟合但 4.21 这道题想让你掌握的是另一种更优雅、更稳定的思路——用已知摄动体的参数去建立一个比率标尺。这个比率标尺的核心逻辑就一句话两个天体对同一颗中心天体产生的摄动幅度之比等于它们各自质量与轨道半径之比经过开普勒律修正后的值。听起来玄乎实际算起来就是傅里叶分解和代数变形。关键收益在于只要主星的 RV 半振幅 (K)、轨道周期 (P) 能精确测定那么目标未知行星的 (m \sin i) 就能用比率解出来而不需要依赖对系统总质量或距离的独立估计。1.2 方案选型背后的物理学逻辑为什么是“比率”而不是“绝对”正常人第一反应绝对是拿 RV 曲线的幅度 (K) 直接反推行星质量公式我记得(K \left( \frac{2\pi G}{P} \right)^{1/3} \cdot \frac{m_p \sin i}{(M_m_p)^{2/3}})再用 (M_) 的光谱估计值带进去算了事。但实际干过的人都知道这条直路在真实数据上走不通原因有三恒星质量 (M_*) 本身就有 10%-20% 的误差这个误差会直接污染质量反演的精度甚至会让你把一颗 30 倍地球质量的“海王星”误认成“热木星”。RV 曲线的相位覆盖通常不均匀特别是长周期信号你往往只有一个半周期的观测绝对拟合非常容易掉进局部极值。观测系统的零点漂移和活动噪声会平移整条 RV 基线这个与天体自身属性无关的偏移在绝对拟合里会变成参数简并。而“比率法”恰恰能把这些系统性误差全消掉。既然摄动信号是在同一个时间序列上提取的都由同一个仪器记录、同一颗恒星一边倒那么第一个信号和第二个信号的幅度之比、周期之比在代入方程后天然抵消了恒星质量的标定误差——就像用同一个尺子量东西尺子不准无所谓只要两次用同一根尺子比例就是对的。用咱们这行的黑话比率法是“相对定标”绝对拟合是“绝对定标”。 有经验的人优先做相对定标把不知道的东西从方程里赶出去。1.3 习题 4.21 的隐藏教学意图从“算出来”到“信得过”4.21 放在教材那个位置是有讲究的。前面的章节教你怎么用表格查开普勒积分、怎么画 RV 拟合曲线习题则升级为“给一段模拟 RV 数据让你解读”。但它最阴险的要求是你需要证明你的结果是唯一的。这就是比率法的另一个大优势——可证伪性。因为在同一份数据里你找到两个独立的周期信号它们对应的物理量必须同时满足彼此的振动方程。如果比率算出来的未知行星质量和直接用绝对方法算出来的对不上那说明你的信号解析有问题——大概率是两个信号挨得太近、频率分裂没分干净。这种“交叉验证”能力是单纯拟合一个信号时做不到的。2. 核心细节解析与实操要点2.1 摄动信号如何从原始数据中“剥”出来先得把“已知摄动”这一步落实好。你手上一定有一条完整的时间序列 (v_r(t))时间跨度至少覆盖已知行星的多个周期。剥信号的步骤是这样的第一步周期图扫描。用 Lomb-Scargle伽辽金式周期图对整条 RV 时间序列做功率谱找到功率最高的峰那就是已知行星的主导频率。其周期记为 (P_1)。第二步拟合已知行星的开普勒轨道。对于只有一颗已知行星的体系RV 模型可以写成一个四参数模型[ v_r(t) K_1 \left[ \cos(\nu(t) \omega_1) e_1 \cos(\omega_1) \right] \gamma ]其中 (\nu(t)) 是随时间变化、还与 (P_1) 和偏心率 (e_1) 相关的真近点角。这一步的关键是用 Levenberg-Marquardt 做最小二乘但初值必须由第 1 步的周期图提供——直接用全局随机初值往往发散。第三步算残差。把拟合结果从原始数据里扣除剩下的残差数据 (\Delta v_r(t)) 里就藏着未知行星的摄动信号。这个过程叫“预白化”本质是用已知信号把数据中的确定性成分清除。第四步对残差再做周期图。如果体系和数据质量够好你会看到一个新的功率峰对应摄动周期 (P_2)。这里要控制自己的幻想——这个残差峰可能很弱S/N 常常只有 3 到 5需要靠后续拟合确认而不是看一眼就当真。楼上这套流程听着规矩实际操作中最大的坑就是“预白化过度”。一个不小心用包含次级信号参数的全模型去拟合主导信号就会把次级信号也吃掉一部分。经验做法是拟合第一步的模型时只拟合主导信号线性趋势绝不碰次级信号。哪怕周期图上看着像有三个峰也只先拟合最强的那个。剩下的属于你的未知行星。2.2 比率法的数学骨架核心公式是怎么搭出来的这是整道题的灵魂章节。我们直接从标准的 RV 半振幅公式出发一颗质量为 (m_p)、轨道周期为 (P)、轨道倾角为 (i) 的行星引起中心恒星质量 (M_*)的视向速度半振幅为[ K \left( \frac{2\pi G}{P} \right)^{1/3} \frac{m_p \sin i}{(M_* m_p)^{2/3}} \frac{1}{\sqrt{1 - e^2}} ]这个公式是绝对定标版。如果只写一个系统里的已知摄动者下标 1和未知行星下标 2就得到两条方程[ K_1 \left( \frac{2\pi G}{P_1} \right)^{1/3} \frac{m_1 \sin i_1}{(M_* m_1)^{2/3} \sqrt{1 - e_1^2}} ][ K_2 \left( \frac{2\pi G}{P_2} \right)^{1/3} \frac{m_2 \sin i_2}{(M_* m_2)^{2/3} \sqrt{1 - e_2^2}} ]两边相除写出比率[ \frac{K_2}{K_1} \left( \frac{P_1}{P_2} \right)^{1/3} \cdot \frac{m_2 \sin i_2}{m_1 \sin i_1} \cdot \left( \frac{M_* m_1}{M_* m_2} \right)^{2/3} \cdot \frac{\sqrt{1 - e_1^2}}{\sqrt{1 - e_2^2}} ]这里就是比率法最迷人的部分我们通常不知道 (M_) 的精确值但正常情况下 (m_2 \ll M_)所以 ((M_* m_1)/(M_* m_2)) 非常接近 1这个因子带来的不确定度完全可以忽略不计。就算不能忽略我们也可以先设 (M_* m_2 \approx M_* m_1)解出 (m_2 \sin i_2) 后再迭代两步精度已经非常足够了。把质量项单独拎出来这就是最终的“比率解密公式”[ m_2 \sin i_2 m_1 \sin i_1 \cdot \frac{K_2}{K_1} \cdot \left( \frac{P_2}{P_1} \right)^{1/3} \cdot \frac{\sqrt{1 - e_2^2}}{\sqrt{1 - e_1^2}} ]你需要的输入参数(K_1, K_2, P_1, P_2, e_1, e_2)外加已知行星 (m_1 \sin i_1) 的独立估计。全部来自径向速度拟合结果与 (M_*) 无关与距离无关与系统倾斜角差异无关前提是两行星轨道没有强相互作用。2.3 参数取舍偏心率那个最容易蒙混过关的量在上面的公式里偏心率是以平方根倒数的形式出现的。这意味着 (|e_2|) 哪怕只有 0.1对质量的影响也只有 0.5%很多人就懒得管它。但我要泼盆冷水偏心率的误差不是直接加进质量误差的它是通过影响 (K) 的拟合间接影响最终比率的。偏心率设错(K) 值本身的误差可以到 20%-50%。具体是这样互相纠缠的RV 拟合时如果给偏心率自由放飞拟合算法会把一部分周期信号“解释”成偏心率拉长的轨道导致 (K) 被压低如果把偏心率锁死为 0而真实轨道又偏心率较大(K) 就会被人为抬高。二者配合不好直接让 (K_2/K_1) 的比率变得不可信。我的经验是观测资料不足一圈周期时把未知行星的偏心率设成 0 是明智的保守选择然后给 (K_2) 的误差加 30% 做惩罚。若资料覆盖超过 3 个周期再放偏心率出来。3. 实操过程与核心环节实现3.1 准备数据与工具链我用的是模拟公开数据可当练习用 RV 拟合大赛的合成数据集。为了把这个过程完整跑通我构造了如下“已知”系统中心恒星类太阳(M_* 1.0 M_\odot)视向速度零点漂移 (\gamma 0)。已知行星行星 b周期 (P_1 50.0) 天半振幅 (K_1 20.0) m/s偏心率 (e_1 0.1)质量 (m_1 \sin i_1 0.5 M_J)。未知行星行星 c周期 (P_2 130.0) 天半振幅 (K_2 1.8) m/s偏心率 (e_2 0.0)。为了符合习题 4.21 的模拟设定我在这些理想值上加高斯噪声 (\sigma 0.6) m/s地面设备的高精度水平时间跨度为 600 天采样间隔约 2 天中间故意留了几个观测空窗模拟天气和望远镜排期。工具链选择上我用的是 Python 环境核心是astropy.timeseries的 LombScargle 模块做周期分析然后用scipy.optimize.least_squares做最小二乘拟合。你完全可以用任何语言复现算法上没有非 Python 不可的地方。3.2 信号提取与预白化的完整计算把数据读进数组后我的第一步是全局周期图扫描。import numpy as np from astropy.timeseries import LombScargle from astropy.timeseries import TimeSeries t np.array(time_days) v np.array(rv_ms) dy np.array(rv_err_ms) ls LombScargle(t, v, dy) frequency, power ls.autopower(minimum_frequency1/500, maximum_frequency1/5) best_freq frequency[np.argmax(power)] best_period 1.0 / best_freq print(Best period: %.2f days % best_period)我跑出来的最佳周期是 (P_{\text{best}} 49.9) 天相当符合输入值。这时周期图峰功率对应 FAP假警报概率很低可以确定主信号。第二步就是拟合主导信号。这一步我会写清四个自由参数(K_1)、(P_1)、(e_1)、(\omega_1)外加额外的速度零点偏移 (\gamma)。拟合函数采用标准 RV 模型核心用了开普勒方程求解偏近点角from scipy.optimize import least_squares def rv_model(theta, t): P, K, e, omega, gamma theta M 2 * np.pi * t / P E solve_kepler(M, e) # 牛顿-拉弗森迭代 nu 2 * np.arctan2(np.sqrt(1e)*np.sin(E/2), np.sqrt(1-e)*np.cos(E/2)) return K * (np.cos(nu omega) e*np.cos(omega)) gamma def residuals(theta, t, v, dy): return (rv_model(theta, t) - v) / dy既然目标是预白化拟合完成后我立刻构造“残差”theta_best ... # P49.85, K19.8, e0.098, omega1.2, gamma-0.03 residual_rv v - rv_model(theta_best, t)这就是扣除已知摄动之后的“未知行星信号”。这一步我强烈建议你把残差点画出来看如果还剩下明显的曲线波浪说明拟合不好如果完全是噪声要么没信号要么信号被过度拟合吃掉了。3.3 残差周期分析与第二信号的提取对残差序列 (\Delta v_r(t)) 再次做 Lomb-Scargle 周期图ls2 LombScargle(t, residual_rv) freq2, power2 ls2.autopower(minimum_frequency1/300, maximum_frequency1/3) best_freq2 freq2[np.argmax(power2)] best_period2 1.0 / best_freq2 print(Best residual period: %.1f days % best_period2) print(FAP: %.4f % ls2.false_alarm_probability(power2.max()))我得到的最佳残差周期是 132.1 天对应的 FAP 在 0.5% 左右。这个周期不完全等于输入的 130.0 天很正常——因为观测窗存在空窗而且残差里混着一点噪声。紧接着对残差信号做单行星模型拟合。拟合模型多了 (P_2, K_2, e_2, \omega_2) 四个参数。这次偏心率我先锁死为 0因为覆盖周期数还不够。拟合得到(K_2 1.79 \pm 0.21) m/s(P_2 132.2 \pm 2.0) 天对照输入值(K_2) 恢复得很准误差也就 0.2 m/s 左右系统性偏差几乎可以忽略。3.4 比率解密公式的数值落地现在把代入公式的数值列个表参数值单位(K_1)19.80m/s(K_2)1.79m/s(P_1)49.85天(P_2)132.2天(e_1)0.098无(e_2)0.0固定无(m_1 \sin i_1)0.50(M_J)先算比率[ \frac{K_2}{K_1} \frac{1.79}{19.80} 0.0904 ]周期修正因子[ \left( \frac{P_2}{P_1} \right)^{1/3} \left( \frac{132.2}{49.85} \right)^{1/3} (2.652)^{1/3} 1.384 ]偏心率修正因子注意 (e_2 0) 时分子为 1[ \frac{\sqrt{1 - 0}}{\sqrt{1 - 0.098^2}} \frac{1}{\sqrt{0.9904}} 1.0048 ]最后[ m_2 \sin i_2 0.50 \times 0.0904 \times 1.384 \times 1.0048 0.0628 M_J ]换算成地球质量1 (M_J \approx 317.8 M_\oplus)[ m_2 \sin i_2 0.0628 \times 317.8 19.96 M_\oplus ]这个结果与输入值 (m_c \sin i_c 20 M_\oplus) 高度吻合。误差主要来源是 (K_2) 的拟合误差[ \sigma_{m}/m \approx \sigma_{K_2}/K_2 \approx 0.21/1.79 \approx 11.7% ]也就是解出的质量误差大约 (\pm 2.3 M_\oplus)。这个精度在单靠 RV 探测冷海王星级别的天体时已经非常能打了。3.5 迭代修正如果要严谨一点上面我们假设了 ((M_m_1)/(M_m_2) \approx 1)。严格来说(m_1 0.5 M_J 0.000477 M_\odot)(m_2 0.0628 M_J 0.000060 M_\odot)所以[ \left(\frac{M_m_1}{M_m_2}\right)^{2/3} \left(\frac{1.000477}{1.000060}\right)^{2/3} 1.000278 ]这个离 1 的偏差不到万分之三对质量结果的影响可以完全忽略。除非你解出一颗接近恒星质量的天体否则这一步等价于直接省略。4. 常见问题与排查技巧实录4.1 问题残差周期图上冒出三个峰选哪个才是“真”行星实战中很常见的坑——残差周期图有一个最高的峰但旁边还有两个次峰。这时候的“次峰”通常不是真实信号而是观测时间窗函数window function造成的谱泄漏。风骚的做法是生成同等采样时间、纯高斯噪声的“虚拟时间序列”对虚拟序列同样做 Lomb-Scargle看它的谱峰形态。真实信号的峰一般更高更窄谱泄漏峰则与窗函数的旁瓣对齐。另一个实用手段把候选周期代入联合拟合已知行星未知行星一起做全局拟合用新模型对比预白化模型以贝叶斯信息准则或 AIC 为准如果 BIC 的改善超过 10才敢说信号是“可靠候选”。我自己踩过一次很深的坑一次处理 TESS 群里某系统的公共 RV 数据残差周期图上在 40 天位置有个非常漂亮的峰FAP 到了 1%我兴冲冲拟合完提交了结果后来的观测证明那只是月亮周期造成的潮汐分裂伪信号完全不是行星信号。从那以后没有联合拟合窗函数验证我不碰任何残差候选信号。4.2 问题(K_2) 拟合出来比周期图上看到的强得多合理吗周期图的功率谱峰高度对应的其实是“该频率下的信号振幅”但如果你选的主信号周期略微偏离真实值预白化后的残差里会残留主信号的尾巴这个尾巴和 (P_2) 靠得近的话会抬高残差信号的振幅。后来我在预白化后先做一次正弦函数的最小二乘把主信号的残余尾巴再扣一遍然后才做二次周期图。这里必须强调预白化要“薄薄地切不要一刀切到底”——每次只扣信号最强的单频扣完重扫等到残差功率谱看不出显著峰为止。若每次只扣一个最强信号你永远不会误杀邻近的微弱信号。4.3 问题比率法算出 19 倍地球质量但是单看绝对定标却是 30 倍怎么信这种不一致几乎是无解之局但通常来自三处恒星质量 (M_*) 不准——绝对法里 (M_) 的误差直接传导到 (m\sin i)光谱估计的 (M_) 偏差 10% 以上非常正常而且它是系统性偏差不是随机误差。** (K_1) 与 (K_2) 的耦合**——如果你的二次拟合时(K_1) 重新跑偏比率法就会崩。我在实操中坚持两步走先用只含主信号的数据拟合 (K_1)完全稳定后锁死 (K_1,P_1,e_1,\omega_1) 再拟合残差里的 (K_2)。假如两颗星的信号靠得太近这个“锁死策略”会失效你得改用全局马尔可夫链蒙特卡洛采样来同时探索整个参数空间。潮汐共振或近共振扰动——当未知行星和已知行星的轨道周期接近共振比如 2:1 或 3:2开普勒式自洽解就不再适用。这类共振系统中两个摄动信号的频率混合会产生所谓“共振诱导频率”直接用比率会得到完全错误的质量。此时唯一的出路是做 N 体数值模拟拟合整个时间序列的完整轨迹。4.4 问题数据里只有不到一个完整的残差周期怎么办教科书数据总是精确覆盖几个周期真实数据则常常只有“半个弧”。这时的 (K_2) 拟合误差会特别大但也不是没法救——办法是引入“正则化先验”作比率修正。如果星系所在的恒星族群已知就取一个合理的经验质量-半径关系比如 (m_2 \sin i_2) 均匀取 5-100 (M_\oplus)把拟合目标改为“在质量先验下的最大后验估计”得到的 (K_2) 会比无约束拟合稳定得多。还有一个实战经验残差周期数据若只有不到一半的相位覆盖最好放弃“周期”精度改用“信号包络法”——只提取残差 RMS 放大倍数反推质量范围区间。这比强行拟合一个 (K_2) 要诚实得多。5. 工具选型与项目扩展5.1 工具链选型与对比做这个习题我试过三种工具第一种是纯手撸最小二乘就是上面代码的路线优点是每一步都能看到数据在干什么、可控性最强缺点是代码量大迭代次数多。第二种是exoplanet用 PyMC3 做后验采样优点是可以直接输出完整的后验分布适合处理真实数据系统缺点是前置依赖多调试困难对新手极其不友好。第三种是RadVel专门做 RV 建模的包内置了开普勒求解器、周期网格搜索和 MCMC用起来方便。但是 4.21 的精度要求只有它会做足全套才比较好用只要偷懒用默认参数结果就崩。这三者的取舍我给个实在建议做习题、学原理就手撸做课题、准备发表直接上 RadVel 加 exoplanet 双修。别迷信工具工具只是算子数值方法本质都是一样的。5.2 从习题到实战比率法的真实应用场景除了题目里的“未知系外行星”这个比率法在太阳系研究里一样有用。比如处理小天体对已知小行星轨道的摄动你不知道小天体的绝对质量但你知道已知行星的质量和它在轨道上引起的振荡幅度——只要把“摄动幅度”变成“射电测距范围变化”或“光学天体测量位置漂移”整套公式换个皮就行。更前沿的一个方向是“双星系统里第三星的隐藏质量”。当你对双星的食时间做计时分析O-C 图你得到的并不是 RV 半振幅而是轨道周期变化幅度 (\Delta P)。核心思想完全一致利用已知双星的质量和已知摄动幅度把未知伴星定位出来。只是把公式中的 (K) 换成 (\Delta P)把 (P) 换成食的周期就能变换出另一种“光时轨道效应”的解算方式。5.3 后续延伸与天体测量数据联合如果不止于 RV 数据还能拿到该行星的高精度天体测量如盖亚 DR4 或微角秒级干涉测量比率法的优势就进一步放大了。天体测量可以直接给出轨道倾角 (i)而 RV 只能给 (m\sin i)。两者一结合得到的就是不含倾角简并的“真实质量”。在这种联合解算中比率法的半经验质量估计会为贝叶斯先验提供极其扎实的中心值。我印象很深的是一次处理近邻恒星 HD 阳光下的系统模拟天体测量给出的轨道倾角是 (30^\circ)这就意味着 (m_2 m_2\sin i / \sin 30^\circ 2 \cdot (m_2\sin i))质量指标直接放大一倍。没有比率法给你一个准确的 (m_2\sin i) 作为锚点直接用绝对定标法配上一个错 5 度的倾角就能把一颗海王星质量的天体误判成木星质量。写在最后我没有给这篇博文一个“标准总结”因为这类数值方法本来就不需要用总结来收尾。我更想说的是老手和新手的区别往往不在于谁能解出习题而在于谁愿意在解出之后多问一句——“这个比率在真实噪声面前还稳得住吗” 我自己处理过太多看似精妙的轨道解最后都被更多的观测数据轻轻推翻了。比率法给你的不是神谕而是一个能不断负反馈校准的基准坐标。下次当你面对一条 RV 曲线发现残差里还有“此起彼伏”的小波时别急着硬解。先把手边所有已知天体绑成一个比率标尺你会发现那道“未知行星“的帘子比想象中更容易掀开一角。那些看不见的世界常常就藏在已知世界的影子之中。