Задание: по набору астрометрических измерений звезды определить параллакс и собственное движение, оценить погрешности и интерпретировать результат.
Данные и предварительные установки
- Эпоха опорного момента: $t_0=2016.0$.
- Единицы: угловые величины в миллисекундах дуги (mas), времена в годах.
- Модель измерений (смещения относительно положения в $t_0$):
$$\begin{array}{rl}\Delta \alpha \cos \delta (t_i)& =\Delta {\alpha }_0+{\mu }_{\alpha }(t_i-t_0)+\pi P_{\alpha }(t_i)+{ϵ}_{\alpha ,i},\\ \Delta \delta (t_i)& =\Delta {\delta }_0+{\mu }_{\delta }(t_i-t_0)+\pi P_{\delta }(t_i)+{ϵ}_{\delta ,i},\end{array}$$ где $\Delta {\alpha }_0,\Delta {\delta }_0$ — смещения в $t_0$, ${\mu }_{\alpha },{\mu }_{\delta }$ — собственные движения (mas/yr), $\pi$ — параллакс (mas), $P_{\alpha },P_{\delta }$ — параллаксные факторы в соответствующих координатах, $ϵ$ — шум (среднее $0$).
- Параллаксные факторы можно вычислить из барицентрических координат Земли $(X(t),Y(t),Z(t))$ в а.е.:
$$\begin{array}{rl}{P}_{\alpha }(t)& =\frac{-X(t)\sin \alpha +Y(t)\cos \alpha }{\mathrm{AU}},\\ P_{\delta }(t)& =\frac{-X(t)\cos \alpha \sin \delta -Y(t)\sin \alpha \sin \delta +Z(t)\cos \delta }{\mathrm{AU}}.\end{array}$$ (Если даётся уже готовая таблица $P_{\alpha ,i},P_{\delta ,i}$, можно её использовать напрямую.)
Пример входных данных (взяты искусственно для задания; единицы mas и годы):
- Таблица измерений (каждая строка: $t_i,\Delta \alpha \cos {\delta }_i,{\sigma }_{\alpha ,i},\Delta {\delta }_i,{\sigma }_{\delta ,i},P_{\alpha ,i},P_{\delta ,i}$):
$$\begin{array}{ccccccc}t& \Delta \alpha \cos \delta & {\sigma }_{\alpha }& \Delta \delta & {\sigma }_{\delta }& P_{\alpha }& P_{\delta }\\ 2015.0& -27.0& 2.0& 33.0& 2.0& -0.8& 0.5\\ 2015.5& -13.5& 2.0& 25.0& 2.0& -0.4& 0.9\\ 2016.0& 0.3& 2.0& 9.3& 2.0& 0.0& 1.0\\ 2016.5& 14.0& 2.0& -8.6& 2.0& 0.5& 0.6\\ 2017.0& 31.0& 2.0& -32.0& 2.0& 0.9& -0.1\\ 2017.5& 34.5& 2.0& -51.2& 2.0& 0.6& -0.7\\ 2018.0& 41.2& 2.0& -69.5& 2.0& 0.1& -0.9\\ 2018.5& 44.4& 2.0& -80.0& 2.0& -0.5& -0.6\end{array}$$
Метод решения (вкратце)
1. Сформируйте вектор наблюдений $y$ размера $2N$ (сначала все строки для RA, затем для Dec) и матрицу дизайна $A$ размера $2N\times 5$ для параметров
x=(Δα0μαπΔδ0μδ). x = \begin{pmatrix}\Delta\alpha_0\\ \mu_\alpha\\ \pi\\ \Delta\delta_0\\ \mu_\delta\end{pmatrix}.
x= Δα0 μα πΔδ0 μδ . Для каждой эпохи $i$ строки матрицы:
$${\text{RA-row}}_i=[1,(t_i-t_0),P_{\alpha ,i},0,0],$$ $${\text{Dec-row}}_i=[0,0,P_{\delta ,i},1,(t_i-t_0)].$$
2. Взвешенное МНК: взвешивание матрицы $W=\mathrm{diag}(1/{\sigma }_{\alpha ,1}^2,\dots ,1/{\sigma }_{\alpha ,N}^2,1/{\sigma }_{\delta ,1}^2,\dots )$. Решение:
$$\hat{x}=(A^{\mathsf{T}}WA{)}^{-1}{A}^{\mathsf{T}}Wy.$$
3. Ковариационная матрица оценок:
$$C=(A^{\mathsf{T}}WA{)}^{-1},$$ оценки погрешностей параметров: ${\sigma }_{x_j}=\sqrt{C_{jj}}$.
4. Оценка качества приближения:
$${\chi }^2=\sum _i{\left(\frac{r_i}{{\sigma }_i}\right)}^2,r=y-A\hat{x},$$ $${\chi }_{\nu }^2=\frac{{\chi }^2}{\nu },\nu =2N-5.$$ Если ${\chi }_{\nu }^2\approx 1$ — модель согласуется с данными в пределах ошибок; если ${\chi }_{\nu }^2\gg 1$ — недооценены ошибки или модель неверна.
Дополнительные проверки и оценки погрешностей
- Бутстрап или Монте-Карло: сгенерировать $M$ реализаций с шумом и повторно оценить распределения параметров.
- Проверить корреляции между параметрами по матрице $C$.
- Для низкого отношения сигнал/шум по параллаксу (${\mathrm{SNR}}_{\pi }=\pi /{\sigma }_{\pi }\lesssim 3$) применять байесовский подход к оценке расстояния (см. Bailer-Jones), т.к. преобразование $d=1/\pi$ будет сильно несимметричным.
Интерпретация результатов
- Параллакс $\pi$ (mas) в расстояние:
$$d[\mathrm{pc}]=\frac{1000}{\pi [\mathrm{mas}]},$$ приближённое погрешность:
$${\sigma }_d=\frac{1000{\sigma }_{\pi }}{{\pi }^2}.$$ Для больших относительных погрешностей используйте байесовскую оценку расстояния.
- Полная величина собственного движения:
$$\mu =\sqrt{{\mu }_{\alpha }^2+{\mu }_{\delta }^2}(mas/yr).$$ Тангенциальная скорость:
$$v_t[km/s]=4.74\frac{\mu [mas/yr]}{\pi [\mathrm{mas}]}.$$
Требуемые результаты (что сдать)
- Оценки параметров и их $1\sigma$ ошибки: $\Delta {\alpha }_0,{\mu }_{\alpha },\pi ,\Delta {\delta }_0,{\mu }_{\delta }$.
- Матрица ковариаций $C$ и значения корреляций.
- Значения ${\chi }^2$ и ${\chi }_{\nu }^2$.
- Графики: наблюдения и модель по времени (RA и Dec), параллаксная эллипса на небе (показать вклад параллакса отдельно), остатки по времени и их гистограмма.
- Интерпретация: расстояние $d$ с погрешностью (и обсуждение корректности преобразования), тангенциальная скорость $v_t$, значимость параллакса (${\mathrm{SNR}}_{\pi }$), возможные систематические ошибки и рекомендации.
Пояснения по часто возникающим проблемам
- Отрицательный оценённый параллакс возможен при шумных данных — интерпретировать через байесовский подход, а не просто отвергать.
- Сильные корреляции между $\pi$ и $\mu$ возникают при плохом покрытии годовой фазы; желательно иметь измерения в разные сезоны.
- Если ошибки точек негауссовы или есть выбросы — применять робастные методы (например, итеративное отсечение по остаткам).
Короткая шпаргалка команд (реализация в Python / numpy)
- Сформировать $A,y,W$, затем:
$$\hat{x}=(A^{\mathsf{T}}WA{)}^{-1}{A}^{\mathsf{T}}Wy,C=(A^{\mathsf{T}}WA{)}^{-1}.$$ - Посчитать ${\chi }^2$ и ${\chi }_{\nu }^2$, построить графики и оценить $v_t$ и $d$ по формулам выше.
Это задание можно выполнять как с приведённым примером, так и с реальными файлами эпhemeris (JPL Horizons) для вычисления $P_{\alpha },P_{\delta }$.