CC BY
69
УДК 521.1:531.1
Вестник СПбГУ. Сер. 10, 2004,
В. С. Королев
ОПРЕДЕЛЕНИЕ ДВИЖЕНИЯ НАВИГАЦИОННЫХ СПУТНИКОВ С УЧЕТОМ ВОЗМУЩЕНИЙ*)
Прогнозирование движения системы навигационных спутников Земли и построение аналитическими или численными методами [1-4] траекторий движения космических аппаратов (КА) в гравитационном поле с учетом основных действующих возмущающих сил можно выполнять с использованием декартовой прямоугольной системы координат (как это принято в российской спутниковой системе ГЛОНАСС) или оску-лирующих кеплеровых элементов орбиты (в американской системе NAVSTAR GPS). Это позволяет по начальным данным, вычисленным заранее для всех спутников системы и передаваемым на аппаратуру пользователя в процессе наблюдения, получить с достаточно высокой точностью таблицу значений в координатах или элементах орбиты (эфемерид) для узловых моментов времени рабочего промежутка. Пользователь навигационной системы должен иметь возможность на основе этой информации определить движение и уточнить положение спутников в любой момент времени внутри заданного промежутка, а также расстояния до некоторых из них [5], наблюдаемых в данный период. Для этого дополнительно вычисляются параметры движения выделенных спутников с учетом основных возмущающих сил численным интегрированием или аппроксимацией на малые промежутки времени [7]. Полученная информация используется в аппаратуре пользователя для нахождения его координат текущего положения на поверхности Земли или в пространстве.
Движение КА в центральном гравитационном поле в декартовой инерциальной системе координат в рамках невозмущенной задачи двух тел описывают уравнения
^Г = -~JX> гдех = (жьж2,жз), т - 1*1» (!)
или с учетом обозначения для силовой функции центрального поля
W = tr =sradt7°' "о = %(*) = £. (2)
В уравнениях (1) и (2) ц = ае2 - гравитационный параметр центрального поля. Уравнения движения для возмущенной задачи двух тел с учетом малых возмущений следующие:
Ui = Ui(t,x), P = P(t,x), а = (01,02,03),
где aj - проекции вектора ускорения КА, вызванного действием возмущающих сил, среди которых можно выделить потенциальные (U) и непотенциальные (Р) силы.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 02-01-01039). © В. С. Королев, 2004
Уравнения (3) можно переписать в виде системы обыкновенных дифференциальных уравнений первого порядка
ШХ{ = Щ' + = 0i' i = 1'2'3- (4)
Решение однородных уравнений (1) невозмущенного движения хорошо известно и выражается через кеплеровы элементы орбиты
Xi = fi(t,e), e = (ej,...,e{j), ej = const. (5)
Решение уравнений возмущенного движения (3) или (4) можно получить численным интегрированием или аппроксимацией. Оскулирующие элементы орбиты е' возмущенного движения позволяют вычислить декартовы координаты х[ по формулам (5) невозмущенного движения через е = e'(i) для произвольного момента времени t при условии
с=1 J
В малой окрестности момента времени Ь имеем
х^ + А^ = + Аг) + ^сц(Аг)2. (7)
Для введения оскулирующих элементов е'(£) можно использовать^замену переменных
о
= Л(*,е), = (8)
Однако прямая подстановка соотношений (8) в уравнения движения (4) приводит к сложным вычислениям. Вместо этого выбирается другой подход с использованием интегралов невозмущенных уравнений движения
у(х', = 0 или <р(х, у,е',£) = 0. (9)
Дифференцируя (9) по времени и сравнивая соответствующие выражения, имеем
0-1 %— 1
Если в качестве первых интегралов движения (9) выбираются интегралы площадей и энергии, то из уравнений (10) можно получить уравнения Эйлера для нахождения оскулирующих кеплеровых элементов е' = (а, е, г, О, ш, Мо) возмущенного движения
^ = 2с? (esm.es' +рг~1Т'), аг
йе
— =рвшвЗ'+р(со*е + со8Е)Т1, аг
(^-=гсоъи\А/г', ~ = гвтисовесгИ^', (11)
аЬ ш
dfjj _lr ч . .dCl
— = e [—pcosflS' + (r + p) sin вТ'] - cos i—,
dM0
(1 - e2)~^e-1 [(pcos0 - 2er)S' - (r + p) sin0Г'],
dt
где правые части системы уравнений .(11) определяются текущими значениями элементов e(í), истинной аномалйи в, аргумента широты и = эксцентрической аномалии Е и проекциями возмущающих ускорений 61,62, 63 на оси орбитальной системы координат (они направлены по радиус-вектору, по трансверсали и по нормали к плоскости орбиты) в следующем виде:
61 = oi (cos и cos П — sin it sin íí cos i)+
+«2 (cos u sin íí + sin и cos íí cos г) + аз sin и sin i,
62 = ai(— sin u cosí) — cos u sin fieos г) + . (12)
+03 (- sin tí sin íí + eos и cos íí cos i) + as cos и sin i, (
Ьз =-oi sin П sin г — ад cos íí sin i + 03 eos г,
5' = as-1p~Hi, T' = ae-1p~¿62) W" = aTVW -
Когда для системы (11) получено решение, можно определить среднее движение и среднюю аномалию М
t
п ~ М = М0 + Jпdt, (13)
to
a затем вычислить декартовы координаты КА по формулам эллиптического движения через кеплеровы элементы орбиты. Последовательно определяются эксцентрическая аномалия Е, истинная аномалия в и полярные координаты г, и
Е = М + esin Е, ^ = y^^e^f' ! (
г = а(1-е2)(1 + еШв)~1, и = в + ш, а затем и декартовы координаты спутника '
х = г (cos и eos íí — sin и sin íí cos i),
у = r(cos u sin íí + sin и cos íí cos г), z — г sin и sin г.
Если возмущения учитываются только от потенциальных сил," т.е.
9R , „ „
то уравнения Эйлера (11) можно представить таким образом:
da 2 dR
dt па дМ0'
de 1-е2 dR v/í~z~e2 dR
FtR
dt ena2 8Mq ena2 du' ctg i OR cosec i dR
dt naV 1-е2 дш na2 у/ 1-е2 9ÍÍ' dfi _ cosec г сШ di" na2Vl
áw _ y/1 - e2 dR _ ctg i dR dt ena2 de na2y/l — e2 di '
dM0 _ 2 dR 1-е2 dR dt na da ena2 de Систему уравнений (15) или эквивалентные им называют уравнениями Лагранжа. Здесь R - пертурбационная функция, выраженная через кеплеровы орбитальные элементы. Производные от пертурбационной функции R по кеплеровым элементам имеют вид
da
= -^L=[-Pcos0S' + (p-+ r) sin9T'], oe vi — e
^^ = na2 y/l — e2r sin и W', (16)
oí
QJ^ _
дтг = na2 y/l — e2r (eos iT' — sinicosuW'),
C/lí
^ - па2у/Г^ё2гТ', дш
ЯП
= na5 (e sind S' +РГ-1 Т'),
дМ0
где правые части содержат проекции ускорений (12).
Можно выразить проекции возмущающих ускорений на оси инерциальной декартовой системы координат через указанные проекции (12) на оси орбитальной системы координат
ai = £>i (cos и eos fí — sin и sin fi cos г)+
+&2 (— sin и cos íí — cos и sin ü cos i) + Ьз sin Г2 sin г,
о2 = bi (cos и sin fí + sin и cos Ü cos i)+ (17)
sin и sin fl — eos и cos O cos г) — 63 cos Q sin i,
a3 — bi sin и sin г + £>2 cos it sin г + Ьз cos г.
Линеаризованные оценки погрешности в определении параметров (в том числе орбитальных кеплеровых элементов, а также декартовых координат в инерциаль-• ной или гринвичской системе) при наличии поправок в значениях вектора скорости V = vo + Avo, где Avo = (Avr,Avg,Av^), по заданным начальным значениям параметров орбиты в момент ¿о на орбите для произвольного близкого момента времени
t = to + At позволяют получить поправки к значениям г, $ в проекции на оси орбитальной системы координат Дг = (Ri,R.2,0) и Ав — (Qi,Q2,0), где
Ri = — (-3íe sin в0 sin 0 +. — (2(г sin 0О - г0 sin 0) + - sin(0 - 0О))),
р \ 36 в /
k2 = - 3íe sin 0 + — (2г - —ecos0 - г0(1 + —) cos(0 - 0О))),
Эти поправки можно использовать для вычисления изменений оскулирующих элементов заданного начального движения на малых промежутках времени.
Учет неизохронности вариации параметров, когда величины г, 0 будут выбираться в момент времени tK, а вариации нужно вычислить для момента t = tK 4- At, дает дополнительные поправки ARtAt, AQtAt, где
ARt=^sm9, = (19)
y/9 г2
Линейные вариации декартовых координат получаются с учетом связи орбитальных и декартовых инерциальных координат:
(Дхг,Дж2,Дх3)т =BA(Avr>Ave,Avc)T- (20)
В векторно-матричном уравнении £20) матрицы А = (A¿j) и В = (By) содержат элементы
Aii = Ri, Ai2 = #2, Ли = 0, A2i = Qu A22 = Qt, \ ^23 = o, ¿31 = 0, Л32 = 0, ¿33 = sin(0- 0O),
B\i — COS U COS U) — sin и sin и sin i, В12 = — sin и coso; + cos и sin и cos i,
1 £?i3 = sin и sin i, B21 = cos tí sin w + sin и cos üj cos г,
B22 = — sin tí sin из — eos tí eos ш eos i y B2z = — eos w sin г,
Bzi = sin и sin i, J?32 = eos «sin г, Б33 = cost.
Зависимость поправок в декартовых координатах от поправок в элементах орбиты (или некоторых функций от элементов) выражается с точностью до членов первого порядка малости
= (21)
Здесь Де* рассматриваются как погрешности начальных значений элементов оскули-рующей орбиты на момент времени to > учет которых позволяет вычислить более точно возмущения. Формулы для производных по элементам можно получить явно аналитически или численно определить коэффициенты разложения.
Вместе с тем вариации координат Sxi можно выразить через вариации параметров Suk,Srk,Sik,Sü, которые определяются прогнозируемыми в системе GPS на момент to амплитудами поправок , зависящих от второй зональной гармоники, и углом Ф = в + ш по формулам
Silk = Сц sin 2Ф + С\2 COS 2Ф,
бгк = С21 sin 2Ф + с22 cos 2Ф, (22)
6ik = С31 sin 2Ф + С32 cos 2Ф,
\
Гк = а(1 -ecosEk) + и* = Ф + Suk,
h = ¿o + Sik, íífc = fío + SÜ, (23)
^ x'k = rk eos uk, y'k= rk sin uk, xk= x'k eos fijt - у к sin Ük eos ik, (24)
у к = x'k sin ük + У к COS ük COS ¿fe, zk = г/fe sin ¿fe, ^
в которых индекс к указывает на значения параметров в момент í* = t — t0. Аргумент широты Ф определяется средней аномалией М* = Mq + ntk для этого момента времени последовательным решением уравнений
Ек-еsin Ек = Мк, tg вк = \Л -e2sinЕк(cos Ек - е)-1, Ф* = + и.
Выпишем подробнее систему условных уравнений
Sxk = duSrk + di2¿Uk + disóik + duSük,
Ьу к = d2iórk + d22^Uk + dizSik + chióük, (25)
6zk = dziSrk + dz-zSuk + dzz$ik + dziSQk*
Коэффициенты dy определяются производными по параметрам при условии 6ик = 0,6гк = 0, 5ik = 0,5Ü = 0.
. Для построения более точного решения можно учитывать возмущения следующего порядка малости от сжатия Земли. Это прежде всего вековые возмущения величин М, oj, fl, которые имеют вид
^М = аеа"§+£>ь = D2, jfü = D3, (26)
где правые части зависят от второй зональной гармоники, а также от а, е, i. Возмущения второго порядка получаются после интегрирования уравнений
drx »жл dDix . dDíx , dDlt.
_(<W = —¿a+ — 5е + ж8г, (27)
-{S2ü) = —óa+—óe+-Br6i,
в которых Sa,Se,Si - возмущения из первого приближения. В результате возмущения второго порядка имеют множителем вторую зональную гармонику и пропорциональны возмущениям первого порядка для элементов 6а, Se, Si.
Размножение эфемеридной информации для потребителя с помощью интерполяции по вычисленным значениям координат для равномерного распределения узлов tk = TQ + kh на интервале po,7i], tn = to + nh = Т0 + Н = Tit при поступлении новой информации для следующего промежутка времени может приводить к появлению скачков или негладкости в аппроксимации координат как функций времени.
Устранить особенности йа границах промежутков можно различными вариантами сглаживания интерполирующих функций [6, 7]. Отметим две основные возможности: гладкие восполнения или сплайн-интерполяция со'сглаживанием.
Процедура построения гладких восполнений более проста, но достаточно эффективна для сглаживания аппроксимирующих полиномов, которые построены на примыкающих (если Т\ = Т2) или пересекающихся (если Г2 < Т{) интервалах времени ßo,Ti] и рг, Тз]. Обозначим Pi(t) полином степени п на первом промежутке, а P2(t) - полином степени тп на втором промежутке. Построим полином Q(t) степени р не выше п+тп+1, который удовлетворяет условиям
Q(T1) = P1(Tl), Q(T2) = Р2(Г2). (28)
При этом также совпадают соответствующие производные этих полиномов к порядка
dkQ. dkP1{ '
•¿¿¡Г\t=Ti = « = • • • fii < n,
^г1*=т3 = k = '' • »mi» mi - ^ (29)
Сглаживающий полином Q(t) степени p = п\ +mi + l однозначно определяется условиями (28), (29) и значениями интерполируемой функции в узлах. Он заменяет один или оба полинома P(t) на соответствующих участках и обладает хорошими аппроксимирующими свойствами. В случае ni — mi =2, что соответствует гладкости кубических сплайнов, получим полином степени р = 5, для построения которого не требуется решать линейную алгебраическую систему.
Построение сглаживающих сплайн-функций q(x) аналогично процедуре построения интерполирующих сплайн-функций д(х), которые на каждом отрезке где
к — 1,... ,п, являются кубическим многочленом вида
з
9(х) = ai (Хк ~ = i = 1
= bk-i(xk - х)3 + Ьк(х - хк-1)3 + dk-i(xk -х)+ dk(x - Sfc-i). (30)
Коэффициенты, b функций д(х) находятся из решения линейной системы
Ab — Hf v (31)
через значения интерполируемой функции / в узлах хк. Для построения сглаживающих функций q(x) из условия минимального среднеквадратичного отклонения от значений функции / в узлах при минимальном изгибании функции
п '
V = f [q'fdx + ]Г Sk[q(xk) - f(xk)]2 min (32)
потребуется решать линейную систему (31) с матрицей
А' = A + HD~lH*.
Здесь матрица £> имеет диагональный вид и состоит из весовых коэффициентов отклонения 8к в узлах промежутка [жо, яп]. Матрицы А размером (п — 1) х (п — 1) и Н размером (п — 1) х (п +' 1) являются трехдиагональными для кубических сплайнов и для постоянного шага Л между узлами промежутка имеют вид
А =
Ш
I 6 О
h
А
i в
о
л
к 3
О
\ О О О
я
/ I _2 1
f 0 1Л
" ь Т
^ 0 0 0 0
о \ о о
ш
3 /
о
Q
/
Предлагаемые алгоритмы могут быть реализованы в современной аппаратуре пользователя навигационной системы и повысят эффективность работы и точность вычислений.
Summary
Korolev V. S. Computation of navigation system motions with perturbations.
Analytical and numerical methods are used for computation of trajectory motion and position of navigation system with perturbations on a small time period.
Литература
1. Субботин М. Ф. Введение в теоретическую астрономию. М., 1968. 800 с.
2. Справочное руководство по небесной механике и астродинамике / Под ред. Г. Н. Дубо-шина. М., 1971. 584 с.
3. Марчук Г. И. Методы вычислительной математики. М., 1980. 535 с.
4. Емельянов Н. В. Методы составления алгоритмов и программ в задачах небесной механики. М., 1983. 128 с.
5. Коваленко А. Н., Королев В. С. Размножение эфемеридной информации потребителем в системе спутниковой навигации // Труды 26-х Чтений, посвященных К. Э. Циолковскому. Калуга, 1992. С. 92-95.
6. Коваленко А. Н., Королев В. С. Задача оптимизации траекторий с учетом ограничений на время или импульс // Управляемые системы в гравитационном поле. Управляемое движение / Под ред. В. С. Королева. СПб., 2000. С. 104-109.
7. Королев В. С. Аппроксимация функций при определении траекторий возмущенного движения навигационных спутников // Управляемые системы в гравитационном поле. Оптимальные траектории / Под ред. В. С. Королева. СПб., 2001. С. 61-69.
Статья поступила в редакцию 19 октября 2004 г.
