/
Текст
КОСМИЧЕСКИЕ ИССЛЕДОВАНИЯ
1 9 в 3
В ы п. 1
УДК 629.197.7
Э. Л. Аким, Т. М. Энеев
ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ ДВИЖЕНИЯ КОСМИЧЕСКОГО
ЛЕТАТЕЛЬНОГО АППАРАТА ПО ДАННЫМ
ТРАЕКТОРНЫХ ИЗМЕРЕНИЙ
Приводится методика определения параметров движения космиче-
ского летатольцого аппарата по дапным траекторных измерений. При-
водится также алгоритм оценки точности прогнозирования движения
аппарата по результатам обработки траекторных измерений. Методика
предполагает ее использование на универсальных быстродействующих
вычислительных машинах.
ВВЕДЕНИЕ
В настоящей статье излагается методика определения орбит косми-
ческих летательных аппаратов по данным траекторных измерений. При
этом под траекторными измерениями понимаются измерения относитель-
ных координат и компонент скорости аппарата и небесного тела. Такие
измерения могут производиться с наземных измерительных пунктов или
с борта космического аппарата. Предполагается, что измеряемые коорди-
наты и скорости могут быть произвольного вида и состава.
Проблема определения орбит небесных тел по данным наблюдений в
течение длительного времени разрабатывалась астрономией и прежде
всего таким важным ее разделом, как небесная механика. Однако до не-
давнего времени в центре внимания находились лишь задачи определения
орбит естественных небесных тел.
С развитием техники космических полетов назрела острая необхо-
димость в создании эффективных методов определения орбит искусствен-
ных небесных тел. Развитые ранее в небесной механике методы требовали
известной доработки, учитывающей специфику траекторий космических
аппаратов, возможности современных средств наблюдения и вычисли-
тельной техники. К настоящему времени в литературе описан ряд методов
определения орбит космических летательных аппаратов [1—3].
Вместе с тем развитие техники космических полетов продолжает предъ-
являть все новые, более высокие требования к эффективности методов оп-
ределения орбит, в первую очередь к такой их важной характеристике,
как быстродействие.
При разработке излагаемого метода главное внимание было уделено
максимальному повышению его быстродействия при сохранении'необ-
ходимой точности расчетов. Это было достигнуто выбором сравнительно
простых алгоритмов для стандартных расчетов, которые многократно
проводятся при массовой обработке измерений. Полученные алгоритмы
и расчетные схемы учитывают специфические особенности силового поля
и динамики движения летательного аппарата в этом поле. К их числу сле-
дует отнести схему разбиения траектории на участки маловозмущенного
движения, выбор универсальных параметров орбиты, пригодных для
любых типов орбит, модель и конечные формулы для вычисления частных
производных и т. д.
6
Э. Л. Аким, Т. М. Энеев
Серьезное внимание было обращено также на повышение надежности
предлагаемого метода, обеспечивающей определение траектории в широ-
ком классе случаев (невысокая точность измерений, недостаточно полный
состав измерений, плохое нулевое приближение для определяемых пара-
метров и т. д.). С этой целью наряду с обобщенным методом Ньютона при-
влечен метод наискорейшего спуска, получена схемная возможность под-
ключения к результатам траекторных измерений дополнительной апри-
орной информации о траектории.
Управление полетом космического аппарата требует надежной оценки
точности прогнозирования его движения по результатам обработки тра-
екторных измерений. Важный круг вопросов, связанных с оценкой точ-
ности прогнозирования, рассматривается в работе. В частности, приво-
дится алгоритм оценки точности прогнозирования для излагаемого метода
обработки траекторных измерений. Алгоритм может быть использован
также в процессе баллистического проектирования полета для выбора
эффективного состава п программы траекторных измерений.
Предлагаемый метод определения параметров движения летательного
аппарата по данным траекторпых измерений предполагает его исполь-
зование на быстродействующей универсальной вычислительной машине.
Метод является достаточно общим, однако в первую очередь он предназ-
начен для решения задач по определению орбит космических аппаратов,
летящих к Луне или к планетам Солнечной системы,
В заключение авторы пользуются возможностью поблагодарить про-
фессора П; Е. Эльясберга за ряд полезных замечаний, высказанных по
поводу настоящей статьи.
1. ПОСТАНОВКА ЗАДАЧИ
Рассматривается движение космического аппарата на пассивном уча-
стке траектории.
В прямоугольной геоцентрической системе координат XYZ с фикси-
рованными относительно звезд направлениями осей это движение полно-
стью характеризуется в каждый момент времени t радиусом-вектором
летательного аппарата
r(x, y,z) (1.1)
и его вектором скорости
V{rx, vu,vz}. (1.2)
Пусть
.....Qm (1.3)
— независимые параметры, однозначно определяющие траекторию дви-
жения Тогда
r = r(Q,t), (1.4)
v = v(^,t), (1.5)
где Q — вектор с компонентами (1.3).
1 В простейшем случае неизвестными параметрами (1.3) движения'являются шесть
независимых величии, однозначно фиксирующих траекторию полета' В более сложных
случаях к числу определяемых параметров относят также отдельные постоянные,
если точность их знания недостаточна для решения задачи и может быть существенно
повышена траекторными измерениями. Примером таких постоянных являются пара-
Определение параметров движения космического летательного аппарата 7
В момент времени tn производится измерение некоторой функции ¥
параметров (1.3) траектории:
Tn = 'F{r(<?,0; (1.6)
со среднеквадратичной ошибкой стп.
Измеренное значение функции обозначим через ('Fn)Ha6.
Совокупность значений (Ч^нав (л = 1, 2, N) различных функций
Y, измерения которых проводятся в процессе полета аппарата, образует
множество траекторных измерений. Набор различных функций W опре-
деляет состав траекторных измерений. В случае однородное состава из-
мерений последний содержит измерения лишь-одной функции ¥ парамет-
ров траектории.
Измеряемыми функциями 'F могут быть, например, сферические, ко-
ординаты космического аппарата в геоцентрической системе координат
(склонение и прямое восхождение), или в топоцентрической системе ко-
ординат (дальность, азимут, угол места), или радиальная скорость аппара-
та относительно пункта измерения и т. п. При измерениях средствами
автономной навигации роль функций Y могут играть углы между направ-
лениями из летательного аппарата на небесные тела с известными коор-
динатами.
Задача обработки траекторных измерепий формулируется следующим
образом: определить неизвестные параметры
Ql, Qi, . . -,Qm, (1.7)
характеризующие движение космического аппарата по информации о его
траектории в виде измерений
('Fi)
наб ’ (Ч'1)иа0.('К*)наб U-8)
заданного состава и известной точности, проведенных в моменты времени
<1, <2, . . ., Z/V.
2. СТАТИСТИЧЕСКИЙ ПОДХОД К РЕШЕНИЮ ЗАДАЧИ
В качестве основной апостериорной информации о траектории мето-
дика определения параметров движения космического аппарата исполь-
зует траекторные измерения. Предполагается, что ошибки этих измерений
носят в основном случайный характер, подчиняясь многомерному нор-
мальному закону распределения с нулевым математическим ожиданием
и известной дисперсией.
Наряду с рассмотренной апостериорной информацией для определения
параметров движения привлекается часто и априорная информация о
траектории. Априорная информация содержит ожидаемые значения от-
метры атмосферы, коэффициенты разложений в ряд гравитационного потенциала
Земли, Луны, планет, скорость света, астрономическая единица и т. д.
Рассмотрим два варианта выбора параметров прямоугольные координаты (1.1)
и компоненты скорости (1.2) космического аппарата, отнесенные к фиксированному
моменту времени То'.
Qi х0, Qa = j/o, Qi = го, = vx>, Qs = »V1, Qe =
и элементы его орбиты qi, q^,..., qt 1
Qi — Qt — <h, Q> — я», Q* = 9*, Qt чъ, Q> ' ?«•
8
Э. Л. Аким, Т. М. Энеев
дельных функций параметров траектории с соответствующими им вероят-
ностными характеристиками возможных ошибок. Эта информация возни-
кает, например, в результате обработки измерений предыдущих участков
траектории и оценки разброса полученных результатов или из каких-
то других источников. Она может быть истолкована как совокупность
некоторых коррелированных измерений с известной для них априорной
статистической характеристикой в виде корреляционной матрицы ошибок.
Избыточность измерений (по сравнению с числом определяемых пара-
метров) и случайный характер их ошибок вызывают необходимость ста-
тистического подхода к решению поставленной задачи. Результатом та-
кого подхода явился выбор в качестве статистического метода одного из
самых эффективных (в смысле минимума дисперсии) методов определе-
ния неизвестных параметров — метода максимального правдоподо-
бия [5].
В общем случае вся привлекаемая информация о траектории распадается
на две группы: некоррелированных и коррелированных измерений. При-
меняя метод максимального правдоподобия, будем учитывать наличие
в информации этих двух групп измерений.
Метод максимального правдоподобия, будучи использован для об-
работки группы некоррелированных измерений, естественно, приводит
к простым формульным схемам метода наименьших квадратов, что. суще-
ственно ускоряет процесс обработки траекторных измерений.
Таким образом, лишь присутствие в информации о траектории группы
коррелированных измерений является основной причиной проведения
всех рассмотрений в рамках единой схемы метода максимального правдо-
подобия (а не метода наименьших квадратов).
3. АЛГОРИТМ РЕШЕНИЯ
1. Нелинейная система уравнений для определяемых параметров.
Метод максимального правдоподобия, применяемый для определения
параметров (1.3) траектории по измерениям
(Ч'Лнаб, (Ч^наб.....(^)паб , . • (ЧГм)иаб (3-1)
(под Тп и при п =/= I понимаются значения, вообще говоря, разных
измеряемых функций в разные моменты времени tn и ti), .приводит к необ-
ходимости минимизации известного функционала [5]
м м
ф = 2 Рп&пЬ = 2 Рш (%. - vn,ia6) (¥, - YiBa6), (3.2)
n,l—l n,l—l
где Pni — элементы «весовой» матрицы Р.
Если обозначить корреляционную матрицу измерений (3.1) через /Gy:
Kv = К?, ' (3;3)
то весовая матрица Р выражается через К соотношением:
Р = Ку\. (3.4)
В результате минимизации функционала (3.2) получим систему не-
линейных относительно параметров Q уравнений
^- = 0 U = l,2.........пг), (3.5)
Определение параметров движения космического летательного аппарата
9
или (ct учетом симметричности матриц К и Р)
М д1
2^n^-=0 (i = 1, 2..............т), (3.G)
П,1 = 1
которая может быть решена методом обобщенных касательных Ньютона
(при наличии достаточно хорошего нулевого приближения и высокой
точности измерений) или методом наискорейшего спуска.
2. Решение нелинейной системы уравнений методом Ньютона. Реше-
ние методом Ньютона сводится, как известно, к серии последовательных
приближений. Число приближений зависит от степени близости выбран-
ного нулевого приближения
<№, Q™, <?£’ (3.7)
к точному решению.
Пусть в результате (s — 1)-го приближения для параметров получены
значения тогда разлагаем в ряд Тейлора в окрестности этих зна-
чений параметров рассогласование £ на s-м приближении
m
В(5) = 6‘-1> + + • •.. (3.8)
где
д<?4(” = <Д” - Qi~l> (3.9)
— поправки на s-м приближении к величинам
Ограничиваясь в (3.8) линейными членами разложения и подставляя
полученное выражение в (3.6), приходим на s-м приближении к нормаль-
ной системе линейных уравнений относительно поправок Д@('!
4AQ(,) = -B, (3.10)
т. о.
2 ЛуД<?/> = -Bi (i = l,2........rn), ' (3.11)
J-i
где элементы матрицы А и компоненты вектора В имеют вид:
Л ЗЕ,
(3.12)
п,(=1 1 i
Bi = 2 Рщ-ЯГ (3.13)
П,<=1
Таким образом, на каждом приближении для определения поправок
к параметрам, полученным на предыдущем приближении, решается нор-
мальная система уравнений (3.11).
Матрица и правые части этой системы формируются для каждого при-
ближения при помощи (3.12) и (3.13) из производных и рассогласований
между наблюденными значениями измеряемых величин и их расчетными
значениями, причем производные и рассогласования считаются для траек-
10 Э. Л. Аким, Т. М. Энеев
----------------------------:-------------------------------!-----
тории, определяемой параметрами предыдущего приближения, для каж-
дого момента наблюдения и каждой наблюдаемой величины.
Система уравнений (3.10) может быть записала и в другом виде. С этой
целью введем прямоугольную матрицу (М X т) производных
(3.14)
с матричными элементами
ат,
(i = 1, . . М;/ = 1, • • rn). (3.15)
Тогда, обозначая через Г матрицу, транспонированную к I, получим
Л =r(v)"(f)' <316>
В-=Г(^-)РЕ, (3.17)
и уравнения (3.10) принимают вид
r (v)w (v)49'’’= - r (v)n- (318)
Для каждого приближения вычисляется величина среднеквадратич-
ной ошибки единицы веса
Где Af0o — общее количество измерений, участвующих в обработке на
данном приближении. Величина сто служит мерой разброса измерений
вокруг траектории, определяемой параметрами, полученными'в результате
предыдущего приближения.
Ввиду возможного присутствия в информации измерений с грубыми
«выбросами», нарушающими статистический характер информации, на
каждом приближении производится отбор измерений, подлежащих об?
работке. Отбор проводится с использованием траектории, полученной на
предыдущем приближении. В результате такого отбора в формировании
матрицы А и правых частей В системы нормальных уравнений (3.10) при-
нимают участие лишь те траекторные измерения, для которых выполняет-
ся неравенство
T-|Yn- Тпнаб|<С, (3.20)
°п
где 5 — константа отбора, значение которой на каждом приближении
определяется соотношением
t = voo. (3.21)
При обработке, информации, содержащей измерения со случайными
ошибками, полагаем v = 3. Наличие систематических ошибок измерений
приводит к необходимости увеличения задаваемого значения v.
Определение параметров движения космического летательного аппарата
11
Эти вычисления завершаются приближением с номером s, для которого
впервые выполняется неравенство
IWI<«i (i = 1,-2,..., zn); (3.22)
£i — набор заданных точностей.
Нормальное течение сходящегося процесса приближепий сопровож-
дается уменьшением среднеквадратичной ошибки ао.
3. Решение нелинейной системы уравнений методом наиекорейшего
спуска. При отсутствии достаточно хорошего нулевого приближения
(3.7) для определяемых параметров траектории или при обработке изме-
рений невысокой точности, когда метод Ньютона пе дает решения задачи,
может применяться метод наискорейшего спуска [Приложешщ I]. ’
Метод наискорейшего спуска приводит к системе дифференциальных
уравнений
А =--------В, , (3.23)
ds Уф
т. е.
m
= <‘ = '..........................................”•> <3-24’
где переменная s служит параметром вдоль кривой спуска в zn-мерпом
пространстве переменных
Qit Qu • • •» (3.25)
так что для каждого значения s
Qt=Qi(s) (i = l,..„ т). (3.2G)
Вычисление матрицы А, вектора В и функционала Ф в системе (3.23)
проводится для траектории с параметрами (3.26) по формулам (3.2),
(3.12), (3.13), используемым для формирования нормальной системы
<3.10) - (3.11).
Значения параметров (3.7) служат начальными данными-при интегри-
ровании системы (3.23)
Qi = Qz = .....Qm = <2^ при s = 0. (3.27)
Система (3.23) интегрируется до значения s, при котором функция
F = В\ (3.28)
становится меньше заданной малой величины.
Полученные при этом значения функций Qi являются искомыми па-
раметрами траектории, обеспечивающими минимум функционала Ф.
Численное интегрирование системы уравнений (3.23)—(3.24) проводит-
ся методом Эйлера.
Участок интегрирования при разумном выборе начальных данных (3.7)
оказывается небольшим.
Метод наиекорейшего спуска в статистической задаче обработки траек-
торных измерений серьезно повышает надежность математической ме-
12
Э. Л. Аким, Т. М. Энеев
тодикп обработки, но уступает по скорости решения задачи методу Пью
тона. Поэтому при решении задачи целесообразно совмещение обоих ме
тодов, в котором методу Ньютона, как более быстрому, отдается предпоч
тение. Метод же паискорейшего’спуска привлекается лишь в случае нару-
шения нормальной сходимости процесса при использовании метода Ньюто-
на, которое проявляется в росте среднеквадратичной ошибки gJ'”’ hi
s-m приближении по сравнению с ее значением на (s — 1)-м приблпжепш
cf‘s>> (3.29)
Последнее означает, что полученная в результате s-го приближения
траектория уводит нас от требуемого минимума функционала Фив этол
смысле является ошибочной. Используя полученную траекторию (s — 1)-го
приближения, привлекаем метод паискорейшего спуска для выявления
направления дальнейшего движения к минимуму. Вычисление парамет-
ров QW. траектории методом наискоройшего спуска можно производить,
например, по формуле:
Q<” = Qt’-1» 4- 0,5AQ(,)- (3.30)
Таким способом делаем один или несколько шагов при помощи метода
паискорейшего спуска и вновь пробуем продолжить процесс методом
Ньютона.
4. Основные особенности алгоритма, определяющие быстродействие
методики обработки. Из приведенных выше рассмотрений следует, что
на каждом приближении при использовании метода Ньютона (или па
каждом шаге при применении метода наискорейшего спуска) производит-
ся вычисление матрицы А и вектора В по формулам (3.12) и (3.13), содер-
жащим производные и разности между расчетными и наблюденными зна-
чениями измеряемых величин для каждого момента наблюдения и каждой
наблюдаемой величины. Поэтому время, необходимое для решения на
быстродействующей вычислительной машине задачи определения пара-
метров траектории по выполненным измерениям, слагается из следующих
компонент:
1) времени, необходимого для вычисления производных на все момен-
ты наблюдения в процессе каждого приближения;
2) времени, необходимого для получения расчетных значений изме-
ряемых величин также на все моменты наблюдения для каждого при-
ближения, и, наконец,
3) общего количества приближений (шагов), необходимого для полного
решения задачи.
Так как общее количество измерений, участвующих в обработке,
может быть велико, а участок траектории, с которого зти измерения сня-
ты, достаточно большим, то уменьшение общего времени обработки траек-
торных измерений является весьма серьезной задачей. Решение этой
задачи реализуется экономией времени для каждой из указанных выше
компонент.
В настоящем разделе кратко сформулирован алгоритм решения по-
ставленной задачи. Его основным стандартным звеном является прибли-
жение в методе Ньютона или шаг интегрирования в методе наискорейшего
спуска. Стандартный характер приводимых в них вычислений, не связан-
ных с номером приближения, позволяет в обозначениях переменных опу-
скать индекс приближения. Без потери общности рассматриваемых ал-
горитмов число определяемых параметров траектории (1.3) ограничим в
дальнейшем значением m = 6.
Определение параметров движения космического летательного аппарата 13
4. РАЗБИЕНИЕ ТРАЕКТОРИИ НА УЧАСТКИ И СХЕМА ФОРМИРОВАНИЯ
НОРМАЛЬНОЙ системы уравнении
Метод касательных Ньютона и метод наискорейшего спуска требуют
для решения системы уравнений (3.6) вычисления изохронных производ-
ных d^JdQ.
Расчет этих производных может быть проведен, например, численным
интегрированием соответствующей системы дифференциальных уравнений
движения в вариациях (36 уравнений первого порядка) или ее конечно-
разностного аналога [1 ].
Оба метода не требуют, однако, высокой относительной точности при вы-
числении производных, и это обстоятельство используется для упрощения
рабочих формул. Упрощение достигается тем, что при вычислении произ-
водных на участке маловозмущенного движения последнее заменяется
невозмущенным движением и производные находятся по конечным фор-
мулам задачи двух тел, простота которых позволяет существенно повы-
сить быстродействие методики.
1. Разбиение траектории на участки маловозмущепного движения.
Вся траектория полота условно разбивается на непересекающиеся по-
следовательные участки, в пределах которых движение происходит в
поле основной центральной силы и малых возмущений
В качестве параметра, упорядочивающего разбиение,
зовать время полета t, значениями которого
__ у(1) у»(2> у,(х)
проводится разбиение на участки (в качестве параметра
жет быть выбрана любая величина, однозначно связанная с временем).
Каждый участок определяется условием
у (*) < у(*+1)
других сил.
можно исноль-
(4-1)
разбиения мо-
(4.2)
и выделяет из множества всех траекторных измерений, подлежащих
обработке, группу измерений, относящихся к данному участку траек-
тории ।
(Ч\)иа0, (Тг)11а0, . . ., (Тм^цао, (4.3)
которая содержит некоррелированные и коррелированные измерения.
Таким образом, разбиение траектории на участки приводит к расна-
дению множества всех траекторных измерений на непересекающиеся
группы. В соответствии с распадением информации на групны матрица А
Л _ V П д1п д^\
2 1 nl dQt dQ- 2 | 2 Pnl dQ dQ- |
Z,n=i 3 k=i 31
и правые части В
М д£ х Iм" д* ]
в' = <2 д&' = 2 2 Р"‘ Й? ч (4-5)
п,1=1 Л=1 п,/=1 '
системы уравнений (3.10) — (3.11) могут быть разбиты на слагаемые
Л=2Л<4), (4.6)
4=1
В = 2 В(4)- (4.7)
*=i
14
Э. Л. Аким, Т. М. Энеев
каждое из которых соответствует своей группе (4.3) траекторных изм
рений н формируется аналогично (3.12), (3.13) в виде
л <*) _ VI n d%n ,> I
Aii ~ 2 PntdQ, dQj ’
n,l=l 1
М* ЭЕ
2 <4«
п.|=1
При этом предполагается, что измерения, принадлежащие различны:
участкам траектории, независимы между собой.
2. Модель для вычисления производных по параметрам участка
Для описания движения точки на каждом участке траектории наряд;
с геоцентрической экваториальной системой координат XYZ, ось X ко
торой направлена в точку весеннего равноденствия, введем систему ко
ординат xwYwZw. Начало О* этой системы совмещено с основным при
тягивающим телом участка, а оси соответственно параллельны осям си
стемы координат XYZ.
Положение точки в момент времени t характеризуется в системе ко*
ординат xwYwZw радиусом-вектором
(4.10)
а ее скорость — вектором V(*> г‘*>, Л В системе координат XYZ соответственно имеем г {х, у, Z), v {vx, vv, vt}. Если теперь радиус-вектор центра Земли и его скорость в ординат XWYWZW обозначить соответственно через г<*’ у<*>, и {<’. <>, <’), то г = г<*> _ г<*>, V = v‘*> - v<*> (4.И) (4.12) (4.13) системе ко- (4.14) (4-15) (4.16) (4.17)
и измеряемая функция V параметров траектории принимает ¥ = ¥ (г, V) = Т (rw, vw; r(B*’, v'*’]. вид: (4-18)
Так как все дальнейшие рассмотрения этого раздела относятся к од-
ному участку, то индекс к у всех величин опускается; при этом исполь-
зуется лишь система координат x^Y^Z^.
Пусть
91» 9г» • • *i Яв
(4.19)
Определение параметров движения космического летательного аппарата 15
— элементы конического сечения с фокусом в О, оскулирующего в момент
времени Т (начало участка), а
q'i, gi. . • gi (4.20>
— элементы конического сечения с фокусом в О, оскулирующего в момент
времени t.
Тогда
q' = ?' (?, Т\ 0, (4.21)
и
т = т (q', 0, (4.22)
v = v (д', 0, (4.23)
так что величины (4.19) однозначно определяют положение и скорость
точки в момент измерения t:
г = г {g, Т; 0, v — v {q, Т', 0- (4.24) (4.25)
В результате из (4.18) и (4.24) — (4.25) получаем
Ч' = Y {г {q, Т; 0, v {д, Т\ t}\ г8, v8). (4.26)
Дифференцируя (4.26) но параметру д, получим, пользуясь (4.21)-
(4.23):
дУ дУ dq, дУ dq' dy dq' (4.27)
— — — = 71 _| 2 . j __L.
дЧ dq dq2 dq dqf dq
Для вычисления производных
dfc _ dy dq dq (4.28)
пренебрегаем малыми возмущениями в пределах участка, полагая при-
ближенно вместо (4.27)
ат _ дУ
dq dq' ’
(4.29)
Таким образом, вычисление производных (4.28) по параметрам
(i = 1, 2, . . ., 6) начала участка заменяется в результате указанного
приближепия вычислением производных в тот же момент t по аналогич-
ным параметрам оскулирующего в этот .момент конического сечения,
которое и является «опорной кривой» для вычисления производных
в момент времени t.
В дальнейшем под обозначением
э9|*>
(4.30)
следует понимать производные (4.29).
3. Формирование матрицы и правых частей нормальной системы урав-
нений для отдельного участка траектории. Продолжая рассмотрения,
связанные с участком, отвлечемся временпо от общей задачи одновремен-
ной обработки информации многих участков. Предположим, что вся об-
16
Э. Л. Аким, Т. М. Энеев
рабатываемая информация принадлежит одному участку и определяются
параметры q этого участка. Выясним в этом простейшем случае механизм
формирования матрицы и правых частей соответствующей системы ли-
нейных уравнений.
Следуя (3.12) — (3.13), построим матрицу G и правые части И нор-
мальной системы уравнений (3.11)
с _ у р
n.TZi a?i
(4-31)
(4-32)
Среди измерений участка содержится в качестве основной группа
некоррелированных измерений <
(^1) наб» (^наб , • • (^)наб,
(4.33)
весовая матрица которой имеет
диагональный вид
п 1 '
Рп = — При
О при
п = I
(l,n = 1...............N),
(4.34)
и группа коррелированных измерений
('F М+1)на0» (М наО» • • •» (Ум) наб
(4.35)
с весовой симметрической матрицей С.
Весовая матрица Р участка имеет, следовательно, клеточно-диаго-
нальный вид. На диагонали весовой матрицы расположены две клетки
(4.36)
Верхняя клетка представляет собой диагональную матрицу весов
для группы некоррелированных измерений, нижняя клетка представляет
собой весовую матрицу С группы коррелированных измерений. Матрицу
G и правые части II можно разбить на два слагаемых:
G = U + R, (4.37)
и = V + S, (4.38)
где U и V соответствуют группе некоррелированных измерений:
р
n a<7j a<ij ’
Vi =
(4.39)
(4.40)
Определение параметров движения космического летательного аппарата 17
а Н и S — группе коррелированных измерений:
п V Г d^N+n d%N+l
(4.41)
M-N .
= n2ic’".'-^7 ^N+l-
(4.42)
Если в информации содержатся лишь независимые измерения, то си-
стема нормальных уравнений формируется- при помощи сортношений
(4.39)—(4.40), тождественно совпадающих с формулами метода наи-
меньших квадратов. При наличии в информации коррелированных из-
мерений обработка ведется методом максимального правдоподобия.
В случае когда группа коррелированных измерений состоит из изме-
рений самих определяемых параметров, для которых задана корреляцион-
ная матрица (С-1), формулы (4.41) и (4.42) принимают особенно простой вид
Пц = Сц> (4.43)
в
2 Cij{g.-(7j.)Ha6). (4.44)
j=i
Таким образом, если априори известны веса Рп независимых измерений
и корреляционная матрица С'1 зависимых измерений, то при определе-
нии параметров q участка по траекторным измерениям, проведенным на
участке, матрица G и правые части Н системы линейных уравнених фор-
мируются при помощи соотношений (4.37)—(4.42).
Предположим теперь, что в задаче обработки измерений определепного
участка траектории производится определение не параметров q данного
участка, а параметров Q какого-то другого участка, связанных с вели-
чинами q зависимостями
q=q{Q,T). (4.45)
Обозначим через I {q/Q} матрицу частных производных с элементами:
(“-f1"1’2......«>• <4«)
Тогда из (4.45) получаем выражение для производных по параметрам
() через производные по параметрам q
В результате приходим к выражению матрицы А и правых частей В
системы нормальных уравнений для величин Q через матрицу G и пра-
вые части Н для параметров q:
Л = ГШС/Н}- <4Л9)
В = /•{-£-} Н, (4.50)
2 Космические исследования, М 1
18
Э. Л. Л кил, Т. М. Энеев
где — матрица, транспонированная к матрице I. Все участвующие
в (4.49) — (4.50) матрицы и правые части .4, В, G, Н относятся, очевид-;
но, к рассматриваемому участку с индексом к.
Параметры Q сами принадлежат какому-то участку разбиения (4.1),
например участку с индексом к = 1, а.каждый предыдущий участок тра-
ектории однозначно определяет последующий участок, так что ।
</*’ = qW (дМ, ГС*-1’; TW} (4.51)
для всех значений к.
Поэтому матрицу проивводных (4.46) можно переписать в виде
/ „(А) \ /о^Х
<4-52>
В результате выражения (4.49) и (4.50) принимают вид
(4.53)
(4.54)
4. Формирование матрицы и правых частей нормальной системы урав-
нений при обработке всей совокупности траекторных измерений. За-
кончив рассмотрения, связанные с обработкой информации одпого уча-
стка, вернемся к общей задаче обработки информации многих участков.
Матрица А и правые части В системы уравнений (3.11) для опреде-
ления неизвестных параметров Q траектории по данным всей совокуп-
ности траекторных измерений были расчленены на слагаемые А^'1, В*<>
(4.6)—(4.7), каждое из которых соответствовало обработке измерений
только одного участка в параметрах Q и формировалось при помощи
(4.8)—(4.9) из измерений на этом участке. Однако обработке траекторных
измерений одного участка и был посвящен предыдущий раздел изложения,
причем для матрицы И<А) и вектора В<А) были получены выражения
(4.53) и (4.54) в предположении, что параметры Q принадлежат участку
с индексом к = 1,
Подставляя выражения (4.53) и (4.54) в (4.6) и (4.7), получаем окон-
чательные выражения для А и В, которые в матричной записи имеют вид:
Таким образом, получаем общую схему формирования матрицы А
и правых частей В нормальной системы уравнений (3.10) — (3.11).
В этой схеме стандартным элементом-ее конструкции являются матри-
ца Gw и вектор Il(fr) участка с индексом к (к = 1, 2, . . ., х), отнесенные
Определение параметров движения космического летательного аппарата 19
к параметрам (i = 1, 2, . . 6). Матрица G(i> и вектор Н(Л> форми-
руются при помощи соотношений (4.37) — (4.42) из производных d^/dqi^
и рассогласований £, вычисляемых для каждого момента измерений и каж-
дой измеряемой величины. Вычислению рассогласований £ по данным
расчета траектории с учетом необходимых возмущающих сил посвящен
раздел 6 этой работы. Вычисление производных д^/дд^, основанное на
модели, изложенной в п. 2 этого раздела, проводится в разделе 7.
Пересчет от параметров А-го участка к параметрам (А — 1)-го участ-
ка производится при помощи матрицы производных I один
раз на каждом участке для уже готовой матрицы G<*) и правых частей
Н**’, полученных по измерениям А-го участка. Матричные элементы матри-
цы I вычисляются по конечным формулам раздела 8.
Все указанные выше вычисления, связанные с общей схемой форми-
рования матрицы А (4.55) и правых частей В (4.56) нормальной системы
уравнений (3.10) — (3.11), производятся па каждом приближении s с ис-
пользованием траектории, определяемой параметрами полученны-
ми в результате предыдущего (s—1)-го приближения.
5. Особенности разбиения на участки траекторий дальнего космиче-
ского полета и орбит искусственных спутников. Разбиение траектории
на участки маловозмущенного движения может быть разумно проведено
в случае дальнего космического полета и при движении по орбите искус-
ственного спутника. Для траектории космического полета такое разбие-
ние возникает в результате деления траектории на участки, соответствую-
щие разным сферам действия. Например, траектория полета к Лупе рас-
падается на геоцентрический и селеноцентрический участки. Траектория
полета к планетам Солнечной системы разбивается па геоцентрический,
гелиоцентрический и планетоцентрический участки.
Траектория движения искусственного спутника Земли разбивается
точками одинаковых значений аргумента широты («синфазными точками»)
на витки. Маловозмущонпый участок траектории спутника может содер-
жать один или несколько его витков (в зависимости от величины и харак-
тера возмущений орбиты спутника). В пределах участка вычисление
производных по элементам орбиты, отнесенным к синфазной точке, осу-
ществляется по формулам задачи двух тел.
Заметим, что для широкого класса орбит искусственных спутников
при совместной обработке не слишком большого количества витков во-
обще не возникает необходимости в разбиении траектории на маловозму-
щенные участки. Однако для спутниковых орбит, подверженных большим
возмущениям (например, для близких орбит), в случае совместной обра-
ботки траекторных измерений, относящихся к большому количеству
витков, появляется необходимость в разбиении траектории на участки.
Это объясняется тем, что при переходе от параметров данного участка
к параметрам предыдущего участка при помощи матрицы производных
I {q^/qik~^} необходимо учесть накопившиеся за время движения па
участке наиболее существенные для вычисления производных возму-
щения орбиты спутника: возмущения, возникающие в результате нецеп-
тральности поля тяготения Земли, и возмущения, обусловленные со-
противлением атмосферы. Учет возмущений может быть произведен весь-
ма грубо, так как он необходим для повышения точности вычисления про-
изводных, а не самих орбит. С этой целью для вычисления матрицы про-:
изводных х>} используется опорная кривая, составленная из
двух оскулирующих конических сечений (см. раздел 8). Параметры кони-
ческих сечений и вычисляются в процессе расчета возмущенной
траектории движения. Так как пересчет от участка к участку произво-
20
Э. Л. Аким, Т. М. Энеев
дится один раз, то это практически не приводит к уменьшению быстро
действия методики обработки, повышая точность вычисления производных
Таким образом, схема формирования (4.55)—(4.56) системы нормаль
ных уравнений, соответствующая разбиению траектории на участки ма
ловозмущенного движения, может быть положена в основу методики об
работки траекторных измерений космических аппаратов или искусствен
ных спутников Земли.
5. РАСЧЕТ ТРАЕКТОРИИ ВОЗМУЩЕННОГО ДВИЖЕНИЯ
КОСМИЧЕСКОГО АППАРАТА
Определение положения и скорости космического аппарата, имепуе
мого в дальнейшем «точкой»:
х, у, z, vx, vu, vz (5.1]
в экваториальной системе координат XYZ, соответствующей рассматри-
ваемому участку полета, осуществляется для каждого момента времени |
при помощи расчета траектории возмущенного движения. Такой расче!
траектории производится численным интегрированием системы дифферен-
циальных уравнений движения с учетом всех необходимых возмущаю-
щих сил.
1. Дифференциальные уравнения движения в прцмоугольных ко-
ординатах. Уравнения движения в прямоугольной системе координат
имеют вид:
— =______L х 4- f
dt* - гз х 1-
= (5-2)
где р. — произведение постоянной тяготения на массу основного притя-
гивающего тела этого участка орбиты, a F {Fx, Fv, Fz} — вектор возму-
щающего ускорения.
Параметры траектории
<?!, <?2, • - - Qm (1-3)
однозначно определяют начальные данные
х = х0, у = у0, z = z0, vx = vx„ vv = vVa, vz = vz„ (5.3)
соответствующие моменту времени t = To:
x0 — x0 (Q, Tq),
.................................... (5.4)
vz, — Vz„ (Qi Fo).
2. Анализ возмущающих сил для случая межпланетного полета.]
Траекторию полета к планетам Солнечной системы можно разбить на три
участка, различающиеся составом возмущающих сил: геоцентрический,
гелиоцентрический и планетоцентрический. В качестве системы коорди-
нат XYZ на каждом участке выбирается соответственно геоцентрическая,
гелиоцентрическая или планетоцентрическая система координат.
На геоцентрическом участке траектории полета необходимо учиты-
вать возмущения, возникающие за счет пецентральности ноля тяготепия
Определение параметров движения космического летательного аппарата 21
Земли F(0(K), а также возмущепия движения точки Луной F„ и Солнцем
Fc. Поэтому вектор возмущающего ускорения принимает вид
F = F(C(H) + F„ + Fc. (5.5)
Ограничиваясь в разложении потенциала земного притяжения по сфе-
рическим функциям только членами первого порядка относительпо сжа-
тия, получаем
Fx(c>i<) = (-7-) {5 ’
= , (5.6)
Fr(0W) = (-у-) {5 ,
где
\ "©DKD /
R0KD — экваториальный радиус Земли, а — сжатие, gBKB— ускорение си-
лы тяготения па экваторе, Q — угловая скорость суточного вращения
Земли.
Ускорения, возникающие за счет возмущений от Луны и Солнца,
вычисляются по формулам:
р ( xj~ X х 1 ( У} —У У; |
|Г_Г.|3 — |77|Г) . fui - Hl||r_r.1S
I Zj — z z- )
= HHi -----is---i—Sri > (5.7)
1 Г’)|Г — Fjl3 I rj I J
где индекс / соответствует номеру возмущающего тела, р; — произведение
гравитационной постоянной на массу возмущающего тела, rj{xj, у,, Zj} —
радиус-вектор возмущающего тела в рассматриваемой системе коор-
динат.
Расчет гелиоцентрического участка траектории полета требует учета
возмущающего действия Лупы, планет, а также светового давления.
Возмущающее действие Земли и Луны при достаточном удалении точки
от Земли (на расстояниях порядка 2—3 млн. км) может быть заменено
с несущественной погрешностью возмущением Fa_a со стороны центра
инерции системы Земля — Лула, которому приписана масса, равная
сумме масс Земли и Лупы. Центр инерции системы Земля — Лупа и пла-
нета назначения Fu являются основными возмущающими телами. Состав
остальных возмущающих тел зависит от характера траектории полета.
Таким образом, вектор возмущающего ускорения в (5.2) принимает вид
_ F = Fa_;, + Fu + ^Fj + F„an; (5.8)
i
Fj — ускорение, возникающее в результате возмущающего действия
/-й планеты, влияние которой на выбранную траекторию полета оказы-
вается существенным. Вычисление возмущающих ускорений Fa_n, F4 и Fj
производится по формулам (5.7).
Так как плотность солнечного излучения и гравитационное действие
одинаково изменяются с расстоянием 1/гг, то учет возмущающего дей-
22
Э. Л. Аким, 'Г. М. Энеев
ствия светового давления может быть произведен формальной заменой
величиной
р = р0 (1 — хс), (5.1
где рс — произведение массы Солнца на гравитационную постояпну!
а х0 — постоянная, характеризующая величину светового давленщ
При расчете планетоцентрического участка траектории необходим
учитывать солнечные возмущения, а также некоторые дополнительны
возмущения F(A0U), обусловленные несферичностью планеты, нрптяж!
нпем ее спутников и т. д. Поэтому вектор возмущающего ускорения
(5.2) имеет вид
F = Fo + F„0Il, (5.10
где Fc вычисляется по формулам (5.7).
3. Дифференциальные уравнения в оскулирующих элементах. В не
которых случаях для расчета возмущенной траектории полета оказы
вается целесообразным интегрировать систему дифференциальных урав
пений в оскулирующих элементах орбиты с фокусом в центре основной
притягивающего тела.
Выбор элементов орбиты в качестве переменных должен обеснечиват)
сравнительную простоту вычисления прямоугольных координат и ком
попоит скорости (5.1) точки по элементам орбиты, а также простоту вы
числения правых чацтеп уравнений. В роли оскулирующих элементе!
орбиты могут использоваться постоянные первых интегралов уравне
ний движения задачи двух тел: компоненты вектора момента количеств
движения
С {Сх, Су, cz} (5.11.1
и вектора Лапласа
Г {/./„, А), (5.И.2J
определяющие геометрию орбиты. Дифференциальные уравнения для эти!
оскулирующих элементов имеют вид:
dc_
-=yFz- zFy,
d-£=zFx- xFz, (5.12.1)
de
dF = - УР"
dfx
— = (VyCz — vzCy) + (FyCz— Fzcv),
df
= {vzcx- Vxcz) + (Fzcx - F^), (5.12.2)
d/r
~dt~ — ( vxPy VyCx) -|- (FxCy FyCx) •
Для симметрии выписаны дифференциальные уравнения для всех
шести элементов (5.11.1) и (5.11.2). Однако в силу условия ортогональ-
ности
(cf) = 0 (5.13)
одно уравнение системы (5.12) может быть исключено (или могут интегри-
роваться все шесть уравнений системы, а условие (5.13) использоваться
для контроля вычислений).
Определение параметров движения космического летатеДного аппарата 23
Система уравнений (5.12) должна быть дополнена уравнением для кине-
матического элемента, в качество которого для гиперболических и эллип-
тических орбит можно использовать время прохождения перицентра ор-
биты т. Недостающее уравнение в этом случае имеет вид
£ = A (vF) — 5(rF), (5.12.3)
где
л = {3 — т) - -уг^2 + Нг) (™г)} ,
=-уг(с2 — К).
П’г = (rv), г =' | г / = | f с = | с |.
Начальные данные для этой системы дифференциальных уравнений
сх = Сх,, Су = Су,, cz = cZl, fx = fx,, fy = fVt, fz = Л,, т = x0, (5.14)
отнесенные к фиксированному моменту времени t = То, однозначно опре-
деляются аналогично (5.4) исходными параметрами (1.3) траектории,
так что
сх. = Cx.(Q, То),
.................................. (5-1.S)
Т0 = Т0 (Q> ^о)-
Для расчета орбиты полета к планетам,Солнечной системы дифферен-
циальные уравнения (5.12) интегрируются численными методами на
каждом из участков траектории.
На каждом участке оскулирующие орбиты имеют фокус, располо-
женный в начале соответствующей системы координат XYZ. Вектор воз-
мущающего ускорения F вычисляется для каждого участка по одной из
формул (5.5), (5.8) или (5.10).
Правые части системы дифференциальных уравнений (5.12) содержат
прямоугольные координаты и скорости точки (5.1). Для вычисления по-
следних на заданный момент времени по оскулирующим элементам орбиты
(5.11) в этот момент времени используются формулы пересчета, указан-
ные в Приложении II. Такой пересчет будем обозначать в дальнейшем
через
{с, f, т; t} —> {г, v; t}. (5.16)
4. Переход от данного участка траектории к следующему участку.
Переход в момент времени Т{к+1} от участка с индексом к к новому участ-
ку с индексом к + 1 связан с пересчетом координат и скоростей точки
в момент 7,(*+1) из старой системы координат X1*-, Y(k\ Z(k} в новую систе-
му координат Х(*+1), У(*+1), Z(*+1). Пересчет требует знания положения
старого центра О* в новой системе координат на момент пересчета или
иоложения нового центра О*+1 в старой системе координат и выполняется
по формулам:
г<*+1> = г(« _ г<*>+1 = r(*) + rg«)t
v<*+1> = v(« _ v{*)+i = v(*> + v<*«). (5Л 7)
24
Э. Л. Аким, Т. М. Энеев |
Полученные в результате пересчета координаты и скорости точь
д;(*+1) y(fc+l) 2<*+1) |>(к+1) 1>(*+1), 1>(*+1) (5.11
’ J ’ » х ’ V ’ х '
могут служить начальными данными (5.3) для интегрирования систем
(5.2) уравнений движения в прямоугольных координатах (к + 1)-г
участка.
Данные (5.18) позволяют получить элементы (5.11) оскулирующв
в момент 7,<*+1) орбиты с фокусом в Ок+1. Формулы для расчета этих ос
купирующих элементов, обозначаемые через
{г, v; t} -> {с, f, т; t}, (5.19
даны в Приложении III.
Оскулирующие элементы орбиты
^(4+1)^ * у(Аг+1)* jlk+l)* (5.20j
могут служить начальными данными (5.14) для системы дифференциаль-
ных уравнений (5.12) в элементах (к + 1)-го участка траектории.
5. Возмущения в движении искусственного спутника Земли. Расче,
траектории искусственного спутника Земли проводится в геоцентрической
экваториальпой системе координат XYZ с фиксированными относительно
звезд направлениями осей. Основные возмущения в движении искусствен-
ного спутника Земли обусловлены нецентральпостью поля тяготения
Земли и сопротивлением атмосферы (для спутников с низким перигеем
орбиты). Наряду с этим в ряде случаев (для спутников с высоким апогеем
орбиты) оказывается необходимым учет лунных и солнечных возмущений;;
Так что в общем случае вектор возмущающего ускорения искусственного
спутника Земли принимает вид
F = FOfK -f- FaTM Fji -f- Fo. (5.21)
Вычисление возмущающих ускорений FCH(, Fn, Fc может быть выпол-
нено по формулам (5.6), (5.7). К этому следует добавить лишь одно заме-
чание.
Использование приближенных формул (5.6) оказывается недостаточ-
ным в ряде задач. В этом случае в разложении потенциала земного при-
тяжения приходится учитывать наряду с гармоникой Ръо ряд последую-
щих гармоник, определяющих аномалии силы тяжести. Так как коэф-
фициенты при гармониках известны довольно грубо, то при высокоточных
траекторных измерениях отдельные коэффициенты должны быть введены
в число определяемых параметров (1.3).
Возмущающее ускорение, обусловленное торможением спутника в
атмосфере, вычисляется по формулам:
Fх (атм) = dp | ^отя | (,VX отн)>
у (атм) = ^P|voth I (^1/отн)> (5.22)
Fz (атм) = dp | V0TH | (l?z отп)>
где Vqth {vx отн> vu отд, v2 отв}—вектор скорости спутника в системе координат,
связанной с вращающейся Землей:
Vx отв = vx + Qy,
Vy отв == Vy — £2х, (5.23)
V2 отв — vt't
Определение параметров движения космического летательного аппарата 25-
d — баллистический коэффициент спутника, определяемый его весовыми,
геометрическими и аэродинамическими характеристиками.
Плотность атмосферы зависит от ряда факторов. Динамическая модель
атмосферы, содержащая учет этих многих факторов, еще не создана. По-
этому ограничиваются построениями приближенных статических моделей
атмосферы, учитывающих зависимость плотности лишь от высоты. При
этом плотность р атмосферы определяется как функция от высоты Л спут-
ника над поверхностью земного эллипсоида:
h = г - Лок„ [1 - асж (-J-)2] (5.24)
по приближенным формулам для принятой модели атмосферы.
Расчет траектории возмущенного движения спутника может произво-
диться интегрированием дифференциальных уравнений движения (5.2)
в прямоугольных координатах или интегрированием системы дифферен-
циальных уравнений (5.12) в оскулирующих элементах.
Переменные (5.11.1), (5.11.2) и соответствующие им уравнения
(5.12.1), (5.12.2) могут быть использованы для всего класса эллиптиче-
ских .орбит, включая круговые орбиты Ч Дифференциальное уравнение
(5.12.3) для кинематического параметра т, как и сам параметр, теряют
смысл для круговых орбит. Выбор вместо т в качестве кинематического
параметра времени прохождения заданного аргумента, широты ио
делает совокупность переменных
с*> cz, fx, fy, jг, (5.25)
и соответствующую ому спстему дифференциальных уравнений вполне
универсальными в смысле их пригодности для орбит любого типа.
6. Выбор переменных для системы дифференциальных уравнении дви-
жения. Мы рассмотрели два возможных метода расчета траектории воз-
мущенного движения: численное интегрирование системы дифференциаль-
ных уравнений в прямоугольных координатах и в оскулирующих эле-
ментах орбиты. Использование оскулирующих элементов в качестве пе-
ременных интегрирования на маловозмущенном участке траектории при-
водит к увеличению шага интегрирования и, следовательно, к уменьше-
нию времени расчетов. (Например, интегрирование уравнений движения
Земли в элементах орбиты приводит к более чем трехкратному увеличению
шага интегрирования по сравнению с его значением при интегрировании
эквивалентной точности в прямоугольных координатах.) Однако слож-
ность вычисления правых частей системы дифференциальных уравнений
в оскулирующих элементах возрастает, что снижает указанный выше
выигрыш во времени. Действительно, в этом случае необходимо, кроме
вычислений, требуемых для расчета правых частей системы уравнений
в прямоугольных координатах (т операций), произвести вычисление на
заданное время прямоугольных координат (/ операций) по оскулирующим
элементам. Обозначая через п количество шагов на данном участке траек-
тории при численном интегрировании системы уравнений (5.2), а через
пв— количество шагов при интегрировании системы уравнений в оскули-
рущщих элементах орбиты, получаем примерное отпошенпе времен рас-
1 Переход к переменным (аналогичным /х, fy, /2)
h = esin (<i> + 62), I = ecos (ш -f- ft)
используется в небесной механике [4] для получения системы дифференциальных урав-
нений в оскулирующих элементах, онисывающей весь класс эллиптических орбит
(включая круговые орбиты).
26
Э. Л. Аким, Т. Л/. Энеев
четов для указанных двух методов
лл (m + Z) Па / t
----------— ----I J —
Н1П п \ ' III
Проведенное сравнение двух методов позволяет попять относите:
ную ценность каждого из них в расчетах. Ускорение расчета при интегр
ровании системы дифференциальных уравнений в оскулирующих э:
ментах (по сравнению с интегрированием в прямоугольных координата
возникает лишь нри выполнении условия
Например, для класса орбит с малыми эксцентриситетами интегр
рование дифференциальных уравнений в оскулирующих элементах ок
зывается более целесообразным, нежели интегрирование в прямоугольна
координатах.: В этом случае решение уравнения Конлера может быть а
меноно непосредственным вычислением прямоугольных координат пр
помощи разложения последних в быстросходящнеся ряды по степенл
эксцентриситета (количество операций I становится сравнительно малым
Если условие (5.27) выполняется, то сравнительный выигрыш в.о времен
расчета при интегрировании в элементах орбиты растет с увеличение
сложности вычисления возмущающих сил (количество операций m ра
тет).
Для орбит, подверженных сильным возмущепиям, интегрирован]
в прямоугольных координатах оказывается по времени расчета боле
экономным.
6. ВЫЧИСЛЕНИЕ РАСЧЕТНЫХ ЗНАЧЕНИИ
ИЗМЕРЯЕМЫХ ФУНКЦИИ
Для вычисления рассогласования
£п = Тп-Ч'ппаб (6.!
в момент измерения tn, участвующего в формировании правых чаете
(4.32) системы уравнений, необходимо знать расчетное значение изме
ряемой функции Y параметров траектории на это время. Вычислепи
расчетного значения производится по данным расчета возмущенно!
движения.
1. Определение положения и скорости космического аппарата отш
сительно небесных тел. Измеряемая функция Т параметров траектори
зависит от иоложения и скорости точки в момент измерения относится]
ио пебесиых тэл. Точный расчет движения должен позволять определи:
в каждый момент времени в системе координат xwYwZw наряду с рад]
усом-вектором и скоростью точки
[г<*> {lW, 2(*>), ,
vW, v«};
радиус-вектор и скорость каждого из интересующих пас тел т
у V Т » J у » у > л •
Определение параметров движения космического летательного аппарата 27
л, следовательно, и их относительные положения и скорости
(р(А) — г^А)) {х^ — уУЧ — — z^} t
(у(А) — у(А)) {vt*} — 1)1® — — V^}. ( )
Ниже в качество примера траекторпых измерений рассматриваются
измерения положения и скорости точки относительно Земли. Только
поэтому в работе отдается предпочтение геоцентрической системе коор-
динат, в которой получаем (6.4) в виде
г = г(*> — jW, V = у(*> — у (А). (6.5)
(Если начало системы координат участка Х(*\ У<л>, Z(fc) выбрано в центре
-Земли, то, очевидно, г = г(*>, v = v(*>).
2. Вычисление положения наземного измерительного пункта в инер-
циальной геоцентрической системе координат. Измерения могут про-
изводиться как относительно центра Земли, так и относительйо измери-
тельного пункта, находящегося на поверхности Земли. В последнем слу-
чае предполагается, что заданы географические координаты измеритель-
ного пункта
ф, X, Н (6.6)
(ф — географическая широта, X — долгота, Н — высота над поверх-
ностью земного эллипсоида). При помощи (6.6) находим:
гии cos Ф, гип sin Ф; sin ф, cos ф; sin X, cos X, (6.7)
где гнп — расстояние измерительного пункта от центра земного эллип-
соида, Ф — геоцентрическая широта измерительного пункта.
Постоянные вычисляются при известных значениях величин (6.6)
по формулам
Пш C°S Ф = (, г "Т- , ,/. + (6.8.1)
\|созг<р4~а2з1п’ф|'• ) ।
г“" sin ф = (гг. V/. + <6-8-2)
\ [cos’<р +a* sin’ф] '• у |
тде 7?экв — экваториальный радиус, а = 1 — аож, аС1К — сжатие.
Координаты измерительного пункта в жестко связанной с Землей
прямоугольной геоцентрической системе координат XrYrZr, ось Zr кото-
рой направлена по оси вращения Земли, а плоскость XfZr совпадаете
.плоскостью гринвичского меридиана, имеют вид:
Г (^яп)г = (Гип cos Ф) COS X; (znn)r = (Пш sin Ф)',
(у ) = (г cos Ф) sin X; г {(х ) , (у ) , (z ) }.•
'^ип'г ' ии ' ’ ип ии'г’ '’иц'г’ ' ип'.г' ।
С каждым измерительным пунктом связывается топоцентрнческая
система координат, характеризующаяся следующими взаимно ортого-
нальными ортами:
ортом, касательным к меридиану (с направлением па север)
с<ип> {(б ц)г, (б„)г, (б13)г),
ортом, касательным к параллели (с направлением па восток)
е*ии) {(би)г, (6И)Г, (аи)г), (6.10.1)
28
Э.Л. Аким, Т. М. 9nei
ортом, нормальным к поверхности земного эллипсоида
е(ип) = [е(пп)хе(цИ)]( е<ип> {(б 31)г, (б32)г, (б33)г).
Через (6i;)r (i, 7=1,2,3) обозначены координаты этих вектор
в системе координат ХгУг2г.
Движение нолюса мира среди звезд вызывает смещение небесно)
экватора и точки весеннего равноденствия. Это перемещение обусловив;
участием оси Земли в двух движениях: вековом (прецессия) и период
ческом (нутация). Пусть система координат XYZ с фиксированными а
правлениями осей отнесена к среднему (истинное положение за вычето
путации) равноденствию и среднему экватору некоторой «начально!
эпохи. Система координат XrYrZr в текущую эпоху связана с положение
истинного полюса мира в эту эпоху. Поэтому переход от вращающей»;
системы координат XTYTZr текущей эпохи к системе координат ХУ,
«начальной» эпохи требует учета наряду с собственным вращением Земл
прецессионного и нутационного движения полюса мира за соответству
ющий интервал времени. Учет может быть произведен, например, дл
момента средней гринвичской полночи заданных суток тремя Bt
личинами:
Wo, g0, Go, (6.11.1
которые обычно имеются в астрономических ежегодниках, где 4% — уго:
между плоскостью гринвичского меридиана в гринвичскую полночь 1
направлением па точку весеннего равноденствия «начальной» эпох1
(т. е. направлением оси X системы координат XYZ), g0, Go — редукциоа
ные величины тригонометрической системы. Отсюда следует, что в снсте
ме координат XYZ координаты и компоненты скорости измерительного
пункта
ГИП {гИП, Упп» Znn), Vlin{VZan> VZBD)t (6.9.2)
а также координаты ортов топоцоптрической системы
е<ип) {б ц, б 12, б 13), е<ип) {б21, б22, б23), е<»“> {б31, б32, б33) (6.10.2)
изменяются с точением времени. Для вычисления их значений в момент
времени t (гринвичское время, отсчитываемое от средней полночи, к ко-
торой отнесены редукционные величины (6.11.1)) производится расчет
матрицы с учетом постоянных (6.11.1) по формулам:
аы — cos Ф cos v’> ai2 — — sin Ф cos ai3 = sin v;
a2L = яшф cos e0; a22 = созф cos e0; = sin e0; (6.11.2)
®31 = ^12^23 ^13^221 ®32 = ^13^21 ^11^231 ®33 = ^11®22 — ®12®21
(формулы не содержат нутационных членов, соответствующих интервалу
времени t), где ф = ф0 + Ш, v = v0 + ut,
Vq = gg cos Go, Bq = —~ gg sin Go.
Вычисление величин (6.9.2) при помощи полученной матрицы про-
изводится но формулам:
Хип — аи (‘Еии)г Н* а12 (Уип)г “I" «13 (2ип)г>
Уип = «21 (*нп)г + (уип)г + а23 (хип)г, (6.11.3)
Хип = a3i (£ип)г + «32 (уип)г+, «33 (zan)r>
Определение параметров движения космического летателфого аппарата 29
Гхш1 — ^а12 (x«n)r а11 (УипМ
гуиц = ^а22 (xuii)r а21 (УИиМ
%п = [азз (*ип)г — а31 (Уиц)г1 П-
(6.11.4)
3. Вычисление расчетных значений и градиентов для различных
измеряемых функций. Ниже в качестве примеров приводятся формулы
для вычисления значений некоторых конкретных измеряемых функций
параметров траектории и векторов-градиептов
L {До Ly, Lt} и M {Mx , My, Mt}\ (6Л2)
Lx = dW dje , 54<* ’ “ “ , d'V L'=-d7' (6.13)
мх= ay &vx Л/Г , My = X— , V dvv az Qvz (6-14)
Состав наземных измерений, для которого выписываются эти формулы,
содержит наиболее употребительные измеряемые функции [1, 3]
Y = Z), ДЛ,т, (6.15)
где дальность/) точки от измерительного пункта, азимут А и угол места у
точки представляют собой сферические координаты точки в топоцентри-
ческой системе координат измерительного пункта, a D — радиальная
скорость точки относительно пункта.
Заметим, что расширение приводимого состава измерений путем до-
бавления новых измеряемых функций Y параметров траектории требует
лишь вывода формул для этих функций, аналогичных приводимым пиже,
без каких бы то ни было изменений в самой методике расчетов.
Переходим к интересующим нас формулам
1. Т = D —дальность точки от измерительного пункта:
D = /(Г — Гцп) (Г — Гцп) = У(х — хИ1,)2 + (у — УипГ + (z — zlin)2, (6.16.1)
Ь = (Л-Янае), (6.16.2)
L = “д' (г гии)1
т. е.
!>» — ~jy з?ип)1
Д = 4-^-2/ип)> (6-16.3)
Lz = ~Q- (z ZHII).
2. XV — D — радиальная скорость точки относительно измеритель-
ного пункта:
В = "д’ (v vBn) (г — Гип) =
= 4" {(v« — %u) (* — *ип) + (Vy — VVhii) (у — уиц) -Ь (г2 — VIhu) (z — Z„n)),
(6.17.1)
£»=(£>-Z)Ba6), (6.17.2)
М = ^-(г — гШ1),
30
Э. Л. Лоди,
Т.
М. Энеев,'
т. е.
Мх — (х ^иа)>
= -^(у- Уип),
Мг = ~(z — z„„),
(6.17,
т. е.
(6.17.!
3. ¥ = A — азимут точки в местной системе координат измерител!
ного пункта (отсчитывается от направления на север по часовой стрелке'
(йи)г = — sin ф cos X, (6 21)г = — sin X,
(б12)г = — sin ф sin X, (6.22)г = cos Л, (6Д8.|
(^1з)г = СОЗф, (б2з)г ~ 0,
6 lj = 2 аЗа (6 la)ri 6 2j= 2 aia № !«)г (/ = 1» 2, 3), (6.18.1
<x=l a=l
<1 = Cj (Г — FBn) = 6 л (x XutI) 4" 6 j2 (у Уип) 4_6 13 Zun)> (6.18.i
b = e2 11 ’ (r — Fun) =6ji (л 3?uu) 4~ б 22 (у Уип) 4~ б 2g (z zH1I),
Г = 4r {b- (sin'XHa6) 4- &}, (6.18Л
A Ob
t. e.
Де = {a621 — 6flu},.
(a’4-62) I
~ i~2 j-Tai {я^зз b^ 12}, (6.18.5
(“’ + b’)
~ ,~s I Tai (a^23 ^13}-
(аг 4- /И)
4. ¥ = т — угол места точки в топоцептрической системе координа'
измерительного пункта (отсчитывается от горизонтальной плоскости)
(6 з1)г = cos ф cos X,
(б зз)г — cos ф sin X,
(б зз)г = sin ф,
з
6sj = 2 аза (^зя)г (/ = it 2, 3),
(6.19.|
(6.19.2
Определение параметров движения космического летательного аппарата ЗГ
С — вз * (г — Гип) — 6 at ^ип) 4“ б за (у Уин) 4“ б зз (2
zUIl), (6.19.3)
(6.19.4}
т. е.
~ <с - (siu Ynao) D},
(6.19.5}
7. ВЫЧИСЛЕНИЕ ПРОИЗВОДНЫХ ОТ ИЗМЕРЯЕМЫХ ФУНКЦИЙ
ПО ПАРАМЕТРАМ МАЛОВОЗМУЩЕННОГО УЧАСТКА ТРАЕКТОРИИ
Мы рассмотрели вопросы определения расчетных значений измеряе-
мых величин на моменты измерений, принадлежащие одному участку
траектории. Комплекс вычислений, связанных с обработкой измерений
одного участка траектории, завершается вычислением изохронных про-
изводных
I'-1.....Ь) (4.30)
по параметрам участка, при помощи которых формируются матрица G*
и вектор И(*> участка, выступающие в роли стандартного элемента кон-
струкции общей схемы формирования (4.55)—(4.56) нормальной системы
уравнений. Формулы-для вычисления производных (4.30) на основании
модели, изложенной в п. 2 раздела 4, приводятся ниже. Выражения для
производных (4.30) получены при помощи теории певозмущенного движе-
ния. Так как формулы для пройзводиых участвуют в массовых расчетах,
то их простоте должно быть уделено самое серьезное внимание.
Все дальнейшие формулы настоящего параграфа относятся к участку
с индексом к, а все содержащиеся в них величины связаны с прямоуголь-
ной системой этого участка. Поэтому без ущерба для пони-
мания индекс участка к, помер измерения п и индекс ^параметра q будут
в дальнейшем для краткости опускаться.
1. Выбор параметров участка. В настоящей работе в качестве пара-
метров участка выбраны элементы орбиты. Выбор продиктован рядом
преимуществ этих параметров по сравнению с прямоугольными коорди-
натами точки и их производными.
Одпо из основных достоинств — сравнительная простота аналитиче-
ских выражений для частных производных, рассчитываемых согласно
принятой модели. Действительно, для вычисления частных производных
по прямоугольным координатам необходимо иметь прямоугольные коор-
динаты и соответствующие компоненты скорости в начальный момент вре-
мени. Вычисление же последних па оскулирующей орбите требует реше-
ния уравнения Кеплера.
Другое серьезное достоинство — простота переноса матрицы нормаль-
ной системы уравнений, отнесенной к элементам орбиты, вдоль маловоз-
мущепного участка траектории. Матрица нормальной системы для эле-
32
Э. Л. Аким, Т. М. Энеев
ментов орбиты (при рассмотренной модели вычисления производны
может быть без изменения отнесена к любой точке маловозмущенно
участка траектории. В необходимом случае учет возмущений элемеш
орбиты при переносе такой матрицы сводится к перемножению ее с псц
единичными матрицами производных, т. е. не может вызвать существу
ной потери точности.
Следует также иметь в виду, что дальнейшее повышение точное
расчета производных (4.30) может быть достигнуто путем учета отдельш
возмущений орбиты. Учет возмущений в аналитических теориях про!
дится в большинстве случаев для элементов орбиты, что позволяет леи
использовать его в матрицах производных (4.30), отнесенных именно
параметрам орбиты'. •
2. Вычисление производных ио параметрам участка. Измеряем)
в момент времени t функция Y параметров траектории имеет вид: <
Y = ¥ {г {д, Т, t}, v {д, Т, t}; r3, v3). (4.2
Поэтому для частной производной по параметру д участка получи
выражение
as _ ат _ г ат ах ,_______. ат az у , г ат ,________, ат dDz ]
dq dq 1 dx dq dz dq J "T” (di>x dq ' di>z dq J’ '
которое в обозначениях (6.12) — (G.14)J принимает форму
L и М — векторы-градиенты, вычисления которых для ряда измеряем!
функций Чг были рассмотрены выше.
Таким образом, для производной (4.30) получаем выражение
а£ _ ( Lw, если ¥ = ¥ {г}, (7
( Z,(Q) + М{<л, если ¥ = VF {г, v). (7
Простота используемых в массовом счете формул для вычислен
А(Ч) и 7И(Ч) иллюстрируется на приводимых ниже примерах, в котор:
компонентами градиентов L и М являются их проекции на радиус, (
нормаль и трансверсаль:
Lr = у (Lr) = у (Lxx + Lvy + L2z),
= 7 (Ч) = «
~ у (Вс) = — (Lxcx -f- Lvcv Ltc^),
(Mr) = -J- (Mx x + Mvy + Mcz),
Mn = 4- (MrK)= 4- (MXXH + MjyH + (7
Mb = 4 (Me) = 4 (Мяся + MyCv +
1. Производные ио элементам p, e, co, cos i, Q,
Пусть параметрами участка g выбраны элементы орбиты
(7
р, е, со, cos I, Q, т;
Определение параметров движения космического летательного аппарата 33
р — параметр орбиты, е — эксцентриситет орбиты, ш — угловое расстоя-
ние перицентра от узла, i — наклонение плоскости орбиты, Q — дол-
гота восходящего узла. Тогда формулы имеют вид:
Lw = А (гЛг - 4 (t - т) {Lrvr + ЛЛ)},
L<e) = (~d^Lr + (г"дг)Ьп’
^)
£(я> = — Lxy + Lvx,
£(?) = — {Lrvr + Лкгк),
м(Р) = £ {у (<-*)£ ^ - у {Mrvr + MHvJ } ,
. М(е) =-£-{мг sine + Мм cos 0 — (-g-)Mr},. (7.8)
= VrMH ~ VnMr<
ctv г
= ~ (с, + ф МЬ,
Mw = — Mxvy + Mvvx,
Мм=±Мг,
(ГЭ = (/«-ц»)' {3 ~ fvn ~ + cV> sin0>’
(-£-) = -?{/("#)9in0 -фсозе}.
Выбор в качестве параметров q элементов орбиты (7.7) особенно удо-
бен в том случае, когда интегрирование уравнений движения проводится
в переменных с, f, т, так как в этом случае в процессе необходимого пере-
счета (см. Приложение II)
(с, f, т; t} -» {г, у; О
в качестве промежуточных величин вычисляются участвующие в форму-
лах (7.8) величины
с2, А г, гн (xn, yn, zH}, vr, vn, sin0, cos0. (7.9)
2. Производные по параметрам h, е, ш, cos i, Q, MQ.
В этом примере в качестве параметров q выбираются элементы орбиты
Л, е, со, cos i, Q, Af0; (7.10)
h — постоянная интеграла энергии, Мо — средняя аномалия в эпоху t0.
Формулы имеют вид:
^<h) = у |у ^о) (^rvr + r^r j 1
L<e) = 1/h {(C2 “ НГ) VuLr ~ (C’ + >) VrU},
£|““’ = _«ТфЛ-
3 Космические исследования, JA 1
£(о; — — yLx,
£'"->=iTk<L'’'+£-’“)’
"« = -i{3 (I - “ (Мл+ МЛ’) ’ <7Л1>
(" - V;
= - <’»Mr - v,Mn),
c’y
Mm" = Mb’
Af(n) = vxMv — vuMx,
Выбор в качестве параметров q элементов орбиты (7.10) особенно
целесообразен в том случае, когда интегрирование уравнений возмущен-
ного движения проводится в прямоугольных координатах. В этом слу-
чае для использования формул (7.11) должны быть дополнительно вычис-
лены следующие величины:
Л = V2 — с = [г х v], г„ = -1- [с х г], (7.12)
Г' = и2 + he2, vr = (rv), vu =
Следует обратить внимание на отсутствие в выражениях для производ-
ных тригонометрических функций от угловых элементов г, 12, со. Чтобы
избежать связанных с ними дополнительных потерь точности и машин-
ного времени при использовании соответствующих стандартных программ
в выражениях для производных фигурируют наряду с прямоугольными
координатами постоянные первых интегралов задачи двух тел с и f, свя-
занные простыми формулами с прямоугольными координатами.
Выбранные элементы участка орбиты (7.7) и (7.10) могут быть исполь-
зованы в задаче обработки траекторных измерений для существенно эл-
липтических и гиперболических орбит, но они но пригодны для круговых
орбит.
Волос универсальными параметрами q являются элементы орбиты вида
h, <Pi, ф2> cos г, 12, v0; (7.13)
h — постоянная интеграла энергии, = е sin со, ср2 = е cos со,
v0 — п (t0 — iyxi), 2у:Ш — время прохождения через узел. Эти элементы
можно использовать для эллиптических и гиперболических орбит
с любым эксцентриситетом.
3. 11 р о и з в о д п ы е и о п а р а м е т р а м h, cpj, сра, cos i, 12, v.
= 4 {4 — ^LrVr+LnvJ - rLr}.
Определение параметров движения космического летательного аппарата 35
Л(СО31,-“71Т^’
^(О) = xLv — yLx,
M<h> = (3 (z ~ z«) £ M' - + M^}.
= ^{(vnMr - vrM„) Л -
= g-a- (^r-
c2r
•fl/(cosi) = a 2 ^bi
(cx+c;)
M (O) = vxMv — vuMx,
(7.14)
1 C . 1 (‘xW1
<P1 1‘ (cxH-<)‘/* l' ‘P2 И (Cx + V'’ ’
7?1 = 0 H! q,,) {fyl C0S U + (J — Фа) sin U>-
51 = '(Г-f <p4 {2(P1 SIU u ~ c°s>h
/?2 = ТГ-U m' 1» <Ф1Ф2 siil u + 1(1 + Фг)” — Ф11 COS uh
s2 = ~(i 4-\p,)» {4)1<p2 cos “ — ((1 + Ф2)2 — q>Jl sill U),
Tl=lT^{2 +ф2) (1 ~ф2) + 51’
7 2 = 0 _|_ tp.j)2 (2 "l" Фг) Ф1Ф2 + ^2'
8. ВЫЧИСЛЕНИЕ МАТРИЦЫ ПРОИЗВОДНЫХ ОТ ПАРАМЕТРОВ ДАННОГО
УЧАСТКА ТРАЕКТОРИИ ПО ПАРАМЕТРАМ ПРЕДЫДУЩЕГО УЧАСТКА
В роли звена, связывающего соседние участки в схеме формирования
(4.55), выступает квадратная (6x6) матрица
(8-1)
производных от параметров данного участка траектории но параметрам
предыдущего участка. Ниже проводится вычисление этой матрицы в об-
щем случае, когда основные притягивающие тела обоих участков не сов-
падают.
Предположим, что в момент Tw перехода с участка (/с — 1) на уча-
сток /с положение в скорость точки заданы в прямоугольной системе
координат
yw, . . у™} (8.2)
36
Э. Л. Аким, Т. М. Энеев
и в системе координат 1)Z<* n
^(*-1) = {X(*-1>, у^\ . . . »<*-»}. (8.3)
Однозначная зависимость между параметрами (/с — 1)-го и /с-ro участ-
ков
Т^. /’(*>} (8.4)
является ношений логическим следствием последовательности очевидных соот- qw = qw {£PW-, TW}, (8.5) Г(Л) = r(A-i) _ г**-!), (8.6) = у(Л—1) y(fc-l) °* ^(*-1) = ^(*-1) ^(*-1), Г(*)}. (8 7)
Соотношения (8.5) — (8.7) приводят к представлению матрицы (8.1)
в виде произведения матриц
I <]«-» J J I qt*-» J
(8.8)
Вычисление матрицы (6 X 6)
(dqi Oqi_ ,
&х г i /I 'r/|) М*171’ \
................ = { L . ,М . ) (8.9)
1 з \ I <"•’ М ("‘) /
dtfa д([$ I \ Ъ IV1 '
~дГ > • • • > щГ- /
производных от параметров q№ участка по параметрам (8.2), где
Р-] (8-9.1)
[ дх ду dz J [ t)vx ’ dvv ovt J
выполняется для значений
г = r(A) v — v(*> (8.9.2)
и величины ц основного притягивающего центра /с-ro участка.
При выборе в качестве параметров переменных
р, е, ш, cos i, Q, т (7.7)
или
Л, е, ш, cos i, Q, Мо (7.10)
расчетные формулы принимают вид (7.12) и далее
f = —цу+Ivxc], f{/x,/У>Л);
= -А-; xs = А; К = -7—,——rv;
1 2 2р/ ’ 3 2/ * 4 (сх/у-%4)
х» ~ (с2 + с2} > хв - р—; — -2ZS-; К - (— sisn Л) •
Ьв = ^3(/-/0)-^]1_-£}хв; Х10 = |з(/ - ta)^ - 2]-U„; (8.9.3)
Определение параметров движения космического летательного аппарата 37
х; = {[з (I - т) - -£-} Х'1О = [з (/ - т) А - 2] А'в;
, _ (с’ —Иг) '
Л11 — --------------------------------/---- л8»
L(h) = 2Х)Г,
М(М = 2v,
L<p) [v X с),
М(р) = -1 [с X г],
L(e) = X2L(h) + X3L(P),
M(e) = X2M(/l) + Х3М(р),
L(ai) = Х4 {XbL(c0Sil - XeL(e) + X,zr - vzv + шг),
где jfcjO.-O.-fr}, шв {0, 0, — rvr},
м(“> = X, {XsM(c03i> - XeM(e) + 2zv — vzr + (!)„},
L(co,0 = yr _ x7L(p), где Yr vu, -1 vx, o} , Y J- Ly, 1I( o} t
M(C08<) = Yv- X7M(P),
L(O) = X12c,
M(Q) = X13c,
L(M,) = X8r -h XuL<0),
M(Mo) = X10v + XuM(e),
L(T) = X'Br -I- X'nL(<”,
m(t) = x;nv + XijM1”.
dx
dgt
J
dvt
9gi
производных от величин (8.3) по параметрам участка проводится
для значений
г = V = v<*-!> (8.10.1)
и величины р. основного притягивающего тела (А — 1)-го участка.
Расчет производится:
*для р, е, о, cos i, Q, т
по формулам (7.3) — (7.6), (7.8), (7.9);
для Л, е, w, cos i, й, Mt
по формулам (7.3) — (7.6), (7.11), (7.12);
для Л, <рц <Ра, cos i, й, v(
по формулам (7.3) — (7.7), (7.12), (7.14).
(7-7)
(7.10)
(7.13)
.38
Э. Л. Аким, Т. М. Энеев
Выбирая в качестве компонент векторов L и М
Lx,Lv,Lt- MX,MV,MZ
элементы строки квадратной (Охб) единичной матрицы н пользуясь од-
ной из указанных выше групп формул, получаем соответствующую стро-
ку матрицы (8.10). Последовательно перебирая все строки единичной
матрицы, вычисляем при помощи этой группы формул все строки матрицы
(8.10).
При записи схемы формирования нормальной системы уравнении в
форме (4.55) вычисление квадратной (0x6) матрицы
(8.И)
производных от параметров </<*> данного участка траектории по уточняе-
мым в задаче параметрам Q проводится по рекуррентной формуле
(Qj !#•<*’][ I 2 J
(8.12)
при помощи матриц (8.9) и (8.10) (при k ~ 1 выражение в квадратных
скобках формулы (8.12) обращается в единичную матрицу).
В случае общего основного притягивающего тела для (к — 1)-го и
к-го участков траектории (например, для орбиты искусственного спутника)
матрица производных (8.1)может быть использована для учета возмущений
в производных, накопившихся за время движения па (к —1)-м участке
(см. н. 5 раздела 4).
9. ОЦЕНКА ТОЧНОСТИ ОПРЕДЕЛЕНИИ ПАРАМЕТРОВ ТРАЕКТОРИИ
Выше была рассмотрена задача определения параметров траектории
Q\, • • •> Qn (9.1)
по данным траекторных измерений.
Ошибки траекторных измерений приводят к погрешностям в величине
определяемых параметров траектории. Как следует из системы уравнений
(3.5), ошибки определяемых параметров связаны с ошибками измерений
нелинейными зависимостями. Если ошибки траекторпых измерений носят
случайный характер (с нормальным многомерным законом распределе-
ния), то определяемые параметры траектории также являются случайными
величинами и решение задачи дает нам в качестве параметров (9.1) мате-
матические ожидания
q\> р.2)
этпх случайных величин.
Для приближенного анализа нелинейных связей подвергнем соотно-
шения, связывающие ошибки траекторных измерений и ошибки парамет-
ров (9,1), линеаризации. В рамках построенного ранее алгоритма опре-
деления математического ожидания параметров траектории получаем
интересующее пас линейное соотношение (3.18).
которое может быть переписано в виде
AQ - I (4) 5’ (9.3)
Определение параметров движения космического летательного аппарата 39
где I — прямоугольная (6 х А) матрица
' (4) _ _ Д-.Д _ _ [Г ( - ) РТ/ (Л)]" р- ( - ) РТ] ; (9.3.1)
5 (£i,..-,£,v) и AQ (А^,..., А<>0}— соответственно вектор ошибок траек-
торных измерений и вектор ошибок параметров траектории. При доста-
точно малых ошибках траекторных измерений линейное соотношение (9.3)
мало отличается от нелинейного и может быть использовано вместо пе-
го в исследованиях.
Алгоритм определения ожидаемых значений параметров траектории
опирается па определенные гипотезы относительно вероятностных харак-
теристик ошибок траекторных измерений. Гипотезы представлены в
алгоритме (а следовательно, и в соотношении (9.3)) корреляционной матри-
цей А'«у. Наряду с такими гипотезами существует реальная статистиче-
ская картина ошибок измерений, характеризуемая корреляционной ма-
трицей, которую обозначим через Ry (в отлично от матрицы Исполь-
зуемые в алгоритме обработки измерений вероятностные характеристики
в большинстве случаев лишь приближенно отражают эту фактическую
картину ошибок измерений. Одна из причин приближенного рассмо-
трения заключается в том, что реальная статистическая картина.ошибок
траекторных измерений почти всегда неизвестна, поэтому приходится
довольствоваться весьма грубыми ео приближениями. Другая серьезная
причина кроется в преднамеренном упрощении алгоритма обработки.
Последнее требует дополнительного разъяснения.
Действительно, целесообразность упрощения алгоритма, используе-
мого в массовой обработке измерений, совершенно очевидна далее и тогда,
когда реальная статистическая картина ошибок измерений известна.
Однако использование приближенных вероятностных характеристик
ошибок измерений в' алгоритме обработки приводит к погрешности в ве-
личине ожидаемых значений параметров траектории. Поэтому прибли-
жение может быть оправдано лишь в тех случаях, когда эта погрешность
оказывается малой по сравнению с величиной разброса параметров вокруг
ожидаемых значений.
Так, например, в большинстве практически важных и интересных слу-
чаев целесообразно применять метод наименьших квадратов для обра-
ботки траекторных измерений, даже если измерения оказываются сильно
коррелированными, так как возникающая при этом методическая по-
грешность весьма незначительна.
Таким образом, положенные соображения приводят к необходимости
различать реальную статистическую картину ошибок измерений
и ту приближенную статистическую картину А'«у, которую мы закладываем
в качестве гипотезы в алгоритм определения ожидаемых значений пара-
метров траектории. Это следует иметь в виду во всех дальнейших рас-
смотрениях.
1. Оценка точности прогнозирования при известной статистической
характеристике ошибок траекторных измерений. В отличие от задачи
определения ожидаемых значений параметров траектории задача оценки
разброса параметров вокруг ожидаемых значений предъявляет более вы-
сокие требования к полноте используемых нами знаний о характере оши-
бок' траекторных измерений. Поэтому для оценки разброса параметров при
рассмотренном выше алгоритме определения их математического ожида-
ния используем для статистической характеристики ошибок измерений
корреляционную матрицу Ку.
40
Э. Л. Лкил1, Т. Л1. Энеев
Разброс параметров Qt, Q2, . . Qa вокруг ожидаемых значений (9.2)
оценивается корреляционной матрицей Kq. Вычисление этой матрицы
при помощи матрицы Ку и линейного соотношения (9.3) может быть про-
ведено обычным в статистике методом 151
В результате получаем выражение для матричного элемента
(i / = 1, . . .,6) матрицы Kq:
N N
ад = 2 2 [/ (-5U \r (4)L • <9-4Л>
а=1 (1=1 L L J1P
Подставляя (9.3.1) в (9.4) перепишем выражение для Kq (с учетом
А = А и Ку = Ку) в виде
Kq = [Л’1/* (-£) Ку [Ку! (9.4.2)
Если в алгоритме определения ожидаемых значений параметров тра-
ектории и в оценке разбросов параметров вокруг них используется одна
и та же гипотеза о статистическом характере ошибок измерений, т. е.
Ку — Ку,
то выражение (9.4.2) принимает простой вид
Kq = Л'» [/’ (X) Ку1! (^-)]Л'1 = Л-4 (9.4.3)
Пусть переменные
W {Жп W2, . . ., Wn) (9.5)
связапы с параметрами
Q{<2d Qt......<?.)
соотношением
W = {Q, То) (9.6)
и задана матрица производных
так что связь между ошибками AQ и AW задается линейным соотношением
AW - / AQ. (9.7)
Тогда разброс параметров W вокруг ожидаемых значений W в ли-
нейном приближении характеризуется корреляционной матрицей К\у:
= (9.8)
Определение параметров движения космического летательного аппарата 41
Л
< Главпая диагональ этой матрицы даот выражения для среднеквадра-
тичных ошибок <iivt (i = 1, 2, . . . п) параметров (9.5):
; = V(^w)h- (9-9)'
У
• 13 /г-мерном пространстве векторов W область возможного разброса
я переменных (9.5) вокруг их ожидаемых значений W с заданной вероят-
ностью Р^ попадания точки в эту область ограничена поверхностью
«эллипсоида ошибок». Эллипсоид ошибок определяется корреляционной
|‘ матрице '' К\у. Центр эллипсоида совпадает с концом вектора W ожидае-
мых значений. 13 системе координат с началом в центре эллипсоида и
осями С/j, . . ., Un, соответственно параллельными осям И^, . . . Wn,
уравнение эллипсоида имеет вид
(tf^U, U) = (9.10)
где U {Ui, . . ., Un), а£ — константа (при £ = 1 уравнение определяет
«среднеквадратичный эллипсоид» ошибок).
Приведение при помощи ортогонального преобразования квадратич-
ной формы
(К$ и, U)
к главным осям позволяет определить геометрические размеры эллипсои-
да (полуоси) и его ориентацию в рассматриваемой системе координат.
В наиболее наглядном случае п = 2 пространство переменных (9.5)
представляет собой плоскость W {Wu W2}, а эллипсоид ошибок вырож-
I дается в эллипс ошибок. Вероятность Рг_ попадания вектора W внутрь
I этого эллипса равна:
Рс = 1 - e-W (9.11).
[ (при £ = 1 Р? ~ 0,39, при £ = 3 Р^ ~ 0,99).
Выражения для параметров эллипса ошибок (большая полуось а,
малая полуось b и угол % между направлением большой полуоси и осью IVj)
содержащие элементы корреляционной матрицы
(9-12)
принимают, как известно, особенно простой вид:
а2 = Ц-{ц7п + ш22 + т]},
= 4* {u’u+и’я — П).
clg * = + т1}’
где Т) = J/ (шп — ш22)2 4- 4ш22.
Очевидно, формулы (9.9)—(9.13) пригодны и для оценки разброса
параметров Q вокруг ожидаемых значений Q*, что соответствует частно-
му случаю соотношения (9.6)
(9.13).
W Q.
(9.14)
Э. Л. Аким, Т. М. Энеев
Полученные параметры траектории космического аппарата можно
использовать для решения ряда самостоятельных задач: дальнейшего
прогнозирования движения аппарата, управления его полетом и т. д.
В этих задачах весьма важным является удачный выбор оцениваемых па-
раметров W. Этот выбор обычно неразрывно связан с физической сущ-
ностью данной конкретной задачи и проводится для нее сугубо индиви-
дуально.
Для примера рассмотрим оценку точности прогнозирования движения
космического аппарата нри сближении с планетой. Предположим, что
в результате обработки траекторных измерений получены ожидаемые
значения (9.2) параметров траектории и корреляционная матрица (9.4).
Движение аппарата в сфере действия планеты происходит почти но ги-
перболической траектории с фокусом в центре масс планеты. Ожидаемое
сближение космического аппарата с планетой без учета ее притяжения
(для устранения возникающих нелинейностей) можно охарактеризовать
вектором W* «прицельной дальности» этой гиперболической траектории.
Вектор W «прицельной дальности» представляет собой перпендикуляр,
опущенный из фокуса гиперболы на ее асимптоту, равный по величине
мнимой полуоси гиперболы.
Построим плоскость, проходящую через центр масс планеты и орто-
гональную вектору иланотоцентрической скорости Voo аппарата при без-
копечном удалении его от планеты («картинная плоскость» планеты).
Введем в плоскости прямоугольную систему координат с началом в центре
масс планеты и единичными векторами l}x), 12х) по осям (х — индекс мало-
возмущепного участка траектории, совпадающего со сферой действия
планеты).
Оси системы координат выбираются при помощи заданного единич-
ным вектором 1цХ) направления следующим образом:
Вектор прицельной дальности расположен в картинной плоскости
планеты. Положение этого вектора задается во введенной системе коор-
динат двумя координатами W (Ц/п РИ2):
= Wl}x) = - — wi'x),
(9.16)
W2
где
W = lc X f]j , W{bx,bu, Ьг). (9.17)
(Здесь п в последующем под с, f, h принимаются элементы орбиты ле-
тательного аппарата в сфере действия планеты.)
Разброс в картинной плоскости планеты вектора прицельной дальности
W вокруг ожидаемого значения W' определяет точность прогнозирования
движения космического аппарата при сближении с планетой без учета ее
притяжения. Величина разброса характеризуется .корреляционной ма-
Определение параметров движения космического летательного аппарата f 43
(ицей (9.12), при помощи которой (по формуле (9.13)) строится в картип-
)й плоскости планеты эллипс ошибок прогпозировапия.
В формуле (9.8) расчета корреляционной матрицы Куу участвует ма-
ица производных (9.6.1), вычисляемая в виде произведения матриц
(9.18)
Матрица I{q^/Q} производных от параметров плапотоцентриче-
ого участка траектории по параметрам Q рассчитывается по формуле
.12).
Прямоугольная матрица (2 хб) производных
(9.19)
: переменных W\, W2 по параметрам q<*> плаиетоцентрического участка
>аектории вычисляется по формулам:
Х6 =
Xg =
^rp- UoVco);
У /*
c
x2 = V1 — x?:
xe = VK2xe; x7 =
х3
х4 =
*1.
X2 ’
с
aiVi
‘>p
de
dWt
дш
dVVj
d cos i
dW3
-ДИ,.Х" = (!;oo)*X4i Xl° ~
- x3VKi:
= -(x^j + J£-x4);
= Xgjy, - ~ x4;
У
= xfl(xBiy1-iXr);
X4
У/Г
arr3
a?
div2
de
dW2
dw
[(^oo)v^ox (^оо)д^01<] •
= х31И2;
= — x7lV2;
= х81У2;
(9.20)
d cos i
хв(ха1У2 — c2) — XsC 1"г
ап
dW2
да
div,
dr
dW2
dt
x2/A ’
= ~ {(bulox — bxlou) — х10ИЛ);
= {77;
= 0;
ox Cxlay — Xlo
= o,
Xj =
х3 =
yj— [с X I] pfl> Vqo {(^oo)x, (l’co)v, (^00)2)1 Iq {^ox> ^02}*
Формулы (9.8)—(9.13) и (9.15)—(9.20) позволяют прогнозировать ожи-
даемое сближение (с оценкой возможного разброса) космического аппа-
рата с планетой.
I 2. Оценка точности прогнозирования при неполной статистической
[арактсрнстике ошибок траекторных измерений. Выше приведена мето-
44
Э. Л. Аким, 'J'. М. Энеев
дика оценки точности прогнозирования, основанная па предположении,
что заданы в виде корреляционной матрицы вероятностные характе-
ристики ошибок траекторных измерений. Одпако основная трудность
оценки точности прогнозирования заключается именно в том, что эти
вероятностные характеристики в большинство случаев известны в явно
недостаточной море. На практике априорная информация о характере
ошибок траекторных измерений исчерпывается, как правило, общими
представлениями о порядке максимальных ошибок последних. Что же
касается корреляционной зависимости между ошибками отдельных изме-
рений, то обычно опа оказывается неизвестной.
В связи с этим возникает необходимость оценки точности прогнози-
рования в тех случаях, когда статистическая характеристика ошибок из-
мерений неполна и содержит лишь звание дисперсии ошибок. Для этой
цели методика оценки, рассмотренная в предыдущем раздело, оказывает-
ся недостаточной и требует известного дополнения. '
В интересующих нас случаях каждому конкретному виду корреля-
ционной зависимости между ошибками измерений соответствует своя,
вполне определенная величина разброса оцениваемых параметров, вы-
числяемая согласно формулам (9.4.2), (9.8). Одпако конкретный вид нор-
мированной корреляционной матрицы Ну (матрицы, составленной из
коэффициентов корреляции)
[7?чНа.р = бааН/Аг1а3 (а, Р = 1---N) (9.21)
неизвестен. Поэтому при заданной дисперсии о® (а = 1, ..., N) ошибок
измерений подбирается наиболее неблагоприятный вид корреляционной
зависимости 7?т, соответствующей максимальному разбросу оцениваемых
параметров. Полученная величина разброса принимается в качество оцен-
ки точности прогнозирования для выбранного состава траекторных из-
мерений.
Допустим, что ошибки измерений, рассматриваемые в зависимости
от времени, представляют собой стационарный случайный процесс. Кор-
реляционная матрица Ку такого стационарного процесса может быть в
большинстве случаев приближенно представлена формулой, содержащей
выражение для матричного элемента (см. [6])
. = 62L-?|T| cos((0lIl + ^)l, т = «я —Гр, (9.22)
I COS ф I
где k, (о, ф — параметры корреляционной матрицы. Матричный элемент
(9.22) характеризует корреляционную зависимость между измерениями
Ya, Yp (а, [1 = 1, ..., ^однородного состава, отнесенными соответствен-
но к моментам времени 1а, tr,. Для любой корреляционной зависимости,
задаваемой набором параметров
1с, СО, ф,
формулы (9.22), (9.4.2), (9.8) позволяют получить соответствующий раз-
брос оцениваемых параметров и, следовательно, провести оценку точности
прогнозирования.
Рассмотренная методика оценки точности прогнозирования может
быть использована также в процессе баллистического проектирования
космического полета для выбора эффективного состава и программы тра-
екторных измерений.
Определение параметров движения космического летательного аппарата 45
ПРИЛОЖЕНИЕ I
Метод наиекорейшего спуска
( 1. Известно, что метод наиекорейшего спуска при решении систем уравнений
может быть реализован различными путями. Пусть пеобходимо решить систему т
уравнений с т неизвестными
Ei(Qb <2>.....Qm)=0 (i = l. 2./П).
(М)
Поставлеппая задача будет решена, если, построив функцию
т
ф = 2 &
(1.2)
найдем значения неизвестных Qi.Qa,..., Qm, для которых функция Ф принимает ми-
нимальное значение, равное пулю.
I Любой способ реализации метода наиекорейшего спуска сводится в конечном
Течете к построению алгоритма, на каждом цикле которого изменение переменных
'Qi, <?2>.-ч Qm ведет к уменьшению функции Ф.
i Наиболее простой и, на первый взгляд, наиболее естественный вариант метода наи-
’скорейшсго спуска заключается в построении алгоритма, на каждом цикле которого
изменение переменных Qi, Q3,..., Qm происходит в направлении, обратном градиен-
ту функции Ф в данной точке. При этом метод наиекорейшего спуска сводится фак-
тически к тому или иному способу интегрирования системы обыкновенных дифферен-
циальных уравнений [7, 8]
dQj
ds
ЭФ
1 (i = 1,2...т)
к г дф \2
(1.3)
с начальными данными, соответствующими первому приближению неизвестных Q(.
Следует, однако, отметить, что практическая реализация указанного варианта метода
наиекорейшего спуска при решении большинства прикладных задач наталкивается, на
большие трудности вычислительного характера, появляющиеся при сложном релье-
фе гиперповерхности Ф((?ь Qt,.--, Qm} и при сильном отличии друг от друга относи-
тельной крутизны спуска по отдельным переменным В настоящей работе приво-
дится вариант метода наиекорейшего спуска, при использовании которого эти труд-
ности легко обходятся. Построение метода основывается па преобразовании перемен-
ных к новым переменным обеспечивающим в некотором смысле «равновесомость»
всем составляющим градиента функции Ф по этим переменным.
Рассматриваемый вариант метода наиекорейшего спуска был разработай в 1957 г.
и использовал при обработке астрономических наблюдений первых искусствепных
спутников Земли [2]. При этом была выявлена его весьма высокая эффективность.
Следует отметить также, что возможность построения подобного метода была незави-
симо показапа М. К. Гавуриным [9], который подошел к идее ого создания с несколько
иной стороны.
2. Введем новые переменные связанные со старыми Qi при помощи формул
(см. (1.1))
(1.4)
«2х. <?>.....Q,„) (i = 1, 2....т).
Очевидно, что решению Q^\ ,..., системы (1.1) в переменных Qi будет
соответствовать решение этой же системы ’= g*** =...=gW = 0 в переменных £(.
Выражения для составляющих градиента функции Ф в подпространстве перемон-
, так же как и уравнения (1.3) движения но градиентным линиям, в новых пере-
силыю упрощаются и принимают вид
ВЫХ £(,
ценных
Л'
ds
(i = l, 2, . . .,т).
/Ф
(Т.5)
46
Э. Л. .1кил(, Т. М. Энеев
Система (1.5) легко интегрируется, и решение ее можно записать в следующем’
виде:
(i = 1,2,. .ш),
(1,6).
где
7П
Sjn) = S4 «?10). Qj0)--------Q™). Ф(о) = 3 [^0)
1=1
При этом величины Q*0), Q^.-.-.Q^ представляют нулевое приближение для иско-
мых Qi, Qa,..., Qm.
Из (1.2) и (1.6) легко видеть, что наискорейший спуск в пространстве Ф, £а,.._
тп
... , £т представляет собой «движение» по образующей параболоида Ф = У| St •
1=1
Ив (1.6) видно также, что соответствующее движение в подпространстве про-
исходит по прямым линиям.
Наискореиший спуск в рассмотренном варианте является в некотором смысле
наиболее быстрым «движением» к искомому решению. Очевидно, что переменные
Qi. Qi,---, Qm должны меняться при этом, согласно соотношениям (1.4), (1.5) или 1.6).
Так как соотношения (1.4) пе могут быть в явном виде разрешены относительно
переменных Qi, Qa,..., Qm, то для фактического вычисления этого изменения целссо* i
образно использовать систему дифференциальных уравнений (1.5), предварительно
преобразовав ее к переменным Qi, Qm. Легко видеть, что после такого преобра-
зования система (1.5) примет вид
dQi dQi dQm Si ,
flil ds + ”’2 ds + • • • + flim ds ~ ~ j/ф = 1. 2, . . ., m) , (1.7) /
\ 1
или же, принимая во внимание (1.2) и (1.6): Т 1
w”| 1
rfQi </Qi dQ,n Siu>
“11 "ST + “12 "ST + • • • + «im . (1-8).
где
Таким образом, рассматриваемый вариант метода наискорейшего спуска сводится
к интегрированию (вообще говоря, численному) системы дифференциальных уравнений-
(1.7) или (1.8) при начальных условиях, соответствующих нулевому приближению
искомых неизвестных Qi, Qa,..., Qm. Интегрирование необходимо вести, согласно (1.6),
до значения независимого переменного s, равного УФ(и\
Более целесообразно, одпако, окончание интегрирования производить, непосред-
ственно контролируя обращение функции Ф в пуль, так как неизбежные ошибки чис-
ленного интегрирования уводят с «кратчайшего» нути и конечное значение s, при
котором Фобращается в нуль, будет отличаться от УФ(0,.По этой же причине при прак-
тическом применении предлагаемого метода целесообразнее использовать систему
уравнений (1.7), а не (1.8), причем систему (1.7) можно интегрировать достаточно
грубо.
3. Выше был рассмотрен случай, когда число уравнений ?{ = 0 равно числу иско-
мых неизвестных Qj. Рассматриваемый вариант метода паискорейшего спуска может
быть, одпако, легко обобщен па случай, когда число уравпоний = 0 превышает число
неизвестных Qj. Последнее будет иметь место, панримор, при статистической обра-
ботке избыточного числа измерений.
Определение параметров движения космического летательного аппарата 47
(" = 1,2.........М). (1.13)
Пусть произведено М измерений параметра Чгп, который связан с искомыми па-
раметрами Qi, Qa Qm зависимостью
'‘'я = «2*. <2»....<2,п) (н = 1,2......М), (1.10>
причем М > т.
Составим разности
Ь^пнаб-Ч'пССь .............Qm). U H)
где представляет собой измеренное значение параметра Уд. Составим далее квад-
ратичную форму
м
® = 2 & (М2)
П=1
* Согласно методу наименьших квадратов, параметры Qi, Qt,..., Qm должны быть
подобраны таким образом, чтобы после подстановки их в формулы (1.11) квадратич-
ная форма (1.12) приняла минимальное значение. Нашей целью является построение
[алгоритма изменения переменных Qi, Q2,..., Qm, на каждом шаге которого это изме-
нение будет приводить к возможно более быстрому уменьшению текущих значений
функции Ф.
Если рассматривать квадратичную форму (1.12) как функцию М независимых
{переменных £п, то наиболее быстрое уменьшение функции Ф будет достигнуто при
[движении по градиентным линиям гиперповерхности (1.12). Движение это будет опи-
сываться дифференциальными уравнениями наиекорейшего спуска (см. также (1.5))
, ds
I
t В малой окрестности данной точки это движение может быть также описано
[конечными формулами
= (114>
где Д« есть, очевидно, малое перемещение вдоль градиентной линии. Однако ввиду
того, что М переменных связаны М соотношениями (1.11) с вд переменными Qi,
причем т < М, движение, удовлетворяющее формулам (1.14), вообще говоря, осущест-
вить невозможно. Можно, однако, попытаться найти такое направленно перемещения
вдоль гиперповерхности (1.12), которое, удовлетворяя соотношениям (1.11), в то же
время было бы наиболее близким к направлению (1.14).
Пусть церемонные Qi,Qs,..., Qm получат малые измолоти
AQi. А<2э....AQ,n. (1.15)
Согласно формулам (1.11), в этом случае будем, очевидно, иметь
Д£п = ап1Д<?1 + ап2Д(2’Н--h °птД<?т. °ni = Jq7 ’ <1Л0)
где ДТ„ — соответствующее AQi, AQ2,. . ., AQm изменение функции £п. Составим
новую квадратичную форму
м
<₽= 2 ЩП-ДГП)’. (1.17)
п=Т
Очевидно, что квадратичная форма (1.17) в некотором смысле характеризует
близость малых перемещений вдоль градиентной линии гиперповерхности
Ф (£ь Ь. • • • > 5д,) и вдоль направления, характеризуемого формулами (I/1G). 11о-
48
Э. Л. Аким, Т. М. Энеев
ставим задачу: для данного As найти такие перемещения AQi, AQI(.... AQm, для
которых квадратичная форма (1.17) принимает минимальное значение. Выписывая
необходимые условия экстремума (минимума) функции m переменных
д<р
~д^=° (‘ = 1.2.....т) . (1.18)
получим m линейных алгебраических уравнений для определения перемещений
AQi, AQa,. . . , AQm. После преобразований и упрощений эти уравнения примут
следующий вид:
Bi
Л1Д(?« + Л2Д<?> +• • •+ AlmbQm = - y=- As, (1.19)
•где
Aij~ 2 dQ{ 0Q- • 2 dQ{ (1-20)
n=l 3 n=l
Из (1.19) и (1-20) легко видеть, что матрица системы (1.19) совпадает с матрицей
системы нормальных уравнений метода наименьших квадратов. Правые части системы
(1.19) с точностью до множителя As/^Ф также совпадают с правыми частями тех
же нормальных уравнений.
Разделим правые и левые части уравнений системы (1.19) на As и совершим про-
дельный переход As-»0. В результате получим систему дифференциальных уравнений
dQi dQi dQm Bi
+ + (‘ = 1.2.....«)• U-20
Система (1.21) интегрируется no s при начальных данных
S = 0, Qi = Q<°>, Qa = Q<°>....Qm =
до значения s= s(/l>, для которого величина F = 0, где
Л = 2 В*.
1=1
Полученные при s = sW значения параметров Qi, Q2,..., и будут искомыми зна-
чениями, наилучшим образом согласующимися с наблюдениями в смысле метода наи-
меньших квадратов.
Если в задаче статистической обработки избыточного числа измерений исполь-
зуется метод максимального правдоподобия (обобщение метода наименьших' квадра-
тов), то вместо (1.12) минимизируется более общая квадратичная форма
м
Ф= 2 (1-22)
п,Ч=1
в которой Рп1 — элементы положительно определенной матрицы (см. раздел 3).
В этом случае, повторяя приведенные выше рассуждения по отношению к (1.12), при-
ходим к системе дифференциальных уравнений, совпадающей по виду с (1.21). Мат-
рица и правые части этой системы формируются при помощи выражений
М -it at М at
2 2 Л..А.(i-23)
П, 1=1 1 n, 1=1 ’
Система дифференциальных уравнений (1.21) является аналитической записью
.алгоритма изменения переменных Qi, Q!t..., Q,n, приводящего к минимуму функции Ф.
Определение параметров движения космического летательного аппарата 49
ПРИЛОЖЕНИЕ II
Вычисление положения и скорости точки
но элементам орбиты
Пусть и инерциальной прямоугольной системе координат XYZ, начало которой
совпадает с притягивающим центром, задана невозмущеппая орбита точки элемен-
тами (см. раздел 5, п. 3)
Н/х. Л), т (П.1)
с фокусом в начале системы координат.
В этом случае вычисление радиуса-вектора
г {X, у, Z} (П.2)
и вектора скорости
v{vx> vir (IL3)
точки в момент времени t производится по формулам:
с« = (сс), /а = (If),
1 - 1
ei = -7-f, е3= — с, е2=[е3Хс1],
/• с
er = cos 0 ci sin 0 е-,
еп = cos 0 02 — sin 0 сц
са 1 н -I- / «’Я О
г р + /cos 0 ’ г с- * ( ‘а)
г = Г ег,
/ с
vr — — sin 0, г>и = — ,
v=₽re, + »BeH,
где р.— произведение гравитационной постоянной на массу притягивающего центра,
О — истинная аномалия точки. Вычисление истинной аномалии точки в заданный
момсит времени t требует численного решения трансцендентного уравнения Кеплера
для эллиптической орбиты и ого аналога для гиперболической орбиты. Решение
этого трансцендентного уравнения может быть проведено, например, методом
Ньютона, в результате чего легко приходим к приведенным ниже формулам
с = Л/ = -А-(1/а-На1)М'Ь-,С) (П.5;
Для случая эллиптической орбиты (е <И) итерационный цикл
^п+1 — Е„ — е sin Е„ — М
(п П 4 9 ..) (1I.G.1)
1 — с cos Zin 4 » » » •
с нулевым приближением Е3 = М Д- е sin Л/. (П.6.2)
Цикл завершается при и = I, для которого выполняется неравенство
i Л { I < б,
где 6 — заданная константа
В результате
sin /4 ^1 — с2 eos e
sinO =—t _ e cos E . cos 0 — i _ e cos Et • (11.0.3)
Для случая гиперболической орбиты (е>1) итерационный цикл
2 (Л/ + IIп) exp (- II п) - е (1 - охр (- 2ЯП))
11n+1 — //п — 2 ехр (— Нп) — е (1 + ехр (— 2//п))
(я = 0,1,2,...) (11.7.1)
с нулевым приближением
//о = 2.
(11.7.2)
4 Космические исследопаиин, |
50
Э. Л. Аким, Т. М. Энеев
Цикл завершается при n = Z, для которого яыполпяется неравенство
|/7[ — где 6' — заданная константа.
В результате
(1 — ехр( —2Л/,)) 1
SIn ° “ е (1 -Н охр ( — 2Я()) — 2ехр (— 7/() ’
2 ехр ( — Нг) е — (1 + ехр (—2/7())
303 ° = е (1 охр ( - 27/,)) - 2ехр( - 77,) ' (1L7'3)
ПРИЛОЖЕНИЕ 111
Вычисление элементов орбиты точки по ее радиусу-вектору
и вектору скорости в заданный момент времени
Пусть в инерциальной прямоугольной системе координат в момент времени I
задал радиус-вектор точки
г {А у, z} (III.1)
и ее вектор скорости
v{rx, vy, Vz} (Ш.2)
Точка движется по невозмущенной орбите с фокусом в начале системы коор-
динат.
Вычисление элементов орбиты
c{cx, Су, сг}, 1{/д., /и, fz}, т (III.3)
производится ч этом случае по формулам:
с — [г х vj,
ц.
f = —г + [v X с), f’ = (ff), (Н1.4)
2u
/ h = vt — -y~.
Для эллиптической орбиты (Л<^0)
sign (cos £) = sign (а2—~г~) ’
т = «-(1МГ’'ЧИ/?-(|М)’л(г>-)).
Для гиперболической орбиты (/г>0)
• Дата поступления
10 мая 1963 г.
ЛИТЕРАТУРА
1. II. Е. Э л ь я с б е р г, В. Д. Ястребов. Сб. «Искусственные спуТнпки Зем-
ли», вып. 4, Изд-во АН СССР, 1960, стр. 18.
2. Т. М. Энеев, А. К. Платонов, Р. К. Казаков а. Сб. «Искусственные
спутники Земли», вып. 4, Изд-во АН СССР, 1960, стр. 43.
3. И. Ш а п и р о. Расчет траекторий баллистических снарядов по данным радио-
локационных наблюдений. ИЛ, М., 1961.
4. М. Ф. Субботин. Курс небесной механики, т. II. Гостохиздат, 1949.
5. Ю. В. Л и н н и к. Метод наименьших квадратов и основы теории обработки па
блюдений. Физматгиз, 1958.
6. В. С. Пугачев. Теория случайных функций. Физматгиз, 1962.
7. Л. В. Канторович. Докл. АН СССР, 48, 7, 483, 1945.
8. А. С. Хаусхолдер. Основы численного анализа. ИЛ, М., 1956.
9. М. К. Г а в у р и н. Изв. высш, учебн. завед., вып. 5, Изд-во Мин-ва высшего об-
разования СССР, 1958, стр. 18.