Код УДК – 004.622 Подготовка данных Напечатано: ВНИИГеосистем, «Геоинформатика», №3, 2008 г. Авторы: Гришин М.Л., д.т.н. Данилкин Ф.А.

Метод быстрой фильтрации потока данных о глобальной позиции наблюдаемого объекта на примере GPS-телеметрии

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

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

Рис.2 Общая схема фильтра GPS-телеметрии      В схеме фильтра, GPS-выборкой обозначен необходимый набор данных, формируемый GPS-приемником (далее просто выборка, в географическом смысле – позиция или точка, рис. 3):

Т – время по Гринвичу (от 1970 года)

j – широта в градусах (от -90 до 90, положительное значение примем в соответствие северному полушарию)

l – долгота в градусах (от -180 до 180, положительное значение примем в соответствие восточному полушарию)

h – эллипсоидная высота в метрах (может иметь отрицательное значение)

PDOP – позиционный  фактор потери точности (в пространстве)

VDOP – вертикальный фактор потери точности (только по высоте)

HDOP – горизонтальный фактор потери точности (на плоскости касательной к рассматриваемой точке поверхности эллипсоида)

            M – режим определения местоположения [1] (учитывая разнообразность протоколов и различный подход к формированию данного параметра, обозначим это как набор состояний: 0 – Сигнал GPS отсутствует, 1 – Доступно только время, 2 – Двухмерный режим позиционирования, 3 – Трехмерный режим позиционирования)

Рис.3. Иллюстрация системы координат и параметров эллипсоида Предварительная фильтрация реализует отсев выборок, которые являются ошибочными, неполными или не соответствуют режиму, для этого анализируется M. Для ряда приложений на данном этапе фильтрации может оказаться полезным, предусмотреть генерацию таких событий, как «потеря сигнала GPS» и «начало работы фильтра». Дополнительно, в зависимости от требований к качеству данных, можно отсеивать выборки, в которых значения DOP выходят за максимально допустимые пределы, указанные в настройках фильтра. Для HDOP рекомендуется брать максимально допустимое значение (HDOPmax) равным 3..6, для VDOP – VDOPmax=4..7. Максимально допустимый PDOP(PDOPmax) можно взять, исходя из того, что все три рассматриваемых DOP связаны отношением: PDOP= [2].

Рис.4 Алгоритм предварительной  фильтрации      Алгоритм предварительной стадии фильтрации, изображен на рисунке 4. При старте предложенный алгоритм проверяет поступающие GPS-выборки (Di) до тех пор, пока GPS-сигнал отсутствует (M=0), для генерации события «начала работы фильтра», затем инициализируются вспомогательные переменные, и начинается основной цикл фильтрации. Для случая, когда достаточно двух координат – широты и долготы, Mmin берется равным 2 (Двухмерный режим позиционирования), а в качестве DOP и DOPmax используем горизонтальный фактор потери точности HDOP и HDOPmax. Для трехмерного случая Mmin = 3 (Трехмерный режим позиционирования), а в качестве DOP и DOPmax используем позиционный фактор потери точности PDOP и PDOPmax. Переменные Mold и DOPold служат для хранения значений M и DOP из предшествующей GPS-выборки. Если основное условие  выполняется, то при подготовке данных для следующего этапа фильтрации вычисляется значение точности позиционирования (ACC) в метрах. В случае если основное условие не выполнилось, а предшествующая выборка удовлетворяла данному условию, формируется событие «потери GPS-сигнала». Точность позиционирования (ACC) необходима для анализа данных на предмет  их избыточности. Аналогично DOP можно выделить 3 типа точности позиционирования: PACC – позиционная точность (в пространстве)  в метрах VACC – вертикальная точность (только по высоте) в метрах HACC – горизонтальная точность (на плоскости касательной к рассматриваемой точке поверхности эллипсоида) в метрах Известно, что PACC, VACC и HACC зависят от псевдо диапазона измерений и прямо пропорциональны значениям PDOP, VDOP и HDOP соответственно [3]. Поскольку точность позиционирования никак не определяется из GPS>-выборки, для получения значений PACC, VACC и HACC необходимо экспериментальным путем подобрать коэффициенты к соответствующим значениям PDOP, VDOP и HDOP, т.е.: PACC=PDOP*KP VACC=VDOP*KV HACC=HDOP*KH где KP,KV и KH – соответствующие коэффициенты для DOP в метрах, далее просто коэффициенты, они указываются в настройках грубого фильтра – рисунок 2.

Рекомендовать строго определенные коэффициенты нельзя, поскольку точность позиционирования зависит не только от качества сигналов со спутников, но так же от типа приемника. Диапазон рекомендуемых значений следующий: 3..5 метров для HDOP, 3..6 метров для VDOP и 3..5 метров для PDOP.

Фильтр избыточных данных в качестве критерия использует расстояние между двумя точками, которое не должно быть меньше суммы их точностей позиционирования. Графически это представлено на рисунке 5, радиус окружности (сферы в трёхмерном случае) берётся равным значению точности (ACC) для данной точки. Новые выборки проходят фильтрацию при условии, что их окружность не пересекается с окружностью предыдущей прошедшей фильтрацию выборки. Если фильтрацию подобного типа не выполнять, то помимо лишних, не несущих полезной информации данных, будет наблюдаться шум, пример которого изображен на рисунке  1.

Рис.5 Иллюстрация к алгоритму фильтрации избыточных данных о позиции наблюдаемого объекта Для вычисления расстояния (L) между двумя позициями на эллипсоиде необходимо использовать формулы [3] перевода эллипсоидных координат в декартовы и расстояния между двумя точками:

где, a и b являются параметрами эллипсоида, используемого в мировой геодезической системе отсчета WGS-84,  используемой в GPS

a = 6378137.00 м. – большая полуось

b = 6356752.31 м. – малая полуось

Рис.6 Алгоритм фильтрации избыточных данных      Алгоритм фильтрации избыточных данных приведен на рисунке 6. Da – это последняя выборка, которая удовлетворила основному условию данного фильтра. Если фильтр только начал свою работу, т.е. Da.M=0, то Di – текущая выборка, прошедшая предварительный фильтр, копируется в Dold и Da, после чего фильтр выдает первую выборку на основе DiDold обозначена копия предшествующей выборки, поступившей на вход рассматриваемого фильтра. В рабочем режиме (Da.M0) вычисляется расстояние (L) между выборками Di и Da, затем проверяется основное условие фильтра на избыточность данных. Если основное условие выполнено, т.е. выборка Di не является избыточной, переходим к формированию выходных выборок (время, позиция) и сообщений. Прежде чем сформировать выборку на основе Di, необходимо проверить совпадает ли время Di и Dold, т.е. прошла ли фильтрацию предыдущая выборка Dold, если нет, на ее основе формируем выходную выборку и событие о «начале движения объекта». Это может понадобиться для дальнейшей обработки выборок, например для более точного вычисления скорости наблюдаемого объекта, подсчета времени его нахождения в подвижном состоянии и других случаев, когда важно точно знать момент начала движения. Если основное условие фильтра избыточных данных не выполняется, но точность у новой выборки лучше, чем у последней прошедшей фильтрацию (Di.ACC < Da.ACC), данные о позиции и значения точности копируются из выборки Di в выборку Da, в некоторых случаях это помогает улучшить качество фильтрации. Например, когда по каким-то причинам точность ухудшилась с нескольких метров до нескольких десятков метров, без данного приема новые данные поступят из фильтра только после того, как объект покинет зону радиусом в эту точность, несмотря на то, что точность новых выборок будет значительно лучше. В двухмерном режиме позиционирования в качестве точности (ACC) используем горизонтальную точность (HACC), а высота (h) не используется. В трехмерном режиме используется позиционная точность (PACC). В алгоритме на рисунке 6 есть критический момент, требующий относительно больших вычислительных ресурсов. Формулы 1, 2, 3 и 4 для вычисления расстояния между двумя точками требуют использования чисел с плавающей запятой двойной точности, что является крайне сложной задачей для большинства недорогих микроконтроллеров.     Для ускорения вычислений необходимо отказаться от перевода координат из эллипсоидной системы в декартову, и проводить фильтрацию в исходной системе координат, но в этом случае не удастся вычислить расстояние между двумя точками.  Поэтому, основное условие фильтра избыточных данных примет вид (для двухмерного случая последнее выражение с участием h исключается):

 и в новом подходе используются в качестве параметров эллипсоида (рис. 2).

Рис.7 Иллюстрация к функции 

 – функция, вычисляющая соотношение градусов широты к одному метру расстояния на указанной в качестве параметра широте. Вычислено, что обратные значения данной функции меняются на протяжении от 0 до 90 градусов широты незначительно, всего на 31 см./сек.Ш. (30,715…31.025 м./сек.Ш.), это значительно меньше максимальной точности позиционирования (2 м.), поэтому  приравнивается константе, равной 1/30.87 = 0,03239391 сек.Ш./м. или °Ш./м. на любой широте. Графически смысл единицы измерения (°Ш./м.) можно увидеть на рисунке 7 – это отношение .

Рис.8. Иллюстрация к функции 

 – вычисляет соотношение градусов долготы к одному метру расстояния на указанной широте. Значение данной функции значительно меняется ближе к высоким широтам, в  переходит в бесконечность, таким образом, данный фильтр нельзя считать надёжным на полюсах, там применение проектируемого устройства маловероятно, но допустимо, если для точек полюсов предусмотреть вместо бесконечности некоторое максимальное значение. Непосредственное  вычисление заменено таблицей, и последующей линейной интерполяцией для промежуточных значений широты, но как показала практика – достаточно таблицы 1 без интерполяции промежуточных данных. Данные для таблицы рассчитаны по формулам 1,2,3,4 в качестве минимального изменения была взята 1 секунда долготы, шаг j по широте равен 1°, в качестве результата в таблицу записаны значения . На рисунке 8 значения можно представить как , очевидно, что оно значительно меняется ближе к высоким широтам.

Таблица 1 Угол по широте Секунд долготы на метр Угол по широте Секунд долготы на метр Угол по широте Секунд долготы на метр 0 0.03233884320 30 0.03731041962 60 0.06451511818 1 0.03234373632 31 0.03769406222 61 0.06653320062 2 0.03235842324 32 0.03809740622 62 0.06870351652 3 0.03238292639 33 0.03852134593 63 0.07104288695 4 0.03241728336 34 0.03896684472 64 0.07357072470 5 0.03246154695 35 0.03943494139 65 0.07630955357 6 0.03251578538 36 0.03992675665 66 0.07928565496 7 0.03258008251 37 0.04044350089 67 0.08252988653 8 0.03265453828 38 0.04098648250 68 0.08607871985 9 0.03273926892 39 0.04155711734 69 0.08997557955 10 0.03283440757 40 0.04215693920 70 0.09427257874 11 0.03294010455 41 0.04278761178 71 0.09903280806 12 >0.03305652830 >42 >0.04345094207 >72 >0.1043333822 >13 >0.03318386566 При расчётах значений иh была принята равной 100м., если говорить о точности, то L из расчёта на 1 секунду широты при изменении высоты от 0 до 1000 м. над уровнем эллипсоида изменяется всего на 5 мм, поэтому высотой можно пренебречь для расчёта этих значений. 

Рис.9 Иллюстрация к быстрому варианту алгоритма фильтрации избыточных данных Графический смысл алгоритма фильтрации избыточных данных в новом подходе примет вид представленный на рисунке 9. Относительно масштабов, фигуру со сторонами  на  можно обозначить как прямоугольник, а если задействовать третье измерение (2*VACC), то как параллелепипед. Таким образом, условием принятия новой позиции (Di) является не пересечение данных фигур, построенных вокруг точек указанных в выборках Da и Di. При таком подходе будет теряться некоторое количество потенциально полезной информации о позиции, если за идеал принять фильтрацию с использованием формул 1-4, графически представленную на рисунке 5. Для движения на плоскости потери можно приблизительно оценить так: , т.е. на 127 отсеянных выборок придется примерно 27 потенциально полезных выборок. Для движения в пространстве:  , аналогично на 191 отсеянную выборку придется около 91-й потенциально полезной выборки. Подобные потери вполне допустимы, для наблюдения за транспортом. Взамен имеем значительное сокращение вычислений, и возможность использования вычислений с фиксированной запятой, т.о. алгоритм фильтрации становится легко реализуемым и требует минимум вычислительных ресурсов.

В таблице 2 представлены результаты обработки NMEA-трека [1] с использованием разработанных алгоритмов, коэффициенты для HDOP и VDOP были взяты равными 4.5 метра, а для PDOP 3.5 метра. Прочерком в поле «тип фильтра» обозначен результат, полученный без применения фильтрации. Звёздочкой (*) отмечены типы фильтров, в которых использовался быстрый алгоритм. Длительность трека 1 час 33 минуты.

Таблица 2 Расчет для движущегося объекта Тип фильтра Путь (м.) Макс. скорость (км/ч) Число выборок – 82627,698 92,530 5597 2D 80328,531 90,918 1743 2D* 80331,885 90,918 1700 3D 80330,503 90,918 1731 3D* 80330,967 90,918 1695 2D/3D 80330,108 90,918 1733 2D/3D* 80331,885 90,918 1700 В таблице 3 приведены результаты обработки 45 минутного NMEA-трека с неподвижного объекта.

Таблица 3. Расчет для неподвижного объекта Тип фильтра Путь (м.) Макс. скорость (км/ч) Число выборок – 448,856 16,294 2719 2D 0 0 1 2D* 0 0 1 3D 0 0 1 3D* 0 0 1 2D/3D 0 0 1 2D/3D* 0 0 1

Фильтр типа 2D рекомендуется использовать в случае, когда необходимости в 3-й координате нет (не нужно очень точно считать пройденный маршрут или объект практически не передвигается по вертикали – хорошо подходит для морского транспорта и большинства наземного), это позволит уменьшить размер выходной выборки и сократить вычисления.

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

Фильтр типа 2D/3D, построен таким образом, что в предварительной фильтрации (рис. 2) используется двухмерное условие. А при фильтрации избыточных данных, если в обеих выборках (Di и DaM=3, то используется трёхмерное условие, иначе двухмерное. Используя данный тип фильтра, в формируемых выборках помимо времени и позиции, необходимо передавать M, или просто логическую переменную, определяющую использование высоты. Для таблиц 2 и 3, перед расчётом расстояния между двумя точками (формулы 1,2,3 и 4), проверялось значение M и если хотя бы в одной из выборок режим не соответствовал трехмерному, то вместо высот использовалась константа. Фильтр типа 2D/3D имеет смысл, когда необходимо знать о местоположении объекта как можно больше.

В качестве основных достоинств разработанного метода фильтрации GPS-телеметрии можно отметить: первое – возможность реализации алгоритма на маломощном и недорогом микроконтроллере. Второе – значительное уменьшение потока телеметрии, в зависимости от характера движения наблюдаемого объекта и того, какие коэффициенты выбраны для DOP. Третье – улучшается точность вычисления производных параметров (путь, скорость, время движения и др.). В ходе разработки метода фильтрации было написано специальное программное обеспечение, демонстрирующее работу GPS-модуля, фильтрацию и осуществляющее расчёты для таблиц, приведенных в статье.

Скачать программу и получить более подробную информацию можно по адресу в Интернете: «http://altmer.arts-.ru/GPSTester.htm». Список литературы: 1) ANTARIS Positioning Engine — Protocol Specification, Copyright 2003 © u-blox AG 2) Richard B. Langley, Dilution of Precision, GPS World 2000. Ссылка на перевод статьи: http://www.navgeocom.ru/gps/dop/ 3) Jean-Marie Zogg, Основы GPS, 26/03/2002, ID документа: GPS-x-02007-R.pdf Список терминов:

DOP – Фактор потери точности (англ. Dilution of precision) — термин, использующийся в области систем глобального позиционирования для описания качества геометрического взаиморасположения спутников друг относительно друга.

GPS – NAVSTAR GPS (англ. Navigation Satellite Time and Ranging, Global Positioning System — Измерение дальности и времени по навигационному спутнику, Глобальная система позиционирования) — спутниковая система навигации, часто именуемая просто GPS. Система разработана, реализована и эксплуатируется Министерством обороны США.

GPS-телеметрия – Подразумеваются данные с GPS приемника удаленного объекта.

NMEA – «National Marine Electronics Association» — полное название «NMEA 0183» — текстовой протокол связи морского (как правило, навигационного) оборудования между собой. Стал особенно популярен в связи с распространением GPS приёмников, использующих этот стандарт.

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