我希望写一篇小短文，来介绍如何回溯F51的历史数据。事实上，这个方法可以适用于COIAS、IASC、TDMP等在内的小行星搜寻项目。

我假设你已经了解了小行星的临时编号信息、MPC是什么这类问题。对于其它数据回溯的方法，比如ADAM，我强烈建议你查看日本同好的相关教程，这些教程上记载的非常详细。

我们将会使用一些必要的工具，比如你需要有一台可以联网的电脑等。在中国的一些区域，你可能需要准备VPN，否则你没有办法访问如JPL之类的工具。

我们需要事先确认一点。与下载FITS等数据自行测量不同，我们并不直接处理原始的FITS图像。我将简单介绍我们的步骤：

首先，我们要输入MPC小行星编号（即在 source = 'XXXX' 这一行填入目标编号）。该脚本会调用CADC（加拿大天文数据中心）的太阳系天体搜索API，以获取在巡天期间可能拍摄到该目标的PS1观测列表。这一步会给出每次曝光的观测时间及相关元数据。

接着，脚本使用astroquery.jpl horizons访问JPL Horizons，并请求目标在每个观测时刻、从PS1台址观测的精确位置，从而生成 (t, t1, t2, …, RA, Dec) 这样的数据。

在获得每个观测时刻的预测坐标后，脚本会将这些位置发送到PS1 CasJobs数据库，并与PS1探测目录进行交叉匹配。随后，再通过PS1图像cutout服务提取汇总结果。

由于望远镜数据已经经过PS1管道处理，因此仪器定标、天体测量解算以及源探测等步骤都由PS1管道完成。该管道还会对每个探测到的源执行PSF光度测量和孔径光度测量，提供诸如psfMag、apMag和magErr等参数。借助JPL Horizons给出的预测坐标和时间，便可以直接从PS1星表中提取相应条目。

由于恒星在短时间尺度上基本不动，而小行星会移动，因此记录到的小行星位置必然会随时间变化。按照 MPC每次提交至少需要2–3次观测的要求，这使我能够区分真实的移动天体与由背景恒星造成的误匹配。换句话说，天体测量归算、位置测量和光度测量都由 PS1管道完成并存储在 PS1 星表中，而我们的工作是从数据库中查询并筛选出相关探测结果。由于这种方法依赖JPL Horizons的高精度反向预报，我只尝试测量那些 U ≤ 2 的目标，因为这类天体的轨道不确定性已经较小，从而进一步降低了误匹配的可能性。

现在让我们正式开始！

1. 注册CasJobs账户

你需要注册一个CasJobs账号，这一步可以在MAST CasJobs网站完成。
打开后你应该可以看到这么一个页面（如图1）

图1

请点击上方的“Create Account”。这样你就会被导向一个新页面。请填写必要的信息并注册，你需要记住你注册时候使用的Userid和Password，这非常重要。

在注册成功后，你可以跳转到这个网页。我们将使用PS1 Public Archive网站上提供的一段脚本（由Rick White编写），紧接着你需要将网页中的代码下载下来整合为一个Python脚本。这一步你可以让AI帮忙，比如ChatGPT、Gemini、Grok......
实际上，脚本在输入CASJOBS_USERID与CASJOBS_PW时候会使用环境变量等的方式。我建议你使用弹窗输入的方式来进行。
比如：
1.1 先检查环境变量是否已经存在.如果存在，就不会进行弹窗输入。

if not os.environ.get('CASJOBS_USERID') or not os.environ.get('CASJOBS_PW'):

1.2 创建 tkinter 窗口，但主窗口隐藏。
root = tk.Tk()
root.withdraw()

1.3 输入用户名

userid = simpledialog.askstring("CasJobs Login", "Enter Casjobs UserID:")
if userid is None or userid.strip() == "":
messagebox.showerror("Login Failed", "Username is required.")
raise SystemExit("User cancelled or did not provide username.")

1.4 输入密码

password = simpledialog.askstring("CasJobs Login", "Enter Casjobs password:", show='*')
if password is None or password.strip() == "":
messagebox.showerror("Login Failed", "Password is required.")
raise SystemExit("User cancelled or did not provide password.")

1.5 把输入结果写回环境变量

os.environ['CASJOBS_USERID'] = userid.strip()
os.environ['CASJOBS_PW'] = password.strip()

1.6 后面拿这些信息去连接CasJobs

jobs = mastcasjobs.MastCasJobs(context="MyDB")

事实上，你可以把这个任务交给AI来做（毕竟AI可以尽可能的帮助你，而且一般不会出错）。

2. 输入小行星的编号信息
接下来，请在代码中找到
source = "rickwhite"

cadc_tab = cadc_ssos_query(source)
print(f"Returned table has {len(cadc_tab)} rows")
cadc_tab[:8]
这一行的内容。以该网页的代码为例。"rickwhite"是要查询的小行星的名称，你需要把该名称替换为你准备查询的小行星的名称。你可以在引号中输入多种形式的信息，如网页所介绍的:
# some other sample objects
# source = "33798" # fast-moving object, 10 arcsec/hour
# source = "134340" # Pluto
# source = "2003 YT1" # asteroid that passes close to North Celestial Pole (only 2 epochs in PS1)
# source = "Cruithne" # fast-moving, unusual Earth-associated orbit
# source = "5261" # Mars-trojan asteroid
# source = "243" # 243 Ida, first asteroid known to have a moon
# source = "2024 YR4"

3. 运行代码

接下来请找到这一行代码
saveplots = False
这是一个开关控制语句，请将False改为True。

在确保万无一失后，请运行这个脚本。脚本会自动完成大部分的工作。但实际上，我建议你使用我附上的这个代码，理论上它应该可以成功运行。该脚本中，默认了一些选项。
saveplots = True # 是否保存生成的图表
save_data = True # 是否将观测表和图像保存到本地
use_jpl = True # 是否尝试从 JPL Horizons 获取高精度星历（设为False跳过）
但实际上，这个脚本并不特别好，因为我的代码水平不高。但是你可以跑通它。我尝试编写了有关下载FITS数据的部分，但我明显意识到这一部分是多余的，甚至是错误的。我会在后续讲解相关的信息。

运行我所提供的脚本

如果你想尽快获得结果。请按照下面的步骤。
1. 配置脚本需要的库

这个脚本导入了astropy、astroquery、mastcasjobs、Pillow、numpy、matplotlib、requests、tkinter等库。你可以尝试
pip install astropy astroquery mastcasjobs pillow numpy matplotlib requests

确认无误后，请修改要查询的小行星信息。示例代码中为：
source = "2016 OR17"
请修改引号中的内容。运行脚本后，你可以看到一个弹窗。请输入CasJobs UserID和CasJobs password。随后点击“OK”。
一段时间后，脚本会创建一个输出目录
output_dir_base = f"PS1_data_{source}"
os.makedirs(output_dir_base, exist_ok=True)

也就是会生成类似：

PS1_data_2016 OR17/

如果 source = "Rickwhite"，则会生成：PS1_data_Rickwhite/
这个目录会出现在运行脚本时的工作目录下。

输出目录会包括以下内容。
1.1 观测表

• observations_目标名.csv
• observations_目标名.fits

1.2 数据文件

• detections_目标名.csv
• detections_目标名.fits

1.3 图像

• images.png
• 可能还有位置偏差图（若前面生成并启用保存）

1.4 缩略图目录

• thumbnails/
  里面是每张 PNG 小图。

1.5 原始 FITS 目录

里面是下载的原始 warp FITS 文件。

2. 信息补充

正如我上面所说的，该脚本在下载FITS数据方面存在很大问题。但是你可以选择无视反馈的下载失败，也可以在脚本反馈进行数据下载的时候结束脚本运行。值得一提的是，该脚本的抗干扰能力不佳，以至于可能出现由于网络波段等导致的代码运行失败。或许我们可以让AI来帮忙修改一下代码的这些问题（划去）。

3. 数据检查与处理

如果你成功运行了代码，你应该可以看到类似于图2中的内容。

图2

你应该首先查看images.png。如图3所示，图中有一些必要的信息。这些是65x65像素的切片（约16x16角秒），每张图的标题都标明了观测所使用的滤镜和观测日期。将物体位置列表与PS1目录匹配，如果该数据较好，它会使用蓝色的圆圈标记出来。如果图像质量较差或者存在问题，它会使用红色圆圈标注。（网页版中使用橙色圆圈标注。但是我觉得红色更加明显。事实上，颜色并不代表该测量结果是否可以使用。有些红色圆圈的结果一样可以用于补充数据。但有些蓝色的结果也可能存在问题）

图3

紧接着，请打开detections_XXXX.csv这个文件夹（在我的示例中是detections_2016 OR17.csv，以下简称csv文件）。你应该会看到一些数据。就像我们之前所说的，我们单个晚上需要最低2到3次的观测。所以你就会发现，上图中只有倒数第二列的两次数据可以使用。也就是hmjd这一列的56283.3798436167和56283.3932838167。

4. 时间转换

我们上文中得到的时间56283.3798436167和56283.3932838167准确来说是TAI时间。也就是TAI 时间 = 2012-12-22 09:06:58.49
随后我们把该时间转为UTC时间：

2012-12-22时，TAI − UTC = 35 s

所以09:06:58.49 − 35 s = 09:06:23.49 UTC

紧接着请查看csv文件的expTime这一列。该测量结果是45，我们将这个时间除以二得到22.5。所以我们需要再一次减去22.5秒
09:06:23.49 − 22.5 s = 09:06:00.99
所以最终结果为2012-12-22 09:06:00.99 UTC
接下来请再一次把这个时间加上22.5秒，紧接着转为MPC的格式。
09:06:00.99 + 00:00:22.5 = 09:06:23.49
所以现在时间是2012-12-22 09:06:23.49 UTC
转换为MPC时间格式，精确到5位小数为2012 12 22.37944

我们再一次尝试另外一个i波段的时间，56283.3932838167。已知MJD 56283 = 2012-12-22

计算小数部分0.3932838167×86400 = 33979.7206 s。换算后也就是TAI = 2012-12-22 09:26:19.7206

紧接着减去 22.5 秒（仍然是TAI）

09:26:19.7206 −22.5s = 09:25:57.2206 (TAI)

转换为UTC。2012-12-22时TAI − UTC = 35 s

所以：09:25:57.2206 − 35s = 09:25:22.2206 (UTC)

再加回 22.5 秒（现在是UTC）

09:25:22.2206 +22.5s = 09:25:44.7206 (UTC)

最后我们把这个时间转为MPC格式，保留5位小数，也就是2012 12 22.39288

实际上在上述流程中你会发现，减去22.5s与加上22.5s相互抵消了。但考虑到为了保证万无一失，我们尽可能的多做一步。

5. 坐标转换

处理完时间后，我们进行坐标的转换。请查看csv文件，找到hra与hdec这一列。56283.3798436167这个时间对应的坐标为86.398836 6.380934。我们需要把这个坐标转为时分秒的格式，也就是：
05h 45m 35.72s +06° 22′ 51.36″
让我们把赤经（RA）精确到0.01，赤纬（Dec）精确到0.1。所以该坐标为

05 45 35.70 +06 22 51.4

按照上面的方法处理hmjd为56283.3932838167的数据。我们可以得到坐标为

05 45 35.02 +06 22 51.9

不管是时间还是坐标，考虑到四舍五入等的情况，理论上可能存在一定的偏差，所以最后一位数字可能不同，但这一般情况下都是可以接受的。

6. 确认视星等

最后，请查看images.png，图片中的圆圈下给出了对应的视星等。56283.3798436167和56283.3932838167这两个时间分别对应的是21.72与21.58。我们需要将视星等精确到0.1，也就是21.7和21.6等。
所以理想情况下，该观测数据应该为：

K16O17R ZC2012 12 22.37944 05 45 35.70 +06 22 51.4 21.7 i F51
K16O17R ZC2012 12 22.39288 05 45 35.02 +06 22 51.8 21.6 i F51

请注意，上面的格式适用于Z代码的情况下。但是目前为止，Z代码已经被MPC禁用，因此你需要进行一些修改。

7. 下载已有数据并使用Find_Orb检查

在处理完所有数据后，我们从MPC下载该小行星的相关数据。并将我们处理的数据添加进去，随后使用Find_Orb进行计算。若轨道拟合后成功有所收敛，代表我们的结果是可信的。随后我们可以将其提交至MPC。

致谢

本研究基于空间望远镜科学研究所（STScI）维护的 Pan-STARRS1（PS1）移动目标查询技术文档（https://outerspace.stsci.edu/spaces/PANSTARRS/pages/298812291/How+to+search+for+moving+targets+in+PS1+images+and+catalogs）及 PS1移动天体数据库代码框架开发完成，感谢 STScI 公开共享的技术资源与代码模板。本研究对原始代码进行了适应性修改，以适配研究需求。同时感谢 Pan-STARRS1 Science Consortium（PS1 科学联盟）提供的观测数据支持，以及加拿大天文数据中心（CADC）、米库尔斯基空间望远镜档案库（MAST）对 PS1 数据的归档与开放服务。此外，感谢 astropy、astroquery 等开源天文数据分析库的开发者团队，为本研究的代码实现提供了基础。