Моделирование комбинированных наблюдений глобальных навигационных спутниковых систем

Texto integral

ВВЕДЕНИЕ

В науке и технике с каждым годом растет интерес к спутниковым системам глобального позиционирования [1]. Данная технология находит свое применение во многих сферах жизнедеятельности человека, в том числе — в научно-технической. Особенно важную роль спутниковые навигационные системы играют в космических проектах, направленных на дистанционное зондирование Земли. Например, опыт работ космических геодезических проектов CHAMP, GRACE и GOCE показал, что для прецизионного определения параметров гравитационного поля Земли необходимо определять абсолютное местоположение космического аппарата (КА) на низкой околоземной орбите с точностью до нескольких сантиметров [2]. Такую точность на орбите можно обеспечить с помощью глобальных навигационных спутниковых систем (ГНСС) [3, 4]. Они позволяют выполнять баллистические и навигационные задачи с высокой надежностью, в том числе и для низкоорбитальных космических аппаратов. К таким системам относят системы GPS, ГЛОНАСС, Galileo и BeiDou. Их характеристики представлены в табл. 1. С их помощью можно с высокой точностью решать задачу навигации не только наземных приемников, но и космических аппаратов, находящихся на низкой околоземной орбите (НОО). С меньшей точностью можно определять средние (СОО), геостационарные (ГСО), наклонные геосинхронные (НГСО), высокоэллиптические (ВЭО), и даже окололунные орбиты [5].

Таблица 1. Характеристики глобальных навигационных спутниковых систем

Для того чтобы корректно обработать информацию, приходящую с различных навигационных спутников, и получить устойчивое и высокоточное навигационное решение, необходимо правильным образом составить модель наблюдений ГНСС [6, 7]. Данная модель содержится во всех методах дифференциального уточнения орбиты и является ключевым фактором (наряду с моделью сил, действующих на движение КА), влияющим на точность определения абсолютных координат бортового приемника низкоорбитального космического аппарата (НКА). В работах, посвященных данной проблематике, обычно приводится перечисление эффектов модели наблюдений, которые были учтены при ее составлении. Для того чтобы понять, какое именно влияние оказывает тот или иной эффект, необходимо провести дополнительные поиски литературы по каждому из них. В зарубежной литературе представлен ряд работ, в которых модель наблюдений ГНСС описана комплексно [8–10]. Однако в этих работах опускаются детали о некоторых эффектах модели наблюдений, приобретающих свою значимость только при высокоточных расчетах.

Кроме того, при выполнении комбинированных наблюдений ГНСС, т.е. с использованием различных навигационных спутниковых систем, возникает необходимость учитывать особенности каждой системы по отдельности. В частности, система ГЛОНАСС, в отличие от систем GPS и Galileo, основана на частотном распределении каналов, что необходимо учитывать при чтении и обработке измерений. В данной работе приводится модель наблюдений ГНСС, учитывающая использование различных систем спутниковой навигации. Данная модель была опробована на имитационных и реальных измерениях низкоорбитального космического аппарата и показала свою состоятельность для получения высокоточного баллистико-навигационного решения.

КООРДИНАТНО-ВРЕМЕННОЕ ОБЕСПЕЧЕНИЕ

Для рассмотрения движения КА и расчета возмущающих ускорений необходимо обозначить системы координат и времени, а также связывающие их преобразования.

Системы времени

На сегодня, для высокоточного моделирования движения искусственных спутников Земли требуется использовать следующие основные шкалы времени.

Международное атомное время TAI (франц. Temps Atomique International) — представляет собой практическую реализацию равномерной шкалы времени, основанной на атомных часах. Отличается от земного времени на постоянное смещение в 32.184 с.

Всемирное координированное время UTC (англ. Coordinated Universal Time) — атомная равномерная шкала времени, отличающаяся от TAI на целое число секунд, называемых високосными секундами.

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

Система времени GPST глобальной навигационной спутниковой системы GPS берет начало своего отсчета с 6.I.1980 и представляет собой равномерную шкалу времени, т.е. не содержит скачков в виде високосных секунд. Эпоха определяется номером недели и номером секунды в неделе. Российская глобальная навигационная система ГЛОНАСС также имеет свою собственную шкалу времени. В пределах 100 нс она привязана к московскому времени UTC+3 и поддерживает эту привязку с помощью водородных мазеров в центральном синхронизаторе ГЛОНАСС. В связи с этим шкала подвержена влиянию високосных секунд.

В эфемеридах навигационных спутников и массиве наблюдений сигналов ГНСС используются, как правило, собственные шкалы времени навигационных систем. Однако при обработке комбинированных измерений (с учетом нескольких ГНСС), эфемериды и наблюдения чаще всего описывают в единой системе времени — шкале GPST.

Системы координат

Для определения местоположения КА в пространстве необходимо задать систему координат, в которой он движется. На сегодня для приложений околоземной космонавтики используются две ключевые системы координат: земная и небесная.

Международная земная система отсчета (англ. International Terrestrial Reference System, ITRS) — это геоцентрическая система координат с началом отсчета в центре масс Земли. Ось Z системы ITRS в пределах 30 мс дуги совпадает с условным международным началом (англ. Conventional International Origin, CIO), что является средним положением земного полюса по измерениям на интервале с 1900 по 1905 гг. [11]. Ось Х лежит в плоскости Гринвичского меридиана и экватора, ось Y лежит в плоскости меридиана с долготой 90° и дополняет систему осей до правой тройки.

Также существует Геоцентрическая небесная система отсчета (GCRS). Начало отсчета в системе координат GCRS находятся в центре масс Земли, однако ее оси не вращаются вместе с Землей. Они «закреплены» относительно удаленных радиоисточников (квазаров), а сама система отсчета движется вокруг барицентра Солнечной системы.

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

МОДЕЛИРОВАНИЕ ИЗМЕРЕНИЙ ГЛОБАЛЬНЫХ НАВИГАЦИОННЫХ СПУТНИКОВЫХ СИСТЕМ

В настоящий момент существует два основных типа наблюдений с навигационных спутников (НС): кодовые и фазовые. Они различаются способом получения измерений и величиной ошибки. Модель наблюдений для кодовых и фазовых измерений, соответственно, выглядит следующим образом:

$\left.\begin{array}{c}{P}_l^h={\varrho }_l^h+c\cdot \delta t_l-c\cdot \delta t^h+a^h+a_l+I+\\ +c\cdot \delta t_{rel}+HC{D}_l+HC{D}^h+IF{B}^h+\\ +ISB+{ϵ}_P,\\ L_l^h={\varrho }_l^h+c\cdot \delta t_l-c\cdot \delta t^h+a^h+a_l-I+\\ +c\cdot \delta t_{rel}+HP{D}_l+HP{D}^h+IF{B}^h+ISB+\\ +\lambda N+\delta {\varrho }_{wind-up}+{ϵ}_L,\end{array}\right\}$ (1)

где $P_l^h$ — кодовое измерение (псевдодальность); $L_l^h$ — фазовое измерение (псевдофаза); ${\varrho }_l^h$ — геометрическое расстояние между центрами масс НС в момент передачи сигнала и НКА в момент приема; $a^h$ — смещение фазового центра антенны НС относительно центра масс; $a_l$ — смещение фазового центра антенны НКА относительно центра масс; $I$ — влияние ионосферы; c — скорость света в вакууме; $\delta t_l$ — влияние ошибки часов НКА; $\delta t^h$ — влияние ошибки часов НС; $\delta t_{rel}$ — влияние релятивизма; $HC{D}_l$ — аппаратурная кодовая задержка НКА (англ. Hardware Code Delay); $HC{D}^h$ — аппаратурная кодовая задержка НС; $HP{D}_l$ — аппаратурная фазовая задержка НКА (Hardware Phase Delay); $HP{D}^h$ — аппаратурная фазовая задержка НС; $IF{B}^h$ — смещение из-за расхождения частот различных НС (англ.Inter-Frequency Bias); $ISB$ — смещение из-за расхождения шкал времени различных ГНСС (Inter-System Bias); N — целочисленная неоднозначность псевдофазы; $\lambda$ — длина волны несущей сигнала; $\delta {\varrho }_{wind-up}$ — поправка, связанная со взаимной ориентацией антенн НС и НКА (эффект набега фазы); $ϵ$ — прочие немоделируемые ошибки, включая эффект многолучевости. Все эти эффекты имеют разный вклад в модель наблюдений, от нескольких сантиметров до сотен километров. Поэтому для получения высокоточных координат приемника необходимо учитывать даже самые незначительные из них. Далее приведены формулы и рекомендации для расчета каждого из этих эффектов.

Ошибка бортовых часов

Наблюдения ГНСС, доступные наземному пользователю, привязываются к шкале навигационной системы, используемой в приемнике НКА. При этом ошибки часов НКА обычно оцениваются заранее и приводятся вместе с массивом наблюдений. В таком случае пользователю остается лишь учесть эту поправку при расчете модели наблюдений, как показано в формуле (1). Если же ошибки часов НКА неизвестны, то следует включить их в число уточняемых параметров, наряду с вектором состояния НКА. Моделирование ошибки часов НКА при этом осуществляется по следующей формуле [12]:

$\begin{array}{c}\delta t_{l,i}=a_0+a_1\left(t_i-t_0\right)+\dots +a_m{\left(t_i-t_0\right)}^m,\end{array}$ (2)

где $t_i$ — эпоха, на которую вычисляется поправка часов относительно эпохи $t_0$; a — полиномиальные коэффициенты. В данном выражении не учитывается периодическое изменение ошибки часов, «вызванное» Общей теорией относительности. Речь об этом пойдет в этом разделе ниже. Коэффициенты $a_0$ и $a_1$ называются отклонением фазы и частоты часов соответственно. Именно эти коэффициенты чаще всего принято уточнять, пренебрегая остальными порядками в выражении (2). Моделирование прогнозирования этих коэффициентов осуществляется с помощью белого гауссовского шума.

Ошибки часов НС содержатся в специальных каталогах, предоставляемых международной службой ГНСС (англ. International GNSS Service — IGS, https://igs.org) и российским информационно-аналитическим центром ГЛОНАСС (https://glonass-iac.ru) в формате SP3 и CLK. Получение ошибки часов из каталога может осуществляться путем интерполирования. Для получения высокоточных результатов следует брать каталоги ошибок часов с наименьшим шагом по времени.

Решение светового уравнения

Величина ${\varrho }_l^h$ в формуле (1) представляет собой геометрическое расстояние между центрами масс НС в момент времени передачи сигнала и НКА в момент времени приема сигнала. На практике известен только момент времени приема, а момент передачи сигнала необходимо вычислять. Отсюда возникает следующее уравнение, названное световым из-за задержки светового сигнала:

$\begin{array}{c}{\varrho }_l^h\left(t_l,t_h\right)={\varrho }_l^h\left(t_l\right)+\underset{0}{\overset{\Delta t}{\int }}{\dot{\varrho }}_l^{h}dt,\end{array}$ (3)

где ${\varrho }_l^h\left(t_l,t_h\right)$ — геометрическое расстояние между центрами масс НС в момент передачи сигнала и НКА в момент приема; ${\varrho }_l^h\left(t_l\right)$ — геометрическое расстояние между центрами масс НС и НКА в момент приема; ${\dot{\varrho }}_l^h$ — производная по времени геометрического расстояния между центрами масс НС и НКА; $∆$t — время полета сигнала от НС к НКА. В левой части уравнения находится геометрическое расстояние, участвующее в модели наблюдений. Однако вычислить его затруднительно, т.к. второе слагаемое правой части остается неизвестным. Поэтому существует специальная процедура, называемая определением световой задержки сигнала, которая позволяет определить момент времени передачи сигнала с НС и его координаты в этот момент. При этом время приема сигнала определяется в шкале времени навигационной системы, так же, как и координаты НС.

Решение осуществляется с помощью итеративного алгоритма, схема которого приведена на рис. 1. На вход алгоритма подается измерение псевдодальности. Рассчитывается задержка сигнала по формуле:

$\begin{array}{c}\begin{array}{r}\delta t_{delay}=\frac{R_{start}}{c},\end{array}\end{array}$ (4)

где $R_{start}$ — наблюдение; c — скорость света в вакууме.

Рис. 1. Процедура вычисления световой задержки ГНСС сигнала для НКА. Наверху изображен ГНСС-спутник в разные моменты времени. Положение НС определяется путем численного интегрирования дифференциальных уравнений движения (propagation) с известными начальными условиями в некоторый момент времени tpast. Внизу изображен низкоорбитальный космический аппарат в момент времени trec, где trec – время приема сигнала. Вычисляется расстояние R0, задержка сигнала $\delta t_{delay}^1$ и определяется положение НС на момент времени $t_{tr}^{est\mathrm{,1}}=\left(t_{rec}-\delta t_{delay}^1\right)$, где $t_{tr}^{est\mathrm{,1}}$ – потенциальное время отправки сигнала. Данная процедура повторяется итеративно, пока не будет соблюдено одно из условий. Либо две последние вычисленные задержки, либо два последних расстояния должны совпадать в пределах малой величины e. После этого можно считать, что последний вычисленный момент времени $t_{tr}^{est}$ есть оценка истинного момента времени отправки сигнала $t_{tr}^{true}$.

Чтобы найти время передачи сигнала, задержка вычитается из времени приема сигнала. На момент передачи сигнала находятся координаты НС по его эфемеридам. После этого вычисляется геометрическое расстояние между НС в момент передачи и НКА в момент приема сигнала. Эта модельная величина сравнивается с измерением, поступившим на вход итеративного алгоритма $R_{start}$. Если эти два расстояния совпадают в пределах достаточно малой величины, то алгоритм завершается. Если же нет, то найденная модельная псевдодальность поступает на вход алгоритма, в уравнение (4), и процедура повторяется до тех пор, пока начальное и конечное геометрические расстояния не совпадут в пределах достаточно малой величины. После завершения алгоритма становится известен истинный момент передачи сигнала и координаты НС на эту эпоху.

Смещение фазового центра антенны спутника глобальной навигационной спутниковой системы

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

$\left.\begin{array}{l}{a}^h=\left(D^h\cdot \left(PC{O}_{GCRF}+PC{V}_{GCRF}\left(\theta \right)\right)\right.,\\ D^h=\frac{r^h-r_l}{\left|r^h-r_l\right|},\end{array}\right\}$ (5)

где $r_l$ и $r^h$ — вектора положения НКА и НС в системе GCRF. Вектор PCO (англ. Phase Center Offset) обозначает смещение номинального центра антенны относительно центра масс космического аппарата. Вектор PCV (англ. Phase Center Variation) представляет собой смещение фактического фазового центра антенны относительно номинального в зависимости от угла $\theta$ — угла между вектором НС – НКА и надиром.

Вектора PCO для НС предоставляются Международной службой ГНСС (IGS) в виде файла формата ANTEX с расширением «.atx». Вектор PCO приведен для каждого существовавшего когда-либо НС и находится в системе координат, жестко связанной со спутником. В этой системе координат ось Z направлена в центр Земли, ось X находится в плоскости, проходящей через ось Z и Солнце, а ось Y дополняет правую тройку векторов. В этом же файле находятся параметры для расчета вектора PCV. Направление вектора PCV совпадает с направлением вектора НС – НКА, а длина вектора рассчитывается в зависимости от угла отклонения от надира.

Поскольку вектора PCO и PCV изначально предоставляются в системе координат спутника SRF (англ. Satellite Reference Frame) и зависят от его ориентации, необходимо преобразовать их в инерциальную систему координат GCRF. Сделать это не составляет труда, поскольку ориентация НС относительно стабильна. Преобразования производятся согласно следующим выражениям.

Вначале вычисляются единичные векторы системы координат спутника в системе GCRF:

$\begin{array}{c}\left.\begin{array}{r}\begin{array}{l}{e}_z=-\frac{r^h}{\left|r^h\right|},\\ s=\frac{r_{\odot }-r^h}{\left|r_{\odot }-r^h\right|},\\ e_y=\frac{\left[e_z\times s\right]}{\left|\left[e_z\times s\right]\right|},\\ e_x=\left[e_y\times e_z\right],\end{array}\end{array} \right\}\end{array}$ (6)

где $r^h$ — вектор положения НС в системе GCRF, $r_{\odot }$ — вектор положения Солнца в системе GCRF. Затем из них составляется матрица поворота из системы SRF в систему GCRF:

$\begin{array}{c}\begin{array}{r}\left.\begin{array}{l}M=\left(\begin{array}{ccc}{e}_x^1& e_y^1& e_z^1\\ e_x^2& e_y^2& e_z^2\\ e_x^3& e_y^3& e_z^3\end{array}\right),\\ PC{O}_{GCRF}=M\cdot PC{O}_{SRF}.\end{array}\right\}\end{array} \end{array}$ (7)

Смещение фазового центра антенны низкоорбитального космического аппарата и эффекты многолучевости

Смещение фазового центра антенны НКА находится схожим образом:

$\begin{array}{c}\begin{array}{r}\left.\begin{array}{l}{a}_l=(D_l\cdot \left(PC{O}_{GCRF}+PC{V}_{GCRF}\left(\theta \right)\right),\\ D_l=\frac{r_l-r^h}{\left|r_l-r^h\right|}.\end{array}\right\}\end{array}\end{array}$ (8)

Однако в случае с НКА ориентация может изменяться непредсказуемым образом. Поэтому при обработке измерений ГНСС необходимо также иметь файлы ориентации НКА. Чаще всего они представляются таблицей по времени кватернионов поворота локальной системы координат НКА относительно системы GCRF. В случае спутников проекта GRACE система SRF определяется следующим образом. Ось X направлена на спутник-компаньон, ось Z направлена в надир, ось Y дополняет систему до правой тройки векторов.

Преобразование из системы координат спутника SRF в систему координат GCRF происходит с помощью файлов ориентации НКА. Перевод в инерциальную систему координат вектора PCO осуществляется следующим образом:

$\begin{array}{c}\begin{array}{c}PC{O}_{GCRF}=Q\otimes PC{O}_{SRF}=\\ =PC{O}_{SRF}+2Q_0\left[Q_{123}\times PC{O}_{SRF}\right]+\\ +2\left[Q_{123}\times \left[Q_{123}\times PC{O}_{SRF}\right]\right],\end{array}\end{array}$ (9)

где $Q=\left\{Q_0,Q_1,Q_2,Q_3\right\}=\left\{Q_0,Q_{123}\right\}$ — кватернион ориентации; $\otimes$ — кватернионное умножение.

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

Также к шуму измерений можно отнести эффект многолучевости. Данный эффект вызван непрямым распространением сигнала между НС и НКА. В космическом пространстве это может происходить только из-за переотражения сигнала от солнечных панелей КА или других его частей. При конструировании космического аппарата, для принимающей антенны, как правило, исключаются такие эффекты.

Как видно, из уравнений (5)–(9), для определения эффектов смещения фазовых центров антенн необходимо знать кватернионы ориентации КА. Поэтому для высокоточного определения параметров движения КА важно совершенствовать методы определения ориентации КА [13].

Задержка распространения сигнала в ионосфере

Ионосферу образуют верхние слои земной атмосферы, в которой газы частично ионизованы под влиянием ультрафиолетового и рентгеновского солнечного излучения. При высоких частотах основную роль в распространении радиоволн играют свободные электроны. Число электронов в кубическом метре (т.е. плотность электронов) изменяется сложным образом с высотой над поверхностью Земли. Распределение плотности электронов зависит от времени суток, времени года, уровня солнечной активности. Ионосферной рефракцией называется разница между измеряемым и истинным расстоянием. Рефракция в ионосфере сводится к изменению длины пути луча в атмосфере. Незнание количества свободных электронов на пути волны в ионосфере определяет ошибки вычисления задержки сигнала. Поскольку орбиты НКА имеют высоты от 160 до 2000 км, а максимум плотности ионосферы приходится на высоты от 250 до 400 км, то для сигналов ГНСС, принимаемых НКА, особенно важно учитывать ионосферную поправку за рефракцию.

При определении орбиты с помощью ГНСС часто используется специальная комбинация измерений, полученных на двух разных частотах. Составить ее можно следующим образом:

$\begin{array}{c}\left.\begin{array}{l}\text{CIF}=\frac{2}{f_1^2-f_2^2}\left(f_1^2{P}_1-f_2^2{P}_2\right),\\ \text{PIF}=\frac{2}{f_1^2-f_2^2}\left(f_1^2{L}_1-f_2^2{L}_2\right),\end{array}\right\}\end{array}$ (10)

где PIF— комбинация фазовых измерений $L_1$ и $L_2$ на частотах $f_1$ и $f_2$, свободная от ионосферы (Phase Ionosphere Free); CIF— аналогичная комбинация кодовых измерений (Code Ionosphere Free). Данная комбинация позволяет избавиться от влияния ионосферного члена первого порядка. Остаются лишь члены более высокого порядка, максимальная величина которых может составлять до 1 см. Если считать дисперсии погрешностей одиночных измерений на обеих частотах равными, то среднее квадратическое отклонение (СКО) погрешности безионосферной комбинации будет в три раза больше СКО погрешности исходных измерений. С учетом малого шума фазовых измерений данное увеличение погрешности комбинации PIF не играет существенной роли. Именно поэтому эту комбинацию чаще всего используют для устранения влияния ионосферы. Однако стоит учесть, что данная комбинация содержит в себе влияние фазовых неоднозначностей. Причем в этой комбинации неоднозначности на двух частотах выражаются единым нецелочисленным членом. Опуская влияние других поправок, выражение для него можно записать следующим образом:

$\begin{array}{c}\begin{array}{r}\left.\begin{array}{l}\text{PIF}={\varrho }_l^h+B,\\ B=\frac{1}{f_1^2-f_2^2}\left(f_1^2{\lambda }_1{n}_1-f_2^2{\lambda }_2{n}_2\right).\end{array}\right\}\end{array}\end{array}$ (11)

Для работы в режиме одиночных измерений, а также для алгоритмов определения фазовых скачков может потребоваться вычисление ионосферной поправки. Она рассчитывается согласно ионосферным картам, полученным центром анализа Международной службы ГНСС. Данные карты представляют общее содержание электронов в ионосфере (TEC, Total Electron Content) для сигнала, испускаемого в зенит из точки на поверхности Земли (ZTEC, Zenith TEC). Данный формат имеет название IOnosphere map Exchange format (IONEX). Описание данного формата позволяет интерполировать значения TEC по сетке карты и по времени для нахождения наиболее близкого значения ZTEC. Однако в реальности требуется находить общее содержание электронов на пути НС – НКА (STEC, Slant TEC). Для этого применяется следующая формула:

$\begin{array}{c}STEC=M\left(z\right)\times ZTEC,\end{array}$ (12)

где $M\left(z\right)$ — картирующая функция, зависящая от зенитного угла z (рис. 2).

Рис. 2. Схема учета ионосферной задержки для системы НКА – НС.

В работе [14] предлагается следующая схема по расчету картирующей функции для случая, при котором приемник находится на низкоорбитальном космическом аппарате. В ней используются две модели картирующих функций для $z<45$ и $z\ge 45$. Модель тонкого ионосферного слоя, которая подходит для $z\ge 45$, описывается формулой:

$\begin{array}{c}\begin{array}{r}{M}_1\left(z\right)={\left(\sqrt{1-{\left(\frac{R_l\text{sin}z}{R_{I\mathrm{,1}}}\right)}^2}\right)}^{-1},\end{array}\end{array}$ (13)

где $R_l$ — расстояние от центра Земли до НКА; $R_{I\mathrm{,1}}$ — расстояние от центра Земли до эффективной высоты ионосферы, согласно выражению (15). В случае $z<45$ картирующая функция выглядит так:

$\begin{array}{c}\begin{array}{r}\begin{array}{c}{M}_2\left(z\right)=\left(1+\frac{R_{I\mathrm{,2}}}{R_l}\right)\times \\ \times {\left(\cos z+\sqrt{{\left(\frac{R_{I\mathrm{,2}}}{R_l}\right)}^2-{\left(\sin z\right)}^2}\right)}^{-1}.\end{array}\end{array}\end{array}$ (14)

Эффективная высота ионосферы для этих двух моделей ($h_{I\mathrm{,1}}$ и $h_{I\mathrm{,2}}$) рассчитывается так:

$\left.\begin{array}{l}{R}_{I\mathrm{,1}}=R_{\oplus }+h_{I\mathrm{,1}},\\ R_{I\mathrm{,2}}=R_{\oplus }+h_{I\mathrm{,2}},\\ h_{I\mathrm{,1}}\left[км\right]=\\ =\left(0.0005F_{107}+1.73\right)h_l\left[км\right]-\\ -3.16F_{107}+\mathrm{458,}\\ h_{I\mathrm{,2}}\left[км\right]=\\ =\left(0.0027F_{107}+1.79\right)h_l\left[км\right]-\\ -5.52F_{107}+\mathrm{1350,}\end{array}\right\}$ (15)

где $R_{\oplus }$ — средний радиус Земли; $h_l$ — высота орбиты НКА (разница радиус-вектора НКА и среднего радиуса Земли); $F_{107}$ — индекс солнечной активности, который находится, например, из файла космической погоды Центра космических стандартов и инноваций (https://ai-solutions.com/_help_Files/cssi_space_weather_file.htm).

Величина $\text{ZTEC}$ определяется с помощью ионосферных карт IONEX. Для ее расчета находятся сферические координаты точки I, т.е. точки пересечения эффективного ионосферного слоя и линии НКА – НС. Стоит также учесть, что с помощью ионосферных карт можно найти общее содержание электронов для эффективной высоты ионосферы в 450 км. Однако на практике оказывается, что орбита НКА зачастую находится выше пика ионосферной концентрации. Чтобы учесть данный эффект, необходимо оценить поправочный коэффициент a, который представляет из себя отношение вертикального содержания электронов, начиная с высоты эффективного ионосферного слоя к общему вертикальному содержанию электронов, начиная с поверхности Земли. Такой коэффициент можно вычислить следующим образом:

$\begin{array}{c}\begin{array}{r}\alpha =\frac{e-\exp \left(1-\exp \left(-\frac{h_I-h_0}{H}\right)\right)}{e-\exp \left(1-\exp \left(\frac{h_0}{H}\right)\right)},\end{array}\end{array}$ (16)

где $h_0$ — высота пика ионосферной концентрации электронов; $h_I$ — эффективная высота ионосферы, которая определяется выражением (15); H — масштабный фактор, характеризующий кривую зависимости концентрации от высоты. Данные величины можно определять из различных атмосферных моделей, например, из модели ионосферы IRI2016 [15]. Однако эти данные не обладают достаточно высокой точностью, и на практике коэффициент a чаще просто оценивается в процессе уточнения орбиты наряду с другими параметрами.

Релятивистский эффект

Релятивистские эффекты можно разделить на слагаемые:

$\begin{array}{c}c\delta t_{rel}=-\delta {\rho }_{h,clock}+\delta {\rho }_{l,clock}+\delta {\rho }_{grav}.\end{array}$ (17)

Первое слагаемое здесь описывает релятивистские эффекты ухода бортовых часов на НС относительно опорной шкалы времени, например, GPST. Общая теория относительности предсказывает эффект ускорения часов на высоте орбит ГНСС, а специальная теория относительности — эффект замедления часов. Оба этих эффекта уже учтены для часов ГНСС путем смещения частоты часов. Однако из-за эллиптичности орбит НС необходимо вводить дополнительную поправку. Эта поправка может быть описана следующей формулой:

$\begin{array}{c}\begin{array}{r}\delta {\rho }_{h,clock}=2c^{-1}\left(\overrightarrow{r},\overrightarrow{v}\right),\end{array}\end{array}$ (18)

где $\overrightarrow{r}$ — радиус-вектор спутника; $\overrightarrow{v}$ — вектор скорости спутника; c — скорость света в вакууме. Данное выражение было получено в рамках Кеплеровской модели сил [16].

Второе слагаемое формулы (17) описывает поправку за уход бротовых часов НКА относительно опорной шкалы времени. Для часов на низкой околоземной орбите оказывается существенным влияние нецентральной части гравитационного поля Земли, поэтому выражение (18) оказалось бы неприменимым. В общем случае поправка для НКА представляется в следующем виде [13]:

$\begin{array}{c}\begin{array}{r}\delta {\rho }_{l,clock}=\int \left(L_G\cdot c-\frac{v^2}{2c}-\frac{U}{c}\right)dt,\end{array}\end{array}$ (19)

где U — гравитационный потенциал Земли в точке КА; v — модуль скорости КА; $L_G$ — безразмерный коэффициент, равный $6.969290134\cdot {10}^{-10}$, характеризующий масштаб шкалы времени. Численное интегрирование уравнения (19) можно проводить вместе с интегрированием уравнений движения НКА на заданный интервал времени. Размер релятивистской поправки может достигать 2 нс. Данная поправка обладает периодическим характером изменения с главной гармоникой, равной периоду орбиты КА, и вспомогательной — в половину периода орбиты.

Из-за эффекта задержки распространения света при прохождении вблизи массивного тела, известного также как эффект Шапиро [17], возникает третье слагаемое в формуле (17):

$\begin{array}{c}\begin{array}{r}\delta {\rho }_{grav}=\frac{2{\mu }_{\oplus }}{c^2}\ln \left(\frac{r_h+r_l+\rho }{r_h+r_l-\rho }\right),\end{array}\end{array}$ (20)

где ${\mu }_{\oplus }$ — гравитационный параметр Земли; $\frac{2{\mu }_{\oplus }}{c^2}$ — гравитационный радиус Шваршильда; $r_h$ — геоцентрическое расстояние до НС; $r_l$ — геоцентрическое расстояние до НКА; p — расстояние между НКА и НС.

Эффект набега фазы

Фазовые измерения зависят от взаимной ориентации антенн приемника и навигационного спутника и от положения в пространстве их линии видимости. Суть эффекта можно описать следующим образом. Если представить, что расстояние между КА и НС зафиксировано, а НС поворачивается вокруг своей оси Z на 360°, то измерение фазы при этом изменится на один цикл, соответствующий 19 см для частоты L1. Это происходит из-за правосторонней поляризации сигнала, который излучают НС. Данный эффект в иностранной литературе принято называть wind-up эффектом. Его можно рассчитать по следующим формулам:

$\begin{array}{c}\begin{array}{r}\left.\begin{array}{l}\delta \rho =\delta \varphi +2\pi N,\\ \delta \varphi =\frac{1}{2\pi }\varphi \lambda ,\\ N=n\mathrm{int}\left(\frac{\delta {\rho }_{prev}-\delta \varphi }{2\pi }\right),\\ \varphi =sign\left(\zeta \right)\arccos \left(\frac{\left(D_h\cdot D_l\right)}{\left|D_h\right|\left|D_l\right|}\right),\\ \zeta =k\left[D_h\times D_l\right],\\ D_h=x_h-k\left(k\cdot x_h\right)-\left[k\times y_h\right],\\ D_l=x_l-k\left(k\cdot x_l\right)+\left[k\times y_l\right],\\ k=\frac{r_l-r_h}{\left|r_l-r_h\right|},\end{array}\right\}\end{array}\end{array}$ (21)

где $r_l$ и $r_h$ — радиус-векторы НКА и НС соответственно; $\left(x_h,y_h,z_h\right)$ — единичные вектора локальной системы координат НС, речь о которой шла в предыдущих разделах; $\left(x_l,y_l,z_l\right)$ — единичные вектора системы координат антенны НКА (ARF); $\lambda$ — длина волны несущего сигнала; $\delta {\rho }_{prev}$ — поправка на предыдущий момент времени. На первую эпоху N полагается равной нулю. Система координат ARF зависит от конкретного проекта НКА. В данной работе, на примере спутника GRACE-A, система координат НКА задается следующим образом:

$\begin{array}{c}\begin{array}{r}\left.\begin{array}{l}{X}_{ARF}=-Y_{SRF},\\ Y_{ARF}=-X_{SRF},\\ Z_{ARF}=-Z_{SRF},\end{array}\right\}\end{array}\end{array}$ (22)

где $\text{SRF}$ — система координат, описанная в разделе «Смещение фазового центра антенны низкоорбитального КА и эффекты многолучевости».

Аппаратурные задержки

Поскольку передача и прием сигнала ГНСС ведутся на двух, трех частотах одновременно, аппаратура передающего и принимающего спутников должна обеспечивать абсолютную синхронизацию в обработке этих сигналов. Однако на практике добиться такой синхронности очень сложно. Отсюда возникают аппаратурные задержки HCD (Hardware Code Delay), также именуемые DCB (Differential Code Bias). Эти задержки исключаются при использовании безионосферной комбинации измерений (10). В режиме одиночных измерений существует ряд методов по их определению [18]. Для навигационных спутников наиболее удобным методом нахождения данных поправок представляется использование модели ионосферы в формате IONEX, предоставляемой IGS. Величины поправок DCB определяются вместе с величинами TEC ионосферы и предоставляются для каждого НС в наносекундах. Значения аппаратурных задержек для двух рабочих частот НС вычисляются по формулам:

$\begin{array}{c}\begin{array}{r}\left.\begin{array}{l}{\text{DCB}}^h\left(L1\right)={\text{HCD}}^h\left(L1\right)=\frac{{\nu }_2^2}{{\nu }_1^2-{\nu }_2^2}c\Delta b,\\ {\text{DCB}}^h\left(L2\right)={\text{HCD}}^h\left(L2\right)=\frac{{\nu }_1^2}{{\nu }_1^2-{\nu }_2^2}c\Delta b,\end{array}\right\}\end{array}\end{array}$ (23)

где $L1={\nu }_1$ и $L2={\nu }_2$ — две разные частоты измерений; $∆$b — величина поправки из файла IONEX в наносекундах. Аналогичные задержки возникают и для фазовых измерений, в этом случае они именуются как HPD (Hardware Phase Delay). В последнее время для обозначения всех видов задержек для всех различных типов наблюдений используется термин OSB (Observed Specific Bias), а служба IGS предоставляет специальные файлы для учета OSBзадержек.

Если аппаратурную задержку приемника требуется оценить отдельно от шума, то это можно сделать по формулам из работы [19]:

$\begin{array}{c}\begin{array}{r}\left.\begin{array}{l}c\Delta b_l={\left(\sum _{}_{i=1}^{narc}\sum _{}_{j=1}^{nobs}{w}_{ij}\right)}^{-1}\sum _{}_{i=1}^{narc}\sum _{}_{j=1}^{nobs}{w}_{ij}\times \\ \times \left[L_{\mathrm{1,}ij}-L_{\mathrm{2,}ij}-\frac{STEC}{F}-\Delta N_i-c\Delta b_{h,i}\right],\\ w_{ij}={\left(\frac{{\sigma }_{TEC,ij}}{F}\right)}^{-2},\\ F=\frac{f_1^2{f}_2^2}{40.3(f_1^2-f_2^2)},\\ \Delta N_i=\frac{1}{n_{obs}}\sum _{}_{j=1}^{nobs}(P_{\mathrm{2,}ij}-P_{\mathrm{1,}ij}-L_{\mathrm{1,}ij}+L_{\mathrm{2,}ij}),\end{array}\right\}\end{array}\end{array}$ (24)

где $∆$b_l — DCB-задержка приемника; $∆$b_h,i — DCB-задержка НС; $n_{arc}$ — количество дуг видимости каждого НС (временной интервал, на котором НС находится в зоне видимости НКА без перерыва); $n_{obs}$ — количество наблюдений в дуге видимости НС; $L_{\mathrm{1,}ij}$ и $L_{2,ij}$ — фазовые наблюдения на двух частотах (в метрах); $P_{\mathrm{1,}ij}$ и $P_{2,ij}$ — кодовые наблюдения на двух частотах (в метрах); ${\sigma }_{\text{TEC},ij}$ — стандартное отклонение величины TEC, которое находится в файле IONEX. Для спутника GRACE-A данная задержка составляет — 19.224 нс [20].

Межчастотное смещение

Сигналы GPS основаны на методе множественного доступа с кодовым разделением (CDMA-сигналы). Каналы передачи в таком случае имеют общую полосу частот, но разные кодирующие последовательности. Однако наблюдения несущей фазы и псевдодальности ГЛОНАСС основаны на методе множественного доступа с разделением каналов по частоте (FDMA-сигналы). Каждый спутник ГЛОНАСС имеет одинаковые кодирующие последовательности, но передает на разных частотах, соответствующих своему частотному каналу. Причем канал может быть одинаковым у двух НС, находящихся с противоположных сторон Земли. Частоты несущих для конкретного спутника определяются следующими выражениями:

$\begin{array}{c}\begin{array}{r}\left.\begin{array}{l}G1\left[МГц\right]=1602+0.5625n,\\ G2\left[МГц\right]=1246+0.4375n,\\ n\in \left[-\mathrm{7,6}\right],n\in \mathbb{Z},\end{array}\right\}\end{array}\end{array}$ (25)

где $G1$ и $G2$ — две основные частоты ГЛОНАСС, на которых излучает каждый спутник; n — номер частотного канала. Каналы с номерами –5, –6 и, в некоторых случаях 7, используются для НС, проходящих тестирование и отладку. На практике между измерениями от спутников с разделением каналов по частоте возникают IFB-смещения (Inter-Frequency Bias), связанные с аппаратурными задержками при обработке различных частот. Эти задержки могут достигать 4 нс [21], поэтому правильное моделирование поправок IFB — это критически важный этап при обработке ГЛОНАСС-измерений. В работе [21] был проведен анализ методов оценки поправки IFB. Среди методов предлагалось:

• пренебрежение поправкой IFB и включение ее в оценку ошибки часов приемника;
• оценка IFB для спутников с одинаковым частотным каналом;
• ценка IFB для каждого спутника ГЛОНАСС.

Анализ показал, что поправки IFB могут не совпадать строго для одного и того же частотного канала. Поэтому предпочтительным оказывается последний метод, в котором выбирается опорный спутник (желательно с нулевым частотным каналом) и в список уточняемых параметров дополнительно вводится $n-1$ (где n — число наблюдаемых ГЛОНАСС спутников) число параметров $∆$IFB, описывающих смещение сигнала во времени относительно опорного спутника. Поправка IFB при этом считается постоянной на всем временном интервале для каждого спутника.

Межсистемное смещение

При обработке данных, полученных с помощью разных глобальных навигационных спутниковых систем, возникает необходимость в учете параметра межсистемного смещения (ISB). Это смещение связано с особенностями приемников различных типов навигационных систем. Но и для одного и того же типа приемника ISB-параметр может незначительно отличаться из-за задержек, вызванных антенной, кабелем или тепловыми эффектами в бортовой аппаратуре. В настоящее время в большинстве ГНСС-приложений поправка ISB оценивается либо как постоянная величина, либо как кусочно-непрерывная на интервале в 30 мин. Также можно моделировать данное смещение белым гауссовским шумом [22]. Поправка ISB для каждой пары спутниковых систем вносится в число уточняемых параметров и оценивается вместе с вектором состояния НКА.

УТОЧНЕНИЕ ОРБИТЫ ПО КОМБИНИРОВАННЫМ НАБЛЮДЕНИЯМ

Ниже приводятся результаты эксперимента определения орбиты по совместным модельным наблюдениям GPS/ГЛОНАСС с целью демонстрации разработанной модели наблюдений. Для этого было разработано специальное программное обеспечение (ПО) LOIS [23], позволяющее уточнеть орбиту космического аппарата по данным глобальных навигационных спутниковых систем. В основе данного ПО лежит метод дифференциального уточнения орбиты — обобщенный фильтр Калмана. В этом методе измерения ГНСС поступают на вход фильтра последовательно, со временем. Исходя из априорной оценки вектора состояния КА составляется модель наблюдений, как это описано в предыдущем разделе. Затем на основе сигнала рассогласования между наблюдениями и моделью вычисляется апостериорная оценка вектора состояния. Эта оценка используется в качестве априорной на следующую эпоху с измерениями ГНСС. Затем цикл повторяется заново, пока не закончатся наблюдения.

Эксперимент по уточнению орбиты проводился на примере синтетических GPS/ГЛОНАСС-измерений геодезического спутника GRACE-A. Проект GRACE 2002–2017 (англ. Gravity Recovery And Climate Experiment) — совместная спутниковая миссия, целью которой было изучение гравитационного поля Земли и его временных вариаций, Национального управления по аэронавтике и исследованию космического пространства (NASA) и Германского центра авиации и космонавтики. Проект состоял из двух спутников, расположенных на практически одной и той же низкой околоземной орбите и разнесенных на ней на расстояние в 220 км. Околополярная орбита спутников GRACE-A и GRACE-B имела наклонение 89.6° и высоту, колеблющуюся от 497 до 521 км. Модель движения НКА GRACE-A, использующаяся в данном эксперименте, представлена в табл. 2.

Поскольку разработанную модель измерений необходимо было проверить и сравнить с аналогами, а получение реальных комбинированных измерений GPS/ГЛОНАСС затруднительно, для моделирования наблюдений было применено зарубежное ПО Bernese [24], которое использует реальные данные космического проекта. Данное ПО позволяет создавать синтетические измерения с учетом различных эффектов модели движения и наблюдений, а также с использованием различных навигационных спутниковых систем. ПО Bernese создавалось группой докторов наук Бернского университета и прошло проверку временем, поэтому сравнение модели измерений из данной работы и метода моделирования наблюдений в ПО Bernese оказывается крайне полезным. При моделировании наблюдений в ПО Bernese использовалась модель движения КА, аналогичная описанной выше. Тем не менее допускается наличие рассогласования в численном интегрировании уравнений движения на уровне 1…2 мм на интервале в 10 с. Для уточнения орбиты использовались фазовые измерения, содержащие шум со среднеквадратическим отклонением 2 мм. Измерения моделировались с частотой 0.1 Гц. Также нужно уточнить, что ПО Bernese не способно моделировать эффект начальных фазовых целочисленных неоднозначностей (только эффект пропуска циклов), поэтому с помощью ПО LOIS в наблюдения были дополнительно введены случайные целочисленные параметры, постоянные для каждой дуги видимости НС. Таким образом, вектор уточняемых параметров в данном эксперименте содержал трехмерные координаты и скорости КА, а также параметры фазовых неоднозначностей. Поправки IFB / ISB считались нулевыми. Вместо определения ионосферного коэффициента использовалась безионосферная комбинация фазовых наблюдений. Ошибка часов НКА GRACE-A считалась нулевой, поскольку доступные пользователю измерения уровня 1B уже скорректированы за эту поправку.

В качестве критерия для сравнения двух моделей измерений можно использовать оценку невязок между синтетическими наблюдениями и непосредственно моделью измерений, так называемых (O – C) невязок. Однако их минимизация не гарантирует получения надежного баллистико-навигационного решения, поскольку одной и той же псевдодальности могут соответствовать разные наборы координат НКА. В связи с этим, следует также обратить внимание на расхождение уточняемых параметров с истинными. В качестве истинной орбиты спутника GRACE-A в настоящем исследовании использовалась оценка орбиты Германского центра авиации и космонавтики. Кроме того, на основе данной оценки моделировались измерения в ПО Bernese.

На рис. 3 изображен график средних невязок (O – C). Под средней невязкой здесь понимается невязка, усредненная на каждую эпоху по всем НС. Таким образом одной эпохе соответствует одно число, характеризующее сигнал рассогласования модели измерений и наблюдений. Поскольку невязки могут иметь и положительные, и отрицательные значения, может случиться так, что отдельные невязки (O – C) будут иметь значительные величины, в то время как средняя невязка будет близка к нулю. Поэтому важно следить и за стандартным отклонением невязок на каждую отдельную эпоху. Средняя величина стандартного отклонения невязок (O – C) на временном интервале 1 сут. составляет 1 см. В некоторые эпохи невязки могут резко возрастать из-за появления новых НС с неизвестными параметрами фазовых неоднозначностей.

Рис. 3. Среднее невязок (O – C) эксперимента по определению орбиты спутника GRACE A по безионосферной комбинации синтетических фазовых измерений систем GPS/ГЛОНАСС.

На рис. 4 приведен график трехмерного отклонения от эталонного местоположения НКА (модуль вектора невязки положения). Средняя невязка местоположения составляет 18 мм. Величина невязок (O – C) и трехмерного положения может объясняться ошибками эфемерид НС и шумом измерений. Также не исключается рассогласование в моделях измерений ПО LOIS (обработка измерений), ПО Bernese (создание синтетических наблюдений) и Германского центра авиации и космонавтики, который предоставлял эталонную орбиту. Таким образом, сравнение разработанной модели с зарубежным программным обеспечением и реальной орбитой НКА показывает совпадение в пределах 1…2 см. Такая точность согласуется с другими зарубежными работами по уточнению орбит НКА с помощью ГНСС [25, 26].

Рис. 4. Невязка трехмерного положения НКА эксперимента по определению орбиты спутника GRACE A по безионосферной комбинации синтетических фазовых измерений систем GPS/ГЛОНАСС.

ЗАКЛЮЧЕНИЕ

В настоящей работе описан процесс построения расчетных значений комбинированных наблюдений глобальных навигационных спутниковых систем. Приводится описание основных формул и выражений для учета всевозможных эффектов модели наблюдений ГНСС. Описанная модель прошла тестирование в разработанном отечественном программном обеспечении для уточнения орбит космических аппаратов. Моделирование наблюдений прошло сравнение с зарубежным программным обеспечением в процессе эксперимента по уточнению орбиты космического аппарата GRACE-A. Среднее трехмерное отклонение местоположения НКА составило 18 мм.

Зарубежные работы по уточнению орбит НКА с помощью данных ГНСС опираются на одни и те же готовые программные обеспечения, не вдаваясь в подробности расчета тех или иных эффектов модели наблюдений. Описанная в данной работе модель позволяет исследователю самому включать или выключать различные эффекты, изучать их вклад в конкретном эксперименте, проводить высокоточное позиционирование космических аппаратов и валидацию результатов. Она будет крайне полезна для различных приложений, экспериментов, научных миссий и программных обеспечений, связанных с навигационными сигналами ГНСС.

КОНФЛИКТ ИНТЕРЕСОВ

Автор данной работы заявляет, что у него нет конфликта интересов.

Sobre autores

П. Запевалин

Autor responsável pela correspondência
Email: pav9981@yandex.ru
Rússia, Москва

Bibliografia

Arquivos suplementares

Arquivos suplementares

Ação

1. JATS XML

Baixar

2. Fig. 1. The procedure for calculating the light delay of the GNSS signal for the NCA. The GNSS satellite is shown at the top at different points in time. The position of the HC is determined by numerically integrating the differential equations of motion (propagation) with known initial conditions at some point in time tpast. The bottom shows a low–orbit spacecraft at time trec, where trec is the signal reception time. The distance R0 is calculated, the signal delay is calculated, and the position of the HC at the time is determined, where is the potential time of sending the signal. This procedure is repeated iteratively until one of the conditions is met. Either the last two calculated delays or the last two distances must coincide within a small value e. After that, we can assume that the last calculated moment in time is an estimate of the true moment when the signal was sent.

Baixar (39KB)

Metadados

3. Fig. 2. Ionospheric delay accounting scheme for the NCA –NS system.

Baixar (20KB)

Metadados

4. Fig. 3. The average of the residuals (O – C) of the experiment to determine the orbit of the GRACE A satellite using an ionospheric combination of synthetic phase measurements of GPS/GLONASS systems.

Baixar (41KB)

Metadados

5. Fig. 4. Discrepancy between the three-dimensional position of the spacecraft of the experiment to determine the orbit of the GRACE A satellite using an ionospheric combination of synthetic phase measurements of GPS/GLONASS systems.

Baixar (42KB)

Metadados