Книга 1. Звезды свидетельствуют — страница 6 из 20

Анализ систематических ошибок звездного каталога

0. Основная идея

0.1. Наглядная аналогия

Необходимость анализа ошибок звездного каталога нами разъяснена ранее. Мы, разумеется, прежде всего, имеем в виду Альмагест, но излагаемый метод будет применяться и к другим реальным и искусственно генерируемым каталогам. В этой главе будет показано, как обнаруживать и компенсировать систематическую ошибку. Идея метода проста и естественна. Более того, он в той или иной форме давно уже используется в математической статистике. Чтобы пояснить его идею, рассмотрим следующий пример. Пусть мы рассматриваем результаты стрельбы в тире, которые показаны на рисунке.

Точками изображены следы от пуль. Спрашивается, какова точность попадания? Ответ очевиден: плохая. Однако, мы видим, что так называемая кучность достаточно хорошая. Это заставляет предположить, что сам стрелок хороший, а то обстоятельство, что все пули попали в сторону от «яблочка» может быть объяснено, например, тем, что прицел установлен неправильно. Разумеется, выяснить, вследствие чего произошло такое смещение, не видя винтовки, нельзя, но определить величину смещения можно. Разумный способ сделать это — определить геометрический центр всех результатов и провести вектор из центра яблочка в найденный центр. На рис. 5.а — это вектор S.

Рис. 5.а. Мишень со следами попадания пуль.

Как формально получить вектор S? Очень просто. Нужно взять векторы x_i соответствующие i-му результату стрельбы, и усреднить их по общему количеству N выстрелов:

Отметим также, что вектор S альтернативно можно вычислить из задачи минимизации среднеквадратичного отклонения: найти вектор S, при котором достигается минимум функции

Здесь мы приняли, что (x_i — S)² = (x_i1 — S₁)² + (x_i2 — S₂)², где x_i1, x_i2 и (S₁, S₂) — координаты векторов x_i и S соответственно.

Тогда точность самогό стрелка можно охарактеризовать разбросом результатов вокруг найденного центра, что существенно выше точности попадания в яблочко. В нахождении вектора S и заключается в данном примере компенсация систематической ошибки. Собственно, S ее и представляет.

Формально, если мы изменим систему координат, сместив ее начало из центра яблочка на вектор S, то результаты стрельбы в новой системе координат будут содержать лишь случайные составляющие (вызванные дрожанием рук и т. п.) и не будут содержать регулярной составляющей.

Теперь вернемся к звездному каталогу. Мы хотим проверить, существует ли систематическая ошибка в какой-то его части и, если существует, определить ее. Пусть сначала не стоит проблема датировки, то есть мы наверняка знаем дату составления каталога t_А (здесь «А», конечно же, означает Альмагест, но все рассуждения верны и для любых других каталогов). Тогда нужно сравнить истинные координаты звезд на момент t_А (известные из современных точных каталогов) со значениями координат из изучаемого каталога, относящиеся к исследуемой его части. При этом сравнении, как и в примере с мишенью, нужно найти среднее отклонение сравниваемых координат. Пусть общее число звезд в избранной области равно N. Обозначим через l_i и L_i соответственно эклиптикальную долготу i-й звезды в изучаемом каталоге и точное значение ее долготы. Аналогично, обозначим b_i и B_i эклиптикальные широты i-й звезды (в каталоге и точное значение). Тогда средняя (систематическая) ошибка по долготе равняется,

а систематическая ошибка по широте

Эти ошибки, как уже сказано, могли произойти из-за неправильного определения плоскости эклиптики и ряда других, вообще говоря, неизвестных нам причин. Причины мы так и не выявим (впрочем, сделаем некоторые предположения на их счет), а вот ошибку, вызванную ими, скомпенсируем. Для этого нужно, попросту, изменить систему координат каталога (аналогично тому, как мы это сделали в примере с мишенью) так, чтобы результирующие средние долготные и широтные ошибки равнялись нулю.

0.2. Реализация метода

В данном пункте мы покажем, как конкретно реализуется общая идея, описанная выше.

Прежде всего, отметим, что мы будем компенсировать только широтную ошибку. Причины тому названы ранее и мы лишь упомянем здесь, что это позволяет уменьшить ошибку расчетов, что существенно ввиду низкой точности древних каталогов.

Итак, у нас есть каталог, из которого мы выбрали большую группу звезд в количестве N с координатами (l_i, b_i)^N_i=1. в результате отождествления мы знаем их «двойников» из современного каталога. Обозначим координаты этих двойников, рассчитанные на момент t, через (L_i(t), B_i(t))^N_i=1. Теперь предположим, что мы хотим проверить, какова была бы систематическая широтная ошибка в предположении, что дата составления каталога есть t_A. Обозначим L^А_i = L_i(t_A), B^A_i = B_i(t_A) и введем широтную невязку ΔB^A_i = B^A_i — b_i.

Наша цель — минимизировать величину

меняя систему координат, то есть, попросту, проводя новую координатную сетку, отличную от принятой в каталоге.

Изменение координатной сетки можно параметризовать двумя величинами, если мы рассматриваем задачу минимизации указанного выше выражения: γ и φ. Объясним их значение. Они показаны на приведенном ниже рис. 5.1. Здесь γ — угол между реальной эклиптикой и эклиптикой каталога, а φ — угол между прямой равноденствия и прямой пересечения реальной эклиптики с эклиптикой каталога.

Итак, минимизируя указанное выше выражение, можно найти значения γ_stat и φ_stat, параметризующие изменение системы координат и дающие исходный минимум. Их явный вид приведен ниже, в формулах (5.5.2) и (5.5.3).

Величина σ_min является остаточной (после компенсации систематической ошибки) среднеквадратичной широтной ошибкой. Явный вид формулы для остаточной дисперсии σ_min см. ниже после формулы (5.5.10). Он получается подстановкой γ_stat и φ_stat в качестве параметров в выражение для среднеквадратичного отклонения. Вывод этих формул приведен ниже.

Однако, мы не можем считать, что нашли систематическую ошибку (вернее, параметры γ_stat и φ_stat ее характеризующие) абсолютно точно. Дело в том, что индивидуальные ошибки измерения, носящие случайный характер, также влияют на значения γ_stat и φ_stat. Следовательно, мы можем лишь утверждать, что истинные значения систематической ошибки лежат где-то неподалеку от γ_stat и φ_stat.

Чтобы сделать утверждение более точным, введем понятие доверительного интервала. Зададимся некоторым уровнем доверия 1 — ε. Если, например, ε = 0,1, то уровень доверия равен 0,9. Уровень доверия представляет собой вероятность, с которой мы гарантируем точность своих результатов. Доверительный же интервал представляет собой отрезок, накрывающий неизвестное нам истинное значение параметра с вероятностью не меньшей 1 — ε. Обозначим

I_γ(ε) = [γ_stat — x_ε, γ_stat + x_ε]

— доверительный интервал для истинного значения параметра γ, и

I_φ(ε) = [φ_stat — y_ε, φ_stat + y_ε]

— доверительный интервал для истинного значения параметра φ. Можно показать (см. ниже), что значения x_ε и y_ε можно вычислить по формулам

x_ε = q_ε, y_ε = q_ε,

где q_ε представляет собой (1 — ε/2) — квантиль стандартного нормального распределения, который находится из таблиц.

Итак, если мы задаемся некоторым уровнем доверия 1 — ε, то с вероятностью не меньшей 1 — ε можно гарантировать принадлежность истинного значения γ интервалу I_γ(ε) и значения φ- интервалу I_φ(ε).

0.3. По значению систематической ошибки нельзя датировать каталог

Дадим теперь несколько иную интерпретацию найденных значений γ_stat и φ_stat. Пользуясь значениями координат звезд (на самом деле, достаточно рассмотреть лишь широты), легко определяются полюса эклиптики Р_А (для изучаемого каталога) и P(t) (для расчетного каталога на момент t), рис. 5.b.

Рис. 5.b. Два полюса — эклиптики и каталога.

Очевидно, что дуговое расстояние между Р_А и P(t) равно в точности γ_stat и компенсация систематической ошибки есть не что иное, как совмещение этих полюсов. Теперь посмотрим, как эта картина меняется со временем. Поскольку движение P(t) происходит в пределах одного градуса, можно воспользоваться плоской картинкой и предположить равномерность движения Р(t), см. рис. 5.b.

Скорость v этого равномерного движения нетрудно подсчитать, зная значения γ_stat в двух различных точках. Тогда легко найти момент t*, когда положение истинного полюса ближе всего отстоит от положения полюса каталога. На первый взгляд может показаться, что следует объявить этот момент искомой датировкой, вычисленной, кстати, обработкой значений координат большого количества звезд. Однако, как мы уже выяснили, логика эта неправильная и датировать каталог моментом t* нельзя. В самом деле, если возможная систематическая ошибка в определении эклиптики Птолемеем может достигать величины δ, то все моменты времени, отвечающие прохождению полюса P(t) через круг радиуса δ с центром в точке Р_А, должны рассматриваться как возможные кандидаты на время датировки. Но величину δ мы не знаем. Конечно, мы можем ее оценить, но лишь при условии, что нам известна датировка каталога. При другой предполагаемой датировке величина оценки будет другой. Таким образом, предполагаемая датировка уже заложена в оценке данной величины.

Поэтому, в зависимости от сделанной Птолемеем систематической ошибки, то есть ошибки в определении эклиптики, момент t* может оказаться либо более ранним, чем истинная дата составления каталога, либо более поздним. В первом случае каталог (вернее, та его часть, для которой ищется γ_stat) «удревняется» — он, попросту, становится похож на каталог, составленный в году t*. Во втором случае, если t* — момент более поздний, чем истинная дата составления, — каталог омолаживается. Мы увидим далее, что обе эти возможности реализованы в Альмагесте. Однако, слова «удревняется» и «омолаживается» относятся к каталогу, в котором систематические ошибки не скомпенсированы. После их компенсации остается «рафинированный каталог», содержащий лишь случайные ошибки, среднеквадратичная величина которых оценивается значением σ_min, но индивидуальные значения оценить невозможно.

Перейдем теперь к более подробной реализации изложенной выше общей идеи.

1. Основные обозначения

Начиная с этой главы, мы считаем, что имеем дело с каталогом, все звезды которого имеют единственное отождествление со звездами из современного каталога. В соответствии с этим, будем идентифицировать звезды индексом i и обозначать через l_i, b_i, соответственно эклиптикальные долготу и широту i-й звезды в Альмагесте. Через L_i(t), B_i(t) мы обозначим истинные долготу и широту i-й звезды в эпоху t. Напомним, что время t мы отсчитываем от 1900 года «назад» и измеряем в столетиях, то есть, например, t = 3,15 соответствует году 1900 — 3,15 × 100 = 1585 году н. э., а t = 22,0 отвечает году 1900 — 22 × 100 = 300 году до н. э.

Пусть t_А — неизвестное нам время составления каталога Альмагеста. Обозначим через L^A_i, B^A_i истинные долготу и широту i-й звезды в год составления каталога, то есть L^A_i = L_i(t_А), B^A_i = B_i(t_А)· Пусть ΔB_i(t) = B_i(t) — b_i — разность между истинной широтой i-й звезды в момент времени t и ее широтой в Альмагесте. Назовем величину ΔB_i(t) широтной невязкой на момент времени t. Эта величина имеет смысл погрешности в определении широты i-й звезды Альмагеста при условии, что он составлен в эпоху t. Естественно, что ΔB_i(t_A) = ΔB^A_i представляет собой истинную погрешность в определении широты.

Как уже отмечалось в главе 3, в случае с Альмагестом приходится анализировать лишь широтные ошибки. Причины этого подробно разъяснены выше.

2. Параметризация групповых и систематических ошибок

Рассмотрим некоторую совокупность звезд, например, созвездие или группу созвездий. Определим групповую ошибку в широтных координатах этих звезд как погрешность в определении широт звезд из этой совокупности, проистекающую из перемещения рассматриваемой звездной конфигурации как единого целого по сфере. Следовательно, — и это обстоятельство мы особо подчеркнем, так как оно существенно используется в дальнейшем, — любое подмножество этой конфигурации также перемешается по сфере как единое целое и на тот же угол, что и вся конфигурация. Такое перемещение имеет три степени свободы, то есть может быть описано тремя параметрами. Определим их.

На рис. 5.1 наглядно изображена соответствующая картина. На звездной сфере с центром в точке O нанесено положение реальной эклиптики на момент времени t_А. На эклиптике обозначены точки Q и R соответственно весеннего и осеннего равноденствий. Точка P отмечает северный полюс эклиптики. Точка E изображает положение некоторой звезды. Как уже говорилось, все групповые ошибки, для фиксированной группы звезд, в эклиптикальной широте, совершаемые составителем каталога, можно без ограничения общности считать следствием неправильного определения полюса эклиптики, то есть, результатом того, что вместо точки P на небесной сфере он принял в качестве полюса точку P_A.

Рис 5.1. Параметры, определяющие систематическую ошибку.

Этой точке соответствует возмущенная эклиптика, которая названа на рис. 5.1 эклиптикой каталога. Ее положение можно однозначно задать следующими двумя параметрами. Во-первых, угол γ между прямыми OP и ОР_А, или, что то же самое, плоский угол между плоскостями реальной эклиптики и эклиптики каталога. Во-вторых, угол φ между прямой равноденствий RQ и прямой CD, образованной пересечением плоскостей реальной эклиптики и эклиптики каталога. Такая параметризация удобна для аналитических целей. Однако далее мы будем иногда наряду с φ использовать величину β, которая интерпретируется следующим образом, рис. 5.1. Смещение эклиптики «разлагается» на два поворота — вокруг оси равноденствия RQ на угол γ и вокруг оси, лежащей также в плоскости эклиптики и перпендикулярной оси RQ, на угол β. Итак, β представляет собой величину дуги Q_AQ, являющейся частью большого круга, проходящего через полюс Р_А и точку Q. Астрономический смысл точки Q_A весьма прозрачен. Это точка весеннего равноденствия на эклиптике каталога. Ясно, что углы γ и φ однозначно определяют углы γ и β, и наоборот. Искомую связь мы находим из сферического прямоугольного треугольника CQ_AQ. Здесь угол при вершине Q_A — прямой, угол при вершине C равен γ, а длина дуги CQ равна β. В результате получаем:

sin β = sin γ × sin φ. (5.2.1)

Третья степень свободы заключается в повороте сферы вокруг оси P_AP'_A, рис. 5.1. Но такой поворот меняет лишь долготы звезд и не меняет их широты. Поэтому эту степень свободы мы не рассматриваем. Отметим, что вместо указанных параметров можно выбрать любые другие базисные параметры, задающие вращение сферы. Ясно, что это не влияет на дальнейшее построение.

Посмотрим теперь, как искажаются истинные координаты звезды i при наличии такой систематической ошибки. Истинные широта и долгота этой звезды равны соответственно длинам дуг ЕЕ' и QE', отсчитываемым по часовой стрелке, если смотреть с полюса Р. Искаженные широта и долгота b_i, l_i равны соответственно длинам дуг EE_A и Q_AE_A. Отметим, что широты звезд, у которых истинная долгота больше долготы точки D и меньше долготы точки C, уменьшаются, а остальные — увеличиваются, рис. 5.1. Этот вывод справедлив, строго говоря, не для всех звезд. Он неверен для звезд, расположенных вблизи полюсов Р и Р' и находящихся на угловом расстоянии γ (или менее) от них. Однако, поскольку γ мало, звезд в такой маленькой области очень немного. Среди звезд Альмагеста их практически нет. Как мы увидим, величина γ составляет около 20′.

Учитывая малость величины у, можно предложить следующую приближенную формулу для широтной невязки:

ΔB^A_i = γ × sin (L^A_i + φ). (5.2.2)

Иначе говоря, систематическую погрешность определения широт звезд можно изобразить синусоидой, показанной на рис. 5.2. Она очень похожа на кривую, обнаруженную ранее К. Петерсом и Е. Кнобелем [1339] при обработке каталога Альмагеста. Погрешность формулы (5.2.2) не превышает 1′ для звезд, у которых |b_A| ≤ 80°. Такая погрешность является для нас несущественной, и мы в дальнейшем говорить о ней не будем, считая формулу (5.2.1) абсолютно точной. Для корректности мы исключим из рассмотрения звезды, имеющие широты более 80 градусов по абсолютной величине. В дальнейшем в этой главе речь будет идти о систематической погрешности, поскольку излагаемые методы работают в предположении, что рассматривается большая совокупность звезд. Проверка того, что найденная погрешность совпадает (или не совпадает) с групповыми ошибками отдельных созвездий, представляет собой самостоятельную задачу. Применительно к Альмагесту она рассматривается ниже, в главе 6.

Рис. 5.2. Зависимость систематической широтной невязки от долготы.

В предположении, что известно время t_A составления каталога, можно определить параметры γ и φ, задающие систематическую ошибку, следующим образом.

1) Вычислим на момент времени t_А истинные ширóты и долгóты для всех звезд из рассматриваемой совокупности.

2) найдем параметры γ* и φ*, дающие решение задачи

σ²(γ*, φ*) → min, (5.2.3)

где

Если бы в каталоге не было других ошибок, кроме систематических, соотношение (5.2.3) превратилось бы в уравнение σ²(γ*, φ*) = 0. Однако наличие случайных ошибок в координатах звезд делает минимум в (5.2.3) отличным от нуля.

В нашем случае момент t_A составления каталога неизвестен. Поэтому мы вынуждены рассчитывать систематические ошибки для всех значений t из рассматриваемого интервала 0 ≤ t ≤ 25, а именно, для каждого t определяются положение реальной эклиптики и ось равноденствия. Затем, как и на рис. 5.1, вводятся параметры γ = γ(t), φ = φ(t) и β = β(t), задающие взаимное положение эклиптики эпохи t и эклиптики каталога. Значения величин γ(t) и φ(t) находятся как решение задачи

Опять-таки, если бы мы имели идеальный случай и каталог не содержал бы других ошибок, кроме систематических, то соотношение (5.2.4) можно было бы, пренебрегая исключительно слабыми эффектами от собственного движения звезд, записать как уравнение σ²(γ(t), φ(t), t) = 0.

По поводу эффектов от собственного движения напомним, что количество заметно движущихся звезд на небе очень мало по сравнению со всеми звездами Альмагеста. Решение последнего уравнения существовало бы при всех t, но дату t_A определить из такого уравнения невозможно. Тем более ее нельзя найти из соотношения (5.2.4), заменяющего указанное уравнение в реальном случае каталога со случайными ошибками. Можно вычислить лишь систематическую ошибку как функцию предполагаемой датировки t. Эта ошибка, естественно, зависит от предполагаемой датировки — благодаря колебанию эклиптики со временем. Именно поэтому мы здесь говорим не о датировке каталога, а о нахождении его систематической ошибки как функции предполагаемой датировки t.

В реальном каталоге, кроме указанной систематической ошибки, присутствуют и случайные ошибки. Поэтому отклонения B_i(t) — b_i являются случайными величинами, значения которых разбросаны вокруг синусоиды их среднего значения, изображенной на рис. 5.2. В предположении, что остальные погрешности каталога, кроме систематических, носят случайный характер, задача определения γ(t) и φ(t) является задачей определения параметров регрессии.

3. Определение параметров γ(t) и φ(t) методом наименьших квадратов

Найдем решение γ(t) и φ(t) задачи минимизации (5.2.4), (5.2.5). Ниже, в конкретных примерах, эта задача рассматривается для совокупностей, состоящих из различного количества звезд. Поэтому в расчетах мы будем использовать следующие нормированные величины, в которых N означает число звезд в изучаемой совокупности:

Отметим, что все эти величины могут быть вычислены для любого момента времени t, исходя из значения современных координат звезд и координат звезд в каталоге Альмагеста.

Очевидно, что задача минимизации (5.2.4) эквивалентна задаче минимизации

σ²₀(γ, φ, t) → min. (5.3.1)

в том смысле, что параметры γ(t) и φ(t), определяемые соотношением (5.3.1), совпадают с параметрами, определяемыми решением задачи (5.2.4).

Как отмечалось, решение задачи (5.3.1) имеет смысл лишь для больших совокупностей звезд, и поскольку ниже мы будем изучать статистические свойства этого решения, то здесь и далее через γ_stat(t) и φ_stat(t) обозначаются величины, удовлетворяющие соотношению (5.3.1). Значение

σ_min(t) = σ₀(γ_stat(t), φ_stat(t), t). (5.3.2)

имеет прозрачный физический смысл. Это среднеквадратичная широтная невязка по рассматриваемой совокупности звезд в момент времени t, получившаяся после компенсации найденной систематической ошибки γ_stat(t), φ_stat(t). Как мы увидим далее, величина σ_min(t) от времени практически не зависит ввиду крайне малой скорости собственного движения подавляющего большинства звезд. Поэтому мы будем использовать также обозначение σ_min. Заметим, что до компенсации этой ошибки среднеквадратичная широтная невязка в момент t равнялась величине

которая, вообще говоря, зависит от t. Таким образом, разность Δσ(t) = σ_init(t) — σ_min(t) оценивает эффект от компенсации систематической ошибки γ_stat(t), φ_stat(t).

Далее при определении величин γ_stat(t) и φ_stat(t) из соотношения (5.3.1) момент времени t будем предполагать фиксированным. Поэтому аргумент t в выкладках мы опускаем, то есть, будем писать L_i вместо L_i(t), и s_b вместо s_b(t) и т. д.

Для нахождения минимума в соотношении (5.3.1), возьмем частные производные функции σ²₀(γ, φ, t) по γ и φ и приравняем их нулю. С учетом формулы sin (L_i + φ) = sin L_i× cos φ + cos L_i × sin φ, получим уравнения

s_bcos φ + c_bsin φ = γ [s₂cos² φ + 2 d cos φ sin φ + c₂ sin² φ], (5.3.4)

— c_bcos φ + s_bsin φ = γ [-d cos² φ + (s₂- c₂) cos φ sin φ + d sin² φ]. (5.3.5)

Разделив уравнение (5.3.4) на (5.3.5), получим

Приведя обе части этого равенства к общему знаменателю, приходим к следующему уравнению относительно tan φ:

(1 + tan² φ)(c_bs₂ — s_bd) + (1 + tan² φ) tan φ (c_bd — s_bc₂) = 0.

Отсюда легко найти тангенс оптимального значения φ_stat:

Равенство (5.3.6) позволяет однозначно определить φ_stat, после чего оптимальная величина γ_stat может быть найдена, например, из (5.3.4):

Формулы (5.3.6) и (5.3.7) дают искомое решение задачи нахождения оценок γ_stat(t) и φ_stat(t) методом наименьших квадратов.

Полезно провести анализ чувствительности в этой задаче. Рассмотрим частные производные второго порядка функции σ²(γ, φ, t) по γ и φ:

Учитывая равенства (5.3.4)-(5.3.7), нетрудно получить для этих частных производных следующие выражения:

a₁₁= 2(s₂cos² φ_stat + 2d cos φ_statsin φ_stat + c₂sin² φ_stat) = (2/γ_stat)(s_bcos φ_stat + c_bsin φ_stat),

a₁₂= 2(c_bcos φ_stat — s_bsin φ_stat), (5.3.8)

a₂₂= 2γ²_stat(s₂sin² φ_stat — 2d sin φ_statcos φ_stat + c₂cos² φ_stat).

Для оценки погрешностей в определении величины среднеквадратичной ошибки σ(γ, φ, t) при отклонении значений γ и φ от найденных оптимальных величин γ_stat(t) и φ_stat(t) воспользуемся следующим разложением функции σ²(γ, φ, t) в окрестности точки (γ(t), φ(t)):

σ²(γ, φ, t) ≈ σ²_min+ a₁₁(t) (γ — γ_stat(t))²+ 2a₁₂(t)(γ — γ_stat(t))(φ — φ_stat(t)) + a₂₂(t)(φ — φ_stat(t))². (5.3.9)

В последней формуле мы пренебрегли членами третьего и более высоких порядков малости по отношению к разностям γ — γ_stat(t) и φ — φ_stat(t).

Формула (5.3.9) позволяет элементарными средствами оценить чувствительность среднеквадратичной ошибки σ(γ, φ, t) к вариациям параметров γ и φ. Для этого достаточно определить величины a₁₁, а₁₂ и а₂₂, входящие в правую часть (5.3.9). После вычисления оценок γ_stat(t) и φ_stat(t) их легко найти по формуле (5.3.8).

Формула (5.3.9) показывает, что «линии уровня» среднеквадратичных ошибок являются эллипсами на плоскости (γ, φ). См. рис. 5.3.

Рис. 5.3. Линии уровня среднеквадратичной ошибки σ(γ, φ, t) при фиксированном t.

Центром эллипса является точка (оценок γ_stat, φ_stat), в которой значение среднеквадратичной ошибки равно σ_min. Направления осей эллипсов и соотношение между ними определяются стандартными формулами аналитической геометрии через величины а₁₁, а₁₂, а₂₂, а именно, угол наклона α одной из осей эллипса определяется соотношением:

Вторая ось перпендикулярна ей. Длины осей относятся друг к другу как λ₁/λ₂, где λ₁ и λ₂ — корни квадратного уравнения λ² — λ(a₁₁ + а₂₂) + (а₁₁а₂₂ — а²₁₂) = 0.

4. Изменение параметров γ_stat(t) и φ_stat(t) с течением времени

Выше мы предполагали, что момент времени t фиксирован. Сейчас мы рассмотрим поведение найденных величин γ_stat(t) и φ_stat(t) в зависимости от времени.

Это поведение можно определить из формул, приведенных в предыдущем разделе. В них входят величины L_i(t) и B_i(t), которые и порождают зависимость γ_stat(t) и φ_stat(t) от времени. Изменение долгот L_i(t) и широт B_i(t) со временем — вещь хорошо изученная. См. главу 1. Соответствующие, достаточно громоздкие, расчеты проделаны нами с помощью компьютера при численном нахождении зависимостей оценок γ_stat(t) и φ_stat(t) от времени. См. главу 6. Здесь мы пока ограничимся анализом лишь качественного поведения этих функций.

Рассмотрим вновь звездную сферу и будем здесь считать для простоты, что все звезды на ней неподвижны. Таким образом, мы возвращаемся сейчас к представлениям Птолемея, хотя делаем это только ради упрощения рассуждений и выкладок. Мы можем так поступить, поскольку количество звезд, имеющих заметную скорость собственного движения, — то есть смещающихся на расстояние нескольких дуговых минут за рассматриваемый промежуток времени в 2500 лет, — сравнительно невелико. Наличие таких звезд практически не влияет на оценки параметров γ_stat(t) и φ_stat(t), которыми мы сейчас занимаемся.

На рис. 5.4 изображена звездная сфера и реальная эклиптика эпохи t_A составления каталога. Полезно сравнить рис. 5.1 и рис. 5.4. В неизвестную нам эпоху t_A полюс эклиптики P(t_A) занимал некоторое, вполне определенное, положение на сфере. Составитель каталога, конечно же, отметил эклиптику на звездной сфере не идеально точно. Поэтому полюс Р_А отмеченной им «эклиптики каталога» занял положение, вообще говоря, отличное от P(t_A).

Рис. 5.4. Геометрия определения углов φ и γ на звездной сфере.

Проведем через полюс P(t_A) дугу большого круга, соединяющую его с точками весеннего равноденствия Q и осеннего равноденствия R. Дополнительно проведем через P(t_A) дугу большого круга D(t_A)D'(t_A), пересекающую только что построенную дугу QP(t_A)R под прямым углом в точке P(t_A). Если бы дата t_A была нам известна, то метод наименьших квадратов, описанный в разделе 3, позволил бы найти параметры γ и φ, определяющие взаимное расположение эклиптики эпохи t_A и эклиптики каталога. Из рис. 5.4 следует, что эти же углы определяют и взаимное положение полюсов P(t_A) и Р_А на звездной сфере, а именно, величина γ равна длине дуги P(t_A)P_A, в дуговых величинах, а угол φ равен углу P_AP(t_A)D'(t_A). С течением времени, как отмечалось в главе 1, положение эклиптики на звездной сфере изменяется. Это — эффект колебания эклиптики. Поэтому полюс эклиптики в момент t, отличный от t_A, окажется в точке P(t), также отличной от P(t_A). Траектория полюса эклиптики на звездной сфере показана на рис. 5.4 пунктирной линией, проходящей через точки P(t) и P(t_A). И тогда, чтобы совместить эклиптику эпохи t и эклиптику каталога, нужно совместить полюса Р_А и P(t). Длина дуги P(t)P_A равна величине γ_stat(t). Положение же оси вращения эклиптики, обеспечивающего данное совмещение, можно параметризовать углом P_AP(t)D'(t), где дуга D(t)D'(t) «параллельна» дуге D(t_A)D'(t_A).

Чтобы разобраться в качественном поведении функций γ_stat(t) и φ_stat(t), обратимся к плоскому рисунку, где изобразим лишь перемещение полюсов эклиптики. Это допустимо, поскольку величины их смешений заведомо лежат в пределах одного градуса. Перенесем с рис. 5.4 на плоскость картину вблизи северного полюса эклиптики, рис. 5.5.

Рис. 5.5. Определение качественной зависимости γ_stat(t) и φ_stat(t) от времени.

Как видно из рис. 5.5, реальный полюс эклиптики с течением времени перемещается вследствие колебания эклиптики. На рассматриваемом интервале времени данное перемещение составляет всего около 25′. Поэтому его можно изобразить отрезком прямой. См. пунктирную прямую на рис. 5.5. Движение полюса эклиптики вдоль этой прямой с большой точностью можно считать равномерным. Поэтому, например, расстояние между полюсами P(t) и P(t_A) равно v(t_A — t), где v — скорость движения полюса эклиптики. Эта скорость равна приблизительно 0,01′ в год. Как уже говорилось, в эпоху наблюдений t_А, из-за ошибки в положении эклиптики, сделанной составителем каталога, полюс эклиптики каталога попал в точку Р_А, отличную от P(t_A). Если при этом перпендикуляр, опущенный из точки Р_А на траекторию движения полюса эклиптики, пересек ее в точке t* > t_А, как изображено на рис. 5.5, то такая ошибка составителя, очевидно, «старит» эклиптику каталога, а именно, эклиптика каталога точнее всего будет отвечать эклиптике года t*. В противном случае, — то есть если указанный перпендикуляр пересек траекторию в точке t*< t_A — ошибка автора, напротив, «омолаживает» каталог. Чтобы дать представление о реальных соотношениях величин, укажем, что для Альмагеста расстояние между полюсом Р(0) эклиптики на 1900 год н. э. и полюсом Р(19) на начало нашей эры составляет около 20′. Приблизительно такое же значение имеет и ошибка γ_stat(t_A).

Как было сказано, величина γ_stat(t_A) равна длине отрезка Р(t_А)Р_А, а φ_stat(t_A) — углу Р_АР(t_А)D'(t_А). Аналогично, γ_stat(t_A) = P(t)Р_А. Здесь чертой сверху обозначена длина отрезка. Однако угол P_AP(t)D'(t) не равен φ_stat(t_A), поскольку к моменту t ось весеннего равноденствия сместилась на величину ω(t_A — t). Здесь ω — угловая скорость прецессии, равная приблизительно 50″ в год. См. главу 1. Это смещение соответствует на рис. 5.5 величине угла D'(t)P(t)S(t). Таким образом, φ_stat(t) = ∠P_AP(t)S(t), причем ∠D'(t)P(t)S(t) = ω(t_A — t).

Чтобы не использовать далее столь громоздкие обозначения, положим

x(t) = P(t)P(t_A),

y = P(t_A)P(t^*),

z = P_AP(t^*),

ψ(t) = ∠ P_AP(t)D'(t),

δ = ∠ D(t_A)P(t_A)P(t).

Величину γ_stat(t_A) можно назвать ошибкой определения эклиптики. Для Альмагеста она имеет порядок 20′. Угол δ не зависит от t и равен углу между направлением движения полюса эклиптики и определенной ранее прямой D(t_A)D'(t_A). Очевидно, что z = γ_stat(t_A) sin (δ — φ_stat(t_A)), y = γ_stat(t_A) cos (δ — φ_stat(t_A)).

Поскольку x(t) = v(t_A — t), то из рис. 5.5 следует, что

Очевидно, что эта функция достигает минимального значения при t = t*. Если же рассматривается случай |t — t_А| << |t_А — t*|, то функция γ_stat(t) ведет себя практически как линейная: γ_stat(t) = γ_stat(t_A) + v cos (δ — φ_stat(t_A))(t_A — t).

Нетрудно найти также функцию φ_stat(t):

И вновь, если |t — t_А| << |t_А — t*|, можно воспользоваться линейным приближением:

Разумеется, найденные формулы дают лишь общее представление о характере функций γ_stat(t) и φ_stat(t). На рис. 5.6 изображен примерный вид этих функций, получаемый из формул (5.4.1) и (5.4.2). Естественно, конкретный их вид зависит от значения ошибки, совершенной составителем каталога, то есть, от величин γ_stat(t_A) и φ_stat(t_A). Формулы (5.4.1) и (5.4.2) определяют также вид зависимости β_stat(t). См. формулу (5.2.1).

Рис. 5.6. Примерный вид функций γ_stat(t) и φ_stat(t).

Обсудим геометрический смысл данных построений. Рассмотрим птолемеевские координаты какой-либо группы звезд, считая что его наблюдения выполнены в момент времени t. В этом предположении, устраним систематическую ошибку γ_stat(t) и φ_stat(t), то есть повернем всю группу звезд на угол γ_stat(t) вокруг оси, отстоящей от оси равноденствия на угол φ_stat(t). Для простоты предположим, что систематическую ошибку мы нашли совершенно точно. Тогда полюс эклиптики каталога Р_А совместится с реальным полюсом P(t). Разумеется, после такого совмещения широтные невязки звезд все равно не станут равными нулю, так как в каталоге присутствуют еще случайные ошибки. Однако случайные ошибки, имея нулевое среднее, не смещают положение полюса эклиптики. Вернее, смещают ее лишь на малую величину, которая тем меньше, чем больше рассматриваемая совокупность звезд.

Из рис. 5.5 видно, что перемещение полюса Р_А в точку P(t) разлагается единственным способом в композицию двух перемещений: Р_А P(t_A) и P(t_A) P(t). Параметры γ_stat(t_A) и φ_stat(t_A), задающие первое перемещение, имеют смысл ошибки наблюдателя, а именно — той ошибки, которую совершил составитель каталога в определении положения плоскости эклиптики. Второе перемещение обусловлено вековым колебанием плоскости эклиптики. Это колебание можно рассчитать по теории Ньюкомба.

Из сказанного вытекает также следующий вывод. Обозначим через ΔB_i(t) широтную невязку i-й звезды, рассчитанную на момент t предполагаемых наблюдений, а через ΔB⁰_i(t) = ΔB_i(t) — γ_stat(t) sin (L_i(t) + φ_stat(t)) — ее широтную невязку на момент t после компенсации систематической ошибки. Тогда для совокупности, состоящей из полностью неподвижных звезд, величины ΔB⁰_i(t) не зависят от t и равны случайным ошибкам, допущенным Птолемеем при определении широт. Ситуация меняется, если в рассматриваемую совокупность входят подвижные звезды. Для них величины ΔB⁰_i(t) будут зависеть от времени t. Характер зависимости определяется как величинами индивидуальных случайных ошибок, так и направлением скоростей собственного движения звезд в совокупности. В частности, в неизвестную нам эпоху t_А величина ΔB⁰_i(t_A) равна случайной широтной ошибке для звезды i. Естественно ожидать, что если эта звезда быстро движется и, кроме того, хорошо измерена, то величина |ΔB⁰_i(t)| достигает минимума в окрестности точки t_А. Величина этой окрестности зависит от величины и направления скорости собственного движения звезды и даже для самых быстрых звезд, например Арктура, составляет сотни лет.

Из приведенного выше рассуждения и, в частности, из рис. 5.5, следует важный вывод. А именно, для определения полюса эклиптики каталога Р_А достаточно знать лишь два значения γ_stat, соответствующих различным значениям моментов времени t₁ и t₂.

В самом деле, из теории Ньюкомба нетрудно найти скорость перемещения полюса эклиптики v. См. главу 1. Зафиксируем два произвольных различных момента времени t₁ и t₂. См. рис. 5.7. С помощью формулы (5.3.7) найдем значения γ_stat(t₁) и φ_stat(t₂). Изобразим прямую, по которой перемешается со временем полюс эклиптики. Отметим на ней точки t₁ и t₂. Выберем такой масштаб, чтобы расстояние между отмеченными точками равнялось v|t₂ — t₁|. Положение полюса эклиптики каталога Р_А определяется как точка пересечения двух окружностей с центрами в точках t_i и радиусами γ_stat(t_i), i = 1,2. Из рис. 5.7 ясно, как при этом определяются величины γ_stat(t) и φ_stat(t) для произвольного значения времени t. Необходимо лишь сказать, что прямая S'S, от которой отсчитывается угол φ_stat(t), пересекает траекторию движения полюса эклиптики под углом δ(t). Угол этот также находится из теории Ньюкомба. Астрономический смысл прямой S'S весьма нагляден. Это «спрямленная» часть большого круга звездной сферы, проходящего через полюс эклиптики Р(t) эпохи t и перпендикулярного в точке Р(t) другому большому кругу, также проходящему через Р(t) и точки равноденствия эпохи t.

Рис. 5.7. Нахождение величин γ_stat(t) и φ_stat(t).

Аналогично, для определения параметров γ_stat(t) и φ_stat(t) при всех t достаточно знать лишь два значения: γ_stat(t₁) и φ_stat(t₂).

Мы, однако, будем работать с углом γ. Он имеет содержательный смысл: это — ошибка в определении угла наклона между плоскостями экватора и эклиптики. Заметим, что этот угол фиксируется, например, в армиллярной сфере. Следовательно, ошибка γ, допущенная в этом угле, может являться инструментальной ошибкой армиллярной сферы. См. главу 1. Таким образом, ошибка γ естественным образом возникает при астрономических измерениях. Кроме того, выбор γ в качестве параметра будет нами далее обоснован и со статистической точки зрения.

5. Статистические свойства оценок γ_stat и φ_stat

Сейчас мы подойдем к задаче оценки параметров γ и φ, задающих систематическую ошибку каталога, как к задаче статистической. Для этого поступим так. Пусть составитель каталога в момент времени t_А совершил систематическую ошибку, задаваемую параметрами γ_А и φ_А. Пусть, кроме того, широта каждой измеренной им звезды подвергалась, — вследствие ошибки наблюдения, — случайному возмущению ξ имеющему нулевое среднее, то есть Eξ_i = 0. Предполагается, что случайные погрешности ξ_i, отвечающие различным звездам, независимы и имеют одно и то же распределение. Пусть σ² = Еξ²_i — дисперсия случайной величины ξ_i. Эта дисперсия нам, вообще говоря, неизвестна.

В этих предположениях широта i-той звезды в каталоге будет иметь вид

b_i = В_i(t_A) — γ_А sin (L_i(t_А) + φ_A) + ξ_i. (5.5.1)

Со статистической точки зрения мы имеем выборку, состоящую из N реализаций случайных величин {b_i}^N_i=1, вида (5.5.1). По этой выборке требуется определить статистические оценки γ' и φ' параметров γ_А и φ_А, а также оценить величину σ, представляющую собой среднеквадратичную ошибку наблюдения. Мы сразу ограничим задачу и будем изучать оценки φ' = φ_stat и γ' = γ_stat, получаемые методом наименьших квадратов. Эти оценки имеют вид (5.3.6), (5.3.7). Основное внимание будет уделено оценке величины γ_А по причинам, объясненным в конце раздела 4.

Равенство (5.5.1) имеет вид, традиционный для регрессионного анализа. В самом деле, это равенство утверждает, что ошибка наблюдения Δb_i = B_i(t_A) — b_i является случайной величиной со средним γ_А sin (L_i(t_A) + φ_А), зависящим от неизвестных параметров γ_А и φ_А, и дисперсией σ². Требуется оценить значения неизвестных параметров методом наименьших квадратов и установить статистические свойства полученных оценок. В такой постановке кривую Y(x) = γ_A sin (x + φ_А) обычно называют линией регрессии.

Определим величины φ и γ с помощью соотношений (5.3.6) и (5.3.7). Отклонения Δb_i случайны по предположению. Поэтому и получаемые из соотношений (5.3.6) и (5.3.7) оценки φ_stat и γ_stat, также являются случайными величинами. Изучим их статистические свойства и рассмотрим, как они связаны с истинными, но неизвестными нам, значениями φ_А и γ_А.

Подставим в приведенные выше формулы для s_b и c_b вместо Δb_i разность γ_А × sin (L_i(t_A) + φ_А) — ξ_i и используем эту подстановку в формулах (5.3.6), (5.3.7). Получим следующие выражения для величин φ_stat и γ_stat:

Введем величину R = (γ_A(d² — s₂ c₂) cos φ_А)^-1

Тогда (5.5.2) можно записать в виде

Из условия Eξ_i = 0 находим, что полученная оценка параметра γ_stat является несмещенной, то есть Eγ_stat = γ_A.

Дисперсия же оценки γ_stat, обозначаемая D_γ, имеет вид

Если ошибки наблюдения ξ_i нормально распределены, то и величина γ_stat является нормально распределенной, и первые два момента (5.5.5) и (5.5.6) полностью определяют ее распределение. Этот факт даст нам возможность построить доверительный интервал для значения γ_А. Анализ оценки φ_stat несколько более сложен. Воспользуемся следующим равенством, получаемым из (5.5.4):

и тем фактом, что при больших N второе слагаемое в знаменателе правой части (5.5.7) — величина малая. В самом деле, эта величина — случайная с нулевым средним и дисперсией

Если ξ_i нормально распределены, то и рассматриваемая величина также нормально распределена. Из этого факта для каталога Альмагеста следует, что уже для N = 30 вероятность p_n того, что знаменатель правой части (5.5.7) будет отрицательным, не превосходит 5 × 10^-3. С ростом N данная вероятность быстро убывает: p₅₀ ≤ 2,5 × 10^-4, p₈₀ ≤ 4 × 10^-6, p₁₀₀ ≤ 3 × 10^-7, p₂₀₀ ≤ 8 × 10^-13, р₃₀₀ ≤ 2,5 × 10^-18.

Из формулы (5.5.7) следует, что, вообще говоря, E tan φ_stat tan φ_A. Однако из этой формулы легко получить функцию распределения F(x) случайной величины tan φ_stat — tan φ_A, необходимую при нахождении доверительного интервала для φ_А. В самом деле, если пренебречь тем маловероятным случаем, что знаменатель в (5.5.7) становится отрицательным, то из этой формулы получается выражение для F(x):

F(x) = P(tan φ_stat — tan φ_A< x) = Р(η_x< x),

где случайная величина η_х имеет вид

Следовательно, если величины ξ_i нормально распределены с дисперсией σ², то и величина η_x имеет гауссовское распределение со средним, равным нулю, и дисперсией

Следовательно,

где

Найденные выше значения γ_stat и φ_stat являются, как говорят, точечными оценками неизвестных параметров γ_А и φ_А. Поскольку найдены функции распределения этих оценок, то можно исследовать вопрос об их возможной погрешности. Дадим ответ на этот вопрос в стандартных терминах доверительных интервалов, основываясь на формулах (5.5.5), (5.5.6), (5.5.8), (5.5.9).

В математической статистике задача нахождения доверительного интервала порождена ситуацией, которую поясним на примере оценки величины γ_А. Эта величина является вполне определенной, детерминированной ошибкой, сделанной составителем каталога. В результате статистической оценки γ_А, — в нашем случае по методу наименьших квадратов — получается случайная величина γ_stat. Возникает вопрос, какие границы можно указать для неизвестной нам величины γ_А, если мы определили γ_stat?

Чтобы границы эти не оказались тривиальными, необходимо задать допустимую вероятность ошибки, то есть, вероятность указать такие границы, в которых истинное значение γ_А не лежит. Обозначим допустимую вероятность ошибки через ε. Тогда уровень доверия будет равен 1 — ε. Случайная величина γ_stat распределена по нормальному закону с параметрами, задаваемыми формулами (5.5.5) и (5.5.6). Поэтому при x > 0 имеем

Определим величину ε/2-квантили нормального распределения x_ε из уравнения

или, что то же, из уравнения Ф(-√D_γx_ε) = ε/2.

Тогда интервал I_γ(ε) = (γ_stat — x_ε, γ_stat + x_ε). (5.5.10)

представляет собой доверительный интервал для γ_А с уровнем доверия 1 — ε. Это следует из того, что P(|γ_stat — γ_А| ≥ x_ε) = ε.

При определении величины x_ε мы, в частности, использовали значение D_γ, которое зависит от неизвестных нам параметров σ² и φ_А. Как это обычно делается в математической статистике, вместо σ² подставим в формулу для D_γ сходящуюся к ней остаточную дисперсию

определяемую формулой (5.3.3), а вместо φ_А — величину φ_stat. Момент t_А составления каталога нам также неизвестен, поэтому все перечисленные выше вычисления необходимо проделать для всех моментов времени t с тем, чтобы оценить систематическую ошибку γ_stat(t), φ_stat(t) при условии, что каталог был составлен в произвольную фиксированную эпоху t. Аналогичным образом можно найти доверительный интервал для φ_А с уровнем доверия 1 — ε. Этот интервал I_φ(ε) будет таким:

где y_ε — решение уравнения F(y_ε) — F(-y_ε) = 1 — ε, в котором функция распределения F задана равенством (5.5.9), то есть ε/2-квантиль соответствующего нормального распределения.

Замечание. Полученные выше оценки истинных значений ошибок γ и φ в каталоге, как функций предполагаемой датировки, важны не только для того, чтобы их скомпенсировать, но и для косвенной проверки правильности предлагаемого подхода. Например, если бы в качестве γ_stat получилась величина, в несколько раз превышающая точность каталога, это указывало бы на какие-то неучтенные нами существенные эффекты.

Однако если речь идет лишь о датировке, то само значение γ_stat в соответствующей процедуре не участвует. Нам необходимо лишь знание длины соответствующего доверительного интервала. Поэтому возможно существенное упрощение вычислений, состоящее в следующем. Вычисляются γ_stat и φ_stat, относящиеся к любому фиксированному моменту времени t₀. Например, к 1900 году, для чего не требуется использования уравнений Ньюкомба. Тогда вместо кривых γ_stat(t) и φ_stat(t) мы получим постоянные значения, соответствующие ошибкам наблюдений, но только не в координатах эпохи наблюдений, а в координатах эпохи 1900 года. Затем вокруг этих постоянных значений откладываются доверительные интервалы, ширина которых от t не зависит. В результате статистической процедуры датировки, описываемой ниже, будет получен тот же интервал возможных датировок каталога, что и при оценивании ошибок γ и φ относительно координат на эпоху предполагаемой датировки t. Единственная информация, которая при этом будет потеряна, — это оценки истинных значений величин γ_stat и φ_stat.

6. Выводы

ВЫВОД 1. Групповая ошибка звездной конфигурации сводится к перемещению этой конфигурации как единого целого по небесной сфере. Данное перемещение, при учете лишь широтных невязок, можно параметризовать двумя параметрами, а именно, γ и φ, либо γ и β.

ВЫВОД 2. Существующие в каталоге широтные невязки могут быть уменьшены за счет компенсации групповых ошибок.

ВЫВОД 3. Если в большой части каталога групповые ошибки совпадают, то эта общая ошибка называется систематической и может быть обнаружена статистическими методами.

При условии, что каталог составлен в эпоху t значения параметров φ(t) и γ(t) оцениваются методом наименьших квадратов. Соответствующие оценки γ_stat(t) и φ_stat(t) имеют вид (5.3.6) и (5.3.7) соответственно.

ВЫВОД 4. Знания значений γ_stat(t₁) и φ_stat(t₂) для двух различных моментов времени достаточно для восстановления функций γ_stat(t) и φ_stat(t).

ВЫВОД 5. В предположении нормального распределения случайных ошибок измерения найдены доверительные интервалы I_φ(ε) и I_γ(ε) для истинных значений параметров φ(t) и γ(t). См. формулы (5.5.11) и (5.5.10) соответственно.

В заключение, на рис. 5.8, приведем страницу из Альмагеста издания 1551 года.

Рис. 5.8. Страница из издания Альмагеста 1551 года.

Глава 6