Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
Кроме самой математической модели, в данной работе описываются такие аспекты численной реализации кодово-фазовой дифференциальной коррекции как численная устойчивость, вычислительная эффективность и эффективность сохранения. Представленный в работе метод является вычислительно эффективным, потому что дает возможность максимально полного использования структуры математической модели и численную устойчивость алгоритмов, использующих ортогональные преобразования. Эффективность хранения также принята во внимание в нижеизложенном методе. Цель этой работы состоит в том, чтобы раскрыть основные моменты и нюансы математической модели, численного метода и учета погрешностей системы дифференциальной коррекции, которые будут необходимы для численной реализации этого подхода на ЭВМ.
В данной работе мы предполагаем, что псевдодальности по фазе несущей и коду измеряются как базовой, так и мобильной станцией. Мы рассматриваем одночастотные приемники, работающие только с частотой входного сигнала L1, потому как многие приемники могут получить только сигнал L1. Но не трудно расширить наш подход и для случая двухчастотных приемников.
Эта работа организована следующим образом. В разделе 2 мы даем математическую модель, которая используется для оценки определяемого положения. В разделе 3 дается описание нашего численного метода, состоящего из ортогональных преобразований методом Хаусхолдера, QR-факторизации и LAMBDA-метода. В разделе 4 рассматривается специфика метода, которую нужно знать при численной реализации.
Повсюду в этой работе мы используем строчные буквы для векторов и прописной регистр для обозначения матриц. Единичная квадратная матрица размерности m будет обозначена и ее i-я колонка - , так что . Вторая норма для векторов будет равна . Символом {·} будем обозначать ожидаемую оценку, cov{·} обозначит матрицу ковариации, которая определяется следующим образом: cov{x}={(x-{x})(x-{x})T}.
Обозначение показывает, что вектор v есть случайный вектор с нормальным распределением, с математическим ожиданием и матрицей ковариации V.
Рисунок 1 – Основная блок-схема алгоритма определения координат
На рис. 1 показана основная блок-схема, наглядно иллюстрирующая алгоритм работы приемника по определению координат потребителя, начиная с включения приемника и заканчивая нахождением оценки координат потребителя. Первичная обработка сигнала производится для каждого из видимых спутников, а во время вторичной обработки до каждого спутника высчитывается псевдодальность по коду и фазе; зная эти расстояния, вычисляют местоположение приемника.
Рисунок 2 - Блок-схема первичной обработки сигнала
Спутники GPS передают сигналы на двух частотах: 1575, 42 МГц и 1227, 60 МГц. На вход приемника поступают сигналы от спутников, находящихся в зоне радиовидимости. Так как для решения навигационной задачи необходимо измерить псевдодальности относительно, как минимум, четырех спутников, то приемник должен быть многоканальным (от 4 до 12 каналов).
После включения питания и начального тестирования навигационная аппаратура потребителей (НАП) входит в режим поиска сигнала навигационных спутников. Суть алгоритма поиска заключается в последовательном переборе всех возможных значений доплеровского смещения частоты при фиксированном значении битов генератора псевдослучайной последовательности (ПСП), в общем случае уникальной для каждого спутника. Для ускорения этого процесса в навигационных системах GPS используется информация долгосрочного прогноза положения спутников, которая называется альманахом.
Современные приемники являются аналого-цифровыми системами, сочетающими аналоговую и цифровую обработку сигналов. Переход на цифровую обработку сигналов осуществляется на определенной промежуточной частоте. С техническо-инженерной точки зрения, приёмник можно разделить на три функциональные части: радиочастотную часть, цифровой многоканальный коррелятор и процессор.
С выхода антенно-фидерного устройства (антенны) сигнал поступает на радиочастотную часть. Основная задача этой части заключается в усилении входного сигнала, фильтрации, преобразовании частоты и аналого-цифровом преобразовании. Помимо этого, с радиочастотной части приёмника поступает тактовая частота для цифровой части приёмника. С выхода радиочастотной части цифровые отсчёты входного сигнала поступают на вход цифрового коррелятора.
В корреляторе в цифровой форме формируются отсчеты синфазных I(k) и квадратурных Q(k) составляющих сигнала, которые являются основой работы алгоритмов поиска сигналов по задержке и частоте, вычисления псевдодальности, слежения за фазой сигнала и выделения навигационного сообщения.
Отсчёты корреляционных интегралов поступают в процессор для дальнейшей обработки и замыкания петель ФАП (фазовая автоподстройка) и ССЗ (схема слежения за задержкой). Измерения параметров сигнала в приёмнике производятся не непосредственно по входному сигналу, а по его точной копии, формируемой системами ФАП и ССЗ. Корреляционные интегралы I и Q позволяют оценить степень «похожести» (коррелированности) опорного и входного сигналов. Задача коррелятора, помимо формирования интегралов I и Q, — формировать опорный сигнал, согласно с управляющими воздействиями (кодами управления), поступающими с процессора. Кроме того, в некоторых приёмниках коррелятор формирует необходимые измерения опорных сигналов и передаёт их в процессор для дальнейшей обработки. В то же время, так как опорные сигналы в корреляторе формируются по управляющим кодам, поступающим с процессора, то необходимые измерения опорных сигналов можно производить непосредственно в процессоре, обрабатывая соответствующим образом управляющие коды, что и делается во многих современных приёмниках.
Дальность при радиотехнических измерениях характеризуется временем распространения сигнала от объекта измерения до измерительного пункта. В СРНС излучение сигналов синхронизировано со шкалой времени системы, точнее, со шкалой времени спутника, излучающего данный сигнал. В то же время, потребитель имеет информацию о расхождении шкалы времени спутника и системы. Цифровая информация, передаваемая со спутника, позволяет установить момент излучения некоторого фрагмента сигнала (метки времени) спутником в системном времени. Момент приёма этого фрагмента определяется по шкале времени приёмника. Шкала времени приёмника (потребителя) формируется с помощью кварцевых стандартов частоты, поэтому наблюдается постоянный «уход» шкалы времени приёмника относительно шкалы времени системы. Разность между моментом приёма фрагмента сигнала, отсчитанным по шкале времени приёмника, и моментом излучения его спутником, отсчитанным по шкале спутника, умноженная на скорость света, является псевдодальностью.
Математическая модель
Здесь мы даем математическую модель оценки позиции.
Постановка задачи. Заданы два одночастотных GPS приемника, которые именуются базовым (БС) s и мобильным (МБ) r, находящиеся на поверхности Земли, расположенные на относительно небольшом расстоянии друг от друга (не более 10 км.). Мы хотим найти вектор базовой линии x, т.е. вектор, направленный от приемника s к приемнику r. По базовой линии однозначно определяются координаты мобильного приемника. Оба приемника постоянно принимают сигналы с некоторого множества спутников GPS. Это множество «видимых» с обоих приемников спутников с течением времени меняется: спутники GPS все время покидают и заходят в «зону видимости».
Координаты базового приемника статичны и считаются известными в международной системе ECEF с достаточно высокой точностью (1-2 мм.). В режиме постобработки координаты мобильного приемника считаются также неизменными в определенный промежуток времени, в которое проводятся измерения.
Каждый из приемников способен высчитывать и выдавать свои измерения расстояния до заданного спутника в определенный момент времени (именуемый эпохой). Причем эти расстояния измеряются аппаратной частью приемника двумя способами: кодовым и фазовым. В данной работе описывается алгоритм решения навигационной задачи посредством кодовых и фазовых измерений. В грубом приближении, расстояние от приемника до спутника по коду рассчитывается как время запаздывания сигнала относительно системного отсчета времени, помноженная на скорость света (берется как скорость сигнала). А дальность по фазе определяется как количество длин волн несущей сигнала, которые помещаются на отрезке между спутником и приемником, помноженное на длину волны сигнала.
Цель задачи: определить координаты мобильного приемника в режиме постобработки, то есть для случая когда в качестве входных данных численного метода используются измерения, получаемые как до так и после эпохи, в которое нужно определить искомые координаты мобильного приемника.
Рисунок 1 - Геометрическое расположение БС, МС и НКА
Пусть в некоторую эпоху в поле зрения базовой станции находятся множество {d} спутников, а мобильной станции некоторое множество {f} спутников. Пересечение множеств {d} и {f} как раз таки является множеством спутников, с которыми будет происходить работа по дифференциальной коррекции. Обозначим это множество {m} (в количестве m), которое обычно называется рабочим созвездием из навигационных космических аппаратов (НКА).
Нам необходимы следующие переменные для описания геометрии относительного положения спутников и GPS приемников (рисунок 1):
- вектор, направленный от БС s к i-му НКА;
- единичный вектор, направленный из середины базовой линии x к i-му НКА;
- дальность в длинах волн от БС до i-го НКА;
- длина волны (для частоты L1, используемой здесь, 19 см)
Верхний индекс определяет спутник, нижний индекс определяет приемник. После того, как введем обозначение:
получим что:
дающее:
|
(1) |
|
|
(2) |
Заметим, что примерно равно 1, то есть вектора, направленные от приемников БС и МС к одному и тому же спутнику почти что коллинеарны, и мы можем заменить этот параметр на 1 для приложений, в которых нет необходимости производить высокоточное определение местоположения. Но для некоторых высокоточных GPS приложений нам может понадобиться неединичный параметр . Данная математическая модель позволяет делать такое упрощение, потому что в дальнейшем, в разделе 4, показано каким образом производится перерасчет x c учетом этого упрощения.
Предположим, что сигнал от i-го спутника приходит на приемник s в момент времени и он доходит от спутника до получения его приемником за период времени . В момент времени измеряется фаза несущей (в длинах волн) и измеряется код (также в длинах волн) от приемника s до i-го спутника (см. [14, раздел 5.2.2] или [15, раздел 3.4]).
|
|
(3) |
|
|
(4) |
где размерности каждого из выше приведенных уравнений определяются как “количество длин волн”,
Ионосферная ошибка , тропосферная ошибка и ошибка часов НКА может быть промоделирована [11, раздел 4.4 и 5.3]. В разделе 4 мы рассмотрим численно реализованные нами тропосферную и ионосферную модели для учете данных ошибок.
Вычтем уравнение измерения фазы несущей, соответствующее приемнику r, из (3) и замечая, что величина:
будет незначительна для коротких базовых линий, мы получим уравнения первых фазовых разностей:
|
(5) |
Для уменьшения числа зависимых неизвестных, мы можем просто объединить ошибки часов приемника, задержки в аппаратуре приемника, и начальные фазы сгенерированных несущих сигналов вместе. Заметив, что приемник s и r встречаются в каждом уравнении, мы можем пропустить эти индексы и обозначить эпоху времени k нижним индексом k. Таким образом, определяя:
из (5) и (1) получим что:
|
(6) |
Аналогично, введя обозначения:
мы получим из (4) уравнения первых кодовых разностей:
|
(7) |
Следуя литературе [15, раздел 3.4]), мы положим, что все в уравнениях измерений фазы несущей (и все в уравнениях кодовых измерений) для различных НКА и различных эпох есть некоррелированные шумы с нормальным распределением, следовательно эти и независимы. Введем:
|
(8) |
мы получим (помня, что ) из (6) и (7):
|
(9) |
|
|
(10) |
где мы полагаем, что среднеквадратические отклонения и являются известными. Уравнения (9) и (10) являются желаемыми (искомыми) уравнениями первых разностей в эпоху k.
3. Численный метод
В общем случае, для решения системы линейных уравнений (9), (10) применяется метод наименьших квадратов, но для повышения вычислительной эффективности и устойчивости алгоритма, используют разноплановые модификации данного метода.
Для решения общей линейной проблемы наименьших квадратов, имеется два основных метода: метод нормальных уравнений и метод QR факторизации (в которой заложены ортогональные преобразования). Плохая обусловленность матрицы коэффициентов (что может быть следствием плохой геометрии расположения НКА GPS) может привести к увеличению вычислительных погрешностей метода. Поэтому, для того чтобы не увеличивать коррелированность входной системы уравнений и повысить численную устойчивость алгоритма, вместо простых вторых разностей в методе используются ортогональные вторые разности [5, раздел 5.3].
Следуя [3], мы распишем рекурсивный алгоритм ортогональных преобразований для оценки позиции, основанный на уравнениях (9) и (10) измерений первых разностей фазы несущей и кода. Мы собираемся полностью использовать структуру задачи (для примера, учтем факт того, что вектор целочисленных неопределенностей идентичен на каждом шаге при отсутствии срывов циклов или при уходе или заходе НКА GPS в рабочее созвездие) для того, чтобы сделать алгоритм как можно более эффективным. Отметим, что Тибериус [15] также использует рекурсивный алгоритм. Но он использует стандартный подход SRIF, основанный на уравнениях измерений вторых разностей фазы несущей и/или кода.
На настоящий момент времени мы предполагаем, что в (9) и (10) известны, и что число видимых НКА не изменяется во временной промежуток наблюдения. Сначала мы устраним из уравнения (9) измерений фазы несущей.
Пусть есть преобразование Хаусхолдера (смотрите также [5, стр. 209]) такое, что:
|
((11) |
Отметим, что матрица симметричная и ортогональная. Представим как , тогда:
|
где имеет размерность {m, m-1}. |
((12) |
Умножая (9) на слева, получим:
|
|
((13) |
Так как верхнее уравнение в (13) включает неизвестное (и не желательное) выражение его можно исключить из алгоритма оценки положения мобильной станций. Таким образом, мы имеем:
|
|
((14) |
Заметим, что мы использовали простое преобразование Хаусхолдера, вместо часто используемой техники двойных разностей для устранения неизвестной переменной . Его главным преимуществом является то, что преобразованные измерения являются все еще некоррелированными. Техника ортонормализации Грама-Шмидта, которая может также устранить и сделать инструмент трансформации некоррелированным, не является таким прямым, как техника преобразования Хаусхолдера.
Так как является матрицей , она не имеет полного ранга по столбцу и мы не в состоянии получить однозначную оценку . Если мы положим как новый вектор, то мы потеряли бы целое натуральное . Поэтому, идем другим путем и введем вектор двойных разностей целочисленных неопределенностей (double difference integer ambiguity - DDIA)
|
((15) |
где без потери общности мы выбрали первый попавшийся НКА как «референсный» и присвоили ему номер 1. Заметим, что все еще вектор целых чисел. Определим:
|
((16) |
где несингулярная квадратная матрица размерности {m-1} и . Это легко проверить, используя (12)
|
((17) |
Таким образом, выражение (14) может быть переписано как:
|
|
((18) |
Подобным образом мы можем также устранить из уравнения (10) измерений кода и при этом получим:
|
|
((19) |
Для того чтобы при объединении и получившийся параметр имел такую же ковариационную матрицу как , мы определим и умножим (19) на , а затем объединим результирующее уравнение и (18):
|
|
((20) |
Последний член в уравнений (20), , по своей сути является неоднородностью в системе линейных уравнений, неточностью, вызванная шумами и переотражениями. Эта неоднородность подчиняется нормальному распределению, то есть, ее пределы изменения значений могут быть оценены, опираясь на результаты ограниченного количества экспериментов с дифференциальной системой. Как показывает практика, данный вид погрешностей на порядок ниже погрешностей, возникающих из-за атмосферной задержки. Поэтому, для упрощения дальнейшего математического анализа, опустим эту неоднородность из уравнения.
Теперь приведем уравнение (20) к виду, пригодному для численной реализации. Правую часть (20) можно переписать как
,
где O – квадратичная матрица, состоящая из нулей, размерности {m-1}.
Таким образом, система уравнений (20) может быть записана в следующей матричной форме:
|
Aс=B, |
((21) |
где
где k проходит целые числа от 1 до n включительно.
Матрица имеет размерность {m, m-1}.
Матрица F имеет размерность {m-1, m-1}.
Обозначения:
n – количество эпох, в которые производятся наблюдения;
m – количество спутников в рабочем созвездий;
- матрица, составленная из единичных векторов, поделенных на длину волны L1 и коллинеарных направлению на спутники в эпоху с номером k (эпохой в литературе по СНС называется момент времени, в который производится измерение);
где - единичный вектор, направленный из середины базовой линии к i-му НКА в эпоху ;
– координаты вектора, соединяющего базовый и мобильный приемники, то есть вектора базовой линии в эпоху ;
Переопределенная система (21) решается методом наименьших квадратов (ищется с, минимизирующая эвклидову норму ) приведением матрицы к верхнетреугольной форме методом Хаусхолдера. Метод Хаусхолдера заключается в разложении матрицы A в виде мультипликации двух матриц (другими словами, происходит QR-факторизация):
A=T*H,
где T – квадратная ортогональная матрица, и H – верхнетреугольная прямоугольная матрица.
После чего условие минимизации можно переписать как:
.
Так как T – ортогональная матрица, данное условие минимума равнозначно минимизации другой нормы:
.
Верхнетреугольную матрицу H представим в виде:
,
где и – верхнетреугольные матрицы, причем имеет размерность 3х3. При практической реализации верхнетреугольная форма вычисляется рекурсивно, по мере поступления данных в последовательные эпохи.
А матрицу, получающуюся из мультипликации , обозначим как , так что
Можно заметить из уравнения, что матрица U не принимает участия в поиске с, поэтому нижняя часть получившейся системы линейных уравнений отбрасывается. Именно такая возможность, предоставляемая методом Хаусхолдера, и приводит к экономии времени и процессорной мощности обрабатывающего ЭВМ. Следовательно:
|
((22) |
Матричное уравнение (22) записывается как система двух уравнений:
|
((23) |
и
|
((24) |
Из (23) получаем вещественнозначное смешанное кодово-фазовое решение:
,
где – вектор вещественных величин (вещественнозначных неопределенностей), полученный решением системы уравнений (24). Такое решение (без использования целочисленной природы неопределенностей) позволяет достичь точности в несколько дециметров.
Для того чтобы получить высокую точность оценок положения для как можно минимального числа эпох, должна быть использована природа целочисленных неопределенностей. Заметим, что в нашем алгоритме мы только оценили как вещественный вектор, но не определили его как целочисленный вектор. Для некоторых приложений этого будет достаточно, и для них не имеет смысла дополнительно тратить время вычислений для разрешения целочисленных неопределенностей. Также, если целочисленные неопределенности определены неправильно, это может вызвать большие ошибки в оценке положения.
Имеется несколько подходов для разрешения целочисленных неопределенностей на базе вторых разностей. Один их хорошо известных подходов называется LAMBDA метод, смотри [14, глава 8]. На вход метода подаются вещественнозначные неопределенности (в качестве исходной оценки), матрица ковариации неопределенностей и требуемое количество векторов-кандидатов (на практике достаточно двух). Он основан на ковариации . Ковариация неопределенностей может быть найдена из (24):
|
((25) |
Дадим теоретическое описание использования LAMBDA-метода при решении нашей задачи. Линейные системы уравнений, используемые в математической модели для решения задачи дифференциальной коррекции, можно представить в виде:
|
((26) |
где y – заданный вектор данных, a и b – неизвестные параметрические вектора и e – вектор шумов. Вектор данных y состоит из «наблюдаемых» минус «вычисленных» первых разностей по коду и фазе за весь временной промежуток наблюдений. a – вектор разностей целочисленных неопределенностей размерности {m-1}, измеряемые в циклах. Значения вектора a являются целыми. В векторе b размерности p заключены оставшиеся неизвестные параметры уравнений, такие как координаты базовой линий и тропосферные (ионосферные) параметры атмосферной задержки.
При использовании принципа наименьших квадратов, система (26) решается посредством минимизации:
|
, |
((27) |
где норма вычисляется по формуле
- матрица ковариации наблюдаемых величин,
- целочисленное пространство размерности m-1.
Для того чтобы решить такого рода задачу минимизации, необходимо сначала провести ортогональную декомпозицию целевой функции. Распишем целевую функцию в виде суммы трех квадратов,
|
((28) |
где
,
где с ортогональными проективными операторами и Матрица является ортогональным проективным оператором, который проецируется ортогонально вдоль направления B (относительно метрики ). Аналогично, матрица является ортогональным проективным оператором, который проецируется ортогонально вдоль направления A.
Ковариационные матрицы и задаются формулами:
Вектора и являются вектором веществозначных неопределенностей и вещественным решением. Они определяются из (24) и (23), то есть когда в задаче минимизации (27) отсутствует ограничение по целочисленности . является остаточным вектором соответствующий веществозначному решению.
Вектор является решением b, в предположении, что a известен. Следовательно, он является условным решением b. Такое условное решение и ее ковариационная матрица могут быть записаны в виде:
где ковариационная матрица равна
Из (27) и (28):
Заметим что первое выражение, не поддается минимизации потому что оно не зависит от a и b. Также последнее выражение в уравнении является неотрицательным (так как под квадратом) и для любого a из пространства целых чисел можно подобрать такое b что оно будет равняться нулю. Следовательно, решение данной задачи минимизации находится посредством формул:
|
((29) |
Вектора и представляют собой вектор разрешенных целочисленных неопределенностей и вектор базовой линии при разрешенных целочисленных неопределенностях.
Поиск не может быть проведен по всему пространству целых чисел, поэтому необходимо определить область поиска, центрированного в точке в смысле метрики
. Такую область можно определить как:
где является оптимально выбранной константой-ограничителем. Для того чтобы область поиска не содержала слишком много целочисленных векторов-кандидатов, должен быть небольшим.
Граница области поиска является эллиптической. В случае с фазовыми неопределенностями эллипсоид сильно вытянут из-за высокой коррелированности компонентов вектора неопределенностей. Такая вытянутость вредит вычислительной эффективности поиска, поэтому область поиска первым шагом приводят к более сферической форме:
|
((30) |
используя приемлемое преобразование неопределенностей: Преобразование неопределенностей считается приемлемым, когда результатом преобразований T и являются целые числа. Такие матрицы сохраняют целочисленную природу неопределенностей. Для того чтобы область поиска имела более сферическую форму, строится такое T-преобразование, сохраняющее объем области, что неопределенности подвергаются как можно более сильной декорреляции. Используя треугольное преобразование левая часть квадратичного неравенства расписывается как сумма квадратов:
С левой стороны в неравенстве участвуют условные оценки неопределенностей , которые находятся как условные случайные целочисленные вектора . Из данного неравенства можно вывести оптимальные интервалы поиска по каждому измерению пространства . Интервалы задаются m-1 формулой:
|
((31) |
Таким образом, определив интервалы возможных изменений вектора неопределенностей через (31), можно заметить, что область поиска будет содержать ограниченное количество векторов-кандидатов. Причем, число кандидатов не будет превышать . То есть, если нам позволяют вычислительные мощности, путем простейшего перебора в режиме постобработки, можно найти оптимального кандидата, удовлетворяющего условию минимизации (29). Поиск начинается с целочисленного вектора , который является простым округлением вещественнозначной оценки. На каждом шагу поиска осуществляется коррекция размера области поиска (изначально размер выбирается в зависимости от количества требуемых векторов-кандидатов). После того, как решение найдено в декоррелированных координатах, производится обратное отображение и целочисленные неопределенности вычисляются в исходных координатах.
На выходе метода получаются разрешенные целочисленные неопределенности (целочисленные вектора-кандидаты, с наименьшими отклонениями от исходного вещественнозначного) и отклонения векторов-кандидатов в смысле метрики, соответствующей матрице ковариации.
4. Специфика
При реализации данного алгоритма на практике следует учесть некоторые нюансы.
Во-первых, в математической модели было сделано допущение, что единичные вектора, направленные от приемников к одному и тому же спутнику практически коллинеарны друг к другу. Чтобы минимизировать ошибку от сделанного допущения, введем вектор коррекции:
- координаты вектора коррекции направления псевдодальности для спутника j (j пробегает значения от 1 до n, где n – количество спутников, видимых в данную эпоху i). Он вычисляется следующим образом:
где Hs – вектор, направленный от базового приемника к спутнику;
Blest – вектор базовой линии.
Во-вторых, приходится учитывать разность тропосферных и ионосферных ошибок для базового и мобильного приемников. Подсчет разности тропосферных ошибок производится через выбор определенной модели тропосферы (UNB3m), в которой дана формула подсчета задержки проходящего сигнала при определенных значениях влажности, температуры, давления, дня года и высоты над уровнем моря. Для расчета ионосферных ошибок также существует модель задержки сигнала во время прохождения ионосферы.
Использованная литература
[1] J. Ashjaee, Ashtech XII GPS technology, Proceedings of IEEE PLANTS’90 Position Location and Navigation Symposium, Las Vegas, IEEE, New York, pp. 184–190.
[2] J. D. Bossler, Clyde C. Goad, Peter L. Bender, Using the Global Positioning System (GPS) for geodetic positioning, Bulletin G.eod.esique, Vol. 54, 1980, pp. 553–563.
[3] X.-W. Chang and C.C. Paige, An orthogonal transformation algorithm for GPS positioning, SIAM J. Sci. Comp., Vol. 24., pp. 1710-1732, 2003.
[4] C.C. Goad, Optimal filtering of pseudoranges and phases from single-frequency GPS receivers, Navigation, Journal of The Institute of Navigation, Vol. 37, Spring 1990, pp. 249–262.
[5] G.H. Golub and C.F. Van Loan, Matrix Computations, 3rd ed., The Johns Hopkins University press, Baltimore, MD, 1996.
[6] R. Hatch, The synergism of GPS code and carrier measurements, Proceedings of the Third International Symposium on Satellite Doppler Positioning, New Mexico State University, New Mexico, February 8–12, 1982, Vol. 2, pp. 1213–1231.
[7] B. Hofmann-Wellenhof, H. Lichtenegger, and J. Collins, GPS Theory and Practice, 5th ed.,
Springer, New York, 2001.
[8] P.Y.C. Hwang and R.G. Brown, GPS navigation; combining pseudo-range with continuous carrierphase using a Kalman filter, Navigation, Journal of The Institute of Navigation, Vol. 37, pp. 181–196, 1990.
[9] X.X. Jin, Algorithm for carrier-adjusted DGPS positioning and some numerical results, Journal of Geodesy, Vol. 71, pp. 411–422, 1997.
[10] A. Kleusberg, Kinematic relative positions using GPS code and carrier beat phase observations, Manuscripta Geodaetica, Vol. 10, 1986, pp. 257–274.
[11] P. Misra and P. Enge, Global Positioning System, Signals, Measurements, and Performance,
Ganga-Jamuna Press, Lincoln, Massachusetts, 2001.
[12] B.W. Remondi, Using the Global Positioning System Phase Observable for Relative Geodesy: Modeling, Processing, and Results, Ph.D. thesis. The University of Texas at Austin, 1984.
[13] P.J. Teunissen, The GPS phase-adjusted pseudo-range, Proceedings of the Second International Workshop on High Precision Navigation, K. Linkwitz and U. Hangleiter (eds), D.ummler, Bonn, 1991, pp. 115-125.
[14] P. J. Teunissen and A. Kleusberg (ed). GPS for Geodesy, 2nd ed. Springer, Berlin, 1998.
[15] C.C.J.M. Tiberius, Recursive Data Processing for Kinematic GPS Surveying, PhD thesis, Department of Mathematical Geodesy and Positioning, Delft University of Technology, 1998.
Для программной части
Chistart.m - вычисляет первоначальный размер искомого эллипсоида. Эта библиотека вычисляет или аппроксимирует первоначальный размер искомого эллипсоида. Если просимое количество кандидатов не больше числа (размерность + 1), это делается посредством вычисления квадратичного расстояния от частично округленных плавающих векторов до плавающего вектора в метрике матрицы ковариации. В противном случае используется аппроксимация.
CorrectedEx.m – функция вычисляет коррекцию псевдодальности, возникающей из за допущения рабочей математической модели. В математической модели было сделано упрощение о том что единичные вектора направленные от базового и мобильного приемников к спутникам коллинеарны.
cubfit.m – скрипт производит кубическую интерполяцию входной функций полиномом.
DD_CPR_303ST.m – этот скрипт является основной подпрограммой модуля, производящего решение навигационной задачи кодово-фазовым методом.
Файлы DD_CPR_DATA1.m и DD_CPR_DATA7Delta.m являются входными данными для основных подпрограмм DD_CPR_303ST.m и PRSD300.m . Оба файла содержат массивы из векторов, компонентами которых являются SVN (номер спутника), Xsx, Xsy, Xsz ( координаты спутника в системе EСEF), PR1, PR2 (дальности от данного спутника до базового и мобильного приемников, измеренные кодовым методом), CP1, CP2 (дальности от данного спутника до базового и мобильного приемников, измеренные фазовым методом), Epoch (момент времени, в которое было произведено измерение).
DDIdr.m - коррекция для двойных разностей вызванная ионосферной ошибкой.
DDTdr_M.m - коррекция для двойных разностей вызванная тропосферной ошибкой.
decorrel.m - декоррелирует матрицу ковариации неоднозначностей. Эта библиотека создает декорелированную матрицу Q, посредством нахождения матрицы Z и производя над ней соответствующее преобразование.
ecef2enu.m - функция делает преобразование координат из системы EСEF в локальную систему координат (с началом координат на базовом приемнике, и ось z направлена вверх).
ecef2llh - библиотека для вычисления геодезических координат. Находит широту, долготу, высоту обьекта при заданных ECEF координатах этого объекта.
ElSin.m - расчет синуса угла места спутника относительно местоположения на поверхности земли базового приемника.
FCS.m - функция производит «сглаживание», то есть удаляет проскальзывания цикла (при расчете дальности фазовым методом). Сглаживание может производиться только для измерений, имеющих незначительный срыв цикла фаз несущей (на целое число волны или полволны). Из практики известно, что первые фазовые разности в местах небольших срывов фаз изменяется незначительно, так что для обнаружения срывов приходиться работать с производными вторых фазовых разностей. Для починки срыва в целую волну, используется стандартная матлабовская функция округления “round”, а для срывов в полволны несущей сигнала - внешняя функция “halfround” («полуокругления»).
halfround.m – подключается к работе функций FCS.m. Производит некое подобие «полуокругления».
ionocorr2.m – применяется ионосферная модель для вычисления ионосферной поправки.
lambda2.m – Оценка целочисленных неоднозначностей используя LAMBDA метод (основная версия). Эта библиотека производит оценку целочисленных неоднозначностей используя LAMBDA метод, разработанный в Дельфтском университете Технологии, Математическая геодезия и позиционирование.
ldldecom.m – Нахождение LtDL-декомпозиции Q-матрицы. Эта библиотека находит LtDL декомпозицию входной матрицы вариантности/ковариантности.
locvert.m - вычисление локальной вертикали.
lsearch.m – процедура поиска в разрешении целочисленных неопределенностей. Эта библиотека находит целочисленный вектор который ближе всего к данному плавающему вектору в смысле наименьших квадратов. Лучше всего производить этот поиск только с неоднозначностями, которые были декорелированы с использованием LAMBDA метода.
PRSD300.m - этот скрипт является основной подпрограммой модуля, производящего решение навигационной задачи кодовым методом.
togeod.m – библиотека для вычисления геодезических координат (широта, долгота и высота) при заданных координатах Картезиана X,Y,Z и значении опорного эллипсоида (большой полуоси (a) и инверсии сглаженности finv).
TOPO_ENU.M – Преобразование вектора dx в топоцентрические координаты. Система с началом координат в X. Оба параметра заданы векторами с тремя компонентами. Выходные данные: E Восточное положение в координатах, N Северное положение в координатах, U Направление на вверх.
TROPO.m – Вычисление тропосферной коррекции. Коррекция дальности ddr должно быть вычтено из псевдодальностей и несущей фазы.
UNB3M.m - Эта библиотека использует UNB3m модель для вычисления наклонной нейтральной атмосферной задержки для заданных широты, долготы, высоты, дня года (DOY) и угла места.
VECEF2ENU.m - преобразование координат.
X2corr.m - коррекция на несинхронность получения данных базовым и мобильным приемниками.
Скрипт DD_CPR_303ST, написанный на языке программирования «Mathlab», является основной подпрограммой программно-математического модуля, разрабатываемого по проекту…
Алгоритм работы скрипта может быть представлен в виде нижеследующей блок-схемы (рис. 1). Дальше будет дано описание каждого блока в отдельности.
Входные данные: файл данных DD_CPR_DATA1, DOY, Xr, Press, temp, hum.
В этом блоке описываются входные параметры и данные подпрограммы:
DD_CPR_DATA1 содержит массивы из векторов, компонентами которых являются SVN (номер спутника), Xsx, Xsy, Xsz ( координаты спутника в системе ESEF), PR1, PR2 (дальности от данного спутника до базового и мобильного приемников, измеренные кодовым методом), CP1, CP2 (дальности от данного спутника до базового и мобильного приемников, измеренные фазовым методом), Epoch (момент времени, в которое было произведено измерение).
DOY – день года (передается в навигационном сообщений)
Xr – координаты базового приемника в EСEF (предполагается что местоположение базового (неподвижного) приемника определены с большой точностью до начала измерений)
Press (значение давления), temp (значение температуры), hum (показатель влажности) – параметры для подсчета поправки в тропосферной модели (передаются в навигационном сообщений)
В этом блоке необходимо выделить такой отрезок времени, в котором не наблюдается значительных срывов измерений фаз несущей. В случае срыва фаз ищется другой отрезок времени.
Вводится новый массив данных {b} с компонентами: SVN (номер спутника), Xsx, Xsy, Xsz (координаты спутника в системе EСEF), PR: = PR1 - PR2 (первые разности для кодового метода), CP: = CP1 - CP2 (первые разности для фазового метода), Epoch (момент времени, в которое было произведено измерение).
В этом блоке удаляются измерения спутников, от которых приемник получает сигналы не на всем промежутке заданного отрезка времени, или же измерения, полученные от переотраженных сигналов.
В блоке создаются два массива fstep и endep (в первом находятся измерения, сделанные в начальные моменты времени, а во втором – измерения в конце заданного временного отрезка). После чего для этих массивов измерений находится синус угла места спутников относительно базового приемника, и удаляются те измерения спутников, для которых синус угла меньше sin = 0.17
В этом блоке создаются три матрицы dpr размерности (m, n), dcp размерности (m, n) и Xs размерности (3, n, m), где m – количество эпох, n – количество спутников. Элементами этих матриц являются следующие измеренные значения:
dpr (i, j) – первая кодовая разность в эпохе i для спутника j.
dcp (i, j) - первая фазовая разность в эпохе i для спутника j.
Xs (k, j, i) – компонента k вектора Xs в системе EСEF в эпохе i для спутника j.
Входными данными для этого блока является dcp – матрица первых фазовых разностей. Эта матрица обрабатывается внешней функцией FCS, которая выполняет следующие операций: вычисляет матрицу вторых фазовых разностей относительно значений первого спутника, находит матрицу значений проскальзывания цикла с помощью «полуокругления» (для этого подключается функция halfround), а затем округления локальных разностей вторых фазовых разностей. Матрица значений проскальзывания цикла используется для корректировки входной матрицы первых фазовых разностей.
Значения вторых кодовых разностей определяются как разность между значениями первых кодовых разностей и значением первой кодовой разности для первого в списке спутника.
Значения вторых фазовых разностей определяются как разность между значениями первых фазовых разностей и значением первой фазовой разности для первого в списке спутника.
В этом блоке находится кодово-фазовое решение навигационной задачи без учета целочисленной природы вторых фазовых разностей (то есть без учета того что вторая фазовая разность дальностей ddcp равна N*λ, где λ – длина волны сигнала со спутника, а N – целое число). Таким образом, одним из результатов первого прохода будет вектор, компонентами которого являются целые числа, полученные из округления вычисленных веществозначных N. Этот вектор является входным параметром (кандидатом) для LAMBDA метода.
В блоке используются нижеследующие внешние функций:
cubfit.m – скрипт производит кубическую интерполяцию входной функций полиномом.
CorrectedEx.m – функция вычисляет коррекцию псевдодальности возникающую из за допущения рабочей математической модели. В математической модели было сделано упрощение о том что единичные вектора направленные от базового и мобильного приемников к спутникам коллинеарны.
ecef2llh - библиотека для вычисления геодезических координат. Находит широту, долготу, высоту обьекта при заданных ECEF координатах этого объекта.
locvert.m - вычисление локальной вертикали.
X2corr.m - коррекция на несинхронность получения данных базовым и мобильным приемниками.
DDTdr_M.m - коррекция для двойных разностей вызванная тропосферной ошибкой.
DDIdr.m - коррекция для двойных разностей вызванная ионосферной ошибкой.
VECEF2ENU.m - преобразование координат.
Входными данными блока являются вещественные значения неопределенностей и матрица ковариаций неопределенностей, а выходными данными – определенные методом целочисленные неопределенности.
К Lambda2 подключаются внешние функций decorrel (декоррелирует матрицу ковариации неоднозначностей), chistart (вычисляет первоначальный размер искомого элипсоида в пространстве значений кандидатов) и lsearch (осуществляет процедуру поиска кандидатов).
С найденными в предыдущем блоке целочисленными неопределенностями осуществляется второй проход. Подключается внешняя функция CorrectedEx.
CorrectedEx.m – эта функция вычисляет коррекцию псевдодальности, возникающую из за допущения рабочей математической модели. В математической модели было сделано упрощение о том что единичные вектора направленные от базового и мобильного приемников к спутникам коллинеарны.
На выходе матлабовского скрипта DD_CPR_303ST получают координаты вектора базовой линий в системах ECEF и ENU.
Рисунок 3 – блок-схема алгоритма работы основного скрипта
Данный документ описывает алгоритм дешифровки навигационных параметров, поступающих с приемника на ЭВМ, используемых для определения координат спутника в системе ECEF. Предлагается навигационные параметры записывать в виде массива векторов, компонентами которых являются номер спутника, момент времени в системе GPS, когда сигнал был получен приемником, эфемериды и поправки для эфемерид спутника.
Первый этап. На первом этапе проводится дешифровка параметров, поступающих в сообщениях U-blox класса RXM, с ID соответствующей RXM-EPH.
Дешифровка проводится только с сообщениями, у которых нулевой байт содержит 0xB5, в первом байте - 0x62, во втором байте - 0x02, в третьем байте – 0x31
Ниже идет последовательность шагов алгоритма на первом этапе дешифровки:
Рисунок 4 - Алгоритм обработки сообщения с эфемеридами
Второй этап. На втором этапе проводится дешифровка навигационных измерений, поступающих в сообщениях U-blox класса RXM, с ID соответствующей RXM-RAW.
Дешифровка проводится только с сообщениями, у которых нулевой байт содержит 0xB5, в первом байте - 0x62, во втором байте - 0x02, в третьем байте – 0x10 Измеренные значения присоединяются к тому вектору, компоненты которого были получены на ЭВМ одновременно с вышеупомянутыми измерениями.
Ниже идет последовательность шагов алгоритма на втором этапе дешифровки:
k := NSV.
|
Если же (k=0) тогда алгоритм завершается.
Рисунок 5 – Алгоритм обработки сообщения с сырыми данными
Третий этап. На третьем этапе проводится подведение компонентов вектора к соответствующим истинным значениям навигационных параметров с помощью фактора масштабирования SF (Scale Factor). Для компонент вектора без указания параметра SF следует полагать что фактор масштабирования равен единице. Таким образом,
a[ I ] := SF*a[ I ] (1)
где a[ I ] –компонента вектора.
В результате работы данного алгоритма на выходе получается массив векторов, имеющих по 24 компоненты. Каждый вектор содержит данные по эфемеридам и поправкам эфемерид для заданного спутника в определенный момент времени в системе отсчета GPS.
Примечание: Байт, который идет в начале сообщения следует считать нулевым, а начальный бит в байте – первым, дешифруемые значения навигационных параметров записанными в память в числовом формате (U4) Unsigned Long, если не оговорено обратное. Все мультибайтовые значения в сообщений записаны в порядковом формате Little Endian ("Прямой порядок байт" в исполняемом коде; порядок от младшего к старшему: запись начинается с младшего и заканчивается старшим байтом).
Для расчёта эфемерид нам необходимы данные из навигационного сообщения для каждого спутника в определенный момент времени, содержащиеся в подкадрах 1, 2 и 3, а именно:
- средняя аномалия в начале отсчёта,
- изменение среднего движения от начала отсчёта времени, отриц. и полож.
- эксцентриситет орбиты
- квадратный корень главной полуоси
- долгота восходящего узла за недельную эпоху
- угол наклона в начале отсчёта
- расстояние перицентра от узла
- прецессия восходящего узла
IDOT – скорость изменения наклона орбиты
- амплитуда косинусной гармонической поправки к аргументу широты
- амплитуда синусной гармонической поправки к аргументу широты
- амплитуда косинусной гармонической поправки к радиусу орбиты спутника
- амплитуда синусной гармонической поправки к радиусу орбиты спутника
- амплитуда косинусной гармонической поправки к радиусу орбиты спутника
- амплитуда синусной гармонической поправки к углу наклонения
- опорное время (опорная эпоха эфемерид или начальное время синхронизации данных)
- коэффициент (постоянной составляющей) аппроксимации временных параметров
- коэффициент (первого порядка составляющей) аппроксимации временных параметров
- коэффициент (второго порядка составляющей) аппроксимации временных параметров
- опорное время
ITOW – время приемника GPS в миллисекундах
PR - псевдодальность
Координаты спутника определяются для каждого вектора, полученного на выходе алгоритма дешифровки навигационных параметров. Для расчёта координат навигационных спутников в геоцентрической фиксированной системе координат (ECEF) используется следующий алгоритм:
(0)
где скорость света с = 2.99792458 x 108 м/с2
, (1)
где:
t - системное время GPS на момент передачи сообщения (время в UTC, на которое производим вычисление координат навигационного спутника. Для получения координат навигационного спутника ежесекундно t необходимо на каждом новом шаге прибавлять одну секунду) Системное время GPS устанавливается управляющим сегментом и сверяется с UTC и поддерживается Морской обсерваторией США. Нулевой отчет времени GPS определен в полночь с 5 на 6 января 1980 г. Наибольшая единица, используемая в формировании времени GPS - одна неделя, определенная как 604 800 секунд. Время GPS может отличаться от UTC, что связано с дрейфом непрерывной шкалы времени GPS и периодически корректируемым на целое число скачков секунд временем UTC Управляющий сегмент контролирует системную временную шкалу GPS и она не должна отличаться от UTC более, чем на 1 микросекунду. В навигационные сообщения включаются данные по связи времени GPS и UTC Точность этих данных на интервале передачи находится в пределах 90 наносекунд.
- опорное время (опорная эпоха эфемерид или начальное время синхронизации данных), извлекаемое из навигационного сообщения.
- поправка для перехода от бортовой шкалы времени к системной шкале. Эта поправка вычисляется по формуле:
(2)
где:
- коэффициенты полинома по расчёту поправки для перехода от бортовой шкалы измерения времени к системной шкале;
- поправка на релятивистский эффект, которая вычисляется в пункте 4 данного алгоритма (в первый проход его значение приравнивается нулю).
Если tk>302400 c, то истинное tk=tk-604800 c, если tk< -302400 c, то истинное tk=tk+604800 c.
(3)
где:
- средняя аномалия в начале отсчёта, извлекаемая из навигационного сообщения, n – среднее движение, рассчитываемое по формуле:
, (4)
где:
- расчётное среднее движение, рассчитываемое по формуле
, (5)
в которой A–большая полуось орбиты навигационного спутника , извлекаемая из навигационного сообщения, , - изменение среднего движения от начала отсчёта времени, передаётся в кадре навигационного сообщения.
(6)
где:
- эксцентриситет, извлекаемый из навигационного сообщения; в качестве начального приближения берут следующее выражение:
(7)
Процесс продолжаем до тех пор, пока | Еk+1 – Еk| < 10−7.
Этот критерий не проверять. Критерием истинности значения полученной эксцентрической аномалии является вычисленное значение средней аномалии М эпохи tk по полученному значению:
(8)
Полученное значение средней аномалии М эпохи tk должно быть равно значению из пункта 2 данного алгоритма.
(9)
где:
Е – эксцентрическая аномалия, вычисленная на шаге 3, A–большая полуось орбиты навигационного спутника , извлекаемая из навигационного сообщения, , с – скорость света c = 2.99792458 x 108 м/сек2
После вычисления повторяем все вычисления, начиная с пункта 1 по 4.
(10)
(11)
где:
- эксцентрическая аномалия, которая ищется на шаге 3 при решении уравнения Кепплера, - эксцентриситет, извлекаемый из навигационного сообщения.
Далее с учётом знака sin и cos определяем номер четверти угла . И с учётом номера четверти для нахождения υ вычисляем либо arcsin (), либо arccos ().
(12)
где:
- истинная аномалия, вычисленная на шаге 5, - расстояние перицентра от узла, извлекаемое из навигационного сообщения.
(13)
где:
, извлекаются из навигационного сообщения.
(14)
где:
, извлекаются из навигационного сообщения.
(15)
где:
,,, IDOT извлекаются из навигационного сообщения.
(16)
(17)
где:
извлекаются из навигационного сообщения, - угловая скорость вращения Земли (стандарт WGS - 84), .
(18)
(19)
где время распространения сигнала от навигационного спутника до потребителя вычисляется по формуле:
(20)
здесь:
D – измеренная псевдодальность навигационный спутник – потребитель, с – скорость света c = 2.99792458 x 108 м/сек2
Базовая станция s
Мобильная станция станция r
Вектор x
x/2
x/2
ei
i-й НКА