Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.sai.msu.ru/neb/pcm/pcm05_51.pdf Дата изменения: Wed Sep 12 09:46:09 2007 Дата индексирования: Sun Apr 10 00:11:17 2016 Кодировка: Windows-1251

Н.В.Емельянов

ПРАКТИЧЕСКАЯ НЕБЕСНАЯ МЕХАНИКА
Оглавление. Глава 5. Построение моделей движения небесных тел на основе

наблюдений.

5.51. Метод дифференциального уточнения параметров движения небесных тел на основе наблюдений. Применение метода наименьших квадратов.
В небесной механике существует ряд задач, различных по изучаемым объектам, но сходных по способу их решения. Эти задачи можно сформулировать как "Уточнение параметров движения небесных тел из наблюдений". Детальное объяснение того, что мы понимаем под параметрами и наблюдениями дано в предыдущих разделах. Здесь опишем эти понятия кратко. Параметрами движения небесных тел называются величины, от которых зависит движение тел и которые по крайней мере на некотором этапе изучения считаются постоянными. Мы рассматриваем три типа параметров. Первый тип параметры, которые входят в уравнения движения. Они существуют еще до решения уравнений. Второй тип параметров появляется в процессе решения дифференциальных уравнений движения. Это либо произвольные постоянные в общем аналитическом решении уравнений, либо начальные условия численного интегрирования, то есть координаты и компоненты скорости тел в начальный момент времени. Параметры третьего типа входят в соотношения, связывающие результаты наблюдений и координаты небесного тела. Они не связаны с изучаемым объектом и называются параметрами условий наблюдений. Примерами параметров первого, второго и третьего типов могут служить гравитационный параметр небесного тела, элементы его орбиты и геоцентрические координаты обсерватории, соответственно. В процессе наблюдений измеряются какие-либо величины, зависящие от положения или скорости небесного тела. Они так и называются измеряемые величины. Наблюдения дают нам значения измеряемых величин на моменты измерений. Одновременно могут измеряться несколько величин. Для простоты изложения без нарушения общности задачи будем полагать, что все величины измеряются независимо, 1


каждая в один момент времени. Моменты измерений разных величин могут совпадать. Примерами измеряемых величин являются угловые топоцентрические экваториальные координаты небесного тела, топоцентрическая дальность небесного тела, разность угловых координат двух тел. Измеряемая величина всегда является реальной физической величиной, получаемой с помощью измерительных приборов в определенном месте. Моменты измерений отсчитываются по часам, расположенным в пункте наблюдений. Задачу уточнения модели движения небесного тела сформулируем следующим образом: даны результаты наблюдений, требуется найти параметры движения. Пусть одна из измеряемых величин, а p1 , p2 , ..., pn истинные, но неизвестные значения параметров движения небесного тела. Теория и модель движения дают нам модельное значение измеряемой величины c , как известную функцию времени t и параметров движения:

c = (t, p1 , p2 , ..., pn ).

(1)

Любая модель содержит ошибки. Обозначим ошибку модели через th . Тогда истинное значение измеряемой величины определится путем исключения ошибки th :

= (t, p1 , p2 , ..., pn ) - th .
В действительности измеряемые величины получаются из наблюдений и поэтому содержат ошибки наблюдений. Пусть o наблюденное значение измеряемой величины, а obs ее ошибка. Вычитая ошибку наблюдения, снова получим истинное значение измеряемой величины

= o - obs .
Приравнивая правые части последних равенств, получим

o = (t, p1 , p2 , ..., pn ) +

obs

- th .
m

После выполнения измерений в моменты времени t1 , t2 , ..., t лучим систему уравнений
o i = (ti , p1 , p2 , ..., pn ) + (i) obs

, по(2)

-

(i) th

(i = 1, 2, 3, ..., m)

относительно истинных значений искомых параметров. Величины



(i) sum

=

(i) obs

-

(i) th

2


называют суммарными ошибками теории и наблюдений. (i) При решении рассматриваемой задачи точные значения ошибок obs (i) и th (i = 1, 2, 3, ..., m) остаются неизвестными. Их обычно рассматривают как случайные величины с заданными вероятностными характеристиками (законами распределения, моментами и т. п.). При этом уравнения (2) заменяют так называемой системой условных уравнений
o i = (ti , p1 , p2 , ..., pn ) (i = 1, 2, 3, ..., m),

(3)

представляющей собой систему из m уравнений с n неизвестными p1 , p2 , ..., pn . Система условных уравнений (3) несовместна, она не имеет решения, поскольку получена путем вычитания из правых частей точ(i) ных уравнений (2) случайных независимых суммарных ошибок sum . Можно пытаться найти некоторую приближенную оценку искомых параметров. При этом получаемые значения должны по возможности мало отличаться от истинных. Алгоритм нахождения приближенной оценки называют алгоритмом фильтрации . Основной задачей этого алгоритма является возможное уменьшение (фильтрация) влияния ошибок теории и ошибок наблюдений. Выбор алгоритма фильтрации неоднозначен, его структура зависит от имеющихся сведений о сум(i) марной ошибке sum . На практике таких сведений очень мало или они вообще отсутствуют. Поэтому приходится довольствоваться теми или иными предположениями о свойствах суммарной ошибки и алгоритмом фильтрации, основанном на этих предположениях. Соотношения (3) можно рассматривать как уравнения относительно искомых параметров p1 , p2 , ..., pn . Решить эти уравнения непосредственно на практике не представляется возможным еще и по другой причине. Дело в том, что (ti , p1 , p2 , ..., pn ) сугубо нелинейная функция своих аргументов. Чаще всего ее бывает невозможно даже записать в явном виде. Тем более нельзя получить в явном виде решение уравнений (3). Решение задачи будем выполнять по схеме, которая уже представлена в Главе 1. Эту схему точнее следует назвать методом дифференциального уточнения параметров движения небесных тел из наблюдений. Практически на любом этапе исследований бывают известны некоторые приближенные значения искомых параметров. Назовем эти зна(0) (0) (0) чения предварительными и обозначим их через p1 , p2 , ..., pn . Пусть точные значения параметров p1 , p2 , ..., pn отличаются от пред3


варительных приближенных на величины поправок

p1 = p1 - p1 ,

(0)

p2 = p2 - p2 , ... , pn = pn - p(0) . n

(0)

Тогда (3) можно записать в виде
o i = (ti , p1 + p1 , p (0) (0) 2

+ p2 , ... , p(0) + pn ). n

(4)

Для большинства небесных тел модели движения развиваются уже давно. Поэтому на очередном этапе уточнения предварительные значения параметров уже близки к истинным. Это позволяет считать поправки p1 , p2 , ... , pn малыми и разложить правую часть соотношения (4) в ряд Тейлора по степеням поправок:
o i

= (ti , p

(0) 1

,p

(0) 2

n

, ... ,

p(0) n

)+
k=1

pk

pk + ... .
i

(5)

Производные в правых частях вычисляются при

t = ti , p 1 = p

(0) 1

, ... , pn = p(0) . n

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

i

c(0)

= (ti , p1 , p2 , ... , p a
(i) k

(0)

(0)

(0) n

),

(6) (7)

=

pk

,
i c(0) i

o i = i -

.

(8)

В результате получим
n

i =
k=1

ak pk ,

(i)

(i = 1, 2, ... , m).

(9)

Приближенные соотношения (9) называются условными уравнениями для определения поправок к параметрам. Они являются линейными неоднородными алгебраическими уравнениями относительно искомых поправок pk , (k = 1, 2, ... , n). Условные уравнения являются приближенными по двум причинам. Во-первых, в левых частях отброшены ошибки наблюдений и ошибки теории. Во-вторых, в правых частях отброшены все члены порядка 4


квадратов поправок и выше. Уравнения (9) иногда называют линеаризованными по отношению к условным уравнениям (3). Делая те или иные допущения относительно ошибок теории и ошибок наблюдений, можно выбрать один из разработанных алгоритмов фильтрации и найти приближенное решение условных уравнений (9). Существующие методы позволяют также оценить погрешность решения. После того, как поправки найдены, прибавляем их к предварительным значениям параметров и получаем новые и, как мы надеемся, более точные их значения. Такой метод определения параметров движения небесных тел называется дифференциальным уточнением. В силу приближенности условных уравнений и приближенности их решения новые значения параметров будут недостаточно точными. Однако уточнение можно провести повторно несколько раз. Если процесс уточнения сходится, то есть поправки от шага к шагу убывают, то вычисления можно прекратить, когда поправки станут существенно меньше их погрешностей. В этом случае мы получим значения параметров движения небесного тела p1 , p2 , ..., pn , соответствующие всем используемым при этом наблюдениям. Это соответствие однозначно определяется заданной моделью движения (1) и выбранным алгоритмом фильтрации. Теоретически сходимость дифференциального уточнения почти не исследована. Можно привести примеры, в которых процесс не сходится или сходится к ложным значениям искомых параметров. При использовании метода дифференциального уточнения следует учитывать, что с некоторого шага уточнения поправки начинают колебаться из-за неизбежных ошибок вычислений. После появления таких колебаний дальнейшие попытки уточнения параметров становятся бесполезными. С другой стороны, когда описанный процесс хорошо сходится, не возникает необходимости в точном вычислении производных pk , i так как в процессе уточнения они используются для определения все уменьшающихся поправок p1 , p2 , ... , pn . В этих условиях лежащие в разумных пределах погрешности вычисления указанных производных могут лишь несколько увеличить число шагов уточнения, практически не отражаясь на точности окончательного результата. На первом и последующих шагах уточнения после вычисления поправок p1 , p2 , ... , pn можно найти так называемые невязки услов-

5


ных уравнений
n

i = i -
k=1

ak pk ,

(i)

(i = 1, 2, ... , m).

(10)

После того, как процесс уточнения завершен, и поправки к параметрам стали пренебрежимо малыми, невязки условных уравнений станут окончательными невязками уточненной теории с наблюдениями ("OC") : o i = i - (ti , p1 , p2 , ..., pn ) (i = 1, 2, ... , m). (11) Совокупность невязок часто используют для оценки качества полученного решения. Однако для установления близости решения к истинному этого оказывается недостаточно. Из всех имеющихся алгоритмов фильтрации в практической небесной механике чаще всего применяется метод наименьших квадратов (МНК). Этот метод обладает рядом преимуществ по сравнению с другими алгоритмами фильтрации. Основное преимущество заключается в его простоте. Для краткости изложения метода введем матричные обозначения p1 1 1 p2 2 2 . . . p = . , = . , = . , . . . pn m m

A=

(1) a1 (2) a1

...

a1

(m)

a a ... (m a2

(1) 2 (2) 2 )

... a ... a ... ... (m) ... an

(1) n (2) n

,
sum

. = . . (m) sum



(1) sum (2) sum

.

Тогда систему условных уравнений (9) и их невязки (10) можно записать в виде = A p , (12)

= - A p .

(13)

6


Метод наименьших квадратов основан на выводах теории вероятностей. Он справедлив при соблюдении следующих условий [1]. 1. Задана модель движения (1). 2. Математическое ожидание E(sum ) равно нулю, т. е.

E(

sum

)=0.

3. Вектор ошибок sum является собственным в смысле теории вероятностей. 4. Ковариационная матрица D(sum ) задана с точностью до некоторого произвольного множителя, т. е.

D(

sum

) = 2 K .

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

S (p) = T K-1 = [ - A p]T K-1 [ - A p].

(14)

Заметим, что на практике проверка соблюдения данных условий оказывается невозможной. В частности, матрица K чаще всего неизвестна. Широкое распространение и привычность метода наименьших квадратов часто приводят к некритическому отношению к получаемым по этому методу результатам. Во многих случаях из результатов делаются неверные выводы. Причиной этого является обычно несоответствие реальных условий, в которых решается данная практическая задача, с условиями, принятыми при обосновании метода. Однако нередко МНК приводит к удовлетворительным результатам даже при невыполнении условий. На практике ковариационная матрица ошибок чаще всего оказывается точно неизвестной. В большинстве задач и механических моделей можно принять, что ошибки наблюдений некоррелированы. Тогда матрица K оказывается диагональной. Если имеется какая-то информация о точности одних наблюдений по отношению к точности других, то матрицу K можно сделать единичной. Для этого каждому наблюдению присваивается некоторый вес wi (i = 1, 2, ..., m). Каждое условное уравнение умножают почленно на установленный вес. При этом ошибки 7


также оказываются умноженными на этот вес. В результате такой операции наблюдения приводятся к равноточным, а матрица K становится единичной. В этом случае соотношение (14) принимает вид
m m

S (p) =
i=1

=
i=1

2 i

i - a pk

(i) k

2

.

(15)

Отыскание минимума функции S (p) сводится к решению системы уравнений S (p) = 0 (k = 1, 2, ..., n). (16) pk Как видно из (15), уравнения (16) содержат только нулевые и первые степени искомых поправок pk (k = 1, 2, ..., n). Поэтому система (16) оказывается системой линейных неоднородных уравнений. Эту систему принято называть системой нормальных уравнений. После выполнения дифференцирования в (16) систему нормальных уравнений можно записать в виде

L p = d ,
где L и d суть матрица и вектор

l1n l2n , ... lnn d=

l11 l L = 21 ... ln1



l12 l22 ... ln2

... ... ... ...

d1 d2 . . . dn

,

(17)

элементы которых вычисляются по формулам
m

lkj =
i=1 m

ak aj
(i)

(i) (i)

(k , j = 1, 2, ... , n),

(18)

dk =
i=1

ak i (k = 1, 2, ... , n).

(19)

Для дальнейших вычислений нам понадобится еще величина:
m

d0 =
i=1

(i )2 .

(20)

8


Рассмотрим матрицу коэффициентов нормальных уравнений L . Как видно из (18), матрица L симметричная, с положительными диагональными элементами. Найдем одним из известных способов матрицу L-1 , обратную матрице L. Тогда решение системы нормальных уравнений получим из матричного уравнения

p = L

-1

d.

(21)

Для того, чтобы существовало решение системы нормальных уравнений и его можно было найти, матрица L должна удовлетворять определенным условиям. В частности, ее ранг должен быть равен n . На практике иногда оказывается, что вычислить обратную матрицу L-1 можно с весьма ограниченной точностью. Если поправки все же найдены, можно определить среднеквадратические ошибки поправок следующим образом. Сначала вычислим так называемую среднеквадратическую ошибку на единицу веса 0 по формуле 1 2 [d0 - (dp)] , 0 = (22) m-n где p) найденный вектор поправок. Умножим теперь все элементы 2 матрицы L-1 на 0 . Полученная таким образом матрица
2 D = 0 L -1

(23)

называется ковариационной матрицей Запишем ее в виде D11 D12 D D22 D = 21 ... ... Dn1 Dn2

или матрицей ошибок поправок.

... ... ... ...

D 1n D 2n . ... Dnn

Согласно методу наименьших квадратов любой диагональный элемент матрицы D с номером k равен квадрату среднеквадратической ошибки k поправки pk , то есть

k =

Dkk .

(24)

Из (22), (23), (24) видно, что ошибки поправок уточняемых параметров убывают с ростом числа наблюдений m. Приближенно они пропорциональны 1 , так как на практике число наблюдений намного больше m числа уточняемых параметров. 9


Cледствием метода наименьших квадратов нулю среднеквадратической ошибки решения с блюдений. Сделанные здесь выводы справедливы при ях, налагаемых не только на ошибки теории и на принятую модель движения небесного тела. рассмотрены в книге [1]. Величины rkj , определенные соотношением

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

rkj =

Dk j , k j

называются коэффициентами корреляции между ошибками поправок, а матрица 1 r12 r13 ... r1n r 1 r23 ... r2n R = 21 ... ... ... ... ... rn1 rn2 rn3 ... 1 называется корреляционной матрицей. После нескольких шагов уточнения, когда поправки становятся достаточно малыми, ошибки поправок характеризуют ошибки улучшенных значений параметров, обусловленные как ошибками теории, так и ошибками наблюдений. Качество согласования теории с наблюдениями после уточнения параметров движения небесного тела будет характеризоваться величиной

d0 , m которая после успешного уточнения параметров является среднеквадратической величиной невязок i (11). Применение МНК на практике часто приводит к неожиданным проблемам. Дело в том, в условиях исследований на основе реальных наблюдений небесных тел с использованием ограниченных по точности теорий не всегда строго выполняются те допущения, при которых правомерно применение МНК. Если ошибки теории превалируют над ошибками наблюдений, то суммарная ошибка не будет случайной величиной. Это приведет к тому, что с увеличением числа наблюдений точность результата не будет улучшаться, а решение в свою очередь приобретет зависимость от состава измерений, то есть от того, в какие моменты делались измерения. Наличие такой зависимости делает результат не вполне достоверным. =
10


В конкретных задачах по уточнению параметров движения небесных тел часто оказывается, что определитель матрицы L близок к нулю. В этих случаях мы имеем дело с так называемыми плохо обусловленными системами нормальных уравнений. При решении таких систем поправки к параметрам могут получаться столь грубыми, что процесс уточнения не будет сходиться. Причина плохой обусловленности может заключаться не в самом методе наименьших квадратов, а свойствах механической модели. Примером случая с плохой обусловленностью является процесс уточнения долготы узла кеплеровской орбиты небесного тела при очень малом ее наклоне. Другой пример совместное уточнение долготы перицентра орбиты и средней аномалии в эпоху при малых эксцентриситетах орбиты. Показателем плохой обусловленности может служить близость к единице одного или нескольких коэффициентов корреляции. Чтобы уменьшить плохую обусловленность, можно исключить из списка уточняемых параметров тот из них, который дает сильную корреляцию, фиксируя его предварительное значение. Более подробно и вместе с обоснованием МНК изложен в книге [1]. Метод дифференциального уточнения параметров движения небесных тел из наблюдений с вопросами его практического применения рассмотрен в книге [2]. Упрощенное описание МНК дано в книге [3]. Существенной частью алгоритма уточнения параметров движения небесных тел из наблюдений является вычисление значений измеряемой величины на заданные моменты времени, а также ее производных по улучшаемым параметрам. Для этой цели можно применять как формулы аналитической теории движения небесного тела, так и методы численного интегрирования уравнений их движения. При использовании аналитической теории время вычислений будет пропорционально количеству используемых наблюдений. При численном интегрировании затраты времени на вычисления пропорциональны интервалу времени, на котором велись наблюдения. Порядок таких вычислений объяснен в следующих параграфах.

11


Литература
1. Эльясберг П.Е. Определение движения по результатам измерений. М.: Наука, 1976. 2. Емельянов Н.В. Методы составления алгоритмов и программ в задачах небесной механики. М.: Наука, 1983. 3. Щиголев Б.М. Математическая обработка наблюдений. 3-е изд. М.: Наука, 1969.

12