МЕТОД ПРИБЛИЖЕННОГО КОРРЕКТИРУЮЩЕГО ОПЕРАТОРА В ЗАДАЧАХ ВАРИАЦИОННОГО ОЦЕНИВАНИЯ ПАРАМЕТРОВ ДВИЖЕНИЯ КОСМИЧЕСКИХ АППАРАТОВ
Метод приближенного корректирующего оператора
11
УДК 629.191
В. И. МИРОНОВ, Ю. В. МИРОНОВ, Р. М. ЮСУПОВ
МЕТОД ПРИБЛИЖЕННОГО КОРРЕКТИРУЮЩЕГО ОПЕРАТОРА В ЗАДАЧАХ ВАРИАЦИОННОГО ОЦЕНИВАНИЯ ПАРАМЕТРОВ ДВИЖЕНИЯ
КОСМИЧЕСКИХ АППАРАТОВ
Рассматриваются особенности применения вариационного подхода к оцениванию параметров состояния нелинейных динамических систем в задачах навигации космических аппаратов. Краевая задача оптимального оценивания параметров движения решается на основе конструктивного метода приближенного корректирующего оператора. Приводится численный пример.
Ключевые слова: навигационное оценивание, краевая задача, метод приближенного корректирующего оператора.
Введение. В настоящее время основными при определении орбит космических аппаратов (КА) являются методы, базирующиеся на совместной обработке результатов наблюдений по полной выборке. Они позволяют успешно решать широкий круг важных и сложных прикладных задач [1—8 и др.]. В динамических задачах оценивания условий непосредственно применяются метод максимального правдоподобия (ММП) и метод наименьших квадратов. По смыслу они представляют необходимые условия оптимальности, характерные для прямых методов оптимизации.
Вместе с тем метод максимального правдоподобия может быть реализован на основе использования условий оптимальности оценок вариационного типа. Вопросы обоснования и разработки соответствующей вариационной технологии рассматривались в работах авторов [9, 10] применительно к оцениванию состояния нелинейных динамических систем.
Настоящая статья посвящена вариационному подходу к решению задач навигационного оценивания параметров движения КА по данным измерений установленной на его борту навигационной аппаратуры потребителя (НАП), работающей по сигналам спутниковой радионавигационной системы. При этом в качестве данных измерений рассматриваются оценки векторов параметров текущего состояния КА, получаемые НАП на заданном мерном интервале времени, которые подлежат совместной обработке. С целью сокращения объема вычислений при решении краевых задач вариационного оценивания применяется конструктивный метод приближенного корректирующего оператора, предложенный в работе [11].
В статье определяются и конкретизируются необходимые условия оптимальности оценок вариационного типа применительно к модели движения КА в нормальном гравитационном поле, учитывающем полярное сжатие Земли. Приводятся результаты численных расчетов.
Постановка задачи. Рассмотрим задачу оценивания параметров движения динамического объекта, которая заключается в наилучшем определении п-мерного вектора его исходного состояния x0 на начальный момент времени t = t0 по результатам измерений, проведенных в N точках ti , заданных на интервале измерений τ = T − t0 .
Задача. Пусть динамика объекта описывается векторным дифференциальным уравнением
x = ϕ(x,t), x(t0 ) = x0 , t ∈[t0 ,T ] ,
измеряется m-мерный вектор
ψ(t) = ψ[x(t)] .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
12 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
Измеренное значение вектора ψ в момент ti обозначим как y(ti ) = yi и представим модель измерений в виде
y(ti ) = ψ[x(ti )] + δi ,
i = 1(1)N ; ti ∈[t0,T ] .
Здесь δi — m-мерный вектор случайных ошибок измерений, стохастическое изменение кото-
рого зададим некоторым многомерным непрерывным дифференцируемым распределением
f (δi , αi ) с параметрами αi , отличающимся в общем случае от нормального распределения.
Требуется найти такую оценку вектора x0 , которая обеспечивает минимальное значение
функционала:
N
I = ∑ ρi {y(ti ), ψ[x(ti )], αi } , i=1
где
ρi = ln fi {y(ti ) − ψ[x(ti ), αi ]} ; i = 1(1)N .
Функции ϕ(x,t) и ψ[x(ti )] будем считать однозначными, ограниченными, непрерыв-
ными и дифференцируемыми по всем аргументам во всей области их определения. Нетрудно
видеть, что функционал I есть не что иное, как логарифмическая функция правдоподобия.
Для решения поставленной задачи в работе [9] были получены условия оптимальности
оценок вариационного типа, которые заключаются в следующем: оптимальная оценка векто-
ра x0 и соответствующая ей оптимальная траектория являются решением краевой задачи для
следующей системы дифференциальных уравнений:
x = ϕ(x, t);
λ
=
−
∂ϕт ∂x
λ
при граничных условиях
λ(ti+
)
=
λ (ti−
)
+
∂ρ ∂x
[yi
,
ψ(xi
),
ti
]
;
λ(t0 ) = 0 ; λ(T ) = 0 ;
i = 1(1)N.
Таким образом, рассматриваемую задачу можно интерпретировать как двухточечную
краевую с промежуточными ограничениями на сопряженный вектор λ(t) .
Приведенные выше условия оптимального оценивания можно конкретизировать приме-
нительно к заданному виду распределения вектора случайных ошибок измерений. Так, если
для вектора δi принимается нормальным распределение N (0, Kδi ) с нулевым вектором ма-
тематического ожидания и корреляционной матрицей Kδi , что, как правило, имеет место на практике, решение задачи оптимального оценивания сводится к решению краевой задачи
x = ϕ(x,t);
λ
=
−
∂ϕТ ∂x
λ
;
λ(t0 ) = 0 ; λ(T ) = 0 ;
λ(ti+
)
=
λ(ti−
)
+
∂ψТ (ti ∂xi
)
Wi
{yi
−
ψ[x(ti
)]};
i = 1(1)N .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
Метод приближенного корректирующего оператора
13
Согласно этим условиям, для получения оптимальной оценки вектора x0 необходимо решить краевое уравнение
λ(x0 ,T ) = 0,
заданное неявно на процедурах интегрирования сопряженных систем дифференциальных уравнений. Для этого можно применить известные численные методы поиска корней нелинейных уравнений, например метод Ньютона, его модификации и другие.
С точки зрения сокращения объема вычислений и повышения оперативности расчетов особый интерес представляет применение конструктивных методов решения краевых задач такого рода, в частности метода приближенного корректирующего оператора [11]. Согласно этому методу, для решения исходной краевой задачи применяется вычислительная процедура
∑q k +1
=
M
⎡ ⎢−
k
λ(qi ,T )⎤⎥ ;
⎢⎣ i=0
⎦⎥
k
= 1, 2,3,... ,
где M [•] — приближенный корректирующий оператор (ПКО), представляет собой алгоритм
приближенного решения приближенного краевого уравнения
λ(q,T ) = 0 .
Конструктивной особенностью метода ПКО является возможность синтеза семейства алгоритмов численного решения краевой задачи, различающихся сложностью и скоростью
сходимости в зависимости от выбора приближенной модели λ(q,T ) и принятого способа
приближенного решения приведенного уравнения. Условия сходимости данного метода определены в [8] на основе принципа сжимающих отображений.
С помощью метода ПКО могут быть выбраны приближенные модели множеством различных способов с учетом специфики исходных зависимостей и условий решаемой задачи. Для этого могут применяться как формальные приемы упрощения исходных моделей, так и методы их аппроксимации. Во всех случаях необходимо стремиться к тому, чтобы приближенный алгоритм решения задачи был сравнительно простым и обеспечивалась достаточно быстрая сходимость итерационного процесса.
Важной особенностью метода ПКО является то обстоятельство, что значение неизвестного вектора q уточняется на каждом итерационном шаге путем однократного интегрирова-
ния дифференциальных уравнений движения объекта. Этим обеспечивается высокая экономичность вычислений. Ниже это будет показано на примере.
Значение оператора M [•] может изменяться или уточняться в ходе вычислительного
процесса по его текущим результатам на каждой итерации или через несколько итераций. В этом случае можно говорить о комбинированных вариантах использования метода ПКО.
В целом метод ПКО выражает достаточно общую идеологию конструирования алгоритмов решения краевых задач и нелинейных уравнений.
Определение параметров орбиты космического аппарата по данным спутниковой навигации. Рассмотрим особенности применения вариационного метода максимального правдоподобия на примере решения задачи статистического оценивания параметров движения КА по результатам текущих навигационных оценок, полученных его бортовой навигационной аппаратурой, работающей по сигналам спутниковой навигационной системы. Движение КА будем рассматривать в нецентральном гравитационном поле Земли. Соответствующие уравнения движения в абсолютной геоцентрической системе координат имеют вид [12]
x = vx ; y = vy ; z = vz ;
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
14 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
vx
=
⎛ ⎜⎝
−
µ r3
+
p
⎞ ⎟⎠
x
;
vy
=
⎛ ⎜⎝
−
µ r3
+
p
⎞ ⎠⎟
y;
vz
=
⎛ ⎝⎜
−
µ r3
+
p+
∆p
⎞ ⎟⎠
z
;
p
=
3 2
µ r3
J 20
⎝⎜⎜⎛1 − 5
z2 r2
⎞⎛ ⎟⎟⎠ ⎝⎜
Re r
⎞2 ⎟⎠
+
15 8
µ r3
J 40
⎛ ⎝⎜⎜
−1
+
14
z2 r2
−
21
z4 r4
⎞⎛ ⎟⎠⎟ ⎝⎜
Re r
⎞4 ⎠⎟
;
∆p
=
3
µ r3
J 20
⎛ ⎝⎜
Re r
⎞2 ⎠⎟
+
15 8
µ r3
J 40
⎛ 28 ⎝⎜⎜ 3
z2 r2
−
4
⎞ ⎠⎟⎟
⎛ ⎝⎜
Re r
⎞4 ⎠⎟
;
r = x2 + y2 + z2 , J20 = −1 082 627 ⋅10−9 ; J40 = 2371⋅10−9 ;
где r = [x, y, z]T и v = [vx , vy , vz ]T — векторы координат и скорости движения КА соответст-
венно; Re = 6378,136 км — экваториальный радиус Земли; µ=398 600,44 км3/с2 — постоянная притяжения Земли.
Проведем прямые полные дискретные измерения элементов векторов r и v, так что модель измерений принимает вид
y1i = r(ti ) + δri ; y2i = v(ti ) + δvi ; i = 1(1)N ,
где δri , δvi — ошибки оценок НАП, полученных в моменты времени ti .
Составим соответствующую сопряженную систему дифференциальных уравнений. Анализ показывает, что при обработке навигационных данных на ограниченных мерных интервалах в модели сопряженной системы можно ограничиться членами, соответствующими ньютоновской части гравитационного поля. В этом случае будем иметь следующую систему сопряженных дифференциальных уравнений:
λr
=
µ0 r3
[λ v
−
3 r2
(rrТ
)λv
];
λv
= −λr ;
λ = [λr , λv ]Т = [λx ,λ y ,λz ,λvx ,λvy ,λvz ]Т . Теперь в соответствии с приведенными выше вариационными условиями оптимально-
сти для оценивания вектора начального состояния КА q = [r, v]T необходимо решить двухто-
чечную краевую задачу для приведенных систем уравнений движения и сопряженных уравнений с учетом следующих граничных условий:
λ(t0 ) = 0 ; λ(T ) = 0 ;
λr (ti+ ) = λr (ti− ) + Kδ−r1[y1 − r(ti )] ; λv (ti+ ) = λv (ti− ) + Kδ−v1[y2 − v(ti )] ; i = 1(1)N ,
где Kδr и Kδv — корреляционные матрицы ошибок измерений векторов координат и скорости движения КА соответственно.
Таким образом, решение задачи сводится к поиску корней краевого уравнения λ(q,T ) = 0.
Применение метода Ньютона приводит к следующему итерационному алгоритму:
q k +1
=
qk
−
⎡ ⎢⎣
∂λ
(q, ∂q
T
)
⎤ ⎥⎦
−1 k
λ
(qk
,T
);
k = 0,1, 2,...
Проиллюстрируем возможности конструктивного метода приближенного корректирую-
щего оператора для конкретизации вида оператора M , рассмотрев приближенное решение
исходной задачи, соответствующее условиям движения спутника по орбитам, близким к кру-
говым. Для этого используем модель динамики в однородном центральном поле и модель
прямых эквивалентных измерений полного вектора состояния. В этом случае краевая задача
оптимального оценивания принимает вид
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
Метод приближенного корректирующего оператора
r = V;
v = −ω2r;
λr = ω2λv ;
λv = −λr ;
λ
(t
+ j
)
=
λ
(t
− j
)
+
K
−1 xj
⎣⎡x
j
−
x(t j
)⎤⎦
;
j = 1(1)N1 ;
λ(t0 ) = λ(T ) = 0 ,
где ω — угловая скорость орбитального движения спутника в начальный момент t0 :
15
ω = ω(x0 ) = µ0r0−3 ; r0 = x02 + y02 + z02 .
Решение этой краевой задачи сводится к определению корня x0 нелинейного уравнения
∑N1
V(T ,t
j
)K
−1 δj
⎡⎣x
j
−
U(t
j ,t0 )x0
⎤⎦
=
0
,
j =1
где
U(t
j
,
t0
)
=
⎡ ⎢
cos ω(t j − t0 )E
⎣⎢−ωsin ω(t j − t0 )E
ω−1
sin
ω(t
j
−
t0
)E⎤ ⎥
;
cos ω(t j − t0 )E ⎥⎦
V(T
,
t
j
)
=
⎡ cos ω(T − t ⎢⎣⎢−ω−1 sin ω(T
j )E − t j )E
ωsin
ω(T
−
t
j
)E⎤ ⎥
;
cos ω(T − t j )E ⎥⎦
E — единичная (3×3)-матрица. Нелинейность этого уравнения определяется зависимостью фундаментальных матриц U и V от параметра ω, который, в свою очередь, зависит от неиз-
вестного вектора x0 . Проанализировав структуру приведенных выражений, приходим к заключению, что для
решения исходной краевой задачи вариационного оценивания параметров движения КА мож-
но применить итерационную вычислительную процедуру метода ПКО с корректирующим
оператором
∑M k
=
⎡ N1 ⎢ V(T ,t j ⎢⎣ j=1
/
x0k
)K
−x1j U(t
j
,
t0
/
x0k
⎤−1 )⎥ ⎥⎦
.
Приведем некоторые результаты вычислений. Расчеты проводились для спутника,
находящегося на орбите с высотой h = 900 км и эксцентриситетом 0,003. С помощью датчика случайных величин по нормальному закону распределения на мерном интервале T = 100 с моделировалась с шагом ∆t = 1 с статистическая выборка прямых измерений вектора текущего состояния КА. При этом предельные ошибки измерений задавались значениями 100 м по эле-
ментам вектора координат и 1 м/с — по элементам вектора скорости. Некоторые результаты
расчетов приведены в табл. 1 и 2.
Результаты оптимального оценивания
Таблица 1
Величина
Точное значение Начальное приближение Оптимальная оценка Ошибка оценивания
x0 , км
0 50 –0,006 –0,006
y0 , км
–7249,926 –7199,926 –7249,928
–0,002
z0 , км
0 50 –0,007 –0,006
vx0 , км/с
0,904 95 0,954 95 0,904 97 –0,000 01
vy0 , км/с
0,005 75 0,055 75 0,005 80 0,000 05
vz0 , км/с
7,370 23 7,820 23 7,370 30 0,000 06
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
16 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
Номер итерации
0 1 2 3
Сходимость вычислительного процесса
Ошибки оценивания
δ x0 , м
50 000 –7,6 –6,3 –6,3
δ y0 , м
50 000 –2,8 –1,9 –1,9
δ z0 , м
50 000 –7,9 –6,5 –6,5
δVx0 , м/с
50 0,05 0,02 0,02
δVy0 , м/с
50 0,04 0,05 0,05
Таблица 2
δVz0 , м/с 50
–0,01 0,06 0,06
В табл. 1 даны точные значения параметров начального фазового состояния КА в абсолютной геоцентрической системе координат, принятое начальное приближение элементов уточняемого вектора, полученные в результате вариационной обработки измерений оптимальные оценки, а также характеристики точности оценивания. В табл. 2 представлены значения ошибок оценивания по итерациям. Приведенные в таблицах данные свидетельствуют о высокой точности и скорости сходимости вычислительного процесса. Расчеты показали, что применение метода ПКО в принятых условиях позволяет сократить объем вычислений примерно в четыре-пять раз.
В заключение отметим, что предлагаемые методические средства могут быть использованы при разработке и модернизации алгоритмов навигационного оценивания КА. Алгоритмы вариационного типа могут применяться самостоятельно или параллельно с традиционными алгоритмами прямого оптимального оценивания для контроля правильности вычислений и обеспечения надежности расчетов. Они могут также найти применение при решении задач тестирования приближенных алгоритмов оценивания и обоснования эффективного состава и программы измерений.
Работа выполнена при поддержке РФФИ (проект № 09-08-00259).
СПИСОК ЛИТЕРАТУРЫ
1. Аким Э. Л., Энеев Т. М. Определение параметров движения космических аппаратов по данным траекторных измерений // Космические исследования. 1963. Т. 1, № 1. C. 5—50.
2. Брандин Н. К., Разоренов Г. Н. Определение траекторий КА. М.: Машиностроение, 1978. 216 с.
3. Космические траекторные измерения / Под ред. П. А. Агаджанова, В. Е. Дулевича, А. А. Коростелева. М.: Сов. радио, 1969. 504 с.
4. Линник Ю. В. Метод наименьших квадратов и основы теории обработки наблюдений. М.: Физматгиз, 1958. 350 с.
5. Основы теории полета космических аппаратов / Под ред. Г. С. Нариманова и М. К. Тихонравова. М.: Машиностроение, 1972. 608 с.
6. Шебшаевич В. С., Дмитриев П. П., Иванцевич Н. В. и др. Сетевые спутниковые радионавигационные системы. М.: Радио и связь, 1993. 408 с.
7. Статистические методы обработки результатов наблюдений / Под ред. Р. М. Юсупова. М.: МО СССР, 1984. 563 с.
8. Эльясберг П. Е. Определение движения по результатам измерений. М.: Наука, 1976. 416 с.
9. Миронов В. И., Миронов Ю. В., Юсупов Р. М. Вариационное оценивание состояния нелинейной динамической системы по критерию максимального правдоподобия // Мехатроника, автоматизация, управление. 2009. № 11. С. 2—6.
10. Миронов В. И., Миронов Ю. В., Юсупов Р. М. Метод наименьших квадратов в задачах вариационного оценивания состояния нелинейных динамических систем // Информационно-управляющие системы. 2009. № 6. С. 2—6.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
Рекуррентное систематическое помехозащитное преобразование кодов
17
11. Миронов В. И., Миронов Ю. В., Юсупов Р. М. Синтез итерационных алгоритмов решения краевых задач и нелинейных уравнений // Изв. вузов. Приборостроение. 2010. Т. 53, № 1. С. 9—14.
12. Мамон П. А., Половников В. И., Слезкинский С. К. Баллистическое обеспечение космических полетов. Л.: ВИКИ им. А. Ф. Можайского, 1990. 622 с.
Вячеслав Иванович Миронов Юрий Вячеславович Миронов Рафаэль Мидхатович Юсупов
Сведения об авторах — д-р техн. наук, профессор; Санкт-Петербургский институт информа-
тики и автоматизации РАН; E-mail: mironuv@yandex.ru — д-р техн. наук, доцент; Санкт-Петербургский институт информатики и
автоматизации РАН; E-mail: mironuv@yandex.ru — д-р техн. наук, профессор; Санкт-Петербургский институт информа-
тики и автоматизации РАН; директор; E-mail: spiiran@iias.spb.su
Рекомендована институтом
Поступила в редакцию 07.07.10 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
11
УДК 629.191
В. И. МИРОНОВ, Ю. В. МИРОНОВ, Р. М. ЮСУПОВ
МЕТОД ПРИБЛИЖЕННОГО КОРРЕКТИРУЮЩЕГО ОПЕРАТОРА В ЗАДАЧАХ ВАРИАЦИОННОГО ОЦЕНИВАНИЯ ПАРАМЕТРОВ ДВИЖЕНИЯ
КОСМИЧЕСКИХ АППАРАТОВ
Рассматриваются особенности применения вариационного подхода к оцениванию параметров состояния нелинейных динамических систем в задачах навигации космических аппаратов. Краевая задача оптимального оценивания параметров движения решается на основе конструктивного метода приближенного корректирующего оператора. Приводится численный пример.
Ключевые слова: навигационное оценивание, краевая задача, метод приближенного корректирующего оператора.
Введение. В настоящее время основными при определении орбит космических аппаратов (КА) являются методы, базирующиеся на совместной обработке результатов наблюдений по полной выборке. Они позволяют успешно решать широкий круг важных и сложных прикладных задач [1—8 и др.]. В динамических задачах оценивания условий непосредственно применяются метод максимального правдоподобия (ММП) и метод наименьших квадратов. По смыслу они представляют необходимые условия оптимальности, характерные для прямых методов оптимизации.
Вместе с тем метод максимального правдоподобия может быть реализован на основе использования условий оптимальности оценок вариационного типа. Вопросы обоснования и разработки соответствующей вариационной технологии рассматривались в работах авторов [9, 10] применительно к оцениванию состояния нелинейных динамических систем.
Настоящая статья посвящена вариационному подходу к решению задач навигационного оценивания параметров движения КА по данным измерений установленной на его борту навигационной аппаратуры потребителя (НАП), работающей по сигналам спутниковой радионавигационной системы. При этом в качестве данных измерений рассматриваются оценки векторов параметров текущего состояния КА, получаемые НАП на заданном мерном интервале времени, которые подлежат совместной обработке. С целью сокращения объема вычислений при решении краевых задач вариационного оценивания применяется конструктивный метод приближенного корректирующего оператора, предложенный в работе [11].
В статье определяются и конкретизируются необходимые условия оптимальности оценок вариационного типа применительно к модели движения КА в нормальном гравитационном поле, учитывающем полярное сжатие Земли. Приводятся результаты численных расчетов.
Постановка задачи. Рассмотрим задачу оценивания параметров движения динамического объекта, которая заключается в наилучшем определении п-мерного вектора его исходного состояния x0 на начальный момент времени t = t0 по результатам измерений, проведенных в N точках ti , заданных на интервале измерений τ = T − t0 .
Задача. Пусть динамика объекта описывается векторным дифференциальным уравнением
x = ϕ(x,t), x(t0 ) = x0 , t ∈[t0 ,T ] ,
измеряется m-мерный вектор
ψ(t) = ψ[x(t)] .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
12 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
Измеренное значение вектора ψ в момент ti обозначим как y(ti ) = yi и представим модель измерений в виде
y(ti ) = ψ[x(ti )] + δi ,
i = 1(1)N ; ti ∈[t0,T ] .
Здесь δi — m-мерный вектор случайных ошибок измерений, стохастическое изменение кото-
рого зададим некоторым многомерным непрерывным дифференцируемым распределением
f (δi , αi ) с параметрами αi , отличающимся в общем случае от нормального распределения.
Требуется найти такую оценку вектора x0 , которая обеспечивает минимальное значение
функционала:
N
I = ∑ ρi {y(ti ), ψ[x(ti )], αi } , i=1
где
ρi = ln fi {y(ti ) − ψ[x(ti ), αi ]} ; i = 1(1)N .
Функции ϕ(x,t) и ψ[x(ti )] будем считать однозначными, ограниченными, непрерыв-
ными и дифференцируемыми по всем аргументам во всей области их определения. Нетрудно
видеть, что функционал I есть не что иное, как логарифмическая функция правдоподобия.
Для решения поставленной задачи в работе [9] были получены условия оптимальности
оценок вариационного типа, которые заключаются в следующем: оптимальная оценка векто-
ра x0 и соответствующая ей оптимальная траектория являются решением краевой задачи для
следующей системы дифференциальных уравнений:
x = ϕ(x, t);
λ
=
−
∂ϕт ∂x
λ
при граничных условиях
λ(ti+
)
=
λ (ti−
)
+
∂ρ ∂x
[yi
,
ψ(xi
),
ti
]
;
λ(t0 ) = 0 ; λ(T ) = 0 ;
i = 1(1)N.
Таким образом, рассматриваемую задачу можно интерпретировать как двухточечную
краевую с промежуточными ограничениями на сопряженный вектор λ(t) .
Приведенные выше условия оптимального оценивания можно конкретизировать приме-
нительно к заданному виду распределения вектора случайных ошибок измерений. Так, если
для вектора δi принимается нормальным распределение N (0, Kδi ) с нулевым вектором ма-
тематического ожидания и корреляционной матрицей Kδi , что, как правило, имеет место на практике, решение задачи оптимального оценивания сводится к решению краевой задачи
x = ϕ(x,t);
λ
=
−
∂ϕТ ∂x
λ
;
λ(t0 ) = 0 ; λ(T ) = 0 ;
λ(ti+
)
=
λ(ti−
)
+
∂ψТ (ti ∂xi
)
Wi
{yi
−
ψ[x(ti
)]};
i = 1(1)N .
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
Метод приближенного корректирующего оператора
13
Согласно этим условиям, для получения оптимальной оценки вектора x0 необходимо решить краевое уравнение
λ(x0 ,T ) = 0,
заданное неявно на процедурах интегрирования сопряженных систем дифференциальных уравнений. Для этого можно применить известные численные методы поиска корней нелинейных уравнений, например метод Ньютона, его модификации и другие.
С точки зрения сокращения объема вычислений и повышения оперативности расчетов особый интерес представляет применение конструктивных методов решения краевых задач такого рода, в частности метода приближенного корректирующего оператора [11]. Согласно этому методу, для решения исходной краевой задачи применяется вычислительная процедура
∑q k +1
=
M
⎡ ⎢−
k
λ(qi ,T )⎤⎥ ;
⎢⎣ i=0
⎦⎥
k
= 1, 2,3,... ,
где M [•] — приближенный корректирующий оператор (ПКО), представляет собой алгоритм
приближенного решения приближенного краевого уравнения
λ(q,T ) = 0 .
Конструктивной особенностью метода ПКО является возможность синтеза семейства алгоритмов численного решения краевой задачи, различающихся сложностью и скоростью
сходимости в зависимости от выбора приближенной модели λ(q,T ) и принятого способа
приближенного решения приведенного уравнения. Условия сходимости данного метода определены в [8] на основе принципа сжимающих отображений.
С помощью метода ПКО могут быть выбраны приближенные модели множеством различных способов с учетом специфики исходных зависимостей и условий решаемой задачи. Для этого могут применяться как формальные приемы упрощения исходных моделей, так и методы их аппроксимации. Во всех случаях необходимо стремиться к тому, чтобы приближенный алгоритм решения задачи был сравнительно простым и обеспечивалась достаточно быстрая сходимость итерационного процесса.
Важной особенностью метода ПКО является то обстоятельство, что значение неизвестного вектора q уточняется на каждом итерационном шаге путем однократного интегрирова-
ния дифференциальных уравнений движения объекта. Этим обеспечивается высокая экономичность вычислений. Ниже это будет показано на примере.
Значение оператора M [•] может изменяться или уточняться в ходе вычислительного
процесса по его текущим результатам на каждой итерации или через несколько итераций. В этом случае можно говорить о комбинированных вариантах использования метода ПКО.
В целом метод ПКО выражает достаточно общую идеологию конструирования алгоритмов решения краевых задач и нелинейных уравнений.
Определение параметров орбиты космического аппарата по данным спутниковой навигации. Рассмотрим особенности применения вариационного метода максимального правдоподобия на примере решения задачи статистического оценивания параметров движения КА по результатам текущих навигационных оценок, полученных его бортовой навигационной аппаратурой, работающей по сигналам спутниковой навигационной системы. Движение КА будем рассматривать в нецентральном гравитационном поле Земли. Соответствующие уравнения движения в абсолютной геоцентрической системе координат имеют вид [12]
x = vx ; y = vy ; z = vz ;
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
14 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
vx
=
⎛ ⎜⎝
−
µ r3
+
p
⎞ ⎟⎠
x
;
vy
=
⎛ ⎜⎝
−
µ r3
+
p
⎞ ⎠⎟
y;
vz
=
⎛ ⎝⎜
−
µ r3
+
p+
∆p
⎞ ⎟⎠
z
;
p
=
3 2
µ r3
J 20
⎝⎜⎜⎛1 − 5
z2 r2
⎞⎛ ⎟⎟⎠ ⎝⎜
Re r
⎞2 ⎟⎠
+
15 8
µ r3
J 40
⎛ ⎝⎜⎜
−1
+
14
z2 r2
−
21
z4 r4
⎞⎛ ⎟⎠⎟ ⎝⎜
Re r
⎞4 ⎠⎟
;
∆p
=
3
µ r3
J 20
⎛ ⎝⎜
Re r
⎞2 ⎠⎟
+
15 8
µ r3
J 40
⎛ 28 ⎝⎜⎜ 3
z2 r2
−
4
⎞ ⎠⎟⎟
⎛ ⎝⎜
Re r
⎞4 ⎠⎟
;
r = x2 + y2 + z2 , J20 = −1 082 627 ⋅10−9 ; J40 = 2371⋅10−9 ;
где r = [x, y, z]T и v = [vx , vy , vz ]T — векторы координат и скорости движения КА соответст-
венно; Re = 6378,136 км — экваториальный радиус Земли; µ=398 600,44 км3/с2 — постоянная притяжения Земли.
Проведем прямые полные дискретные измерения элементов векторов r и v, так что модель измерений принимает вид
y1i = r(ti ) + δri ; y2i = v(ti ) + δvi ; i = 1(1)N ,
где δri , δvi — ошибки оценок НАП, полученных в моменты времени ti .
Составим соответствующую сопряженную систему дифференциальных уравнений. Анализ показывает, что при обработке навигационных данных на ограниченных мерных интервалах в модели сопряженной системы можно ограничиться членами, соответствующими ньютоновской части гравитационного поля. В этом случае будем иметь следующую систему сопряженных дифференциальных уравнений:
λr
=
µ0 r3
[λ v
−
3 r2
(rrТ
)λv
];
λv
= −λr ;
λ = [λr , λv ]Т = [λx ,λ y ,λz ,λvx ,λvy ,λvz ]Т . Теперь в соответствии с приведенными выше вариационными условиями оптимально-
сти для оценивания вектора начального состояния КА q = [r, v]T необходимо решить двухто-
чечную краевую задачу для приведенных систем уравнений движения и сопряженных уравнений с учетом следующих граничных условий:
λ(t0 ) = 0 ; λ(T ) = 0 ;
λr (ti+ ) = λr (ti− ) + Kδ−r1[y1 − r(ti )] ; λv (ti+ ) = λv (ti− ) + Kδ−v1[y2 − v(ti )] ; i = 1(1)N ,
где Kδr и Kδv — корреляционные матрицы ошибок измерений векторов координат и скорости движения КА соответственно.
Таким образом, решение задачи сводится к поиску корней краевого уравнения λ(q,T ) = 0.
Применение метода Ньютона приводит к следующему итерационному алгоритму:
q k +1
=
qk
−
⎡ ⎢⎣
∂λ
(q, ∂q
T
)
⎤ ⎥⎦
−1 k
λ
(qk
,T
);
k = 0,1, 2,...
Проиллюстрируем возможности конструктивного метода приближенного корректирую-
щего оператора для конкретизации вида оператора M , рассмотрев приближенное решение
исходной задачи, соответствующее условиям движения спутника по орбитам, близким к кру-
говым. Для этого используем модель динамики в однородном центральном поле и модель
прямых эквивалентных измерений полного вектора состояния. В этом случае краевая задача
оптимального оценивания принимает вид
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
Метод приближенного корректирующего оператора
r = V;
v = −ω2r;
λr = ω2λv ;
λv = −λr ;
λ
(t
+ j
)
=
λ
(t
− j
)
+
K
−1 xj
⎣⎡x
j
−
x(t j
)⎤⎦
;
j = 1(1)N1 ;
λ(t0 ) = λ(T ) = 0 ,
где ω — угловая скорость орбитального движения спутника в начальный момент t0 :
15
ω = ω(x0 ) = µ0r0−3 ; r0 = x02 + y02 + z02 .
Решение этой краевой задачи сводится к определению корня x0 нелинейного уравнения
∑N1
V(T ,t
j
)K
−1 δj
⎡⎣x
j
−
U(t
j ,t0 )x0
⎤⎦
=
0
,
j =1
где
U(t
j
,
t0
)
=
⎡ ⎢
cos ω(t j − t0 )E
⎣⎢−ωsin ω(t j − t0 )E
ω−1
sin
ω(t
j
−
t0
)E⎤ ⎥
;
cos ω(t j − t0 )E ⎥⎦
V(T
,
t
j
)
=
⎡ cos ω(T − t ⎢⎣⎢−ω−1 sin ω(T
j )E − t j )E
ωsin
ω(T
−
t
j
)E⎤ ⎥
;
cos ω(T − t j )E ⎥⎦
E — единичная (3×3)-матрица. Нелинейность этого уравнения определяется зависимостью фундаментальных матриц U и V от параметра ω, который, в свою очередь, зависит от неиз-
вестного вектора x0 . Проанализировав структуру приведенных выражений, приходим к заключению, что для
решения исходной краевой задачи вариационного оценивания параметров движения КА мож-
но применить итерационную вычислительную процедуру метода ПКО с корректирующим
оператором
∑M k
=
⎡ N1 ⎢ V(T ,t j ⎢⎣ j=1
/
x0k
)K
−x1j U(t
j
,
t0
/
x0k
⎤−1 )⎥ ⎥⎦
.
Приведем некоторые результаты вычислений. Расчеты проводились для спутника,
находящегося на орбите с высотой h = 900 км и эксцентриситетом 0,003. С помощью датчика случайных величин по нормальному закону распределения на мерном интервале T = 100 с моделировалась с шагом ∆t = 1 с статистическая выборка прямых измерений вектора текущего состояния КА. При этом предельные ошибки измерений задавались значениями 100 м по эле-
ментам вектора координат и 1 м/с — по элементам вектора скорости. Некоторые результаты
расчетов приведены в табл. 1 и 2.
Результаты оптимального оценивания
Таблица 1
Величина
Точное значение Начальное приближение Оптимальная оценка Ошибка оценивания
x0 , км
0 50 –0,006 –0,006
y0 , км
–7249,926 –7199,926 –7249,928
–0,002
z0 , км
0 50 –0,007 –0,006
vx0 , км/с
0,904 95 0,954 95 0,904 97 –0,000 01
vy0 , км/с
0,005 75 0,055 75 0,005 80 0,000 05
vz0 , км/с
7,370 23 7,820 23 7,370 30 0,000 06
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
16 В. И. Миронов, Ю. В. Миронов, Р. М. Юсупов
Номер итерации
0 1 2 3
Сходимость вычислительного процесса
Ошибки оценивания
δ x0 , м
50 000 –7,6 –6,3 –6,3
δ y0 , м
50 000 –2,8 –1,9 –1,9
δ z0 , м
50 000 –7,9 –6,5 –6,5
δVx0 , м/с
50 0,05 0,02 0,02
δVy0 , м/с
50 0,04 0,05 0,05
Таблица 2
δVz0 , м/с 50
–0,01 0,06 0,06
В табл. 1 даны точные значения параметров начального фазового состояния КА в абсолютной геоцентрической системе координат, принятое начальное приближение элементов уточняемого вектора, полученные в результате вариационной обработки измерений оптимальные оценки, а также характеристики точности оценивания. В табл. 2 представлены значения ошибок оценивания по итерациям. Приведенные в таблицах данные свидетельствуют о высокой точности и скорости сходимости вычислительного процесса. Расчеты показали, что применение метода ПКО в принятых условиях позволяет сократить объем вычислений примерно в четыре-пять раз.
В заключение отметим, что предлагаемые методические средства могут быть использованы при разработке и модернизации алгоритмов навигационного оценивания КА. Алгоритмы вариационного типа могут применяться самостоятельно или параллельно с традиционными алгоритмами прямого оптимального оценивания для контроля правильности вычислений и обеспечения надежности расчетов. Они могут также найти применение при решении задач тестирования приближенных алгоритмов оценивания и обоснования эффективного состава и программы измерений.
Работа выполнена при поддержке РФФИ (проект № 09-08-00259).
СПИСОК ЛИТЕРАТУРЫ
1. Аким Э. Л., Энеев Т. М. Определение параметров движения космических аппаратов по данным траекторных измерений // Космические исследования. 1963. Т. 1, № 1. C. 5—50.
2. Брандин Н. К., Разоренов Г. Н. Определение траекторий КА. М.: Машиностроение, 1978. 216 с.
3. Космические траекторные измерения / Под ред. П. А. Агаджанова, В. Е. Дулевича, А. А. Коростелева. М.: Сов. радио, 1969. 504 с.
4. Линник Ю. В. Метод наименьших квадратов и основы теории обработки наблюдений. М.: Физматгиз, 1958. 350 с.
5. Основы теории полета космических аппаратов / Под ред. Г. С. Нариманова и М. К. Тихонравова. М.: Машиностроение, 1972. 608 с.
6. Шебшаевич В. С., Дмитриев П. П., Иванцевич Н. В. и др. Сетевые спутниковые радионавигационные системы. М.: Радио и связь, 1993. 408 с.
7. Статистические методы обработки результатов наблюдений / Под ред. Р. М. Юсупова. М.: МО СССР, 1984. 563 с.
8. Эльясберг П. Е. Определение движения по результатам измерений. М.: Наука, 1976. 416 с.
9. Миронов В. И., Миронов Ю. В., Юсупов Р. М. Вариационное оценивание состояния нелинейной динамической системы по критерию максимального правдоподобия // Мехатроника, автоматизация, управление. 2009. № 11. С. 2—6.
10. Миронов В. И., Миронов Ю. В., Юсупов Р. М. Метод наименьших квадратов в задачах вариационного оценивания состояния нелинейных динамических систем // Информационно-управляющие системы. 2009. № 6. С. 2—6.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3
Рекуррентное систематическое помехозащитное преобразование кодов
17
11. Миронов В. И., Миронов Ю. В., Юсупов Р. М. Синтез итерационных алгоритмов решения краевых задач и нелинейных уравнений // Изв. вузов. Приборостроение. 2010. Т. 53, № 1. С. 9—14.
12. Мамон П. А., Половников В. И., Слезкинский С. К. Баллистическое обеспечение космических полетов. Л.: ВИКИ им. А. Ф. Можайского, 1990. 622 с.
Вячеслав Иванович Миронов Юрий Вячеславович Миронов Рафаэль Мидхатович Юсупов
Сведения об авторах — д-р техн. наук, профессор; Санкт-Петербургский институт информа-
тики и автоматизации РАН; E-mail: mironuv@yandex.ru — д-р техн. наук, доцент; Санкт-Петербургский институт информатики и
автоматизации РАН; E-mail: mironuv@yandex.ru — д-р техн. наук, профессор; Санкт-Петербургский институт информа-
тики и автоматизации РАН; директор; E-mail: spiiran@iias.spb.su
Рекомендована институтом
Поступила в редакцию 07.07.10 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2011. Т. 54, № 3