close

Вход

Забыли?

вход по аккаунту

?

Моделирование полей видимости в среде ArcGIS средствами приложения «Картометрия».

код для вставкиСкачать
Вестник СПбГУ. Сер. 7. 2012. Вып. 1
УДК 911.2:528
Е. А. Паниди
МОДЕЛИРОВАНИЕ ПОЛЕЙ ВИДИМОСТИ В СРЕДЕ ARCGIS
СРЕДСТВАМИ ПРИЛОЖЕНИЯ «КАРТОМЕТРИЯ»
Введение
Моделирование полей видимости — одна из классических задач в геоинформатике.
Оно может быть применено во всех случаях, когда необходимо рассчитать обеспечена
ли прямая видимость между двумя точками, находящимися на земной поверхности
или над ней, либо выделить окружающую территорию, для которой будет обеспечена
прямая видимость из заданной точки пространства. Наиболее часто задача построения
полей видимости возникает при проектировании и моделировании различных систем
связи и управления (сотовых сетей, радиолокационных систем и т. д.).
На наличие видимости между парой точек оказывают влияние четыре группы факторов:
— рельеф земной поверхности, который непосредственно влияет на наличие прямой видимости между двумя точками, расположенными на земной поверхности или
над ней;
— кривизна земной поверхности, влияние которой сказывается на больших расстояниях и сходно с влиянием рельефа;
— физические факторы, связанные с характером среды, такие как атмосферная рефракция и характеристики прозрачности воздуха;
— физические факторы, связанные с «наблюдателем», для которого моделируется
область видимости, такие как мощность сигнала и его затухание, уровень помех и т. д.
Факторы, входящие в последние две группы, в сущности, не являются постоянными, их влияние зависит от состояния атмосферы в конкретный промежуток времени
и от характеристик применяемой для наблюдения аппаратуры и может быть уменьшено. Оба вида физических факторов могут лишь ограничивать видимость, но не влиять
на её принципиальное наличие или отсутствие. Таким образом, влияние этих факторов в каждом конкретном случае будет меняться и должно быть учтено индивидуально
в зависимости от предметной области, в рамках которой производится моделирование.
В связи с этим, алгоритм моделирования полей видимости, описанный в данной работе, позволяет учесть только влияние рельефа и кривизны земной поверхности.
Однако, при необходимости учёта в модели физических факторов, данный алгоритм
может быть достаточно просто модифицирован и дополнен соответствующими операциями. Использованный автором подход основан на пространственном геометрическом моделировании и допускает использование дополнительных правил и условий
при расчёте геометрических характеристик зрительного луча.
Данная статья является продолжением серии статей, посвящённых разрабатываемой автором библиотеки функций для проведения картометрических изысканий в среде ГИС. [1, 2, 3].
© Е. А. Паниди, 2012
121
7-1-2012.indd 121
28.02.2012 14:41:12
Постановка задачи
Для построения поля видимости использована модель рельефа, представленная
в виде регулярной сети точек (DEM — Digital Elevation Model — Цифровая Модель
Рельефа), что связано с удобствами описания и анализа данной модели. Очевидно, что
при использовании такой модели, корректность описания рельефа будет возрастать
с уменьшением шага сетки.
Рассчитывать видимость будем в пространственной системе координат с учётом
кривизны поверхности земного эллипсоида. А высота рельефа вычисляется с использованием DEM. При этом для описания модели рельефа используем геодезическую систему координат.
Высоты точек, находящихся на земной поверхности, могут быть определены по
DEM неоднозначно. Проблема состоит в том, каким образом будет интерполирована
искомая высота по известным высотам узлов сетки рельефа. Задача интерполяции в
данном случае сводится к вычислению высоты рельефа как функции двух переменных — координат искомой точки. Решена эта задача может быть, по крайней мере, тремя способами:
— метод ближайшего соседа [4],
— билинейная интерполяция [5],
— бикубическая интерполяция [6].
Первый способ наиболее прост и требует минимальной затраты вычислительных
ресурсов. Он заключается в том, что за высоту определяемой точки принимается высота ближайшего к ней узла сетки рельефа.
Второй метод заключается в проведении линейной интерполяции сначала вдоль одной координатной оси, затем вдоль другой. Высота искомой точки (Н) в данном случае
будет вычислена по высотам четырёх соседних узлов сетки рельефа. Если координаты
этих узлов принять за (0,0), (0,1), (1,0) и (1,1), то формула вычисления выглядит следующим образом [10]:
(1)
Бикубическая интерполяция аналогично билинейной интерполяции сводится
к проведению кубической интерполяции последовательно вдоль каждой из координатных осей. В общем виде вычислительная функция в данном случае может быть записана следующим образом [6]:
3
3
H ( x, y ) = ∑∑ ai , j x i y j ,
(2)
i =0 j =0
где i, j — номера в пределах блока пикселей.
При этом интерполяция производится по 16 соседним узлам сетки рельефа. Вычисленная функция (2) является гладкой, в отличие от первых двух способов.
Для нахождения коэффициентов ai , j в уравнение необходимо подставить координаты соответствующих 16 узлов сетки и решить получившуюся систему линейных алгебраических уравнений.
122
7-1-2012.indd 122
28.02.2012 14:41:12
Для интерполяции высот рельефа наиболее желательным представляется третий
метод, так как в целом рельеф является достаточно гладкой поверхностью.
В целом, несмотря на то, что задача построения полей видимости вполне стандартна в геоинформатике и подобная функциональность реализована практически во всех
универсальных геоинформационных системах, зачастую реализация подобных функций имеет особенности, которые в отдельных случаях не позволяют произвести вычисления требуемым образом.
Так, например, алгоритм [7], реализованный в ГИС «Карта» [8] не позволяет учитывать кривизну земной поверхности при построении поля видимости.
В геоинформационной системе ArcGIS [9] в свою очередь, нет возможности задать
при расчетах способ интерполяции высоты рельефа.
Алгоритм
В основе вычислительного алгоритма лежит геометрический подход, основная идея
которого состоит в следующем: между наблюдаемой точкой и точкой наблюдения строится отрезок (визирный луч), для которого контролируется наличие или отсутствие
пересечений с поверхностью рельефа. При этом все построения производятся в пространственной декартовой системе координат, что позволяет одновременно учитывать
кривизну земной поверхности и влияние рельефа.
Исходными данными, помимо модели рельефа, служат следующие значения:
— геодезическая система координат, на основе которой будут производиться вычисления;
— геодезические координаты точки наблюдения;
— высота точки наблюдения над земной поверхностью;
— высота наблюдаемых точек (геодезическая или относительная над поверхностью
рельефа) в случае если она отлична от высоты поверхности рельефа;
— границы сектора видимости.
Последний пункт необходимо пояснить. Дело в том, что поле зрения наблюдателя
в вертикальной плоскости может быть ограничено определённым сектором. Например, радиолокационные станции часто имеют диаграмму направленности, которая
ограничена 0° снизу и 45° сверху (углы измеряются от горизонта).
На основании исходных данных на первом этапе вычислений рассчитываются размеры матрицы (количество ячеек), в которую будут записаны вычисленные значения.
Каждой ячейке будет присвоено значение 0, если видимость для этой точки отсутствует
и значение 1, если видимость имеется.
Так как матрица видимости строится на основе матрицы рельефа, и за шаг сетки
принимается шаг сетки рельефа, то, по сути, задача сводится к выделению в матрице
рельефа блока, в пределах которого будет располагаться линия горизонта и созданию
аналогичной этому блоку матрицы видимости.
Положение линии горизонта в данном случае рассчитывается для поверхности выбранного земного эллипсоида, на основе положения точки наблюдения и высоты наблюдаемой точки. При этом вычисляется экстремальное положение линии горизонта, при котором она будет находиться на максимально возможном удалении от точки
наблюдения с учётом поверхности эллипсоида. На следующих этапах это позволит
исключить несоответствие размеров поля видимости и матрицы видимости, так как
с учётом рельефа, зона видимости всегда будет находиться в пределах вычисленного
таким образом положения линии горизонта.
123
7-1-2012.indd 123
28.02.2012 14:41:12
На основе геодезических координат точки наблюдения и модели рельефа находим
геодезическую высоту точки, при этом учитываем её высоту над земной поверхностью,
которую прибавляем к полученной геодезической высоте. После чего, переходим от
геодезических координат точки наблюдения к её пространственным декартовым координатам [10]:
X = ( N + H ) ⋅ cos B ⋅ cos L
Y = ( N + H ) ⋅ cos B ⋅ sin L
(
)
Z = (1 − e 2 ) ⋅ N + H ⋅ sin B,
где
N=
a
2
1 − e sin 2 B
— радиус кривизны первого вертикала в данной точке;
e = 2α − α 2 e — первый эксцентриситет эллипсоида; a — большая полуось эллипсоида; α — сжатие эллипсоида; B, L, H — геодезические координаты точки.
Затем, начиная от положения точки наблюдения на модели рельефа в четырёх направлениях (на север, юг, запад и восток) для каждой ячейки матрицы, в пределах соответствующей строки или столбца вычисляем геодезические координаты наблюдаемой
точки.
При этом за широту и долготу точки принимаем широту и долготу соответствующей ячейки матрицы, а высоту, в зависимости от исходных данных либо интерполируем по модели рельефа, либо берём постоянную, заданную ранее. От полученных таким
образом геодезических координат, переходим к пространственным декартовым координатам наблюдаемой точки.
В результате для каждой выбранной ячейки мы получим положение наблюдаемой
точки и точки наблюдения в декартовой системе координат, после чего останется проконтролировать отсутствие пересечений отрезка, соединяющего эти точки, с поверхностью рельефа, а следовательно, и нахождение наблюдаемой точки в пределах горизонта. Сделать это можно следующим образом: вычислим координаты центра отрезка,
воспользовавшись стандартными формулами:
X=
X í à÷ + X êî í
Y +Y
Z + Z êî í
; Y = í à÷ êî í ; Z = í à÷
.
2
2
2
Здесь X íà÷ , Yíà÷ , Z íà÷ и X êîí , Yêîí , Z êîí — соответственно координаты начальной и
конечной точек отрезка.
И от полученных пространственных координат перейдём к геодезическим координатам [10]:
вычислим вспомогательную величину D :
D=
X 2 + Y 2 , значение которой обусловит дальнейший ход решения.
124
7-1-2012.indd 124
28.02.2012 14:41:12
Если D = 0, то В и Н вычислены по следующим формулам:
Если D > 0 , то L = arcsin
Y
. При этом:
D
если Y < 0, X > 0, то L = 2π − L,
если Y < 0, X < 0, то L = 2π + L,
если Y > 0, X < 0, то L = π − L,
если Y > 0, X > 0, то L = L.
Если Z = 0, то B = 0, H = D − a,
Если : Z ≠ 0, то В определяем способом итераций:
⎛
p sin ( 2b ) ⎞
⎟,
Bi = c + Si , Si = 0 = 0, Si = arcsin ⎜
⎜ 1 − e 2 sin 2 B ⎟
i
−
1
⎝
⎠
где
r=
X 2 +Y 2 + Z2
,
⎛Z⎞
c = arcsin⎜ ⎟ ,
⎝r⎠
2
e a
.
p=
2r
Затем вычисляем H :
(3)
В случае если высота визирного луча будет положительной, повторим аналогичные вычисления для центральных точек, полученных на предыдущем этапе частей
визирного луча. Затем поделим пополам получившиеся четверти и так далее. В случае
если все вновь получаемые точки будут находиться над поверхностью эллипсоида,
остановить вычисления целесообразно в тот момент, когда визирный луч будет разбит на части, линейные размеры которых не превышают шаг матрицы рельефа. В таком случае можно считать, что визирный луч не пересекает поверхность эллипсоида
и наблюдаемая точка находится в пределах горизонта. Если же на одном из этапов
высота визирного луча окажется отрицательной, это будет означать, что наблюдаемая
точка и все точки, удалённые на большее расстояние, находятся за пределами горизонта и для оставшихся ячеек матрицы рельефа в данном направлении вычисления
можно не производить.
Повторив аналогичные вычисления для северного, южного, восточного и западного
направлений, определим соответствующие границы матрицы видимости. Таким образом выделим по матрице рельефа ячейки, которые войдут в матрицу видимости.
На заключительном этапе проверим, не попадают ли в пределы горизонта какиелибо ячейки, непосредственно прилегающие к границам матрицы видимости с внешней
стороны. В случае если такие ячейки будут обнаружены, расширим пределы матрицы
125
7-1-2012.indd 125
28.02.2012 14:41:12
видимости. Если же таких ячеек не обнаружится, остановим вычисления, построим
матрицу видимости и сохраним её в памяти компьютера.
Производя вычисления для ограниченного количества ячеек матрицы рельефа и не
прибегая к их полному перебору, мы получаем возможность построить матрицу видимости, которую затем сможем заполнить значениями.
Дальнейший вычислительный процесс можно разделить на следующие этапы:
— построение визирного луча из точки наблюдения в точку, для которой производится расчёт видимости;
— контроль нахождения визирного луча в пределах сектора видимости;
— контроль пересечений визирного луча с поверхностью рельефа;
— заполнение поля видимости значениями.
В качестве наблюдаемой точки будем последовательно выбирать каждую ячейку
матрицы видимости. Указанные четыре этапа вычислений повторим для каждой ячейки матрицы видимости.
Визирный луч построим в пространственной декартовой системе координат из точки наблюдения в наблюдаемую точку уже известным способом.
Попадает ли построенный луч в сектор видимости — определим, воспользовавшись
формулой для вычисления угла между вектором, направленным из точки наблюдения
в наблюдаемую точку, и вектором, направленным вдоль нормали к поверхности эллипсоида из точки наблюдения от поверхности эллипсоида [11]:
Здесь ( x1 , y1 , z1 ) — координаты первого вектора, ( x 2 , y 2 , z 2 ) — координаты второго вектора.
Нахождение дополнения этого угла до 90º в пределах заданного сектора видимости
будет означать, что и визирный луч находится в пределах сектора видимости.
Если вычисленный угол места не будет находиться в пределах заданных границ
сектора видимости, то вычисляемой ячейке матрицы видимости присвоим значение 0
и перейдём к следующей ячейке. Если же визирный луч пройдёт в пределах сектора
видимости, перейдём к следующему этапу вычислений.
Поиск пересечений визирного луча с поверхностью рельефа будем производить аналогично поиску пересечений визирного луча и поверхности эллипсоида. Воспользовавшись широтой и долготой вычисляемой ячейки матрицы видимости, определим по модели рельефа высоту поверхности рельефа и сравним с полученной геодезической высотой центральной точки визирного луча. Если высота визирного луча окажется больше
высоты рельефа, повторим аналогичные вычисления для центральных точек половин
визирного луча, затем поделим пополам получившиеся четверти и так далее. Вычисления остановим в тот момент, когда визирный луч будет разбит на части, линейные
размеры которых не превышают шаг матрицы рельефа. Если на одном из этапов высота
визирного луча окажется меньше высоты рельефа, это будет означать, что визирный
луч пересекает поверхность рельефа и дальнейшие вычисления можно не производить,
а вычисляемой ячейке матрицы видимости следует присвоить значение 0. Если же все
полученные точки визирного луча будут находиться над поверхностью рельефа, это
126
7-1-2012.indd 126
28.02.2012 14:41:12
будет означать наличие прямой видимости на наблюдаемую точку. В таком случае вычисляемой ячейке матрицы видимости следует присвоить значение 1.
Пользовательский интерфейс
Описанный алгоритм реализован автором при разработке библиотеки функций для
компьютерной картометрии «Картометрия». Данная библиотека создана в виде приложения к геоинформационной системе ArcGIS [9] и является полностью интегрированной в состав её картографического редактора — ArcMap. Написана библиотека на
языке программирования C#.
Интерфейс модуля, выполняющего моделирование полей видимости, представлен
на рис. 1.
Рис. 1. Пользовательский интерфейс
Модуль позволяет производить как единичные вычисления с заполнением исходных данных в соответствующие поля ввода, так и пакетные вычисления. В пакетном
режиме программа считывает исходные данные из текстового файла.
На рис. 2 представлен результат вычисления поля видимости аэродромного радиолокатора аэропорта Мурманска.
Выводы
Представленный алгоритм не учитывает физические характеристики атмосферы
и особенности «наблюдателя», но позволяет корректно учесть геометрические условия,
возникающие при расчёте поля видимости. Использование в качестве поверхности относимости поверхности земного эллипсоида позволяет повысить точность моделирования по сравнению с вычислениями на поверхности сферы или в плоскости картографической проекции. Использование различных способов интерполяции высот по
модели рельефа даёт возможность учесть характер рельефа исследуемой территории.
127
7-1-2012.indd 127
28.02.2012 14:41:12
128
7-1-2012.indd 128
28.02.2012 14:41:12
Рис. 2. Отображение рассчитанного поля видимости в окне ArcMap
Данный алгоритм может быть достаточно просто реализован программно как в среде современных ГИС, так и в самостоятельных программных продуктах.
Литература
1. Афанасьев В. А., Павлова О. А., Паниди Е. А, Щербаков В. М. Автоматизированная система
генерации номенклатуры российских топографических карт и пересчёта координат между различными системами // Статья в сборнике Теория и практика эколого-географических исследований. СПб.: ТИН, 2005. 604 с.
2. Паниди Е. А. Алгоритм перепроецирования растровых изображений средствами программного комплекса «Топографический калькулятор» // Вестник С.-Петербург. ун-та. 2008.
Сер. 7. Вып. 3. С. 141–146.
3. Паниди Е. А., Голубков С. Н., Павлова О. А., Щербаков В. М. Автоматизированная система
для анализа основных метрических свойств картографического изображения // Вестник С.-Петербург. ун-та. 2008. Сер. 7. Вып. 4. С. 188–193.
4. Интернет-энциклопедия Wikipedia URL: http://en.wikipedia.org/wiki/Nearest-neighbor_
interpolation (Дата обращения: 24.12.2010)
5. Интернет-энциклопедия Wikipedia URL: http://en.wikipedia.org/wiki/Bilinear_interpolation
(Дата обращения: 24.12.2010)
6. Интернет-энциклопедия Wikipedia URL: http://en.wikipedia.org/wiki/Bicubic_interpolation
(Дата обращения: 24.12.2010)
7. Форум сайта ЗАО КБ “Панорама” (разработчик ГИС “Карта”) URL: http://www.gisinfo.ru/
forum/phpBB2/viewtopic.php?t=6946&highlight=%E2%E8%E4%E8%EC%EE%F1%F2%E8 (Дата обращения: 24.12.2010)
8. Сайт ЗАО КБ “Панорама” URL: http://www.gisinfo.ru (Дата обращения: 24.12.2010)
9. Сайт Environmental Systems Research Institute (разработчик ArcGIS) URL: www.esri.com
(Дата обращения: 24.12.2010)
10. Государственный Стандарт Российской Федерации ГОСТ Р 51794-2001 Системы координат.
11. Выгодский М. Я. Справочник по высшей математике. М.: Наука, 1977. 872 с.
Статья поступила в редакцию 26 сентября 2011 г.
129
7-1-2012.indd 129
28.02.2012 14:41:13
Документ
Категория
Без категории
Просмотров
6
Размер файла
761 Кб
Теги
среды, поле, моделирование, средства, arcgis, картометрия, видимости, приложение
1/--страниц
Пожаловаться на содержимое документа