CC BY
99
2008
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Управление, вычислительная техника и информатика
№ 3(4)
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
УДК 681.03 513
К.Г. Гульбин
МОДЕЛИРОВАНИЕ ПРОЦЕССА ПОЛУЧЕНИЯ ЛАЗЕРНО-ЛОКАЦИОННЫХ ДАННЫХ
В статье рассмотрены различные аспекты получения модельных лазернолокационных данных: задание и внутреннее представление модели местности, моделирование процесса лазерной локации и визуализация его результатов. Предложены практические способы моделирования отражений лазерного луча от частично проницаемых объектов, учета ошибок лазерной локации и оптимизации работы моделирующего алгоритма.
Ключевые слова: лазерная локация, цифровая модель рельефа, 3В-объекты.
Лазерная локация - один из новейших методов воздушной стереотопографи-ческой съемки, основанный на совместном использовании систем глобального позиционирования и лазерной дальнометрии, что позволяет получать массивы («облака») точек поверхности Земли и наземных объектов. Эти массивы содержат миллионы точек с плотностью до нескольких точек на один кв. м. Метод лазерной локации имеет множество приложений: топографические, геодезические, гидрографические (специализированные лазерные локаторы позволяют проводить съемку поверхности и дна водоемов), инженерно-изыскательские, лесоустроительные, экологические, электроэнергетические и другие [1]. Информация может быть получена на основе первичных данных лазерной локации разными способами: применением статистических методов, выделением групп точек разной морфологической принадлежности, построением цифровой модели рельефа, ортофотомозаики, распознаванием топографических объектов и т.д. В связи с широтой спектра приложений лазерной локации, необходима моделирующая процесс получения лазерно-локационных данных система. Наличие такой системы позволяет проводить детальное тестирование алгоритмов обработки данных, отрабатывать частные случаи. В данной статье рассмотрены различные аспекты создания такой системы. Также описаны принятые автором для её создания технические решения.
В рамках поставленной задачи выделены три основные подзадачи: задание сцены (рельефа местности с находящимися на ней объектами), моделирование процесса съемки лазерно-локационных данных на заданной сцене и визуализация сцены, смоделированных данных и промежуточных результатов работы тестируемых алгоритмов.
1. Внутреннее представление и визуализация сцены в моделирующей системе
Один из вариантов задания модели местности - представить рельеф триангуляцией с привязанными к ней топографическими объектами, состоящими из наборов простых геометрических фигур в трехмерном евклидовом пространстве. Таким образом, здание можно представить в виде набора плоских пространственных фигур (скаты крыши, стены, окна - плоскости). Столб, ствол дерева - цилиндры, крона дерева - параболический конус и т.д. Поверхности разделим на абсолютно непроницаемые (здания, поверхность Земли) и частично проницаемые с разными коэффициентами (кроны деревьев, травяная подстилка и т.п.). Также припишем поверхностям коэффициенты яркости, определяющие мощность отражаемых от нее сигналов. Для визуализации предлагается конвертировать сцену в файл или набор файлов общепринятого BD-формата. Это позволит переложить на стороннее программное обеспечение задачи отображения сцены и навигации по ней. Одним из наиболее подходящих BD-форматов является формат языка Virtual Reality Modeling Language (VRML), описывающий наборы BD-объектов в трехмерной декартовой системе координат. Пример визуализации простой модельной сцены представлен на рис. 1.
Рис.1. VRML-визуализация модельной сцены
Для внутреннего структурирования сцены в моделирующей системе удобно представить ее в виде древовидной структуры часть - целое, в которой индивидуальные и составные объекты трактуются единообразно (такой способ представления называется Composite - компоновщик). Объекты наследуются от одного класса, в интерфейсе которого определены такие методы, как привязка к данному экземпляру дочернего объекта (например, привязка здания к рельефу местности), получение типа данного объекта, нахождение точек пересечения с лучом и т.п. При вызове одного из этих методов для составного объекта он рекурсивно применяется ко всем его подчастям.
2. Процесс лазерно-локационной съемки
В зависимости от назначения лазерно-локационных данных, к характеристикам распределения точек (продольному и поперечному расстояниям между точками, равномерности распределения и т.п.) предъявляются определенные требования. Для получения удовлетворяющего этим требованиям облака точек перед выполнением лазерно-локационной съемки подбирают значения параметров полета и выполнения съемки [1].
Рассмотрим эти параметры и основные геометрические соотношения, характеризующие процесс лазерно-локационной съемки в трехмерной декартовой системе координат XYZ. Пусть летательный аппарат движется вдоль оси X на высоте, равной H, со скоростью V (рис. 2). Лазерный локатор наклонен в плоскости XZ на фиксированный угол ю, в перпендикулярной XZ плоскости 0 его угол наклона (у ) меняется от -ф до +ф и обратно с частотой F. С частотой f (f >> F) генерируются импульсы лазерного излучения и регистрируются координаты их отражений. При этом группы точек, соответствующих одному полному изменению угла наклона в плоскости 0 (линии сканирования), располагаются зигзагообразно, а расстояние между начальной и конечной точками каждой линии сканирования -ширина захвата (W) - зависит от параметров съемки следующим образом:
W = 2 .
cos W
Расстояние между соответствующими точками смежных линий сканирования колеблется от нуля до величины Sx, называемой продольным интервалом сканирования:
Sx = V / f.
Поперечный же интервал сканирования
S, = 2 Wf .
J F
А количество точек на единицу поверхности (плотность облака) D вычисляется как
D =- 1
Z
Рис. 2. Геометрические соотношения, характеризующие процесс лазерно-локационной съемки
На рис. 3 представлено облако точек отражения от сцены (см. рис.1), полученное в моделирующей системе. Для моделирования использовались следующие значения параметров: ю = 15°, Н = 300 м, высота здания - 30 м. Траектория дви-
жения сканирующего летательного аппарата проходит справа от здания по направлению к наблюдателю. При таком сканировании область сцены слева от здания недоступна для сканирующих импульсов. Этой области не принадлежит ни одна из точек полученного облака (рис. 3).
Рис. 3. Визуализация полученного на основе модельной сцены облака точек
Схематически процесс получения облака точек рис. 3 представлен на рис.4.
Рис. 4. Схема процесса получения облака точек лазерной локации
3. Моделирование процесса лазерно-локационного сканирования
Смоделировав сцену и задав значения параметров лазерно-локационной съемки, перейдем к моделированию процесса лазерно-локационного сканирования. В начальный момент времени сканерный блок находится в точке {X=Х0, У = У0, 2 = Н}, полет должен завершиться в точке {X = Хь У = Уь 2 = Н }. Курсовой угол носителя в начальный момент времени
С \
а = агссо8
X! - X
4іX - Х0 )2 + (^ - У0 )2
два других угла его ориентации (тангаж в и крен у ) равны 0, а вектор движения
соответственно таков: т = {X -Х0, У\ - У0, 0}. Сканирующий блок посылает импульсы через промежутки времени, равные 1 // За это время:
- координаты носителя изменятся на расстояние V// по направлению вектора движения;
- угол наклона у - на 4 -ф- (Р / /).
Используя значения наклона сканирующего блока и ориентации носителя, произведем серию поворотов вертикального вектора (0, 0, -1) и получим вектор, по которому направлен соответствующий данному импульсу лазерный луч:
Г 0 ї
I = АВ
0
-1
где
А =
008 а -8Іп а о 1 (1 0 0 1 ' 008 у 0 - 8ІП уЛ
8ІП а 008 а 0 0 008 в Іп 8І 1 0 1 0
0 0 1, V 0 8ІП в 008 в , ч 8ІП у 0 008 у ,
' сое ю 0 - вій юЛ (1 0 0 >
В = 0 1 0 0 сов у - вій у
ч вій ю 0 сов ю , V 0 вій у сов у )
4. Нахождение точек пересечения лазерного луча с элементами сцены
Рассмотрим, как определить набор координат всех точек пересечения луча с пространственными многоугольниками, составляющими элементы сцены. Для ускорения работы алгоритма не будем искать точки пересечения луча с гранями сцены, с которыми он заведомо не может пересечься. Например, воспользуемся для этого индексацией пространства с использованием Я-дерева. Поскольку элементы сцены имеют небольшой разброс координат по оси 2, ограничимся индексацией по координатам X и У.
Шаг 1. Вычислим экстенты (охватывающие прямоугольники) всех граней сцены.
Шаг 2. Корневому узлу Я-дерева соответствует вся область сцены. Ему приписан список всех граней сцены. Выберем его в качестве текущего и перейдем на шаг 4.
Шаг 3. В список граней текущего узла заносятся те, экстент которых пересекает экстент его области. Эти грани - подмножество граней из списка родительского узла.
Шаг 4. Проверим текущий узел на необходимость более глубокой индексации. Критерием этого является пересечение с экстентами более чем к граней сцены. Если этот критерий удовлетворяется, переходим на шаг 5. Иначе делаем текущим следующий узел и повторяем шаг 4.
Шаг 5. К текущему узлу добавляется 4 узла-потомка. Область, соответствующая текущему узлу, разбивается на 4 равных части, по одной для каждого потомка.
Шаг 6. Каждый из потомков поочередно обрабатывается (начиная с шага 3).
Покажем, как вычисляется точка пересечения луча с выпуклой пространственной гранью. Общая точка содержащей эту грань плоскости
Ах + Ву + Сг + Б = 0
и прямой х = х0 +1, у = у0 + т?, г = г0 + ы
находится по формуле
(А1 + Вт + Сп)г + Ах0 + Ву0 + С20 = 0.
Для установления принадлежности точки пересечения выпуклой грани преобразуем координаты ее вершин и точки пересечения с лучом в координаты плоскости, содержащей эту грань. Критерием нахождения точки внутри плоской выпуклой фигуры является ее нахождение по одну сторону от всех векторов, полученных последовательным обходом ребер этой фигуры.
При замене системы координат ХУ2 новой системой ХнУн2н с тем же началом старые координаты точки {х, у, z} выражаются через новые {хн, ун, гн} формулами
{х, у, г} = {хн, ун, гн}А,
где
' ^(4, I) СОвОН, I) ^(4, I) ^
А = ^(4, I) С0в( /, I) ^(4,1)
чс^(4, ?) со^^.7н > *) с08(4 > *),
а (гн, г) - угол между новой и старой осью абсцисс, (/н, г) - угол между новой осью
ординат и старой осью абсцисс и т.д.[2]. Таким образом, новые координаты точки выражаются через старые:
{хн > Ун, 2н } = {Х У, г)АГ- .
Косинус угла между двумя векторами а1 {хьу1, z1}и а2 {х2,у2, Z2} находится из их скалярного произведения:
Х1 х2 + У1У 2 + 2122
008 ф =
I 2 . 2 I 2 . 2
Vх: + у і + гі V х2 + У2 + г2
Осями базовой системы координат являются вектора (1, 0, 0), (0, 1, 0) и (0, 0, 1), а в качестве осей новой системы координат примем: вектор, образованный одной из сторон треугольника, вектор нормали к поверхности и вектор, дополняющий их до прямоугольной системы координат. Этот третий вектор является вектором нормали по отношению к плоскости, образованной первыми двумя. Пусть заданы два вектора а1 {хьу\,гх} и а2 {х2,у2,г2}, лежащих в одной плоскости. Они перпендикулярны вектору нормали п {хп,уп, гп}, то есть соответствующие скалярные произведения равны нулю. Имеем
Х1 Хп + УіУп + 212п = 0
Ш- ()
Х2 Хп + У2 Уп + 22 2п = 0]
Помножим верхнее уравнение на х2, и вычтем из него нижнее, помноженное на х1. Получим
(УХ2 - У 2 Х1) Уп = ( 22 Х1 - 21Х2 ) *п'
()
Х2 Хп + У Уп + 22 2п = 0
Приняв zn = у1х2 - у2х1, получим уп = z2x1 - z1x2 из (2). Подставив теперь zn и уп во второе из уравнений системы (1), найдем и хп:
_ У2 (22 х1 - 21х2 ) + 22 (У1х2 - У2 х1)
хп _ .
х2
Воспользовавшись вышеприведенными выкладками, мы можем вычислить вектор нормали к плоскости треугольника, взяв в качестве векторов а1 и а2 его стороны, а затем вычислить и третий осевой вектор новой системы координат.
5. Моделирование проницаемости
Лазерный луч расширяется по мере распространения. В связи с этим лазерный луч отражается не точечно, имеет место пятно отражения. Поэтому при наличии неполных препятствий лазерно-локационные данные могут содержать по нескольку точек отражения одного лазерного импульса.
Будем последовательно рассматривать точки пересечения смоделированного луча с элементами сцены по возрастанию дальности от сканирующего устройства (точки 1,2,3 и 4 на рис.5). Добавляем очередную точку в результирующий набор, если соответствующая поверхность является непроницаемой. В противном случае она сигнализирует о переходе в частично проницаемую среду, тогда между ней и следующей по отдаленности точкой могут произойти одно или несколько частичных отражений лазерного луча. Смоделируем их координаты и добавим в набор, и как только число точек в нем превысит 3 - 5, завершим выполнение алгоритма (такая ситуация соответствует практически полному отражению лазерного луча от нескольких неполных препятствий).
Для описания в моделирующей системе неполных препятствий поверхностям сцены приписан коэффициент проницаемости. Будем интерпретировать его как вероятность неотраженного распространения лазерного луча на 1 м в глубь расположенной за поверхностью среды. Пусть кроне дерева, точками пересечения лазерного луча с которой являются точки 1 и 2 на рис. 5, соответствует коэффициент проницаемости к, а расстояние между ними п метров. Вероятность отражения лазерного луча при распространении на 1 метр в глубь кроны дерева равна р = 1 - к, а при распространении на т метров р1,т. Поделим отрезок между точками 1 и 2 на п/т подотрезков длиной т метров. На каждом из них с вероятностью р1,т произойдет случайное событие - частичное отражение лазерного луча. За координату такого отражения примем координаты случайной, равномерно распределенной точки соответствующего отрезка.
6. Моделирование погрешностей
При проведении лазерно-локационной съемки неизбежно возникают погрешности. Рассмотрим, как учесть этот фактор при моделировании лазернолокационного измерения.
В большинстве случаев точность измерений гарантируется производителями лазерных локаторов и составляет 15 - 20 см в абсолютных геодезических координатах. Поэтому простейшим способом учета погрешности является отклонение координат смоделированной точки отражения на вектор погрешности. Направление этого вектора определяется двумя углами поворота случайной величины, распределенной по равномерному закону. Его длина также случайная, распределенная по нормальному закону с СКО ~ 6 - 7 см.
Основными источниками погрешности лазерно-локационного измерения является совокупность векторов ошибок GPS-приемника, лазерного дальномера и инерциальной системы (таблица). Поэтому для более точного учета такой погрешности эти ошибки необходимо промоделировать отдельно. Ошибки GPS-приемника и инерциальной системы могут возникать при определениях ими координаты и углов ориентации носителя, производимых дискретно с частотами Г1 = 1 - 20 Гц и Г2 = 50 - 200 Гц соответственно. Примем за ошибку навигационного блока вектор ошибки со случайной длиной, нормально распределенной с СКО ~ 3-4 см. За ошибку инерциальной системы - вектор ошибки с длиной, также нормально распределенной с СКО ~ 10 см. Значения этих векторов будут переопределяться с частотами ^ и _Р2. Ошибку в определении наклонной дальности до поверхности смоделируем для каждого импульса как вектор случайной нормально распределенной с СКО ~ 3 - 5 см длины. Направлением этого вектора является направление соответствующего конкретному импульсу лазерного луча.
Основные источники погрешности лазерно-локационного измерения
Определяемый параметр Источник Точность
Пространственные координаты носителя ОР8 8 - 10 см
Наклонная дальность Лазерный дальномер 10 - 15 см
Ориентация носителя Инерциальная система 1 - 2 мрад (ошибка позиционирования 15 - 30 см при высоте съемки 300м)
Таким образом, для получения модели лазерно-локационных данных с учетом погрешностей прибавим к каждой из смоделированных точек соответствующие ей векторы ошибок.
ЛИТЕРАТУРА
1. Данилин И.М., Медведев Е.М., Мельников С.Р. Лазерная локация земли и леса: Учеб. пособие. Красноярск: Инст. леса им. В.Н.Сукачева СО РАН, 2005.
2. Выгодский М.Я. Справочник по высшей математике. М.: Наука, 1973.
Статья представлена кафедрой теоретических основ информатики факультета информатики Томского государственного университета, поступила в научную редакцию 25 января 2008 года.
