[Все статьи]

Сглаживание GPS треков сплайнами

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

Выбор вида функции аппроксимации

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

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

Возьмем для примера трек представленный на рисунке 1.

Рисунок 1

Аппроксимируем эту кривую полиномом третьей степени вида:

Результат аппроксимации представлен пурпурной линией на рисунке 2.

Рисунок 2

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

Рассмотрим для сравнения другой трек, вместе с аппроксимирующей его линией (см. Рисунок 3).

Рисунок 3

Очевидно, что в этом случае результат аппроксимации значительно хуже, чем на рисунке 2.

Внимательный читатель должен заметить, что трек на рисунке 2 и трек на рисунке 3 на самом деле одинаковы… Просто трек на рисунке 3 повернут зеракльно относительно прямой Y = X. Если бы мы искали аппроксимирующую функцию, как зависимость X от Y, то и аппроксимирующая кривая развернулась так же и результат остался бы удовлетворительным.

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

Первое что приходит в голову: попробовать аппроксимировать двумя способами: y=f(x) и x=f(y) затем из двух выбрать наилучшее приближение. Однако, это требует двойных вычислительных расходов.

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

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

Наилучшим выходом является переход к параметрическому виду кривой:

Существует несколько общепринятых способов аппроксимации, соответствующих данной формуле. Мы рассмотрим два из них: cardinal splines и B-Splines.

Особенностью обоих способов является то, что кривая составляется из кусков, которые стыкуются в точках сопряжения.

Cardinal spline и B-spline различаются по своим свойствам. Особенность Cardinal spline является то, что кривая проходит точно через все точки трека, а в точках сопряжения соседних кусков имеет связность С1 (непрерывную первую производную). B-spline же проходит примерно через точки трека, но взамен на такую неточность предлагает связность типа С2 (непрерывную вторую производну).

Оба сплайна по сути являются параметрическими кривыми третьего порядка, т.е. в конечном счете вычисляются по формулам вида:

где t изменяется от 0 до 1 для каждого куска сплайна.

Существует более удобная с практической точки зрения запись формул для обоих сплайнов:

здесь:  Pi – точки трека, представленные как вектор-строка;

t- параметр;

P(t) – результирующая функция, возвращающая точку на сплайне;

M – главная матрица, которой собственно и отличаются данные два способа интерполяции/аппроксимации.

Как видно из формулы, сплайн строится по 4-м точкам. Результирующая кривая проходит при этом где-то в районе точек Pi – Pi+1. Кстати, прелестью данной формулы является то, что она может без изменений применяться как на плоскости так и в трех и более мерном пространстве.

Матрица M для B-Spline имеет вид:

Для Cardinal spline:

Параметр s в последней матрице влияет на закругленность линий, при s=0 сплайн вырождается в ломаную линию. Рекомендованное значение s=0.5. Хотя, можно поиграть с этим значением и посмотреть, какой вариант будет лучше.

Пример приближения нашего трека обоими способами представлен на рисунках 4 и 5.

Рисунок 4. B-Spline

Рисунок 5. Cardinal spline при s=0.5

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

Практическое применение

Первое применение, которое приходит в голову, это отрисовка трека плавной линией. Тут, полагаю, особых объяснений не требуется. Для каждого куска пробегаем по t от 0 до 1 с маленьким шагом и соединяем на экране получившиеся точки.

Второе, что можно определить имея трек выраженный сплайном, это азимут, которые может быть вычислен как вектор, с координатами dx/dt и dy/dt.

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

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

Радиус кривизны для параметрически заданной кривой (а мы помним, что рассматриваемые сплайны после перемножения матриц превращаются в параметрические кривые) вычисляется по формуле:

 

Производные берутся по параметру t. Формула выглядит страшно, но с вычислительной точки зрения не все так плохо: первая производная от x и y используется два раза, так что ее можно посчитать один раз, а потом подставить, а вторая производная от кубического полинома вообще линейная функция, так что, здесь вычисления еще не такие страшные.

Здесь следует сделать важное замечание, Cardinal spline не пригоден для расчета радиусов кривизны, так как его вторая производная в точках сопряжения кусочков разрывна, что приведет к тому, что радиус кривизны посчитанный в конце одного кусочка и в начале следующего за ним будут сильно отличаться. Для расчетов радиуса кривизны нужно применять B-Spline

До сих пор вычисления были относительно простые, зубодробительная математика начинается, когда мы вспомним, что параметр сплайна это обезличенное t, которое просто изменяется от 0 до 1 в пределах каждого кусочка, в то время, как для практических целей целесообразно адресовать точки на дороге, как расстояние в метрах от начала. Таким образом необходимо иметь механизм нахождения кусочка сплайна и параметра t в пределах него, куда попадет точка, скажем, отстоящая на 120 метров от начала трека.

Вспомним формулу длины дуги для параметрически заданной кривой:

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

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

Тогда остается определить значение параметра t, которое соответствует искомой точке.

Для этого нужно решить следующее уравнение

Здесь l – это остаток от искомого расстояния от начала трека. Т.е. допустим что нам нужно найти точку расположенную на расстоянии 120 метров от начала. Первые три кусочка имеют длины соответственно 35, 40, 50. Очевидно, что искомая точка должна находиться в пределах третьего кусочка, на расстоянии l=120-35-40=45 метров от его начала.

Поговорим о подынтегральном выражении из двух последних формул.

Функции x(t) и y(t) являются полиномами третьего порядка относительно t, следовательно их первые производные по t будут квадратными полиномами.

Под знаком корня стоят квадраты производных, а квадрат полинома второй степени будет полиномом четвертой степени. Сумма полиномов четвертой степени так же является полиномом четвертой степени.

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

Здесь коэффициенты аi вычисляются после вычисления производных и раскрытия всех скобок, как множители перед соответствующими степенями t.

Проблема в том, что нет аналитической формулы интеграла от данного выражения (по крайней мере, если верить MathCAD и Mapple).

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

Для решения это проблемы, автор предлагает разложить последнее выражение в ряд Маклорена (Тейлора), до 4-го-5-го члена (энтузиасты могут попробовать больше).

Пример разложения до 4-х членов приведен ниже:

Опять вроде бы страшно, но если вспомнить что ai это константы, вычисленные на предыдущих этапах, то эта страшная формула преобразуется к виду:

Здесь bi равны множителям перед соответствующими степенями t из предыдущего выражения.

Таким образом мы нашли приблизительное выражение для подынтегральной функции в интеграле длины дуги.

Интеграл от полинома находится тривиально. Интегральное уравнение:

придет к виду:

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

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

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

[На верх страницы]
Loading... Загрузка...