CC BY
147
УДК 629.78
КИНЕМАТИЧЕСКИЙ СИНТЕЗ ПРОГРАММНОГО УГЛОВОГО ДВИЖЕНИЯ КОСМИЧЕСКОГО АППАРАТА ПРИ ОПТИКО-ЭЛЕКТРОННОЙ СЪЕМКЕ ЗЕМЛИ
© 2007 С.А. Бутырин
Самарский научный центр РАН
Представляются методики расчета программного углового движения космического аппарата для трех вариантов оптико-электронной съемки поверхности Земли: съемка с выравниванием продольной скорости движения изображения в фокальной плоскости орбитального телескопа, съемка при совмещении следа линии визирования на поверхности Земли с геодезической линией и съемка с постоянной ориентацией оси визирования телескопа в орбитальной системе координат. Приводятся некоторые результаты, полученные с помощью созданного программного обеспечения.
Введение
Оптико-электронная съемка с выравнивания продольной скорости движения изображения (СДИ) обладает рядом важных преимуществ и применяется для получения изображения с высоким качеством. Недостатком этого метода является отклонение линии маршрута от первоначально заданного направления сканирования, причем скорость такого отклонения увеличивается пропорционально угловому расстоянию от трассы полета. Метод съемки по геодезической линии имеет практическое преимущество, связанное с возможностью обеспечения прохождения маршрута через заданные объекты на поверхности Земли. Метод съемки с постоянной ориентацией линии визирования (ЛВ) в орбитальной системы координат (СК) не требует существенных угловых эволюций КА. В отличие от известных [2,3] работ здесь развиваются [4] методы синтеза маршрута съемки в кватернионах относительно инер-циальной СК. Это позволяет повысить точность расчетов и использовать синтезированное программное движение в контуре управления космическим аппаратом (КА) без дополнительных преобразований. Орбитальное движение центра масс КА считается известным.
Модели и постановка задач
Вводятся следующие правые ортогональные системы координат: экваториальная инер-циальная СК (ИСК) Оех*у^ эпохи т, имею-
щая начало в центре масс Ое Земли: ось Оех: направлена в точку весеннего равноденствия эпохи т, а ось О^1 направлена по нормали к плоскости экватора к Северному полюсу мира; гринвичская СК (ГСК) Оех8у^8, связанная с эллипсоидальной, вращающейся в ИСК относительно оси О^:(О^8) моделью Земли [2]; орбитальная СК (ОСК) Осх°у^°, имеющая начало в центре масс КА Ос: ось Осу° которой направлена по радиали орбиты, а ось - по отрицательному направлению нормали орбиты; горизонтная СК (ГорСК) Орхьу^ь, имеющая начало в точке Р, на заданном превышении Ир от поверхности Земли: ось ОРхь направлена по орту внешней нормали к поверхности эллипсоида в точке Р; ось - по касательной к геодезическому меридиану точки Р на север, а ось Оруь - на местный восток; маршрутная СК Орхшу^ш получается из ГорСК разворотом последней на угол ^ /2 - А вокруг оси Орхь, где А - текущий азимут маршрута съемки; визирная СК (ВСК) ОД^^ связанная с корпусом КА, началом в центре проектирования S и визирной осью О5ху, направленной по оптической оси ОТ в сторону наблюдаемого объекта. Для простоты принимается, что центр масс КА (С) и центр проектирования совпадают. Начало координат фокальной плоскости (ФП) орбитального телескопа Ог, в которой расположены оптико-электронные преобразователи (ОЭП), сдвинуто от начала ВСК в отрицательном направлении оси Окху на величину фокусного расстояния орбитального телескопа (ОТ) F,
оси ФП совпадают по направлению
с одноименными осями ВСК.
Поверхность Земли описывается в виде сжатого в направлении ее полюсов эллипсоида вращения с полуосями ае, се. Программное положение точки Р маршрута на земной поверхности может определяться как декартовыми координатами х?, у?^?, так и эллипсоидальными геодезическими координатами: широтой Вр, долготой Ьр и высотой Нр.
Математическая модель движения изображения позволяет определить координаты и скорости точек изображения наблюдаемого участка земной поверхности в фокальной плоскости ОТ по заданному поступательному и угловому движению КА в ИСК.
Угловое движение КА определяется кватернионом ориентации ВСК относительно ИСК Л[ = Л = (X0Д}, X = (X., i = 1:3}, получаемым из решения дифференциального кинематического уравнения
Л 0) = | Л 0) о ш(1), (1)
с начальным условием Л (1;0), при известном в функции времени 1 векторе абсолютной угловой скорости ш = шу = шу(1;) = (юД*), i = 1:3}, заданном в разложении на оси ВСК.
Пусть вектор г = гу из центра проектирования S до точек изображения задан в разложении на оси ВСК, а вектор R = R1 из центра проектирования до точек наблюдаемой поверхности - в ИСК (рис. 1).
В силу коллинеарности векторов г и Я связь между их проекциями в разложении на оси ВСК будет иметь вид
г = к • Я, или г = к • Л о Я о Л (2)
где г = (х,у^}, к = -БШ; к = к• 3-проекция Я на визирную ось.
Задача определения скорости движения изображения точек поверхности Земли в фокальной плоскости сводится к вычислению полной производной по времени от выражения (2) для координат. Получим
г = V = к • (Л о Я о Л+Л о Ко л+Л о Я о Л)+к • Л о Я о л . (3)
Подставляя выражение для производных, преобразуем сумму первого и третьего слагаемого в (4) к виду
Ло Я о Л + Л о Я о Л = (Л о Я о Л) х ш .
Для вектора Я можно записать очевидное равенство (рис. 1)
К = К © - к с,
(4)
где К - вектор, направленный из центра Земли в центр проектирования S, а Яф - вектор, направленный из центра ГСК в наблюдаемую точку земной поверхности (см. рис. 1).
Производная по времени вектора Я (4) в разложении на оси ИСК будет равна
я=к©-Rs, (5)
где К© = ш©хк©, ш© =(0,0,ю©}, -
модуль вектора ш© угловой скорости вращения Земли, а К s - известный вектор линейной скорости движения центра масс КА.
Вектор Я определяется как пересечение направления, определяемого заданными координатами точки фокальной плоскости, с поверхностью земного референс-эллипсоида. Из (4) находим компоненты вектора наблюдаемой точки земной поверхности К © в проекции на оси ИСК в виде
(6)
К © = т • е + К,,
Рис. 1. Схема оптической проекции
где е - орт вектора наблюдаемой точки (вектора к) в проекции на оси ИСК, т > 0 -неизвестный постоянный множитель. Орт е в проекции на оси ВСК равен еу = (Б,у^}, где у^ - координаты заданной точки фокальной плоскости, Б - фокусное расстояние.
Точка конца вектора К © лежит на поверхности референс-эллипсоида, и, следовательно, его координаты удовлетворяют уравнению поверхности эллипсоида
(X?)2 + (У?)2 + (^)2(1 + 8 2) = (ае + Н?)2, (7) где 82 = (а2 - с2 )/с2 - квадрат второго эксцентриситета земного эллипсоида, Нр - заданное превышение.
Подставляя (6) в (7) приходим к квадратному уравнению относительно множителя m > 0, решение которого позволяет из (6) определить вектор R в ИСК. Для определения эллипсоидальных геодезических координат B,L наблюдаемой точки земной поверхности, достаточно перепроектировать вектор Rs в ГСК, используя матрицу перехода из ИСК в ГСК Gg, определяемую углом вращения Земли на момент расчета.
Предполагается, что ОЭП расположены на отрезке прямой (линейке ОЭП) по оси Oz СК ФП симметрично относительно центра ФП. Продольные Vy(y,z) и поперечные Vz(y,z) составляющие вектора скорости движения точек изображения направлены соответственно перпендикулярно и вдоль линейки ОЭП.
Съемка с выравниванием продольной составляющей СДИ
Этот метод накладывает на движение КА следующие условия:
1. Продольная составляющая СДИ в центре ФП должна быть равна номинальному значению Vy(0,0) = V .
2. Поперечная составляющая СДИ в центре ФП Vz(0,0) = 0 .
3. Продольная составляющая СДИ должна быть постоянной
Vz(0,z) = V* = const, -1/2 < z < l/2.
Последнее требование не может быть выполнено строго, так как продольная составляющая скорости, нелинейно зависит от координат y,z точек фокальной плоскости. Однако для узкоугольных оптических систем отклонения от линейности невелики.
Решение задачи разбивается на 2 части:
1 - определение начальной ориентации ВСК относительно ИСК в момент начала съемки;
2 - определение вектора угловой скорости в разложении на оси ВСК.
Задача определения начальной ориентации КА состоит в определении кватерниона ориентации ВСК относительно ИСК Л 0 = Л (t0) по следующим данным: начальному азимуту маршрута А0; геодезическим координатам и превышению наблюдаемой точки земной поверхности Bg,Lg,H; параметрам орбитального движения Rs, Vs = Rs. Решение этой задачи содержит 2 этапа:
1. Вычисление матрицы MV(t0) ориентации ВСК относительно ИСК.
2. Вычисление кватерниона ориентации ВСК относительно ИСК Л [(i) по элементам матрицы M[(t) .
Матрица M [ (t0) определяется из условий: 1) проекция орта оси Oy ВСК должна быть направлена по оси Oy маршрутной СК; 2) орт линейки ОЭП (оси Oz ВСК) определяется как перпендикуляр одновременно визирной оси (оси Ox ВСК) и оси Oy маршрутной СК
MV = MV(t0) = [ed I eJz x ed | e^]T. (8)
Здесь ed - орт наклонной дальности (оси Ох ВСК) D:(t) в проекции на оси ИСК; ez = (ey x ed)/1| ey x ed ||, а e^ - орт оси Oym маршрутной СК в ИСК.
Орты ed и e:y определяется по известным координатам точки начала R0(t) и азимуту А маршрута в ГСК
D:(t) = Gg(t)Rg(t)-Ri(t), ed = ed(t) = D:(t)/||Dd ||, t = ^
(9)
ey = Mhehy, ehy = {0;sin(A);cos(A)}, (10)
где ey - его проекция в ГорСК, определяемая заданным азимутом маршрута А, а матрица M h перехода из ГорСК в ИСК определяется через вектор R0 (см. рис. 1) в ИСК
Mh = Mh(t0) = [eR I ey j eR xey], (11)
где
ey = e^x eR/^x eR |; eR = R0/|R J, e^ = [0,0,1]T-орт оси Oz ГСК в ИСК.
Кватернион ориентации ВСК относительно ИСК Л[(t) по элементам матрицы M[ по известному алгоритму Стенли (S.W. Stanley)
ЛV(t) = fst(M V). (12)
Требуемая угловая скорость движения КА определяется из условий: 1) скорость движения изображения в центре ФП Vy(0,0) = Vz(0,0) = 0 , 2) производная скорости движения изображения в центре ФП по поперечной переменной 5Vy (y, z) / dz = 0 . В результате получается соотношение
<(t) "0 c31 /c11 0 " " 0 "
®V(t) = 0 0 -1 - v — s 0
<(t) 0 1 0 V*/F
где вектор
~s = Mh-[Л о (R s - ©0 x Rs) о Л]^ (14)
представляет СДИ, нормированную к фокусному расстоянию телескопа, а с , с - элементы матрицы
С£(1) = (] = 1:3,] = 1:3} = М ?(1)М ¿(1).
Для получения маршрута движения требуется интегрировать кинематическое уравнение (1) на интервале дискретности Д1 для определения новых значений кватерниона Л в момент времени 1 + Дь. Далее по (13)-(14) вновь вычисляются угловые скорости. Процедура циклически повторяется до завершения интервала времени маршрута.
Съемка маршрута по геодезической линии
Здесь маршрут съемки задается двумя точками на земной поверхности: начальной, с геодезическими координатами: широтой В0 и долготой Ь0 и конечной - с координатами Вк , Ьк причем след линии визирования представляет геодезическую линию, т.е. линию пересечения земного эллипсоида плоскостью, проходящей через эти точки и центр Земли. Компоненты вектора угловой скорости КА относительно ГСК юу, шz вычисляются при таких требований: в центре ФП продольная СДИ V должна иметь требуемую величину V, а поперечная СДИ Vz должна равна нулю. Эти условия дают соотношение
(15)
где V ж = V ж ^) = (уя. , I = 1:3} - нормированная к наклонной дальности 3 линейная скорость движения центра масс КА относительно ГСК в проекции на оси ВСК, определяемая по (14). Методика определения компоненты угловой скорости КА, направленной по визирной оси ш^ , основана на численном дифференцировании кватерниона ориентации ВСК относительно ИСК.
Пусть Л(10) - кватернион ориентации ВСК относительно ИСК) в некоторый момент времени t0, определенный в начальной точке маршрута В0 = В0(^),Ь0 = Ц^) и Лр(1 + Д1) требуемое (программное) значение кватерниона ориентации КА по истечении малого промежутка времени Д1, при условии, что точка пересечения линии визирования (ЛВ) с земной поверхностью
4(t)~ Vs3
®v(t)_ _- Vs2(t) - V*/F_
принадлежит геодезической линии. Тогда из кинематического уравнения (1) имеем
ro(t) = 2Ä(t) о A(t), (16)
или приближенно
ю (t) и 2k(t) о AA(t)/ At, АЛ(t) = Лp(t + At) - Л(t).
(17)
Выражение для определения шХ следует из развернутой записи формулы (17) raV(t) = 2(X0(t)A^j (t) - \- + / At.
(18)
Кватернион А p (t + At) определяется по алгоритму S.W. Stanley как функция от матрицы Mp = M V(t + At)
А p (t + At) = fst(Mp), (19) по соотношениям (8)-(12), в которых формула (10) для определения e:y заменяется на
ey = N1 х e), (20)
где N1 - орт нормали к плоскости геодезической линии, e) - орт вектора, направленного из центра Земли в точку на маршруте съемки. Орт нормали равен
N1 = (el0 хe)k)/| e)0 хe)k|, (21) где e)0 х e)k - вектора, направленные из центра Земли в начальную и конечную точки на маршруте съемки.
Если определить начальный кватернион A(t0) по соотношениям (19)-(21), то при значениях угловой скорости вида (15) и произвольной проекции угловой скорости на визирную ось точка пересечения ЛВ с земной поверхностью будет двигаться по касательной к геодезической линии. Используя этот факт, можно определить вектор угловой скорости КА относительно ИСК в виде
ю = {°;©; X} +
ю! = {<} = AV(t) о юе о AV(t), (22) юе = {0;0;©е}
Считая, на малом промежутке времени [t, t + At], угловую скорость КА приближенно постоянной, для момента времени t + At определяются положение ЦМ
Rl(t + At) = Rl(t) + VsI(t)At, (23) и кватернион ориентации КА
A(t + At) = A(t) о [cos(| ю | At / 2), e • sin(| ю | At /2)],
e = ю/|ю |
По известным R1(t + Л0 и Л^ + Л0, используя соотношения, приведенные выше, определяется радиус-вектор R 0т^ + Лt), направленный из центра Земли в новую точку пересечения ЛВ с поверхностью Земли, и вектор наклонной дальности Б1^ + Л0 . Далее, по формулам (8),(9),(20) определяется новое значение программного кватерниона ориентации Л р ^ + Лt), а по формуле (18) значение проекции вектора угловой скорости КА на визирную ось. После чего окончательно формируется вектор угловой скорости КА
ш(0 = К; ш^} + {0; ш;; ®12}. (25)
Интегрируя уравнение (1) при начальных условиях Л(^), ю(t0) с определением угловой скорости ю(0 на каждом шаге интегрирования по соотношениям (18)-(25) получим программное движение КА Л(0,ю(0^ е[0,Т] при съемке маршрута по геодезической линии. Момент времени окончания интегрирования определяется условием прохождения линией маршрута траверса заданной конечной точки.
Алгоритм имеет методические погрешности, связанные с неточностью численного дифференцирования (16)-(18), которые будут приводить к отклонению расчетного маршрута от заданной геодезической линии. Численные расчеты показывают, что для маршрутов с угловой длиной 1о дуги геодезической линии (111 км на местности) погрешности метода не превышают 1 метра. При удалении точек маршрута от трассы, соответствующих углам отклонения визирной оси от надира до 60о, погрешности увеличиваются, но не превышают нескольких десятков метров, что вполне достаточно для практического применения метода.
Съемка с постоянной ориентацией визирной оси в ОСК
При этом методе съемки траектория следа центральной ЛВ по поверхности Земли проходит на некотором расстоянии от трассы, величина которого, при заданной ориентации ЛВ в ОСК, зависит от высоты полета. Ориентация ЛВ в ОСК задается двумя углами, причем первый поворот выполняется на угол а вокруг оси О2 ВСК (ОСК), второй -на угол Р вокруг нового положения оси Оу ВСК.
Угловое движение КА определяется кватернионом Л(0 ориентации ВСК относительно ИСК и вектором абсолютной угловой скорости © = юу(0, заданной в разложении на оси ВСК для произвольного момента времени 1
Кватернион ориентации ВСК относительно ИСК равен
Л [(0 = Л 0(0 о ЛХ). (26)
где Л 0(0, Л 0(0 кватернионы ориентации ИСК относительно ОСК и ОСК относительно ВСК соответственно. Кватернион Л 0 определим через кватернион Л 0/¥=0(0 при нулевом угле V поворота вокруг визирной оси и кватернион поворота вокруг визирной оси Л ¥ = {^(у /2);е„ вт( V/2)}, е¥ = {1;0;0},
Л0(0 = Л 0/¥=0(0 о Л¥. (27)
Кватернион Л 0/¥=0(0 определяется углами ориентации ЛВ в ОСК а и Р в виде
Л о (0 — {cacß; CaSß; SaSß; SaCß};
co — cos(^/2), sv — sinO^), p — a, ß
, (28)
а кватернион Л 0(t) — fs(M 0) перехода из ИСК в ОСК - как функция fs от матрицы M0 по алгоритму S.W. Stanley. Здесь матрица M0 составлена из ортов столбцов R'(t), Vs'(0 и зависит от параметров движения центра масс КА, именно
M0(t) — [е, хез I е, | ез], (29)
где
е, — RI(t)/|Rs(t) |; ез — (V'(t) х RI(t))/1V'(t) х RI(t) |.
Из условия неизменной ориентации ЛВ в ОСК следует равенство угловой скорости КА в проекции на оси Oz и Oy ВСК соответствующим проекциям орбитальной угловой скорости КА ю0 на эти оси. Спроектируем ю0 на оси ВСК ш! (t) — К; ш;; ш!} — Л! (t) о ю0 о Л] (t). (30) Искомый вектор угловой скорости КА шv(t) будет включать ®V и ш!! из (29), а проекция угловой скорости КА на визирную ось определяется из условия, что вектор скорости движения изображения (СДИ) в центре фокальной плоскости ОТ должен быть перпендикулярен линейке ОЭП, т.е. оси Oz ВСК. Угол поворота V равен
уф =
I- аг^У^/У^О), У* * 0,
л-Мвп(У^(1))/2, У* = 0,
(31)
где У^) - вектор линейной скорости точки следа ЛВ относительно поверхности Земли в проекции на оси ВСК,
У; (t) = VI (fЩ(г); V; (г)} =
=Л;^ а) ° (У а) - У с)) ° л]1г_0 ^)+у: а). (32)
Здесь кватернион Л V ¥=0 определяется по формулам (26)-(29) при V =0, а векторы У1, УШ вычисляются по формулам
УШ(1) = Dv(t) X ш!(1); (33) У^) = ше х кР(1); (34)
(t) = Л^ (t) ° О1 (0 о ЛV/v=o (t), (35) в которых обозначено: юе = {0;0; ше}- вектор угловой скорости вращения Земли в ИСК; Rp - вектор, направленный из центра Земли в наблюдаемую точку Р в ИСК; О1^) = Rp(t)-Rs(t)- вектор наклонной дальности в ИСК.
Соотношения (26) - (35) позволяют для каждого дискретного момента времени tг, г = 0 : N из временного интервала съемки маршрута [0, Т], Т = Д1 - ^ -1), где Д1 - шаг
расчета по времени, последовательно определить: угол V кватернион Л и угловые скорости шЩу(1) и ш^) . Искомая угловая скорость по визирной оси определяется численно из условия, что на малом интервале времени Д1 ее можно считать постоянной Ш1 О,-,) = (V (^ ) - V О,- )) / Д1, г = 1: N.
Некоторые результаты расчетов
На основе представленных методик разработан комплекс программ визуализации маршрутов съемки поверхности Земли. На рис. 2. приведен пример съемки трех маршрутов, последовательно проходящих через Лиссабон (метод с постоянной ориентацией визирной оси в ОСК а = 0°, Р = 20°, длительность 40 с), Мадрид (метод с выравниванием продольной СДИ с азимутом 45о длительностью 120 с) и Лион (метод съемки по геодезической линии, когда геодезические координаты начальной точки маршрута 45о с.ш. и 3о.5 в.д., а координаты конечной точки: 46о с.ш. и 5о в.д.). Расчет выполнен для круговой орбиты высотой 680 км, фокусного расстояния телескопа F = 6 м, длины линейки ОЭП 500 мм при требуемой продольной СДИ в центре ФП у*= 30 мм/с.
Рис. 2. Маршруты съемки
Заключение
Представлены методики кинематического синтеза программного углового движения КА для трех практически применяемых вариантов оптико-электронной съемки поверхности Земли и некоторые численные результаты, полученные с помощью созданного программного обеспечения. Контрольный расчет скоростей движения получаемого изображения в фокальной плоскости телескопа полностью подтвердил достоверность предложенных методов синтеза маршрутных движений КА.
Работа поддержана РФФИ (04-01-96501, 07-08-97611), Президиумом РАН (программа фундаментальных исследований 22) и Отделением энергетики, механики, машиностроения и процессов управления РАН (программа 15)
СПИСОК ЛИТЕРАТУРЫ
1. В.Н. Баранов, И.И. Краснорылов и др. Космическая геодезия. М.: Недра, 1986.
2. Батраков А.С. Общая модель расчета и анализа скорости сдвига оптического изображения при съемке Земли из космоса // Исследование Земли из космоса. 1989. № 4.
3. Анатольев А.Ю, Батраков А.С. Математическое моделирование космической оптико-электронной съемки с программным сканированием. СПб.: ВИКА им. А.Ф. Можайского, 1995.
4. Бутырин С.А. Синтез маршрута космического аппарата наблюдения Земли с выравниванием продольной составляющей скорости движения изображения // Вестник СамГТУ, сер. Технические науки. Самара: СамГТУ 2005. Вып.37.
KINEMATIC SYNTHESIS OF THE SPASECRAFT PROGRAMMED ATTITUTE MOTION AT THE EARHT OPTIC-ELECTRONIC SURVEY
© 2007 S.A. Butyrin
Samara Science Centre of Russian Academy of Sciences
Methods of calculating the spacecraft programmed attitude motion for three variants of optic-electronic survey by the Earth surface, are presented, namely: survey with alignment of an image motion lengthwise velocity into focal plane of orbital telescope, survey with coinciding the line-of-sight trace on the Earth surface with geodesic line and survey with fixed position of the telescope line-of-sight in orbital reference frame. Some results are represented which were obtained by elaborated software.
