Апостериорная оценка точности определения вектора состояния земного наблюдателя по измерениям дальности и скорости системы космической навигации GPS

(The posteriori precision evaluation of the definition of user's state vector on measurements of pseudorange and dopler frequency for GPS satellite navigation system
Preprint, Inst. Appl. Math., the Russian Academy of Science)

Оглавление

Введение

1.      Назначение, состав и общая характеристика системы GPS

2.      Интерпретация измерений

3.      Вычисление вектора состояния НКА GPS по эфемеридным данным

4.      Алгоритм определения вектора состояния наблюдателя по измерениям
псевдодальности и псевдоскорости

5.      Статистическая оценка точности определения вектора  состояния

Литература

Введение

Широко распространенные в настоящее время GPS-приемники ориентированны в основном на информационное обеспечение потребителей, расположенных на поверхности Земли. Применение средств космической навигации в сложных технических системах, к которым относятся, например, бортовые автономные системы управления космических аппаратов, предъявляет более высокие требования к надежности работы аппаратуры потребителя. В Институте прикладной математики проводятся исследования, направленные на изучение возможности применения средств космической навигации в сложных  технических системах, включая борт КА. В рамках этих исследований выполнен первый этап работы, посвященный оценке точности определения вектора состояния земного наблюдателя по результатам обработки псевдодальности и псевдоскорости.

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

Для апостериорной оценки точности определения вектора состояния были использованы ежесекундные измерения на месячном интервале. Определены статистические характеристики для каждого дня этого интервала.

В первой части работы описана общая структура глобальной навигационной системы GPS, ее назначение и состав.

Во второй части настоящей работы описана интерпретация измерений псевдодальности и псевдоскорости, полученных с использованием  отечественного одночастотного GPS-приемника, который был изготовлен в РНИИ Космического приборостроения.

Третья часть посвящена расчету вектора состояния навигационных КА по данным, получаемых приемником.

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

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

Авторы благодарят В.А. Степаньянца за помощь в работе и В.В. Тюбалина за содействие в проведении измерений.

1. Назначение, состав и общая характеристика системы GPS

Глобальная спутниковая система GPS (Global Position System) Navstar предназначена для высокоточного определения пространственно-временных координат и составляющих скорости объектов-потребителей. Система разработана по заказу и находится под управлением МО (ВВС) США. Система GPS состоит из космического сегмента,  сегмента управления и сегмента потребителей.

Космический сегмент образован орбитальной группировкой, состоящей из 24 основных и 3 резервных навигационных КА (далее просто НКА). НКА распределены по шести плоскостям, которые разнесены по долготе на 60°. В каждой плоскости находится четыре, и, возможно, один резервный НКА, которые двигаются по круговым орбитам с наклонением 55° и с полуосью около 26,5 тыс. км. Период обращения НКА составляет 12 часов. Основными разработчиками и создателями космического сегмента являются Rockwell International Space System Division и Martin Marietta Astro Space Division.

      Рис. 1 Общий вид НКА GPS           Рис. 2 Развернутая навигационная система GPS

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

Рис. 3 Наземные станции слежения

 С помощью наземного сегмента управления осуществляются высокоточные измерения параметров орбит НКА, которые собираются и обрабатываются в Фальконе. Результатом обработки является информация об орбите, частотно временные поправки, ионосферные поправки. Полученная информация передается на борт НКА для последующей ретрансляции потребителю. Частота обновления ретрансляционной информации приблизительно раз в два часа.

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

2. Интерпретация измерений

Метод определения координат и скорости потребителя основан на измерениях дальности и радиальной скорости относительно НКА.

В настоящей работе были использованы измерения, полученные с использованием  отечественного одночастотного GPS-приемника, который был изготовлен в РНИИ Космического приборостроения. Отметим, что была накоплена месячная база ежесекундных измерений для неподвижного приемника.

Передатчики НКА GPS излучают два непрерывных сигнала на частотах L1 и L2. Одночастотный 12-канальный GPS-приемник может одновременно принимать сигналы  L1 на частоте   и модулировать общедоступный псевдослучайный C/A код (Coarse Acquisition Code)  не более чем от 12 НКА.

Несущая частота L1 состоит из двух компонентов, сдвинутых на  для удобства их разделения. Первая компонента модулируется двумя двоичными последовательностями (дальномерный псевдослучайный P-код  и информационная последовательность), складывающимися по модулю 2. Вторая - также модулируется двоичными последовательностями (дальномерным псевдослучайным C/A-кодом  и информационной последовательностью), складывающимися по модулю 2. Информационная последовательность передается со скоростью 50 бит/с и содержит навигационные сообщения, включающие данные об орбите НКА,  частотно-временные поправки по которым расчитывают сдвиг шкалы времени  i-го НКА от шкалы времени GPS-системы и поправки, связанные с задержкой распространения сигнала.

Принимая C/A-код, приемник измеряет временную задержку  прохождения сигналов синхронизированных между собой НКА и доплеровское смещение  частоты сигнала .

Так как один цикл передачи C/A-кода состоит из 1023 бит и повторяется 1000 раз в секунду, зона однозначного измерения составляет , где с - скорость света. Под псевдодальностью до i-го НКА называют величину  удовлетворяющую следующему соотношению:

                                           ,                                                  (1)

где  - целое число,  - фазовый сдвиг псевдошумовой последовательности C/A-кода, вызванный различием в синхронизации часов i-го НКА и GPS-системы.

Задача раскрытия обычно решается в приемнике с использованием регистрации фиксированных меток времени в цикле передачи информационной последовательности.

В качестве псевдоскорости для i-го НКА приемник выдает следующее измерение

                                                 .                                                                  (2)

В связи с нестабильностью кварцевого генератора  приемника, модуляция псевдошумового C/A-кода начинается каждый раз заново для каждой группы одновременных измерений.

Отметим, что измерения псевдодальности и псевдоскорости относятся, вообще говоря, к разным временам. Для измерения смещения доплеровской частоты   требуется некоторый интервал времени, называемый временем накопления сигнала. Обычно время регистрации сигнала привязывают к середине этого интервала. В используемом приемнике время накопления сигнала фиксировано и составляет 1.4 секунды. Следовательно, измеренное значение псевдодальности относится к более позднему моменту времени и отличается от времени измерения псевдоскорости на величину  [с].

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

Т.к. на регистрацию одного бита измерительной информации C/A-кода требуется ,то если бы псевдодальность измерялась с точностью до 1 бита псевдошумовой последовательности, точность измерения псевдодальности составила бы  , что соответствует ошибке  в определении координат наблюдателя. При измерении псевдодальности с точностью половины бита, точность в определении координат составила бы  . В используемом GPS-приемнике заложена более сложная схема для разделения бита псевдошумовой последовательности.

3.      Вычисление вектора состояния НКА GPS по эфемеридным данным

Для нахождения вектора состояния наблюдателя необходимо уметь вычислять вектор состояния НКА. Исходные данные для построения алгоритма содержатся в интерфейсном документе GPS-системы [1].

Для описания движения НКА системы GPS используют геоцентрическую вращающуюся систему координат WGS-84 ("World Geodetic System") [6]. Время в системе представляется номером недели и смещением от начала недели в секундах. Отсчет системного времени GPS ведется непрерывно начиная с  0 часов 00 секунд 5 января 1980 года по Гринвичу. 

Расчет вектора состояния НКА производится с помощью эфемеридных данных или эфемерид,  передаваемых с борта. Эфемеридные данные обновляются раз в два часа и относятся к времени эпохи (опорному, исходному моменту времени). Согласно [1], срок действия эфемерид составляет 604800 секунд или неделю. Эфемериды состоят из стандартных элементов кеплеровской орбиты и некоторых параметров:

 - эпоха,

 - квадратный корень из полуоси орбиты,

  - отклонение значения среднего движения от ,

 - средняя аномалия на время эпохи,

- эксцентриситет орбиты,

- угол перигея,

- наклонение орбиты на время эпохи,

- скорость изменения наклонения орбиты,

- долгота восходящего узла плоскости орбиты на время эпохи,

- скорость изменения долготы восходящего узла плоскости орбиты,

- амплитуда синусоидальной поправки аргумента широты,

- амплитуда косинусоидальной поправки аргумента широты,

- амплитуда синусоидальной поправки радиуса орбиты,

- амплитуда косинусоидальной поправки радиуса орбиты,

- амплитуда синусоидальной поправки наклонения орбиты,

- амплитуда косинусоидальной поправки наклонения орбиты,

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

   -  время с начала эпохи;

 - среднее движение ,

где  - универсальная гравитационная постоянная Земли в системе WGS-84;

 - скорректированное среднее движение;

 -  средняя аномалия;

 - эксцентрическая аномалия (уравнение Кеплера);

- истинная аномалия;

 - аргумент широты;

 - поправка аргумента широты;                                     

- поправка радиуса-вектора;                                            

- поправка к наклонению;                                             

- скорректированный аргумент широты;                                               

- скорректированный радиус орбиты;                                                 

 - скорректированное наклонение.                                                    

 - положение долготы восходящего узла в инерциальной системе координат оси которой совпадают с подвижной WGS-84 на рассматриваемый момент времени,

где - угловая скорость вращения земли в системе WGS-84.

Вектор положения в плоскости орбиты НКА описывается следующими соотношениями:

.

Компоненты вектора положения НКА в системе координат WGS-84, согласно [1],  удовлетворяют следующим соотношениям:

.

В настоящей работе были получены соотношения для определения компонент вектора скорости. Для получения этих соотношений необходимо вычислить следующие производные:

,

,

,

,

,

,

,

,

.

Следовательно, вектор скорости в плоскости орбиты НКА описывается следующими соотношениями:

.

Дифференцируя соотношения для компонент вектора положения, получаем компоненты вектора состояния i-го НКА во вращающейся системе координат WGS-84:

                       (3)

4. Алгоритм определения вектора состояния наблюдателя по данным измерений псевдодальности и псевдоскорости

Обозначим через  фактическое время приема сигнала от i-го НКА, которое отличается от шкалы времени системы GPS  на величину :

                                                          .                                                                  (4)

Связь истинной дальности  и измеренной псевдодальности  до i-го НКА в момент излучения сигнала  определяется соотношением:

                                                ,                                               (5)

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

Для измеренного значения псевдоскорости  справедливо следующее соотношение:

                                                   ,                                                           (6)

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

Рассмотрим инерциальную систему координат, направление осей которой совпадает с направлением осей подвижной системы координат WGS-84 на момент времени . Будем считать, что в некоторой окрестности момента регистрации сигнала связь вектора положения в системе WGS-84 и инерциальной системе координат описывается следующем соотношением:

,                                                                           (7)

где ,.

Легко видеть, что соотношение для связи скоростей имеет следующий вид:

.                          (8)

Вывод соотношений для обратного перевода очевиден, учитывая, что .

Предположим , что в течение половины интервала  времени накопления сигнала скорость наблюдателя постоянна. Пусть ,   - положение и скорость     i-го НКА в момент времени  соответственно, а ,  - положение и скорость наблюдателя в инерциальной системе координат. Тогда подставив в (5) и (6) выражение истинной дальности и радиальной скорости с учетом времени распространения сигнала и используя (4) получаем следующие соотношения:

,                                                                (9)

,       (10)

где  - время распространения сигнала.

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

 , .                                                  (11)

Зная модель движения наблюдателя, и решив систему (11), можно найти вектор положения , вектор скорости  в инерциальной системе и  неизвестные параметры: .

Система (11) решается в два этапа. Сначала находится приближенное решение системы, а потом при избыточности измерений оно уточняется методом наименьших квадратов. Необходимо отметить, что система однозначно разрешима при количестве измерений не менее чем от пяти КА.

Предположим, что наблюдатель неподвижен. С учетом этого предположения положим и . Определим положение наблюдателя и уточним ошибку в определении скорости за счет  - разности по времени в регистрации измерений псевдодальности и псевдоскорости.

Пренебрегая учетом времени распространения сигнала и ионосферной ошибкой  в уравнениях для псевдодальности системы (11), запишем систему уравнений в виде:

, , .               (12)

Система (12) явно решается при наличии измерений от 4 НКА. Для ее решения используем алгоритм, предложенный в [4]. Необходимо отметить, что первое уравнение системы (10) решается неоднозначно и в итоге мы имеем два решения системы:

, , , , .          (13)

Рассмотрим алгоритм устранения неоднозначности решения. Для этого используем измерения от пяти НКА следующим образом.

Для каждой четверки НКА из пяти имеющихся ищется два решения (13). Таких пар решений будет . Далее применим следующий алгоритм разделения решений на два множества. На первом шаге алгоритма первая пара решений (13) распределяется по двум множествам произвольным образом.

На k-ом шаге , получаем пару решений (13). Среди этих двух решений ищем то, у которого вектор положения ,  наиболее приближен к среднему одного из множеств. Найденное решение относим к этому множеству, а второе решение к другому  множеству.

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

В качестве начального приближения решения системы (11) берутся средние значения компонент вектора положения, скорости,  и  выбранного множества:

       , , , .            (14)

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

Итак, мы получили приближенное решение системы (11) по одновременным измерениям не менее чем от 5 НКА.

Уточним найденное приближенное решение системы (11) методом наименьших квадратов [2]. Как и в случае с системой (12), пренебрежем ионосферной ошибкой измерений и составим следующую систему:

.                                                                                                                          (15)

Обозначим вектор неизвестных параметров через  и

                         ,                                  (16)

,                                    

                                                                                                                                                       (17)

где .

Тогда вектор  является решением системы (15) в смысле метода наименьших квадратов в том случае, если функция

                                                                                                (18)

достигает своего минимума. Здесь , , .

Минимум функции (18) будем искать итерационным методом. На каждой итерации будем уточнять решение с использованием линейной модели. Найдем параметры этой модели.

Пусть

             ,                              (19)

                 , .                                      (20)

где , ,                                                                                                  

        , .       (21) 

Легко видеть, что:

 ,    (22)

                                         ,                                       (23)

                 ,                           (24)

                 ,                           (25)

                .                            (26)

                                                                                                       (27)

                ,                  (28)

                ,                  (29)

                .                  (30)

В качестве нулевого приближения используем приближенное решение (14)  и                     .

Для начала итерационного процесса необходимо определить значения следующих параметров: , , и .

Для каждой итерации выполняются следующие действия.

1. Последовательно решаются системы уравнений:

                                                (31)

                                                 (32)

2. , ,                (33)    

    , ,                  (34)

3. Если и , где - наперед заданные точности, то алгоритм завершается на шаге k, иначе происходит переход к следующему шагу.

В конце итераций мы получаем решение системы уравнений (15):              

    ,   .      (35)

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

5. Статистическая оценка точности определения вектора состояния

Для апостериорной оценки точности определения вектора состояния были использованы ежесекундные измерения на месячном интервале. Статистические оценки проводились для каждого дня мерной базы. Отметим, что разделение измерений по неделям в точности соответствует разделению эфемеридной информации. Под днем недели (1-7) подразумевается интервал времени внутри недели длиной 86400 секунд.

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

Необходимо отметить, что обработка группы ежесекундных измерений проводилась в том случае, если "возраст" эфемеридных  данных всех НКА, от которых был принят сигнал, не более чем два часа. Это сделано для того, чтобы уменьшить ошибку определения вектора состояния наблюдателя, вызванную ошибкой определения вектора состояния НКА. Следствием этого является достаточно малое, по сравнению с 86400 (количество секунд в дне),  количество дневных измерений. Во втором столбце таблицы 2 показано количество обработанных измерений.

Из-за различных сбоев приемника возникали аномальные измерения. Они были отсеяны по принципу выхода за пределы априорно известной точности определения вектора состояния. Поиск аномальных измерений производился апосториорно для каждого интервала измерений. В третьем столбце таблицы 2 представлено количество аномальных  измерений. В четвертом столбце показано количество точек, которые использованы для получения статистических параметров. В пятом столбце показан процент решений выходящий за пределы  хотя бы по одной из компонент вектора состояния относительно среднего за дневной период. В шестом столбце показан процент точек, лежащих в пределах .

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

Для визуального представления геометрии разброса построены графики проекций на плоскости определения компонент вектора состояния для 3-го дня 1101 недели измерений (рис. 4). Каждая полуось эллипса соответствует  разброса соответствующих компонент.

На рис. 5 изображена динамика фазового сдвига псевдошумовой последовательности. На рис. 6 и 7 изображены график ухода частоты генератора приемника от номинала GPS-системы и график точности времени регистрации сигнала соответственно.