太阳系柯伊伯带天体轨道通常在海王星轨道以外，由于日心距离远，温度低，是行星吸积盘原始遗迹之一，因此携带了有关行星形成的物理和化学的宝贵信息。目前已知的柯伊伯带天体超过4 000颗(http://www.minorplanetcenter.net)。2004~2005年，Brown团队(https://archive.org/details/howikilledplutow00mike)和Ortiz团队(http://www.infoastro.com/200809/26ataecina-haumea.html)声称发现了矮行星Haumea。它被临时命名为2003 EL61，并于2006年9月7日由官方小行星收录，编号为136108。Haumea是最大的柯伊伯带天体之一，并且是继冥王星和Makemake(136472)后最亮的海外天体，其轨道周期为284地球年，近日点35 AU，远日点51 AU，绝对星等0.358 mag^[1]，半长轴43.3 AU，偏心率0.19，轨道倾角28°。

光变曲线是获取天体物理性质的主要来源，为小天体及行星起源和演化的研究提供了重要线索。例如，文[2]推导出Haumea自转光变曲线的周期为3.915 4 ± 0.000 2 h，峰峰值幅度为0.28 ± 0.02 mag，颜色B-V=0.97 ± 0.03 mag。文[3]使用夏威夷大学2.2 m望远镜观测，观测时使用Kron-Cousins测光系统的B，V，R和I滤光片，导出周期为3.915 5 ± 0.000 1 h的双峰自转光变曲线，其峰峰值幅度为0.29 ± 0.02 mag。因自转周期短，Haumea呈旋转椭球体(Jacobi ellipsoid)，假设其已经达到流体静力平衡，文[2-3]结合这两个特性推测其质量密度ρ~2 500 kg·m^-3。

2005年，文[4]发现了Haumea的两颗卫星Hi′iaka和Namaka，这两颗卫星的星等较Haumea分别暗2.98 ± 0.03 mag和4.6 mag。由于地面望远镜观测受到大气的影响，观测时通常不能分辨Haumea及其卫星，但是天体测量结果可以用于确定光心的轨道。根据文[5]，光心轨道是双星系统相对轨道的缩放版本，因此可以通过光心的运动轨道初步推断双星系统的相对轨道。于是，文[5]结合径向速度测量结果、历史测量结果以及其他所有可用的天体测量结果(包括相对位置测量值、Hipparcos中间数据以及追踪光心移动的测光观测结果)，对双星系统GJ67AB进行轨道和质量求解，但是质量测量的精度仍然主要受到相对位置测量的影响。然而，对于光变幅度较大的天体，位置测量受到光度变化的影响。参考文[5]的工作，在求解Haumea与卫星的光心轨道时也会受到相同的影响，我们需要先估计影响范围。

基于紫金山天文台姚安观测站80 cm望远镜观测的Haumea CCD图像资料，本文推导Haumea的光度变化，为后续改进相对位置测量奠定基础。

1 观测资料和数据处理

我们使用紫金山天文台姚安观测站80 cm望远镜(详细参数如表 1)对矮行星Haumea进行两轮测光观测，观测时使用Johnson-Cousins Ⅰ滤光片，具体的观测资料情况如表 2。

我们利用MaxIm DL(https://diffractionlimited.com/maxim-dl/)软件对获得的原始图像进行预处理，主要包括本底、暗场和平场校正^[6]。由于反射光在CCD硅片薄层间来回反射，形成干涉条纹，导致图像背景不均匀。尤其对于Johnson-Cousins Ⅰ滤光片，CCD硅片的吸收率随着波长增加迅速下降，有相当一部分低能光子没有被完全吸收之前，在CCD硅片的薄层来回反射好几次，产生干涉条纹，不利于精确的光度测量。因此需要对拍摄的图像进行处理，剔除干涉条纹，具体的步骤如下。首先，从每晚观测的CCD图像中选择1或2幅图像。接着，对筛选的CCD图像(剔除了3σ)使用“3-2”公式(3倍中值-2倍均值)进行天空背景估计，得到每幅图像的背景估计值。所有筛选的CCD图像减去各自的背景估计值后，将这些图像每一个对应的灰度值取中值运算，最后的中值结果作为干涉条纹模板对应像素的灰度值。对于构造的干涉条纹模板，需要乘上一个比例因子以获得归一化干涉条纹模板，该比例因子为360 s除以待去除干涉条纹图像的曝光时间。最后，从待去除干涉条纹的原始图像中减去归一化干涉条纹模板，得到剔除干涉条纹后的图像。典型的处理结果如图 1。

2 光度测量结果与分析 2.1 光度测量数据

本文使用相对测光技术获取目标的相对仪器星等。由于两轮观测时视场中的星像较少，而且在每轮观测期间，目标仅运动了一小段距离，因此，我们可以使用相同的参考星进行测光^[7]。具体地，在第1轮观测时，选择3颗参考星，在第2轮观测时，选择5颗参考星。所有的参考星都是恒星且光度测量稳定。两轮观测中典型的目标和参考星的图像如图 2，其中，目标Haumea使用绿色圆圈标记，参考星使用红色圆圈标记。

我们使用二维高斯定心算法获取目标星和参考星的像素坐标，采用“3-2”公式对目标星和参考星周围环形区域天空背景进行估计，使用photutils(https://pypi.org/project/photutils/)测光库提供的孔径测光程序对目标星和参考星进行光度测量。photutils孔径测光程序允许我们使用一系列不同的孔径，选择光度测量精度最高的孔径作为最优孔径，得到每个夜晚目标星相对仪器星等。误差主要来源于孔径的目标信噪比，同时我们也考虑了目标与参考星之间因颜色不同带来的误差。

由于没有进行标准测光，为了将两轮观测的测光结果联合起来进行周期分析，根据文[8]，我们分别拟合两轮观测求出的目标相对仪器星等，拟合结果之间存在一个常数偏差，最后使用这个常数改正求出的目标相对仪器星等。每次拟合采用一轮多晚观测的数据，而不是每个夜晚的观测数据。

2.2 自转周期

为了测定准确的Haumea自转周期，本文使用两种周期求解方法，分别是相位弥散最小化方法以及Lomb-Scargle周期图法。两种方法求出目标Haumea自转周期结果和讨论分析如下。

相位弥散最小化方法由文[9]在1978年提出。相位弥散最小化方法通过求统计量θ的最小值获取周期，θ定义为

$ \theta=\frac{s^2}{\sigma^2}, $

(1)

$ \sigma^2=\frac{\sum\limits_{i=1}^N\left(x_i-\bar{x}\right)^2}{N-1}, s^2=\frac{\sum\limits_{j=1}^M\left(n_j-1\right) s_j^2}{\sum\limits_{j=1}^M\left(n_j-M\right)}, $

(2)

(1) 式中σ²为观测的N个离散数据点的方差。(2)式中x_i是指在t_i时刻观测的目标星等，x为N个数据点的均值；从N个数据点中选择M个有着相同相位Φ_i的样本，每个样本中包含数据点的个数为n_j；s_j²为第j个样本的方差。通常，θ值在0到1之间。如果求出的是正确的周期，那么θ会达到一个局部最小值，接近0。使用相位弥散最小化方法计算得到的周期如图 3。在频率为12.259 39(cycles/day)处达到局部最小值，最小值约为0.17。根据自转频率Ω与自转周期Ρ之间的关系Ω=24/Ρ，可以求出自转周期Ρ=1.957 68 h。

Lomb-Scargle周期图法^[10-11]适用于测定非均匀采样数据的周期。该方法基于傅里叶谱分析，考虑非均匀采样数据，并赋予每个数据点权重而不考虑时间间隔。假设有N个数据点h_j=h(t_j)(j=0, …, N-1)，t_j为实际测量时刻，按照

$ P_N(\omega)=\frac{1}{2 \sigma^2}\left\{\frac{\left[\sum\nolimits_j\left(h_j-\bar{h}\right) \cos \omega\left(t_j-\tau\right)\right]^2}{\sum\nolimits_j \cos ^2 \omega\left(t_j-\tau\right)}+\frac{\left[\sum\nolimits_j\left(h_j-\bar{h}\right) \sin \omega\left(t_j-\tau\right)\right]^2}{\sum\nolimits_j \sin ^2 \omega\left(t_j-\tau\right)}\right\} $

(3)

定义Lomb归一化周期图，其中，ω为角频率，σ²为数据的方差，τ定义为

$ \tan (2 \omega \tau)=\frac{\sum\nolimits_j \sin 2 \omega t_j}{\sum\nolimits_j \cos 2 \omega t_j} . $

(4)

Lomb归一化功率谱最大值对应的值为求解的最佳周期。使用Lomb-Scargle方法计算得到的周期如图 4。图 4展示了一个明显的最大峰值在自转频率为12.259 40 ± 0.000 04(cycles/day)，可以求出自转周期Ρ=1.957 7 ± 0.000 1 h。与相位弥散最小化方法求出的结果仅相差0.072 s，表明两种方法求出的结果具有良好的一致性。

2.3 自转光变曲线分析

根据文[12]，对于像Haumea这种因快速旋转而变形拉长成椭球体的TNO(Trans-Neptunian Object)，其完整的光变曲线应包含两个极大值和两个极小值，极大值代表天体的前后两个侧视面，而极小值是两端的视图。因此，将2.2节求解的Haumea自转周期拉长至两倍，求得Haumea的自转周期T=3.915 4 ± 0.000 2 h，该结果与文[2-3]的结果吻合。

我们使用一般的线性最小二乘法拟合测光数据(剔除了3σ)，拟合模型为

$ m_i=c+\sum\limits_{j=1}^3\left[a_j \sin \left(j \frac{2 \pi}{T} t_i\right)+b_j \cos \left(j \frac{2 \pi}{T} t_i\right)\right], $

(5)

其中，t_i为观测时刻的JD时间(校正了光行时)。各系数拟合结果和误差如表 3。观测数据点和拟合模型光变曲线如图 5，图中显示两个极大值分别为0.12 mag和0.10 mag，两个极小值分别为0.14 mag和0.09 mag，对应的χ²=252，拟合结果服从自由度为205的χ²分布。χ²定义为

$ \chi^2=\sum\limits_{i=1}^N\left[\frac{y_i-m_i}{\sigma_i}\right]^2, $

(6)

其中，N为拟合数据点总个数；y_i为第i个数据点对应的相对仪器星等；m_i为(5)式计算的结果；σ_i为第i个数据点的测量误差。

通过图 5可以看出，模型光变曲线拟合了大部分的观测数据点。少部分数据点残差较大，可能是因为时间间隔内采样稀疏或者数据点本身误差较大，也可能由于受到Haumea旁边两颗卫星Hi′iaka和Nāmaka光变的影响，该数据点弥散较大。

3 光变对位置测量的影响

对于像Haumea这种光度变化幅度较大的天体，光度会对位置测量造成一定的影响。由于Nāmaka较暗，参考文[1]的工作，只考虑Haumea与卫星Hi′iaka构成的双星系统(如图 6)。在天空平面上，设主星与卫星之间的角距离为d；质心与光心之间的角距离为Δ；P和S分别代表主星和卫星；C_M和C_L分别是质心和光心。根据杠杆原理，

$ \varDelta=\left[\frac{1}{1+r_{\mathrm{L}}}-\frac{1}{1+r_{\mathrm{M}}}\right] d, $

(7)

其中，r_L=L_p/L_s；r_M=m_p/m_s；L_p和L_s分别是主星和卫星的亮度；m_p和m_s分别是主星和卫星的质量。假设卫星Hi′iaka的光变幅度不变，Δ的变化主要由Haumea光变引起，可以根据误差估计公式估计Δ受光变影响的变化范围。令

$ k=\frac{1}{1+r_{\mathrm{L}}}-\frac{1}{1+r_{\mathrm{M}}}, $

(8)

则有Δ=kd，那么Δ受光度的变化影响为

$ \varDelta^{\prime}=(k+\delta k) d, $

(9)

k的变化只与主星和卫星的亮度变化有关，因此有

$ \delta k \leqslant-\frac{L_{\mathrm{s}}}{L_{\mathrm{\Sigma}}} \frac{\delta L_{\mathrm{\Sigma}}}{L_{\mathrm{\Sigma}}}. $

(10)

其中，L_Σ是主星与卫星的亮度总和；根据文[4]，卫星Hi′iaka的亮度是主星Hamuea的5.9%，因此r_L=16.9，L_s/L_Σ=1/17.9。根据文[14]，r_M=223.80。根据星等转换关系式有δL_Σ/L_Σ=10^-0.4×Δm-1，其中，Δm为Haumea光变曲线上对应的星等值。

从IMCCE历表(https://www.imcce.fr/)可以获得任意时刻主星与卫星之间的距离d。在整个观测期间，选取Haumea与卫星Hi′iaka的距离最远(d≈1.446″)的时间节点，每隔1 min取数据点，一共300个数据点，覆盖了Haumea一个自转周期(T=3.915 4 h)。根据2.3节求得的Haumea光变曲线图(图 5)，我们可以获得Haumea光度的变化对位置测量的影响曲线图(如图 7)。其中，横坐标为Haumea的自转相位，与图 5的横坐标相对应，纵坐标为该时刻因光度变化，Haumea位置测量受到的影响。可以看出，Haumea光度变化对位置测量的影响最大可达到-9~9 mas。因此，利用Haumea的光变特征改正其位置测量是有必要的。

4 总结与展望

我们对紫金山天文台姚安观测站80 cm望远镜拍摄的矮行星Haumea CCD图像剔除了干涉条纹，使用相对测光技术获得Haumea的相对仪器星等，通过联合两轮观测资料探究了其光变情况。本文使用相位弥散最小化方法以及Lomb-Scargle周期图法计算Haumea的自转周期，结果显示一条明显的双峰自转光变曲线，周期为3.915 4 ± 0.000 2 h，峰峰值幅度为0.26 ± 0.01 mag；最终估计了Haumea由于自转产生的光度变化对位置测量的影响最大为-9~9 mas。在后续的工作中，我们将结合Haumea光度测量与位置测量，利用其光度变化特征改进位置测量精度。

致谢: 感谢紫金山天文台姚安观测站80 cm望远镜运行组全体成员。感谢郑中杰在本工作期间提供的帮助。