基于随机森林的TESS变星自动分类：从特征工程到14092颗新发现

1. 项目概述当机器学习遇见百万恒星如果你也处理过TESS凌星系外行星巡天卫星的数据那你一定对那种“幸福的烦恼”深有体会。动辄上百万条光变曲线每一条都代表着一颗恒星在数十天里亮度变化的“心跳”记录。我的工作就是从这片数据的海洋里把那些特殊的“心跳”——变星——给找出来并且分门别类。变星可不是普通的星星。它们的亮度会随时间发生周期性变化原因五花八门可能是两颗星相互绕转、彼此遮挡食双星可能是星体自身像心脏一样有规律地膨胀收缩脉动变星也可能是表面星斑随着自转导致亮度变化自转变星。研究它们就像拿到了恒星内部结构和演化过程的“诊断报告”意义重大。但问题是面对TESS前57个扇区产生的约105万条光变曲线传统的人工“看脸”目视检查分类法彻底行不通了。你算算就算每条曲线只看10秒不眠不休也得干上3000多个小时这还不算分类决策的时间。所以我们转向了机器学习更具体地说是随机森林Random Forest。这个选择不是拍脑袋定的。相比需要海量数据、像个“黑盒子”的深度学习方法随机森林有几个对我们天文数据处理者非常友好的特质它训练速度快能直接处理我们手工设计的物理特征比如傅里叶分解参数、偏度等并且结果相对可解释。在算力有限、且希望理解分类背后物理原因的场景下它往往是更务实、更高效的选择。这个项目的核心目标就是构建一套基于随机森林的自动化流水线对TESS海量数据进行可靠的变星分类并挖掘出尚未被记录的新变星。最终我们不仅系统性地给出了数万颗食双星EA型、EW型、脉动变星造父变星CEP、盾牌座δ型变星DSCT、天琴座RR型变星RRab/RRc和自转变星ROT的名单更重要的是我们从中新发现了14092颗变星。这个数字背后是方法、耐心和无数次调优的结果。接下来我就把这套从数据准备、特征工程、模型训练到后处理筛选的完整流程以及其中踩过的坑和总结的经验毫无保留地分享出来。2. 数据基石TESS光变曲线与Gaia标签工欲善其事必先利其器。整个项目的根基在于高质量的数据和可靠的“教师”标签。2.1 TESS数据105万条恒星“心电图”我们使用的数据全部来自TESS卫星的2分钟短曝光观测。选择2分钟数据而非30分钟的全帧图像FFI主要是为了时间分辨率。对于多数变星尤其是短周期脉动变星2分钟的采样率足以捕捉其亮度变化的细节。我们从MAST数据中心下载了第1至57扇区的所有2分钟目标数据提取其中的PDCSAP_FLUX经过系统误差校正后的流量。这一步很直接但数据量巨大最终我们获得了大约105万条独立的光变曲线。一个重要的决策点在于如何处理不同扇区的重复观测。同一颗星可能被多个扇区观测到。我们的做法是不进行去重而是将不同扇区的观测视为独立的样本。这相当于一种天然的数据增强。因为不同扇区的观测条件、噪声水平、甚至数据间隙都不同同一个目标在不同扇区可能呈现出略有差异的光变曲线这能帮助模型学习到更鲁棒的特征而不是过拟合于某一次特定的观测。2.2 标签来源与Gaia DR3变星表的交叉匹配监督学习需要标注数据。我们的“教师”是Gaia DR3发布的变星表。我们以3角秒的匹配半径将我们的TESS目标与Gaia变星表进行交叉匹配。匹配半径不能太大否则会引入错误匹配也不能太小因为不同星表的坐标系统存在微小差异。3角秒是一个在天文交叉匹配中常用的、较为保守的折中值。匹配后我们获得了9022颗食双星EBs、1528颗脉动变星和3744颗自转变星ROT的初步标签。这里要清醒认识到一点Gaia的标签并非100%准确。它的分类本身也存在误差尤其是对于光变曲线形态相似的类别比如某些EW型食双星和脉动变星。因此我们后续的模型训练和结果验证都必须包含对抗和修正这些标签错误的机制。注意数据一致性检查在开始任何分析前务必对匹配结果做一次快速抽查。随机选取几十个匹配目标在Aladin或类似的星图软件中叠加TESS和Gaia的星像确认是同一颗星。我曾遇到过因为星表坐标系统差或自行运动导致错误匹配的情况这一步能提前避免系统性偏差。3. 特征工程从光变曲线中提取物理指纹光变曲线是一条随时间变化的亮度序列。直接把它扔给模型比如用一维卷积神经网络是一种思路但对于随机森林我们需要手动提取有物理意义的特征。这一步是决定分类成败的关键也是整个项目中最体现天文“手艺”的部分。3.1 周期搜索使用广义Lomb-Scargle方法变星的核心特征是周期性。我们使用Astropy库中的广义Lomb-Scargle (GLS) 方法进行周期搜索。GLS方法对非均匀采样、存在数据缺失的天文时间序列特别友好。我们的策略分两步走粗搜索在0.01到33.33天的范围内进行搜索。上限设为33天略长于一个TESS扇区的27.4天是为了确保能探测到在一个扇区内至少出现一次的周期信号。这一步使用单阶傅里叶项nterms1快速定位功率谱峰值。精搜索在粗搜索找到的最大功率频率f0附近缩小频率范围[f0/1.9, f0*2]并使用5阶傅里叶项nterms5进行精细搜索以得到更精确的周期值。这里有一个针对食双星的重要修正对于食双星GLS周期图经常返回其真实周期的一半因为光变曲线在一个周期内有两个相似的极小值。因此在后续处理所有被初步判定为食双星的目标时我们会将其周期直接乘以2。3.2 傅里叶分解与特征提取得到周期后我们将光变曲线按此周期进行相位折叠以最小流量点为相位0点并归一化到中值流量。接下来进行四阶傅里叶拟合y_fit a0 Σ [ai * sin(2πiφ) bi * cos(2πiφ)]其中 i1 到 4。从拟合系数中我们可以计算出振幅Ai sqrt(ai² bi²)和相位φi arctan(bi/ai)。这些是构建特征的基础砖石。我们主要提取了以下几类特征振幅比与相位差例如A21 A2/A1,φ21 φ2 - 2φ1。不同类别的变星其谐波分量之间的关系不同。比如天琴座RRab型变星的光变曲线不对称其A21通常比对称的RRc型星更大φ21也有显著差异。拟合优度r²衡量傅里叶模型对原始数据的拟合程度。规则变星如脉动变星、部分自转变星的r²值接近1而非变星或噪声主导的曲线r²值很低。统计特征偏度Skewness和峰度Kurtosis。偏度描述光变曲线分布的不对称性食双星的偏度通常较大峰度描述分布的陡峭程度。GLS功率值我们记录了粗搜索和精搜索得到的最大功率值power1和power2。功率值越接近1信号越可能是真实的周期性信号而非噪声。我们最初也尝试引入颜色色指数、有效温度等物理参数但发现TESS数据中这些信息的缺失率很高强行使用会损失大量样本因此最终放弃了。操心得相位折叠的零点选择有些研究建议在折叠前裁剪掉流量最高和最低的1%的数据点以去除离群值。但我们实测发现这对于EA型食双星是灾难性的因为其食甚阶段的数据点很少裁剪会直接抹掉关键信息因此我们选择不裁剪保留完整数据。虽然这会引入更多噪声但通过后续的r²和功率值过滤可以有效控制。4. 随机森林分类器的构建与训练有了特征和标签就可以搭建分类器了。我们选择使用Scikit-learn库中的RandomForestClassifier。4.1 为什么是随机森林在项目初期我们对比过几种方案。深度学习如CNN确实在图像和序列分类上强大但它需要极大的数据量且像一个“黑箱”我们很难理解它到底根据光的什么形状做出了分类决策。这对于追求物理可解释性的天文研究有时是个障碍。而随机森林的优势非常突出训练速度快在同样的硬件上训练一个随机森林模型的时间通常比训练一个深度神经网络少一个数量级。处理高维特征能直接处理我们提取的十几个特征无需复杂的降维。提供特征重要性评估训练完成后我们可以知道哪些特征如φ21、A21、r²对分类决策贡献最大这反过来可以验证我们的特征设计是否合理。对缺失值相对不敏感虽然我们最终去掉了有缺失值的特征但这一特性在处理不完美的天文数据时总是一个加分项。不易过拟合通过构建多棵决策树并集成其泛化能力通常很强。4.2 训练策略分层分类与样本平衡直接用一个模型区分所有7个子类EA, EW, CEP, DSCT, RRab, RRc, ROT效果很差。主要原因是样本不平衡和类间混淆。例如EW和DSCT的光变曲线都是近正弦形极易混淆而ROT的数量远多于其他类别。因此我们设计了一个分层分类流程第一层四大主类分类目标将目标分为食双星EBs、脉动变星Pulsations、自转变星ROT和非变星Non-variable。训练集构建从Gaia匹配的标签中每类选取1000个样本共4000条曲线。对于非变星我们设定一个经验阈值power2 0.01且r² 0.01从大量数据中筛选出最“平静”的曲线作为负样本。关键调整由于Gaia中脉动变星的r²普遍较高我们限定训练集中脉动变星的r² 0.64以确保模型学到的是清晰的脉动信号。但对EBs和ROT不做此限制以保留多样性。模型与结果设置n_estimators700700棵树max_features5。使用袋外估计OOB评估最终模型OOB分数达到0.9178。在测试集上各类别的精确度和召回率均在85%以上见下表。这个模型作为“粗筛”将105万条曲线分到了四个大篮子里。类别精确度召回率非变星92%89%食双星96%94%脉动变星98%98%自转变星85%89%第二层子类精细分类将第一层分出的“篮子”分别进行下一步处理。4.3 食双星EBs子分类EA vs. EW食双星主要分为EA测光双星有平缓的亮度变化和陡峭的食和EW相接双星光变曲线呈连续平滑的正弦状。我们设计了一个简单的统计特征c_bin来区分它们将相位折叠并归一化到[0,1]的光变曲线分成100个相位箱。统计流量值低于0.5的数据点数量。这个值就是c_bin。原理EW型的光变曲线近似正弦波数据点均匀分布c_bin应接近50。EA型的光变曲线大部分时间处于“平顶”的高亮度状态只有食发生时亮度骤降因此流量低于0.5的点很少c_bin值很小。我们统计了所有EBs的c_bin分布发现其在20处存在一个明显的极小值因此选择c_bin 20作为阈值。c_bin 20的判为EA反之为EW。最后再用r² 0.1过滤掉混入的非变星非变星的拟合通常很差r²很低。通过这个方法我们初步得到了6770颗EA和2971颗EW。4.4 脉动变星Pulsations子分类CEP、DSCT、RRab、RRc这是最复杂的一步因为子类多、形态相似。我们的策略是“先分长短再辨形状”。初步按周期划分造父变星CEP周期 1天。TESS一个扇区仅27天所以我们主要发现的是短周期造父变星。盾牌座δ型变星DSCT周期 0.3天。天琴座RR型变星周期介于0.2-1天之间其中RRab型周期通常为0.45-1天RRc型为0.2-0.45天。训练小模型区分DSCT和RR Lyrae由于Gaia标签中DSCT样本很少我们首先用周期阈值初筛DSCT (P0.14d), RRc (P在[0.32, 0.37]d), RRab (P在[0.45, 1]d)各取200颗作为训练集。注意周期本身不能作为训练特征否则就是数据泄露。我们只用傅里叶参数A21,φ21等来训练。这个小模型的OOB分数为0.76875效果一般说明单靠傅里叶参数区分它们有难度。引入周期约束与形态判据根据文献我们强制规定周期小于0.12天的目标只要不是EA一律归为DSCT。这是因为在ASAS-SN和ZTF等大型巡天数据中周期小于0.12天的变星绝大多数都是DSCT。这是一个强有力的先验知识。对于RRab和RRc我们利用其光变曲线形状的差异RRab不对称陡升缓降或缓升陡降RRc近似对称的正弦波。我们对分箱后的光变曲线进行单正弦拟合计算拟合优度。RRc的拟合优度通常更高。我们通过多项式拟合找到分布谷值设定阈值为0.868来区分两者。经过这一系列操作我们得到了980颗CEP、8347颗DSCT、457颗RRab和404颗RRc。4.5 自转变星ROT的净化ROT的挑战在于其光滑、近正弦的光变曲线很容易与非变星或某些脉动变星混淆。我们的策略不是进一步细分ROT而是尽可能地从初步分类为ROT的样本中剔除非变星。我们利用了ROT信号质量高的特点设定了严格的过滤条件power2 0.7且r² 0.7。同时应用上述的“周期小于0.12天归DSCT”规则。最终我们得到了12348颗高置信度的ROT。5. 结果验证与新发现目视检查的不可替代性机器学习给出了分类结果但天文发现必须严谨。我们坚持对所有机器分类出的变星进行了人工目视检查Visual Inspection, VI。这是最耗时但价值无可替代的一步。5.1 目视检查发现了什么分类错误纠正机器会把部分光变曲线形状相似的CEP误判为EW也会把一些噪声模式有规律的源误判为变星。人工检查修正了这些错误。周期准确性确认对于长周期变星尤其是EAGLS方法有时会找到谐波周期。我们在目录中专门增加了“周期是否正确”的标注列P_correct。特殊目标的识别多周期变星我们发现了一些目标在一个扇区内就表现出两个显著的周期例如既有食双星的掩食周又有脉动周期。对于这些我们目前只记录了最显著的周期但单独建表注明。类型变化的目标更有趣的是我们发现同一颗星如TIC 238853963在不同扇区表现出完全不同的变星类型如一个扇区是CEP另一个是ROT。这可能是恒星活动本身发生变化也可能是不同观测条件下主导的信号不同。这类目标也被单独列出。经过目视检查的最终“净化”我们获得的可靠变星数目为6046颗EA3859颗EW2058颗CEP8434颗DSCT482颗RRab416颗RRc以及9694颗ROT。5.2 交叉验证与14092颗新变星的诞生为了评估我们星表的可靠性并确认新发现我们做了两件事与权威星表交叉匹配将我们的目录与Gaia DR3、VSX、TESS EB Catalog、ZTF、ASAS-SN等已知变星星表进行交叉匹配。共有16298颗目标被至少一个已有星表收录。这证明了我们方法发现已知变星的能力。新发现统计用我们星表的总数减去上述匹配数剩下的便是新发现的变星。最终的数字是14092颗。这是一个非常可观的增量它们大多可能是亮度较低、信号较弱、或位于之前巡天未深入覆盖天区的目标。我们将最终星表绘制在赫罗图H-R Diagram和颜色-周期图上其分布与各类变星的已知理论区域基本吻合从侧面验证了分类的物理合理性。6. 踩坑实录与经验总结回顾整个项目从数据下载到最终星表发布每一步都充满了挑战。这里分享几个让我印象最深的“坑”和对应的解决方案。6.1 特征工程中的“过犹不及”问题最初我们试图引入尽可能多的特征包括从Gaia匹配的颜色、视星等、自行等。我们天真地认为更多信息总会带来更好结果。教训这导致了严重的数据缺失问题。很多TESS目标尤其是暗星在Gaia中没有对应的天体测量或测光参数。如果强行使用这些特征要么损失大量训练样本要么需要用各种方法填充缺失值这又会引入难以评估的偏差。解决方案做减法。最终我们只使用了完全从TESS光变曲线自身推导出的特征周期、傅里叶参数、统计矩、GLS功率。这些特征对于所有目标都是完备的。模型的表现证明这些“内生”特征已经包含了足够强的分类信息。6.2 样本不平衡与模型欺骗问题在尝试直接用7个子类训练一个“全能”随机森林时模型准确率看起来很高90%但仔细看混淆矩阵发现模型几乎把所有样本都预测成了数量最多的类别如ROT。教训在不平衡数据集上准确率是一个具有欺骗性的指标。模型学会了“偷懒”。解决方案采用分层分类和人工平衡训练集。先分四大类再在子类内部处理。在每一层都确保训练集中各个类别的样本数量大致相等通过下采样实现。这迫使模型去学习区分各类别的真实边界。6.3 周期“幽灵信号”与TESS系统噪声问题在周期搜索中我们发现不少目标在7天或13.7天TESS轨道周期的一半附近出现功率峰值。这显然是卫星轨道运动引起的系统噪声而非天体物理信号。排查我们随机抽取了一批在这些周期上有峰值的目标仔细检查其光变曲线和周期图。发现其光变曲线并无明显的7天周期性变化功率峰值往往很宽或伴随许多其他峰值。解决方案我们没有简单地在代码中硬性屏蔽这些周期。而是通过设置GLS功率power2和拟合优度r²的高阈值如ROT筛选中的0.7来过滤。真正的天体物理信号通常在这两个指标上都非常高而系统噪声产生的假信号很难同时满足高功率和高拟合优度。6.4 目视检查的规模化管理问题要对数万颗变星进行人工检查如何保证效率、一致性和可追溯性解决方案我们开发了一个简单的网页交互工具。前端以统一格式展示目标的TIC编号、光变曲线、相位折叠曲线、GLS周期图、机器分类结果和关键特征值。检查者只需点击下拉菜单选择确认的分类或修正的分类并勾选周期是否正确。所有结果自动记录到后台数据库。这比在终端看图、手动改Excel表格快了不止一个数量级也避免了版本混乱。这个项目让我深刻体会到在大数据天文时代机器学习不是要取代天文学家而是要成为天文学家手中更强大的望远镜和过滤器。它将我们从重复性的初筛劳动中解放出来让我们能更专注于那些真正有趣、复杂的案例和背后的物理。这套基于随机森林的TESS变星分类流程已经形成了一个相对稳定的 pipeline。如果你手头也有类似的时序测光数据不妨以此为基础框架进行尝试和调整。最重要的永远是理解你的数据设计贴合物理意义的特征并且永远不要完全相信机器的判断最后那一道人工审视的目光依然是科学发现最可靠的守门人。