Известия РАН. Теория и системы управления, 2021, № 2, стр. 14-24

Введение. Предполагается, что наземный наблюдатель в дискретные моменты времени может измерять только угловые координаты объекта. Объект находится в свободном движении со скоростью, меньшей первой космической, на высотах не более 1000 км. Решается вопрос о том, можно ли по угловым измерениям оценить параметры движения с практически интересной точностью. Для ответа учитываем лишь основное, пренебрегая многочисленными влияющими на точность факторами: считаем, что движение кеплеровское, а измерения производятся с дискретностью Δt и со случайной погрешностью, дисперсия которой – σ². Цель работы – описать для инженера-разработчика приближенное, но достаточно простое средство анализа.

Задача оценивания параметров орбит космических тел (КТ) по угловым измерениям не является новой [1–9], однако с изменением технических средств измерения, требований к точности и вычислительных возможностей задача приобретает различные особенности. В частности, актуальна ситуация пассивной радиолокации, высокой частоты измерений, мощных вычислительных средств и высоких требований к точности определения расстояния от измерителя до точки падения.

Простые алгоритмы оценивания [1, 5, 6, 10] основаны на построении точечной оценки состояния КТ в фиксированный момент времени. Начальную оценку угловых величин получают регрессией угловых измерений, используя полиномы; оценки радиальных величин находят как корни трансцендентных уравнений орбитального движения, связывающих радиальные и угловые величины. Полученную оценку применяют в качестве начального приближения в различных методах уточнения, [1, 3–6], связанных с минимизацией квадратичной формы отклонения. В [7] излагается альтернативный подход, связанный с аппроксимацией апостериорной плотности вероятности состояния КТ, получаемой при аппроксимации априорной плотности вероятности состояния КТ и функции преобразования состояния [11–13].

Анализ эффективности методов оценивания обычно связывают с границами дисперсий многомерного неравенства Рао–Крамера, [1, 14–19]. В [1] теоретическая граница Рао–Крамера используется как эталонная мера ошибки оценивания, а в фундаментальной отечественной работе [5] авторы приводят универсальные соотношения для приближенного вычисления ошибок оценивания параметров орбит, которые в частном случае независимых гауссовых ошибок измерений совпадают с границей Рао–Крамера. Вычисление границ в ряде случаев представляет отдельную нетривиальную задачу, [14–16], например, в [15, 16] для вычисления границ Рао–Крамера в задаче фильтрации двумерного движения используется специализированная система координат угловых измерений. В данной работе при вычислении матрицы Фишера численно решается дифференциальное уравнение, представляющее второй закон Кеплера.

Разработчикам алгоритмов необходимо знать, с какой частотой и точностью нужно измерять и какую точность на каких орбитах можно получить. Эти вопросы могут прояснить информационные оценки точности, которая зависит от параметров орбиты. На практике особо выделяют класс траекторий, направленных в сторону наблюдателя, для которых весьма характерными являются малые изменения азимута. Малые изменения азимута имеют место в тех случаях, когда движение наблюдателя, связанное с вращением Земли, происходит вблизи плоскости орбиты. В этих случаях задачу можно приближенно считать двумерной, проецируя движение наблюдателя на плоскость орбиты, а единственным источником информации – измерения угла места. Ситуации с малыми изменениями азимута интересны также тем, что они наименее благоприятны для оценки параметров движения, и потому анализ их точности может служить ориентиром для использования угловых измерительных систем.

Цель данной работы состоит в том, чтобы предложить инженерам-разработчикам алгоритмов достаточно простые инструменты для приближенной оценки возможной точности при угловых измерениях. Задача оценки точности решается, опираясь на два момента, отличающих ее от других. 1. Выделено семейство неблагоприятных для наблюдателя орбит; это орбиты, плоскость которых касается орбиты движущегося наблюдателя, и точка касания является точкой пересечения соответствующих траекторий. 2. Обоснована простая приближенная плоская модель движения–измерения, в которой истинное движение наблюдателя (вне плоскости орбиты КТ) без изменения его скорости заменяется движением в плоскости орбиты КТ по окружности земной сферы.

Далее необходимо сделать следующие пояснения. Во-первых, нас интересуют орбиты, которые могут закончиться в окрестности наблюдателя. Во-вторых, семейство орбит, которые заканчиваются в точке наблюдателя (и в этой точке окружность наблюдателя касается плоскости орбиты), являются наиболее неблагоприятными в том смысле, что иначе информации будет больше, и тогда и качество оценок должно быть выше. Действительно, если в момент падения в эту точку вектор скорости наблюдателя не лежит в плоскости орбиты, то измерительной информации по орбите будет больше (чем в случае касания) за счет азимутальных измерений. Если же отодвинем точку падения от точки наблюдателя, то диапазон изменения угла места увеличится и мы опять имеем увеличение информации.

1. Двумерная модель измерения. 1.1. Предварительная схема: неподвижный наблюдатель. Рассматривается ситуация, в которой точка наблюдения Н (рис. 1) находится в плоскости орбиты и она неподвижна. В дальнейшем учтено влияние движения вне этой плоскости. Объект движется по эллипсу, главная ось которого составляет угол φ относительно направления на точку Н наблюдения. Орбита характеризуется трехмерным параметром а:

φ – угол между наблюдателем Н и главной осью эллипса, р – фокальный параметр, е – эксцентриситет орбиты. Уравнение орбиты в полярных координатах где r(θ) – расстояние до центра Земли, θ – истинная аномалия – угол направления на объект, отсчитываемый от направления на апогей. Измерения угла места начинаются после момента t₀ входа в зону видимости, т.е. когда объект выходит из-за горизонта G. Угол места ε = ε(t) в любой момент t ≥ t₀ определяется очевидным соотношением: тангенс угла есть отношение катетов (рис. 1): $\theta = \theta (t)$ – аномалия в момент t измерения. После умножения (1.2) на имеем

Это соотношение определяет ε = ${\varepsilon }({\theta };{\varphi },p,e)$, как функцию θ и параметров орбиты, но $\theta = \theta (t)$ – функция времени, и потому для любого фиксированного момента t

– функция трех переменных (заметим, что$\theta (t) = \theta (t;{\varphi ,}p,e)$ – тоже функция (φ, р, е)). По значению ε и фиксированным параметрам (φ, р, е) орбиты можно определить аномалию $\theta $, решая уравнение (1.3), сводимое к квадратному: ${\theta } = h(\varepsilon ;{\varphi ,}p{\text{,}}e)$ (см. Приложение). Для произвольного момента t значение $\theta (t)$ при известных параметрах орбиты р и е определяется по дифференциальному уравнению (второй закон Кеплера о постоянстве секториальной скорости) с начальным условием: при t = t₀ ${\theta = }{{{\theta }}_{0}}$; ${\omega }(\theta ;p,e) = {{\sqrt {p{\mu }} } \mathord{\left/ {\vphantom {{\sqrt {p{\mu }} } {{{r}^{2}}}}} \right. \kern-0em} {{{r}^{2}}}} = (\sqrt {\mu } {\text{/}}{{p}^{{3/2}}}){{(1 - e\cos {\theta })}^{2}}$ – угловая скорость; $\mu = 3.9860044 \times {{10}^{{14}}}$ м³/с² – геоцентрическая гравитационная постоянная.

В дискретные моменты времени t₀, t₁, …, t_n, t_i = t_{i – 1} + Δt, i = 1, …, n, имеем измерения

δ_i – независимые случайные ошибки измерений со средним 0 и дисперсией $D\delta _{i}^{2} = {{\sigma }^{2}}$; x = = (x₀, …, x_n) – исходная совокупность измерений.

В некоторый известный момент времени t₀ имеется начальное измерение

По значению $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\varepsilon } {}_{0}$ по (1.3) определим приближенно ${{{\tilde {\theta }}}_{0}}(a)$ – аномалию ${{{\theta }}_{0}}$ в момент t₀ как функцию параметра а (начальное условие при t = t₀):

${{{\tilde {\theta }}}_{0}}(a)$ – оценка начальной аномалии ${{{\theta }}_{0}}(a)$, δ_θ – погрешность оценки, с.к.о. (среднее квадрати́ческое отклоне́ние) которой ${\sigma /|\varepsilon }_{{\theta }}^{'}{\text{|}}$, ${\varepsilon }_{{\theta }}^{'}$ – производная функции ε по θ. При известной орбите $({\varphi ,}p{\text{,}}e)$ погрешность в ${{{\tilde {\theta }}}_{0}}(a)$ определяется погрешностью измерения угла в момент t₀. Заметим, что уравнение (1.4) определяется двумя параметрами (p, e), а зависимость решений от параметра φ происходит только через начальное условие.

1.2. Учет движения наблюдателя вследствие вращения Земли. Наблюдатель вследствие вращения Земли движется по окружности вдоль своей параллели. Будем считать траекторию КТ неблагоприятной для измерителя, если плоскость орбиты КТ касается орбиты движущегося измерителя, и точка касания является общей точкой (точкой пересечения) для соответствующих траекторий. Эти траектории характеризуются, во-первых, настолько малым диапазоном измеряемых азимутальных углов, что измерениями азимута можно пренебречь, и, во-вторых, диапазон изменения угла места минимален при последних (наиболее информативных) наблюдениях (так как траектория входит в точку О наблюдения). Следовательно, траектории неблагоприятны, если в момент падения вектор скорости наблюдателя лежит в плоскости орбиты.

Если реверсировать время, то можно сформулировать так: неблагоприятны для измерителя те траектории, которые исходят из точки наблюдения, и орбитальная плоскость проходит через вектор скорости наблюдателя. В этом случае взаимное движение выглядит следующим образом: из точки касания О начинаются одновременно два движения: наблюдатель движется по окружности – параллели с широтой φ_Н, а КТ движется в орбитальной плоскости.

Простые соображения показывают, что существует фиктивное движение измерителя (по окружности сферы S_З Земли в орбитальной плоскости), эквивалентное (в смысле получения углов места) истинному движению (вне орбитальной плоскости). Для этого достаточно, чтобы в любой момент времени два расстояния (от КТ до фиктивного движения и до истинного движения) были одинаковы. С помощью такого эквивалентного фиктивного движения Δφ(t) задача сводится к плоской. Однако фиктивное эквивалентное движение зависит от точки расположения КТ.

Построим приближенно эквивалентное движение, независимое от КТ. Заменим равномерное движение наблюдателя вдоль своей параллели движением по окружности S₃ сферы в орбитальной плоскости, сохранив его линейную скорость ${{{\omega }}_{{\text{З}}}}{{R}_{{\text{З}}}}\cos {{{\varphi }}_{{\text{Н}}}}$. Заметим, что для экватора ${{{\varphi }}_{{\text{Н}}}} = 0$ и на северном полюсе ${{{\varphi }}_{{\text{Н}}}} = {\pi /}2$, это фиктивное движение является эквивалентным.

Оценим это приближение. Оба движения, истинное и фиктивное, исходят (вернее, приходят) одновременно в общую точку О. Определим разность длин двух отрезков, исходящих из точки КТ на истинную траекторию и на фиктивную. Если длины в любой момент одинаковы, $\Delta d = 0$, то углы места равны. Если они не одинаковы, $\Delta d \ne 0$, то оценим величину погрешности Δε(t; d, измеряемого угла места или величину сдвига $\Delta l(t;d,{\varepsilon }) \approx {{\Delta d} \mathord{\left/ {\vphantom {{\Delta d} {\cos {\varepsilon }}}} \right. \kern-0em} {\cos {\varepsilon }}}$ фиктивного наблюдателя, при котором угол места равен истинному. Необходимые формулы даны в Приложении. В табл. 1 в каждой клетке приведены расчетные значения $\Delta l$ и $\Delta {\varepsilon }$, зависящие от расположения КТ (угла ε и дальности d) и времени t движения до точки О; значение широты принято φ_Н = 60°.

Видно, что смещение $\Delta l(t,\varepsilon ,d)$ относительно эквивалентных точек меньше 5 м в диапазоне до 300 с и до 72° при дальностях от 100 до 1000 км. Если $\Delta l$ переводить в погрешность $\Delta {\varepsilon }$ по углу места, то, например, за 300 с до точки О при дальности 300 км и ε = 48° имеем $\Delta l$ = 3.84 м, что означает $\Delta {\varepsilon } \approx 0.033$ угл. мин., что меньше с.к.о. измерения в сотни раз.

Фиктивное движение при вычислениях учитывается добавлением к неизвестному параметру φ в (1.2) линейно изменяющееся во времени (со скоростью ${{{\omega }}_{{\text{З}}}}\cos {{{\varphi }}_{{\text{Н}}}}$) слагаемое Δφ(t) = = $ \pm {\kern 1pt} {{{\omega }}_{{\text{З}}}}{\text{cos}}{{{\varphi }}_{{\text{Н}}}}$(T – t), T – момент падения, причем $\Delta {\varphi (}T) = 0$. Заметим, что формулы производных по параметрам, необходимые для матрицы Фишера, остаются без изменения; изменяются лишь их значения при вычислениях.

2. Вычисление матрицы Фишера. Информационная матрица Фишера, как известно, есть

где M, Т – символы математического ожидания и транспонирования. Если измерения в (1.5) принять независимыми и нормальными, то для матрицы Фишера справедливо выражение: где слагаемое под знаком суммы есть матрица Фишера для момента t_i, (∂ε_i/∂a) – вектор-строка. В соответствии с неравенством Рао–Крамера диагональные элементы матрицы являются нижними границами для дисперсий несмещенных оценок.

Производные (∂ε_i/∂a) функций определяются из уравнения (1.2), и при фиксированном $t = {{t}_{i}}$ их можно выписать в явном виде (см. Приложение). Эти функции зависят от аномалии ${\theta = \theta (}t{\text{)}}$ и от ее производных по параметрам:

где обозначено ${{{{\theta }_{{\varphi }}} = \partial \theta } \mathord{\left/ {\vphantom {{{{\theta }_{{\varphi }}} = \partial \theta } {\partial {\varphi }}}} \right. \kern-0em} {\partial {\varphi }}}$, ${{\theta }_{p}} = {{\partial \theta } \mathord{\left/ {\vphantom {{\partial \theta } {\partial p}}} \right. \kern-0em} {\partial p}}$, ${{\theta }_{e}} = {{\partial \theta } \mathord{\left/ {\vphantom {{\partial \theta } {\partial e}}} \right. \kern-0em} {\partial e}}$. Заметим, что значения производных аномалии являются функциями времени. Для вычисления $\theta = \theta (t)$ и ее производных используем дифференциальное уравнение (1.4) с начальным условием (1.6). Формула для ${{{\theta }}_{0}}(a)$ выписана в Приложении. Продифференцировав (1.4) и (1.6) по параметрам, получим три дифференциальных уравнения: где угловая скорость ${\omega = \omega }(p,e,\theta )$ зависит от трех переменных; а справа стоят ее производные по параметрам; начальные условия: Решения (1.4) и (2.3) для всех четырех функций $\theta = \theta (t)$, ${{\theta }_{p}}(t)$, ${{\theta }_{e}}(t)$, ${{\theta }_{{\varphi }}}(t)$ в моменты ${{t}_{i}} = i\Delta t$, $i = 1$, ..., n, получаем простым численным методом. Вычислительные формулы производных скорости по параметрам даны в Приложении. Производные ${{\partial {{\varepsilon }_{i}}} \mathord{\left/ {\vphantom {{\partial {{\varepsilon }_{i}}} {\partial a}}} \right. \kern-0em} {\partial a}}$ по параметрам для матрицы Фишера с изменением i вычисляются последовательно.

3. Некоторые расчетные результаты. 3.1. Неподвижный наблюдатель. Исходные параметры орбит: α и V – угол падения и скорость КТ (формулы связи $(V,{\alpha })$ с $(p,e)$ даны в Приложении); L_а = ${\varphi }{{R}_{З}}$ – расстояние до точки апогея; H_a – высота в апогее; φ – угол главной оси выбирался так, чтобы орбита была “неблагоприятной”, Δt = 5с – дискретность измерений, σ = 6' – стандартная ошибка измерений; N – количество наблюдений (доля времени наблюдения составляет 0.9).

Результаты нижних границ для с.к.о. приведены в табл. 2 для следующих параметров: σ_Va – скорости V_a в апогее, σ_Vy – вертикальной V_y скорости, σ_Δ – расстояния Δ до точки падения. По результатам делаем следующие выводы:

– горизонтальную составляющую V_a можно оценить с точностью до единиц и десятков м/с,

– точность по месту падения для высоких и с малой скоростью орбит можно получить порядка сотен метров (выделено 571, 645, 719, …),

– вертикальную составляющую вектора скорости оценить невозможно (погрешности порядка единиц километров в секунду).

3.2. Подвижный наблюдатель. Широта наблюдателя принята равной ${{{\varphi }}_{{\text{н}}}}$ = 60°. При встречном движении наблюдателя и КТ результаты показаны в табл. 3. Видим увеличение всех значений (понятна причина: увеличение скорости взаимного приближения). При движении в одну сторону наблюдателя и КТ результаты показаны в табл. 4. Наблюдаем уменьшение всех значений (из-за уменьшения скорости взаимного приближения).

Заключение. Выделено семейство неблагоприятных для наблюдателя орбит, плоскость которых касается орбиты движущегося наблюдателя, и точка касания является точкой пересечения соответствующих траекторий. Обоснована простая приближенная плоская модель движения–измерения, в которой истинное движение наблюдателя (вне плоскости орбиты КТ) без изменения его скорости заменяется движением в плоскости орбиты КТ по окружности земной сферы. Обоснование проведено вычислением погрешностей.

На этой модели проанализирована потенциальная точность с помощью информационной матрицы Фишера. Вычисление матрицы Фишера основано на численном решении дифференциального уравнения, выражающего второй закон Кеплера.

Рассчитаны границы для точностей параметров для семейства неблагоприятных орбит. Примеры показывают, что по угловым измерениям возможно получение практически интересных точностей для некоторых параметров (например, горизонтальной составляющей скорости).

В Приложении приведены все формулы, необходимые для проведения аналогичных расчетов.