В.В.Калашников, Г.В.Носовский, А.Т.Фоменко
ЗВЕЗДЫ ЗОДИАКА
Первоисточник - "Астрономический анализ хронологии. Альмагест. Зодиаки", 2000 г.

Глава 5

АНАЛИЗ СИСТЕМАТИЧЕСКИХ ОШИБОК ЗВЕЗДНЫХ КОНФИГУРАЦИЙ В ПРОИЗВОЛЬНОМ КАТАЛОГЕ, НАПРИМЕР В АЛЬМАГЕСТЕ
5. Основные обозначения

Начиная с этой главы, мы считаем, что имеем дело с каталогом, все звезды которого имеют ЕДИНСТВЕННОЕ отождествление со звездами из современного каталога. В соответствии с этим, будем идентифицировать звезды индексом i и обозначать через li,bi соответственно эклиптикальные долготу и широту i-й звезды в Альмагесте. Через Li(t),Bi(t) мы обозначим истинные долготу и широту i-й звезды в эпоху t. Напомним, что время t мы отсчитываем от 1900 года "назад"и измеряем в столетиях. То есть, например, t = 3, 15 соответствует году 1900 - 3, 15 ⋅ 100 = 1585 году н.э., а t = 22, 0 отвечает году 1900 - 22 ⋅ 100 = 300 году до н.э.

Пусть tA – неизвестное нам время составления каталога Альмагеста. Обозначим через LiA,B iA истинные долготу и широту i-й звезды в год составления каталога, то есть LiA = L i(tA),BiA = B i(tA). Пусть Bi(t) = Bi(t) - bi - разность между истинной широтой i-й звезды в момент времени t и ее широтой в Альмагесте. Назовем величину Bi(t) широтной невязкой на момент времени t. Эта величина имеет смысл погрешности в определении широты i-й звезды Альмагеста ПРИ УСЛОВИИ, ЧТО ОН БЫЛ СОСТАВЛЕН В ЭПОХУ t. Естественно, что Bi(tA) = BiA представляет собой истинную погрешность в определении широты.

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

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

Рассмотрим некоторую совокупность звезд, например, созвездие или группу созвездий. Определим групповую ошибку в широтных координатах этих звезд как погрешность в определении широт звезд из этой совокупности, проистекающую из перемещения рассматриваемой звездной конфигурации КАК ЕДИНОГО ЦЕЛОГО по сфере. Следовательно, – и это обстоятельство мы особо подчеркнем, так как оно существенно используется в дальнейшем, – ЛЮБОЕ ПОДМНОЖЕСТВО ЭТОЙ КОНФИГУРАЦИИ ТАКЖЕ ПЕРЕМЕЩАЕТСЯ ПО СФЕРЕ КАК ЕДИНОЕ ЦЕЛОЕ И НА ТОТ ЖЕ УГОЛ, ЧТО И ВСЯ КОНФИГУРАЦИЯ. Такое перемещение имеет три степени свободы, то есть может быть описано путем задания трех параметров. Определим их.

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

Этой точке соответствует возмущенная эклиптика, которая названа наРис. 5.1 ЭКЛИПТИКОЙ КАТАЛОГА. Ее положение можно однозначно задать, определив следующие два параметра. Во-первых, угол между прямыми OP и OPA, или, что то же самое, плоский угол между плоскостями реальной эклиптики и эклиптики каталога. Во-вторых, угол между прямой равноденствий RQ и прямой CD, образованной пересечением плоскостей реальной эклиптики и эклиптики каталога. Такая параметризация удобна для аналитических целей. Однако далее мы будем иногда наряду с использовать величину b, которая интерпретируется следующим образом,Рис. 5.1. Смещение эклиптики "разлагается"на два поворота – вокруг оси равноденствия RQ на угол и вокруг оси, лежащей также в плоскости эклиптики и перпендикулярной оси RQ, на угол b. Итак, b представляет собой величину дуги QAQ, являющейся частью большого круга, проходящего через полюс PA и точку Q. Астрономический смысл точки QA весьма прозрачен. Это точка весеннего равноденствия на эклиптике каталога. Ясно, что углы и однозначно определяют углы и b. Как и наоборот. Искомую связь мы находим из рассмотрения сферического прямоугольного треугольника CQAQ. Здесь угол при вершине QA – прямой, угол при вершине C равен , а длина дуги CQ равна b. В результате получаем:
sin b = sin g sin f.
(5. 1)


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

Посмотрим теперь, как искажаются истинные координаты звезды i при наличии такой систематической ошибки. Истинные широта BiA и долгота этой звезды LiA равны соответственно длинам дуг EE и QE, отсчитываемым по часовой стрелке, если смотреть с полюса P. Искаженные широта и долгота bi,li равны соответственно длинам дуг EEA и QAEA. Отметим, что широты звезд, у которых истинная долгота больше долготы точки D и меньше долготы точки C, уменьшаются. А остальные – увеличиваются,Рис. 5.1. Этот вывод справедлив, строго говоря, не для всех звезд. Он неверен для звезд, расположенных вблизи полюсов P и P и находящихся на угловом расстоянии (или менее) от них. Однако, учитывая малость искажений, вносимых составителем каталога, таких звезд практически не существует. Как мы увидим, величина составляет около 20.

Учитывая малость величины , можно предложить следующую приближенную формулу для широтной невязки:
ΔBAi = g ⋅ sin(LAi + f).
(5. 2)

Иначе говоря, систематическая погрешность определения широт звезд может быть изображена в виде синусоиды, показанной наРис. 5.2. Полезно сравнить ее с кривой, полученной К.Петерсом и Е.Кнобелем [92]. Погрешность формулы (5.2.2) не превышает 1 для звезд, у которых ∣bA∣≤ 80o. Такая погрешность является для нас несущественной, и мы в дальнейшем говорить о ней не будем, считая формулу (5.2.1) абсолютно точной. Для корректности мы исключим из рассмотрения звезды, имеющие широты более 80 градусов по абсолютной величине. В дальнейшем в этой главе речь будет идти о СИСТЕМАТИЧЕСКОЙ ПОГРЕШНОСТИ, поскольку излагаемые методы работают в предположении, что рассматривается БОЛЬШАЯ СОВОКУПНОСТЬ звезд. Проверка того, что найденная погрешность совпадает (или не совпадает) с групповыми ошибками отдельных созвездий, представляет собой самостоятельную задачу. Применительно к Альмагесту она рассматривается ниже, в главе 6.
Рис. 5. 2: Зависимость систематической широтной невязки от долготы

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

1) Вычислим на момент времени tA истинные широты BiA и долготы L iA для всех звезд из рассматриваемой совокупности.

2) найдем параметры * и *, дающие решение задачи
s2(g*, f*) → min,
(5. 3)


где

sum s2(g,f) = (BAi - bi- g sin(LAi + f))2. i


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

В нашей ситуации момент tA составления каталога неизвестен. Поэтому мы вынуждены рассчитывать систематические ошибки для всех значений t из рассматриваемого интервала 0 ≤ t ≤ 25. А именно, для каждого t определяются положение реальной эклиптики и ось равноденствия. Затем, как и наРис. 5.1, вводятся параметры = (t), = (t) и b = b(t), задающие взаимное положение эклиптики эпохи t и эклиптики каталога. Значения величин (t) и (t) находятся как решение задачи
s2(g(t),f(t),t) → min,
(5. 4)


где
sum s2(g, f,t) = (ΔBi(t) - g sin(Li(t) + f))2. i
(5. 5)


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


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

В реальном каталоге, кроме указанной систематической ошибки, присутствуют и случайные ошибки. Поэтому отклонения Bi(t) -bi являются случайными величинами, значения которых разбросаны вокруг синусоиды их среднего значения, изображенной наРис. 5.2. В предположении, что остальные погрешности каталога, кроме систематических, носят случайный характер, задача определения (t) и (t) является задачей определения параметров регрессии. Решением этой задачи являются статистические оценки параметров (t) и (t).

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

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

s2(g,f, t) = 1-s2(g,f, t), s (t) = -1 sum N sin2 L (t), 0 N 2 N i=1 i 1 sum N sum sb(t) = --- ΔBi(t) sin Li(t), c2(t) = N1 Ni=1cos2 Li(t), N i=1 1 N sum sum cb(t) = --- ΔBi(t) cos Li(t), d(t) = N1 Ni=1 sin Li(t)cosLi(t). N i=1


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

Очевидно, что задача минимизации (5.2.4) эквивалентна задаче минимизации
2 s0(g,f,t) → min
(5. 6)
в том смысле, что параметры (t) и (t), определяемые соотношением (5.3.1), совпадают с параметрами, определяемыми решением задачи (5.2.4).

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

Значение
smin(t) = s0(gstat(t),fstat(t),t)
(5. 7)


имеет прозрачный физический смысл. Это среднеквадратичная широтная невязка по рассматриваемой совокупности звезд в момент времени t, получившаяся ПОСЛЕ КОМПЕНСАЦИИ НАЙДЕННОЙ СИСТЕМАТИЧЕСКОЙ ОШИБКИ stat(t),stat(t). Как мы увидим далее, величина smin(t) от времени практически не зависит ввиду крайне малой скорости собственного движения подавляющего большинства звезд. Поэтому мы будем использовать также обозначение smin. Заметим, что ДО КОМПЕНСАЦИИ ЭТОЙ ОШИБКИ среднеквадратичная широтная невязка в момент t равнялась величине
---------------- N sum sinit(t) = s0(0,0,t) = V~ 1- (ΔBi(t))2 , N i=1
(5. 8)
которая, вообще говоря, зависит от t. Таким образом, разность s(t) = sinit(t) - smin(t) оценивает эффект от компенсации систематической ошибки stat(t), stat(t).

Далее при определении величин stat(t) и stat(t) из соотношения (5.3.1) момент времени t будем предполагать фиксированным. Поэтому аргумент t в выкладках мы опускаем. То есть, будем писать Li вместо Li(t), sb вместо sb(t) и т.д.

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

s cos f + c sinf = g[s cos2f + 2d cosf sin f + c sin2 f], (5. 9) b b 2 2 2 2 - cbcos f + sbsinf = g[- d cos f + (s2- c2)cosf sin f + d sin f(]5. .10)


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

2 -sb +-cb-tan-f-= ---s2-+-2d-tan-f-+-c2-tan-f---- . - cb + sbtan f - d + (s2- c2)tan f + d tan2f


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

(1 + tan2 f)(cbs2- sbd) + (1 + tan2 f)tan f(cbd - sbc2) = 0.


Отсюда легко найти тангенс оптимального значения stat:
sbd---cbs2 tan fstat = cbd - sbc2 .
(5. 11)


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



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

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



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


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



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

Формула (5.3.9) позволяет элементарными средствами оценить чувствительность среднеквадратичной ошибки s(,,t) к вариациям параметров и . Для этого достаточно определить величины a11,a12 и a22, входящие в правую часть (5.3.9). После вычисления оценок stat(t) и stat(t) их легко найти по формуле (5.3.8).

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


Формула (5.3.9) показывает, что "линии уровня"среднеквадратичных ошибок являются эллипсами на плоскости (,). См. рис. 5.3. Центром эллипса является точка (stat,stat), в которой значение среднеквадратичной ошибки равно smin. Направления осей эллипсов и соотношение между ними определяются стандартными формулами аналитической геометрии через величины a11,a12,a22. А именно, угол наклона a одной из осей эллипса определяется соотношением:

(2a12) tan 2a = a-----a--. 11 22


Вторая ось перпендикулярна ей. Длины осей относятся друг к другу как c1/c2, где c1 и c2 – корни квадратного уравнения

c2 - c(a11 + a22) + (a11a22 - a212) = 0.

5. 4 Изменение параметров stat(t) и stat(t) с течением времени

Выше мы предполагали, что момент времени t является фиксированным. Сейчас мы рассмотрим поведение найденных величин stat и stat в зависимости от времени.

Это поведение можно определить из формул, приведенных в предыдущем разделе. В эти формулы входят величины Li(t) и Bi(t), которые и порождают зависимость stat и stat от времени. Изменение долгот Li(t) и широт Bi(t) со временем - вещь хорошо изученная. См. главу 1. Соответствующие, достаточно громоздкие, расчеты были фактически проделаны нами с помощью компьютера при численном нахождении зависимостей оценок stat(t) и stat(t) от времени. См.главу 6. Здесь мы пока ограничимся анализом лишь качественного поведения этих функций.

Рассмотрим вновь звездную сферу и будем здесь считать для простоты, что все звезды на ней НЕПОДВИЖНЫ. Таким образом, мы возвращаемся сейчас к представлениям Птолемея, хотя делаем это только ради упрощения рассуждений и выкладок. Мы можем так поступить, поскольку количество звезд, имеющих заметную скорость собственного движения, - то есть смещающихся на расстояние нескольких дуговых минут за рассматриваемый промежуток времени в 2500 лет, - сравнительно невелико. Наличие таких звезд практически не влияет на оценки параметров stat(t) и stat(t), которыми мы сейчас занимаемся.

Рис. 5. 4: Геометрия определения углов и на звездной сфере


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

Проведем через полюс P(tA) дугу большого круга, соединяющую его с точками весеннего равноденствия Q и осеннего равноденствия R. Дополнительно проведем через P(tA) дугу большого круга D(tA)D(tA), пересекающую только что построенную дугу QP(tA)R под прямым углом в точке P(tA). Если бы дата tA была нам известна, то метод наименьших квадратов, описанный в разделе 3, позволил бы найти параметры и , определяющие взаимное расположение эклиптики эпохи tA и эклиптики каталога. Из рис. 5.4 следует, что эти же углы определяют и взаимное положение полюсов P(tA) и PA на звездной сфере. А именно, величина равна длине дуги P(tA)PA, в дуговых величинах, а угол равен углу PAP(tA)D(tA). С течением времени, как отмечалось в главе 1, положение эклиптики на звездной сфере изменяется. Это – эффект колебания эклиптики. Поэтому полюс эклиптики в момент t, отличный от tA, окажется в точке P(t), также отличной от P(tA). Траектория полюса эклиптики на звездной сфере показана на рис. 5.4 пунктирной линией, проходящей через точки P(t) и P(tA). И тогда, чтобы совместить эклиптику эпохи t и эклиптику каталога, нужно совместить полюса PA и P(t). Длина дуги P(t)PA равна величине stat(t). Положение же оси вращения эклиптики, обеспечивающего данное совмещение, можно параметризовать углом PAP(t)D(t), где дуга D(t)D(t) "параллельна"дуге D(tA)D(tA).

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

Рис. 5. 5: Определение качественной зависимости stat(t) и stat(t) от времени


Как видно из рис. 5.5, реальный полюс эклиптики с течением времени перемещается вследствие колебания эклиптики. На рассматриваемом интервале времени данное перемещение составляет всего около 25. Поэтому его можно изобразить отрезком прямой. См. пунктирную прямую на рис. 5.5. Движение полюса эклиптики вдоль этой прямой с большой точностью можно считать равномерным. Поэтому, например, расстояние между полюсами P(t) и P(tA) равно v(tA - t), где v – скорость движения полюса эклиптики. Эта скорость равна приблизительно 0, 01 в год. Как уже говорилось, в эпоху наблюдений tA, из-за ошибки в положении эклиптики, сделанной составителем каталога, полюс эклиптики каталога попал в точку PA, отличную от P(tA). Если при этом перпендикуляр, опущенный из точки PA на траекторию движения полюса эклиптики, пересек ее в точке t* > t A, как изображено на рис. 5.5, то такая ошибка составителя очевидно "старит"эклиптику каталога. А именно, эклиптика каталога точнее всего будет отвечать эклиптике года t*. В противном случае, – то есть если указанный перпендикуляр пересек траекторию в точке t* < t A, – ошибка автора, напротив, "омолаживает"каталог. Чтобы дать представление о реальных соотношениях величин, укажем, что для Альмагеста расстояние между полюсом P(0) эклиптики на 1900 год н.э. и полюсом P(19) на начало нашей эры составляет около 20. Приблизительно такое же значение имеет и ошибка stat(tA).

Как было сказано, величина stat(tA) равна длине отрезка P(tA)PA, а stat(tA) – углу PAP(tA)D(tA). Аналогично, stat(t) = P(t)PA. Здесь чертой сверху обозначена длина отрезка. Однако угол PAP(t)D(t) не равен stat(t), поскольку к моменту t ось весеннего равноденствия сместилась на величину w(tA -t). Здесь w – угловая скорость прецессии, равная приблизительно 50 в год. См. главу 1. Это смещение соответствует на рис. 5.5 величине угла D(t)P(t)S(t). Таким образом, stat(t) = ⁄PAP(t)S(t), причем ⁄D(t)P(t)S(t) = w(tA - t).

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

---------- -------- x(t) = P (t)P (tA), y = P(tA)P (t*), z = PAP (t*),
′ y(t) = ⁄-PAP (t)D (t), d = ⁄ D(tA)P (tA)P (t).


Величину stat(tA) можно назвать ОШИБКОЙ ОПРЕДЕЛЕНИЯ ЭКЛИПТИКИ. Для Альмагеста она имеет порядок 20. Угол d не зависит от t и равен углу между направлением движения полюса эклиптики и определенной ранее прямой D(tA)D(tA). Очевидно, что

z = gstat(tA) sin(d- fstat(tA)), y = gstat(tA) cos(d- fstat(tA)).


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

V~ ---------------2---2- gstat(t) = (v(tA - t) + y) + z V~ --2-----------------------2-------2 = g stat(tA) + 2yv(tA - t) + v (tA - t) . (5. 15)


Очевидно, что эта функция достигает минимального значения при t = t*. Если же рассматривается случай ∣t-tA∣≪∣tA -t*∣, то функция stat(t) ведет себя практически как линейная:

gstat(t) ≈ gstat(tA) + vcos(d - fstat(tA))(tA - t).


Нетрудно найти также функцию stat(t):
( ) z fstat(t) = d + w(tA - t)- arctg ------------- . y + v(tA - t)
(5. 16)
 

И вновь, если ∣t - tA∣≪∣tA - t*∣, можно воспользоваться линейным приближением:

[ ] fstat(t) = fstat(tA) + w + v-sin(d---fstat(tA))- (tA - t). gstat(tA)


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

Рис. 5. 6: Примерный вид функций stat(t) и stat(t)


Обсудим геометрический смысл данных построений. Рассмотрим птолемеевские координаты какой-либо группы звезд, считая что его наблюдения были выполнены в момент времени t. В этом предположении, устраним систематическую ошибку stat(t), stat(t), то есть повернем всю группу звезд на угол stat(t) вокруг оси, отстоящей от оси равноденствия на угол stat(t)). Для простоты предположим, что систематическую ошибку мы нашли совершенно точно. Тогда полюс эклиптики каталога PA совместится с реальным полюсом P(t). Разумеется, после такого совмещения широтные невязки звезд все равно не станут равными нулю, так как в каталоге присутствуют еще случайные ошибки. Однако случайные ошибки, имея нулевое среднее, НЕ СМЕЩАЮТ ПОЛОЖЕНИЕ ПОЛЮСА ЭКЛИПТИКИ. Вернее, смещают ее лишь на малую величину, которая тем меньше, чем больше рассматриваемая совокупность звезд.

Из рис. 5.5 видно, что перемещение полюса PA в точку P(t) разлагается ЕДИНСТВЕННЫМ способом в композицию двух перемещений: PA → P(tA) и P(tA) → P(t). Параметры stat(tA) и stat(tA), задающие первое перемещение, имеют смысл ошибки наблюдателя. А именно – той ошибки, которую совершил составитель каталога в определении положения плоскости эклиптики. Второе перемещение обусловлено вековым колебанием плоскости эклиптики, Это колебание можно рассчитать по теории Ньюкомба.

Из сказанного вытекает также следующий вывод. Обозначим через Bi(t) широтную невязку i-й звезды, рассчитанную на момент t предполагаемых наблюдений, а через Bi0(t) = B i(t) - stat(t) sin(Li(t) + stat(t)) – ее широтную невязку на момент t после компенсации систематической ошибки. Тогда ДЛЯ СОВОКУПНОСТИ, СОСТОЯЩИЙ ИЗ ПОЛНОСТЬЮ НЕПОДВИЖНЫХ ЗВЕЗД, ВЕЛИЧИНЫ Bi0(t) НЕ ЗАВИСЯТ ОТ t И РАВНЫ СЛУЧАЙНЫМ ОШИБКАМ, ДОПУЩЕННЫМ СОСТАВИТЕЛЕМ КАТАЛОГА ПРИ ОПРЕДЕЛЕНИИ ШИРОТ. Ситуация меняется, если в рассматриваемую совокупность входят ПОДВИЖНЫЕ звезды. Для них величины Bi0(t) будут зависеть от времени t. Характер зависимости определяется как величинами индивидуальных случайных ошибок, так и направлением скоростей собственного движения звезд в совокупности. В частности, в неизвестную нам эпоху tA величина Bi0(t A) равна случайной широтной ошибке для звезды i. Естественно ожидать, что если эта звезда быстро движется и, кроме того, хорошо измерена, то величина ∣Bi0(t)∣ достигает минимума в окрестности точки tA. Величина этой окрестности зависит от величины и направления скорости собственного движения звезды и даже для самых быстрых звезд, например Арктура, составляет сотни лет.

Из приведенного выше рассуждения и, в частности, из рис. 5.5, следует важный вывод. А именно: ДЛЯ ОПРЕДЕЛЕНИЯ ПОЛЮСА ЭКЛИПТИКИ КАТАЛОГА PA ДОСТАТОЧНО ЗНАНИЯ ЛИШЬ ДВУХ ЗНАЧЕНИЙ stat, СООТВЕТСТВУЮЩИХ РАЗЛИЧНЫМ ЗНАЧЕНИЯМ МОМЕНТОВ ВРЕМЕНИ t1 И t2.

Рис. 5. 7: Нахождение величин stat(t) и stat(t)


В самом деле, из теории Ньюкомба нетрудно найти скорость перемещения полюса эклиптики v. См. главу 1. Зафиксируем два произвольных различных момента времени t1 и t2. См. рис. 5.7. С помощью формулы (5.3.7) найдем значения stat(t1) и stat(t2). Изобразим прямую, по которой перемещается во времени полюс эклиптики. Отметим на ней точки t1 и t2. Выберем такой масштаб, чтобы расстояние между отмеченными точками равнялось v∣t2 - t1∣. Положение полюса эклиптики каталога PA определяется как точка пересечения двух окружностей с центрами в точках ti и радиусами stat(ti),i = 1, 2. Из рис. 5.7 ясно, как при этом определяются величины stat(t) и stat(t) для произвольного значения времени t. Необходимо лишь сказать, что прямая SS, от которой отсчитывается угол stat(t), пересекает траекторию движения полюса эклиптики под углом d(t). Угол этот также находится из теории Ньюкомба. Астрономический смысл прямой SS весьма нагляден. Это "спрямленная"часть большого круга звездной сферы, проходящего через полюс эклиптики P(t) эпохи t и перпендикулярного в точке P(t) другому большому кругу, также проходящему через P(t) и точки равноденствия эпохи t.

Аналогично, для определения параметров stat(t) и stat(t) при всех t достаточно знать лишь два значения: stat(t1) и stat(t2).

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

5. 5 Статистические свойства оценок stat и stat

Сейчас мы подойдем к задаче оценки параметров и , задающих систематическую ошибку каталога, как к задаче статистической. Для этого предположим следующее. Пусть составитель каталога в момент времени tA совершил систематическую ошибку, задаваемую параметрами A и A. Пусть, кроме того, широта каждой измеренной им звезды подвергалась, – вследствие ошибки наблюдения, – случайному возмущению q, имеющему нулевое среднее, то есть Eqi = 0. Предполагается, что случайные погрешности qi, отвечающие различным звездам, независимы и имеют одно и то же распределение. Пусть s2 = Eq i2 – дисперсия случайной величины q i. Эта дисперсия нам, вообще говоря, неизвестна.

В этих предположениях широта i-той звезды в каталоге будет иметь вид
bi = Bi(tA) - gA sin(Li(tA) + fA) + qi.
(5. 17)
 

Со статистической точки зрения мы имеем выборку, состоящую из N реализаций случайных величин {bi}i=1N вида (5.5.1). По этой выборке требуется определить статистические оценки gˆ и fˆ параметров A и A, а также оценить величину s, представляющую собой среднеквадратичную ошибку наблюдения. Мы сразу ограничим задачу и будем изучать оценки fˆ = stat и ˆg = stat, получаемые методом наименьших квадратов. Эти оценки имеют вид (5.3.6), (5.3.7). Основное внимание будет уделено оценке величины A по причинам, объясненным в конце раздела 4.

Равенство (5.5.1) имеет вид, традиционный для регрессионного анализа. В самом деле, это равенство утверждает, что ошибка наблюдения bi = Bi(tA) -bi является случайной величиной со средним A sin(Li(tA) + A), зависящим от неизвестных параметров A и A, и дисперсией s2. Требуется оценить значения неизвестных параметров методом наименьших квадратов и установить статистические свойства полученных оценок. В такой постановке кривую Y (x) = A sin(x + A) обычно называют ЛИНИЕЙ РЕГРЕССИИ.

Определим величины и с помощью соотношений (5.3.6) и (5.3.7). Отклонения bi случайны по предположению. Поэтому и получаемые из соотношений (5.3.6) и (5.3.7) оценки stat и stat также являются случайными величинами. Изучим их статистические свойства и рассмотрим, как они связаны с истинными, но неизвестными нам, значениями A и A.

Подставим в приведенные выше формулы для sb и cb вместо bi разность A sin(Li(tA) + A) - qi и используем эту подстановку в формулах (5.3.6), (5.3.7). Получим следующие выражения для величин stat и stat.

1- sum N qi(s2cosLi(tA)- dsinLi(tA)) tan-fA-+--N---i=1gA(d2-s2c2)cosfA-------- tanfstat = -1 sum N qi(c2sinLi(tA)- dcosLi(tA)) (5. 18) 1 + N---i=1gA(d2-s2c2)cosfA------- 1- sum N gstat = gA - N---i=1qi(sin-Li(tA) +-tan-fA-cos-Li(tA))-. (5. 19) (s2 + 2dtan fA + c2 tan2 fA) cosfA


Введем величину

R = (gA(d2 - s2c2) cosfA) -1.


Тогда (5.5.2) можно записать в виде
R- sum N tan f = tan-fA--+ sum N--i=1qi(s2cosLi(tA)---d-sin-Li(tA)), stat 1 + NR Ni=1 qi(c2sinLi(tA) - d cosLi(tA))
(5. 20)
а (5.5.3) – в виде

1 sum N -N---i=1-qi(sinLi(tA)-+-tan-fA-cosLi(tA)) gstat = gA - (s2 + 2d tan fA + c2tan2 fA) cosfA .


Из условия Eqi = 0 находим, что полученная оценка параметра stat является несмещенной, то есть
Egstat = gA.
(5. 21)
 

Дисперсия же оценки stat, обозначаемая Dg, имеет вид
2 Dg = ---------------------s---------------------. N (cos2 fAs2 + 2d cosfA sinfA + c2sin2fA)
(5. 22)
 

Если ошибки наблюдения qi нормально распределены, то и величина stat является нормально распределенной, и первые два момента (5.5.5) и (5.5.6) полностью определяют ее распределение. Этот факт даст нам возможность построить доверительный интервал для значения A.

Анализ оценки stat несколько более сложен. Воспользуемся следующим равенством, получаемым из (5.5.4):

tanf - tan f = stat R sum NA N---i=1-qi((s2-+-d-tanfA)-cos-Li(tA)--(d-+-c2-tan-fA)-sinLi(tA))- = 1 + R- sum Ni=1 qi(c2 sin Li(tA)- d cosLi(tA)) (5. 23,) N
и тем фактом, что при больших N второе слагаемое в знаменателе правой части (5.5.7) – величина малая. В самом деле, эта величина – случайная с нулевым средним и дисперсией

2 ---------s-c2---------. N g2A(s2c2- d2)cos2fA


Если qi нормально распределены, то и рассматриваемая величина также нормально распределена. Из этого факта для каталога Альмагеста следует, что уже для N = 30 вероятность pN того, что знаменатель правой части (5.5.7) будет отрицательным, не превосходит 5 ⋅ 10-3. С ростом N данная вероятность быстро убывает: p50 ≤ 2, 5 ⋅ 10-4, p 80 ≤ 4 ⋅ 10-6, p 100 ≤ 3 ⋅ 10-7, p 200 ≤ 8 ⋅ 10-13, p 300 ≤ 2, 5 ⋅ 10-18.

Из формулы (5.5.7) следует, что, вообще говоря, E tan stat⁄= tan A. Однако из этой формулы легко получить функцию распределения F(x) случайной величины tan stat - tan A, необходимую при нахождении доверительного интервала для A. В самом деле, если пренебречь тем маловероятным случаем, что знаменатель в (5.5.7) становится отрицательным, то из этой формулы получается выражение для F(x):

F (x) = P(tan fstat- tan fA < x) = P(jx < x),
где случайная величина jx имеет вид

R sum jx = --- qi((s2 + d(tan fA + x) cosLi(tA))- (d + c2(tan fA + x) sin Li(tA))). N i


Следовательно, если величины qi нормально распределены с дисперсией s2, то и величина j x имеет гауссовское распределение со средним, равным нулю, и дисперсией
2 2 D(j ) = R--s-(c s - d2)(s + 2d(x + tan f ) + c (x + tan f )2). x N 2 2 2 A 2 A
(5. 24)
 

Следовательно,
V ~ ----- F (x) = Φ(x/ D(jx)),
(5. 25)

где (x) = (2)-1/2 -∞x exp(-u2/2) du.

Найденные выше значения stat и stat являются, как говорят, ТОЧЕЧНЫМИ ОЦЕНКАМИ неизвестных параметров A и A. Поскольку найдены функции распределения этих оценок, то можно исследовать вопрос об их возможной погрешности. Дадим ответ на этот вопрос в стандартных терминах ДОВЕРИТЕЛЬНЫХ ИНТЕРВАЛОВ, основываясь на формулах (5.5.5), (5.5.6), (5.5.8), (5.5.9).

В математической статистике задача нахождения доверительного интервала порождена следующей ситуацией, которую поясним на примере оценки величины A. Эта величина является вполне определенной, детерминированной ошибкой, сделанной составителем каталога. В результате статистической оценки A, – в нашем случае по методу наименьших квадратов – получается СЛУЧАЙНАЯ величина stat. Возникает вопрос, какие границы можно указать для неизвестной нам величины A, если мы определили stat?

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

V ~ --- V~ --- P (∣g - g ∣ < x) = Φ( D x) - Φ(- D x). stat A g g


Определим величину x из уравнения

V~ --- V~ --- Φ( Dg xe) - Φ( - Dg xe) = 1- e,
или, что то же, из уравнения Φ(- V~ --- Dg x) = /2.

Тогда интервал
Ig(e) = (gstat- xe, gstat + xe)
(5. 26)
представляет собой доверительный интервал для A с уровнем доверия 1 -. Это следует из того, что P(stat - A∣≥ x) = .

При определении величины x мы, в частности, использовали значение Dg, которое зависит от неизвестных нам параметров s2 и A. Как это обычно делается в математической статистике, вместо s2 подставим в формулу для Dg сходящуюся к ней ОСТАТОЧНУЮ дисперсию. См. (5.3.3). Получим:

N 2 -1- sum 2 s 0(gstat,fstat,tA) = N (ΔBi(tA) - gstatsin(Li(tA) + fstat)) , i=1
а вместо A – величину stat. Момент tA составления каталога нам также неизвестен, поэтому все перечисленные выше вычисления необходимо проделать для всех моментов времени t с тем, чтобы оценить систематическую ошибку stat(t),stat(t) при условии, что каталог был составлен в произвольную фиксированную эпоху t.

Аналогичным образом можно найти доверительный интервал для A с уровнем доверия 1 - . Этот интервал If() будет таким:
( ) fstat- ------------ye------------,fstat + ------------ye------------ , 1 + tan2 fstat- ye tan fstat 1 + tan2 fstat + yetan fstat
где y – решение уравнения F(y) - F(-y) = 1 - , в котором функция распределения F задана равенством (5.5.9).

ЗАМЕЧАНИЕ. Определение параметров stat(t) и stat(t) нацелено на последующую компенсацию этих ошибок в эпоху t, рассматриваемую как возможное время составления каталога. При этом зависимость данных параметров от t обязана лишь вековому движению полюса эклиптики. Полученные выше оценки истинных значений ошибок и в каталоге, как функций априорной датировки, важны не только с чисто познавательной точки зрения, но и для косвенной проверки правильности предлагаемого подхода. Например, если бы в качестве stat была получена величина, в несколько раз превышающая точность каталога, это указывало бы на какие-то неучтенные нами существенные эффекты.

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

5. 6 Выводы

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

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

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

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

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

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

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

Рис. 5. 8: Страница из издания Альмагеста якобы 1551 года
Главная страница Оглавление Продолжение