CC BY
46
Интерполяция распределённых данных горизонталей для получения
цифровой модели рельефа
В. С. Тарасян, Н. В. Дмитриев Уральский государственный университет путей сообщения, Екатеринбург
Аннотация: В работе рассматривается интерполяция плана горизонталей, полученного в рамках преобразования изображения топографической карты в геоинформационную систему местности, для получения цифровой модели рельефа. Авторами предложен метод, с помощью которого можно отразить рельефность местности, учитывая распределённый характер данных. Метод сравнивался с методами линейной интерполяции, интерполяции кубическими сплайнами и интерполяцией Сибсона по критериям максимального отклонения, среднеквадратического отклонения и скорости обработки. Экспериментальное исследование проводилось на сгенерированных с помощью метода diamond-square местностях. Разработанный метод показал улучшение показателей качества по сравнению с существующими методами и рекомендован к использованию при интерполяции плана горизонталей.
Ключевые слова: интерполяция, геоинформационная система, топографическая карта, цифровая модель рельефа, распределённые данные.
Введение
В настоящее время всё чаще встаёт вопрос об оптимальном использовании местности. Окружающее пространство может использоваться для решения задач различной проблематики, при этом обычно создаются системы поддержки принятия решений различного характера, имеющие в своей основе геоинформационную систему (ГИС), осуществляющие интеллектуальную обработку пространственно распределённых данных. С их помощью решаются задачи экологического [1, 2], управленческого [3], проектно-транспортного [4, 5] характера. В частности, для оптимального проектирования объектов транспортной и логистической инфраструктуры ГИС должна содержать цифровую модель рельефа (ЦМР).
Чтобы получить ГИС необходимо проанализировать источники информации, содержащие данные об определённой местности. Такими источниками могут быть данные со спутника, кадры аэрофотосъёмки и сканированные топографические карты (ТК). Для получения ЦМР наиболее
оптимальным способом видится использование именно топографических карт, так как они содержат полную информацию о рельефе местности, изображённую с помощью горизонталей с некоторым шагом. Относительно данных со спутника, с помощью ТК можно получить более точную модель. Сами изображения можно получить в открытом доступе в сети Интернет, что обуславливает экономическую эффективность относительно получения аэрофотоснимков, получаемых с помощью запуска беспилотных летательных аппаратов.
К настоящему моменту созданы алгоритмы и программное обеспечение, выполняющее сегментацию изображения ТК, комплексное распознавание объектов, выполняющее восстановление плана горизонталей [4, 5]. Последняя операция необходима, так как не у всех горизонталей существуют подписи со значениями высот. В итоге по текущей топографической карте может быть получена матрица, соответствующая пикселям карты, у которой в местах расположения горизонталей располагаются величины соответствующей высоты, аналогично проставлены значения геодезических опорных пунктов, а в остальных точках данные отсутствуют.
В случае проектирования транспортных линий коммуникации приходится использовать несколько соседних ТК, но существующие программные продукты (ЛгсОК, Мар1п&, Панорама) не
позволяют работать с данными большого объёма, так как требуют большого объёма оперативной памяти.
Постановка проблемы
Требуется по существующим данным получить цифровую модель рельефа, то есть интерполировать существующие данные на всю область карты. Сложность заключается в том, что данные с одной стороны имеют распределённую структуру (расстояние между двумя соседними
горизонталями может составлять 100 и более пикселей), а с другой стороны -компактную, так как данные по горизонталям располагаются непрерывно.
Кроме того в реальных картах горизонтали могут в некоторой степени не соответствовать реальной высоте местности, для более точной передачи рельефности [6]. Поэтому требуется создать алгоритм, который бы лучше всего передавал именно рельефность, то есть степень соответствия формы горизонталей с промежуточным шагом форме горизонталей с основным шагом. Кроме того алгоритм не должен вызывать такие артефакты, как «бычий глаз» - рисунка концентрических контурных линий вокруг точек выборки. В качестве пространства используется тот же самый растр, а не триангуляционная нерегулярная сеть (TIN), для более точного дальнейшего анализа.
Существуют различные методы интерполяции пространственно распределённых данных, отличающиеся принципом выбора точек и функцией интерполирования [7, 8]. По принципу выбора точек различают алгоритмы N-ближайших соседей и выбор соседей из окрестности определённого радиуса. Большинство существующих алгоритмов предлагают использовать не все существующие точки, а только некоторые, например, расположенные друг от друга на некотором минимальном расстоянии, при этом становится невозможным передать рельеф. Такой современный алгоритм, как кригинг не может быть применён, так как количество точек, которые нужно интерполировать, слишком большое.
Разработанный алгоритм
Так как пространство карты разделено непрерывными горизонталями, то они могут быть либо замкнуты сами на себя, либо на границы карты, при этом получается множество непересекающихся подпространств. Для выполнения интерполяции данные подпространства были разделены на три типа: находящиеся между двумя горизонталями, внутри замкнутой
горизонтали локального экстремума, между горизонталью и границей карты. Для определения типа подпространства достаточно определить тип хотя бы одной его точки, для этого воспользуемся поиском в ширину ближайших точек с данными. Если найденные точки относятся к двум различным горизонталям, то подпространство точек принадлежит к первому типу. Если найденные точки относятся к одной горизонтали, то подпространство точек принадлежит ко второму типу. Если найденные точки относятся к одной горизонтали и границам карты, то подпространство точек принадлежит к третьему типу.
Учитывая принцип сохранения рельефности для интерполяции точек первого типа, выберем некоторое количество ближайших соседних точек из ближайшей горизонтали, а от дальней горизонтали возьмём одну ближайшую точку, но в соответствующем количестве. Для вычисления интерполяционной высоты воспользуемся методов обратных взвешенных расстояний. При этом около первой горизонтали всегда будут выбраны ближайшие точки, что будет способствовать сохранению рельефности, а вторая горизонталь будет способствовать изменению высоты в целом.
Если считать, что точка находится вблизи первой горизонтали, то формула средневзвешенной функции будет принимать значения (для простейшего случая 0 = 1/ dJ 0):
2 о=
2 *уп +^ г 1*У" —2,—+г2 1 ^'= —,0 —2,0 1 ^'= 2 г *р+г 2
1 + П п —2,0 р + 1
V" —+— уп ^+1
г ,0 2,0 —
где 2\, 22 - значения высот первой и второй горизонталей; - расстояние от г-ной точки первой горизонтали до текущей, —2,0 - среднее расстояние от
текущей точки до второй горизонтали, р - среднее отношение расстояния ко второй горизонтали к расстоянию до точек первой.
Полученная формула с одной стороны отражает влияние формы рельефа вблизи горизонталей (при ^ 0: р ^ 0, 20 ^ 2{), а с другой стороны показывает, что посередине между двумя горизонталями (р = 1) будет получаться среднее значение высот.
Для точек второго типа важно наличие отметки высоты в точке локального экстремума. Остальные точки можно интерполировать, зная крутизну склона к нормалям в точках ограничивающей горизонтали, вычисленные после обработки внешних точек (они принадлежат к точкам первого типа) и факт, что в точке экстремума крутизна равна нулю. При интерполяции данного типа точек считаем, что крутизна склона изменятся равномерно:
а(Г) = а 0 — Ц * Г,
где г - расстояние от точки на горизонтали в сторону точки экстремума, а0 -значение крутизны склона в точке на горизонтали, ц - линейный коэффициент.
Тогда, с учётом того, что угол крутизны очень маленькая величина, относительная высота выразится по формуле:
* V2
АН = |Бта(х)ёг « |а(=а
ц * х о-х 2
—™ * х —
Г ,с 1!г 2
* Ц * Гс
= а0* Г,с--е-
0 2
где г- расстояние от точки на горизонтали до точки экстремума. Из данной формулы можно вывести значение линейного коэффициента.
Фактически данная интерполяция является квадратичной, так как возвышенность или низина приближается полиномом второй степени по направлению от каждой точки горизонтали к точке экстремума. Так как площади таких участков невелики, что получаемую точность считаем
достаточной, но её нужно привести в растровый формат, так как при данном алгоритме сетка получается непрямоугольная, а полярная.
Точки третьего типа можно интерполировать, если рядом находятся отметки высот (алгоритм аналогичен интерполяции по второму типу), иначе возникает задача экстраполяции, для которой сложно будет вычислить ожидаемую точность. Данную проблему можно решить, если объединить несколько соседних топографических карт: тогда в этом типе останутся только подпространства, лежащие у границ объединённой карты.
Для сравнения с разработанным алгоритмом будет использоваться метод линейной интерполяции, кубическая интерполяция сплайнами, интерполяция Сибсона [9].
Для того чтобы проверить качество и быстродействие методов будут созданы искусственные местности различных размеров, построенные с помощью широко использующегося алгоритма для создания псевдопланетарных поверхностей diamond-square [10], затем из них будут получены карты горизонталей с точками экстремумов. Интервал высот будет нормализован к [0; 1], шаг горизонталей зададим равным 0.05, что примерно соответствует такой местности, как Уральские горы для карт масштаба 1: 50000. Пример карты местности и полученный план горизонталей представлен на рис. 1: для карты чёрный цвет соответствует минимальной высоте, белый - максимальной.
В качестве критериев оценивания алгоритмов будем использовать максимальное отклонение, среднеквадратическое отклонение (в процентах относительно шага горизонталей) и скорость обработки, равную времени обработки на количество пикселей. Все методы будут запускаться на одной ЭВМ несколько раз, чтобы нивелировать возможные отклонения по длительности выполнения, вызванные использованием процессора фоновыми программами.
Рис. 1. - Полутоновая карта местности с соответствующим планом горизонталей
Результаты экспериментального исследования
Результаты экспериментального исследования представлены в таблице №1. На рис. 2 представлены гистограммы отклонений предсказанной высоты от действительной в процентах от шага горизонталей для методов в порядке, указанном в таблице №1.
Таблица № 1
Критерии оценивания методов интерполяции
Метод интерполяции Максимальное отклонение, % Среднеквадратичное отклонение, % Скорость обработки, с/Мп
Линейный 92 13.2 4.2
Сибсона 83 12.4 11.8
Кубический 91 14.7 2.7
Предложенный 76 12.0 13.2
Рис. 2. - Гистограммы отклонений Заключение
Предложенный алгоритм показал улучшение по критериям качества при некотором снижении быстродействия. Стоит учитывать, что в условиях преобразования значений точек рельефа в горизонтали происходит неизбежная потеря данных, поэтому получить меньшие значения среднеквадратичной ошибки видится сложнодостижимым. Планируется дальнейшее исследование на возможность использования в алгоритме данных соседних горизонталей. Максимальные отклонения возникают в областях, описанных в предложенном алгоритме как точки второго типа: их количество незначительно в общем объёме карты.
Литература
1. Давлетбакова З. Л., Абдуллин А. Х., Павлов С. В. Обработка пространственной информации о границах санитарных зон полигонов отходов на основе методов нечеткой логики // Инженерный вестник Дона, 2013, № 4. URL: ivdon.ru/ru/magazine/archive/n4y2013/1887.
2. Сладкова Ю. М., Лычагин А. А., Сурков Ф. А. Комплексный картографический анализ экологического состояния г. Ростова-на-Дону //
Инженерный
вестник
Дона, 2007, №2. URL:
ivdon.ru/ru/magazine/archive/n2y2007/34.
3. Ефремова О.А. Применение системного подхода к анализу проблемы использования пространственной информации для поддержки принятия решений региональными органами исполнительной власти//Инженерный вестник Дона, 2014, №2. URL: ivdon.ru/ru/magazine/archive/n2y2014/2371.
4. Тарасян В. С., Дмитриев Н. В. Восстановление плана горизонталей при обработке топографических карт//Современная наука: актуальные проблемы теории и практики. Серия: Естественные и технические науки. М.: Научные технологии. 2017. № 1. с. 56-60.
5. Тарасян В. С., Дмитриев Н. В. Алгоритм автоматического восстановления значений высот горизонталей//Молодежь в науке: новые аргументы: сб. научн. работ V Международного молодежного конкурса. Ч. I. 2016. с. 176-179.
6. Кошкарев А.В., Тикунов В.С. Геоинформатика. Под ред. Д. В. Лисицкого. М.: Картгеоцентр - Геодезиздат. 1993. 213 с.
7. Saveliev A. A., Mukharamova S.S., Chizhikova N.A., Budgey R., Zuur A.F. Spatially continuous data analysis and modelling. Analysing ecological data. Springer-Verlag, 2007. Chapter 19. pp. 341-372.
8. Геостатистика: теория и практика / Савельева Е.А., Демьянов В.В.; под ред. Р.В. Арутюняна; Ин-т проблем безопасного развития атомной энергетики РАН. М.: Наука, 2010. 327 с.
9. Sibson R. A Brief description of natural neighbor interpolation// Interpolating multivariate data. New York: John Wiley & Sons. 1981. Chapter 2. pp. 21-36.
10. Jacob Olsen Realtime Procedural Terrain Generation. Technical Report, University of Southern Denmark. 2004. 20 p.
References
1. Z. L. Davletbakova, A. Kh. Abdullin, S. V. Pavlov Inzenernyj vestnik Dona (Rus), 2013, № 4. URL: ivdon.ru/ru/magazine/archive/n4y2013/1887.
2. Sladkova Yu. M., Lychagin A. A., Surkov F. A. Inzenernyj vestnik Dona (Rus), 2007, №2. URL: ivdon.ru/ru/magazine/archive/n2y2007/34.
3. Efremova O.A. Inzenernyj vestnik Dona (Rus), 2014, №2. URL: ivdon.ru/ru/magazine/archive/n2y2014/2371.
4. Tarasyan V. S., Dmitriev N. V. Sovremennaya nauka: aktual'nye problemy teorii i praktiki. Seriya: Estestvennye i tekhnicheskie nauki. M.: Nauchnye tekhnologii. 2017. № 1. pp. 56-60.
5. Tarasyan V. S., Dmitriev N. V. Molodezh' v nauke: novye argumenty: sb. nauchn. rabot V Mezhdunarodnogo molodezhnogo konkursa. Ch. I. 2016. pp. 176179.
6. Koshkarev A.V., Tikunov V.S. Geoinformatika [Geoinformatics]. M.: Kartgeotsentr-Geodezizdat. 1993. 213 p.
7. Saveliev A.A., Mukharamova S.S., Chizhikova N.A., Budgey R., Zuur A.F. Spatially continuous data analysis and modelling. Analysing ecological data. Springer-Verlag, 2007. Chapter 19. pp. 341-372.
8. Savel'eva E.A., Dem'yanov V.V. Geostatistika: teoriya i praktika [Geostatistics: theory and practice]. M.: Nauka, 2010. 327 p.
9. Sibson R. A Brief description of natural neighbor interpolation. Interpolating multivariate data. New York: John Wiley & Sons. 1981. Chapter 2. pp. 21-36.
10. Jacob Olsen Realtime Procedural Terrain Generation. Technical Report, University of Southern Denmark. 2004. 20 p.
