Ю.Л. Костюк, А.Л. Фукс
ПОСТРОЕНИЕ ЦИФРОВОЙ МОДЕЛИ РЕЛЬЕФА МЕСТНОСТИ НА ОСНОВЕ СТРУКТУРНЫХ ЛИНИЙ И ВЫСОТНЫХ ОТМЕТОК
Для цифровой модели рельефа, заданного структурными линиями и высотными отметками, вводится новый тип ограничений, позволяющий более полно учитывать влияние горизонталей на форму рельефа. Предлагается эффективный алгоритм построения цифровой модели на основе триангуляции с ограничениями данного типа.
Важным объектом исследования современных геоин-формационных систем является земной рельеф. Как правило, рельеф задается нерегулярными наборами высотных отметок и структурных линий, которые получаются с помощью методов дистанционного зондирования или векторизации картографических материалов. Высотные отметки обычно представляют локальные экстремумы и другие характерные точки рельефа. Структурные линии определяют множества точек с резким изменением наклона рельефа (границы оврагов, обрывов, береговые линии) или одинаковыми высотами (горизонтали, изолинии), т.е. накладывают дополнительные ограничения на форму рельефа. Все линии задаются наборами узловых точек и для упрощения дальнейшей обработки считаются ломаными.
Важнейшая задача при работе с рельефом заключается в построении его цифровой модели (ЦМР), т.е. цифрового представления с помощью прямоугольной или треугольной сетки, в узлах которой заданы высоты. Цифровая модель позволяет получать производные данные как для анализа, так и для построения ортофотопланов местности на основе космо- и аэроснимков.
Качество цифровой модели определяется, прежде всего, качеством исходных данных. Типичными недостатками исходных структурных линий являются:
1. Наличие топологических нарушений (пересечение горизонталей разных уровней, самопересечение линий и др.). Такие ошибки легко обнаружить, однако корректно исправить их можно только интерактивно. Далее мы будем предполагать, что они уже исправлены.
2. Осцилляции горизонталей (изолиний), которые трудно измерить с высокой точностью. Устранению этих недостатков посвящена работа [4], в которой предлагаются методы сглаживания этих линий.
3. Некорректное определение высот и неполное задание изолиний. Такие ошибки можно обнаружить только в процессе получения цифровой модели.
При построении ЦМР на основе прямоугольной сетки необходимо рассчитать высоты в ее узлах. Для этого обычно применяются различные алгоритмы локальной аппроксимации и интерполяции по значениям высот в нескольких соседних исходных точках. Данная сетка благодаря своей простоте позволяет упростить алгоритмы анализа и обработки рельефа. Однако она имеет два существенных недостатка:
1. Узлы сетки являются расчетными, поэтому значения высот в них всегда приближенные. В то же время исходные точки - узлы структурных линий, локальные экстремумы и другие характерные точки рельефа - обычно определяются с высокой точностью, но в получаемой цифровой модели они утрачиваются.
2. В цифровой модели, использующей прямоугольную сетку, практически невозможно учитывать ограничения, налагаемые структурными линиями.
Модель рельефа на основе треугольной сетки строится с помощью триангуляции с ограничениями. Она свободна от указанных недостатков, но известные алгоритмы не учитывают в полной мере влияние горизонталей на форму рельефа. В результате получаемая модель искусственно спрямляет рельеф, игнорируя априорно известные изломы, и поэтому является недопустимо грубым приближением. Для уст-
ранения этих недостатков в большинстве геоинформацион-ных систем используются интерактивные методы корректировки исходных линий и/или полученной треугольной сетки. Однако последнее средство практически бесполезно, если сетка содержит многие тысячи треугольников.
В настоящей работе предлагается способ построения кондиционной цифровой модели рельефа в автоматическом режиме.
ТИПЫ ОГРАНИЧЕНИЙ ТРИАНГУЛЯЦИИ
При задании рельефа наборами точек и линий построение его цифровой модели на основе треугольной сетки включает два шага [1].
На первом шаге по проекциям всех исходных высотных отметок и узлов структурных линий на плоскости ХОУ строится триангуляция Делоне. Задание высот в вершинах триангуляции определяет систему пространственных треугольников - простейшую кусочнолинейную поверхность, интерполирующую рельеф.
На втором шаге производится перестроение триангуляции, позволяющее включить в модель рельефа все структурные линии. Каждая линия должна целиком располагаться на интерполяционной поверхности, поэтому необходимые условия перестроения (ограничения триангуляции) формулируются очень просто: каждый отрезок структурной линии должен совпадать с одним из ребер сетки пространственных треугольников. Далее ограничения этого типа будем называть слабыми.
Способ включения отрезка структурной линии, не вошедшего в начальную триангуляцию Делоне (рис. 1, а), зависит от дополнительных требований к ЦМР. Если необходимо, чтобы сетка на ХОУ всегда оставалась триангуляцией Делоне, то отрезок нужно разбить на несколько частей. Точки разбиения станут новыми вершинами триангуляции, а каждая часть отрезка - ребром треугольника (рис. 1, б). Если в качестве вершин треугольной сетки можно использовать только исходные точки, то начальную триангуляцию Делоне нужно перестроить так, чтобы каждый отрезок исходной линии стал ребром сетки (рис. 1, в).
а б в
Рис. 1. Включение отрезка: а - начальная триангуляция, б -отрезок разбивается на части, в - отрезок становится ребром
Учет слабых ограничений при построении модели рельефа является необходимым, но недостаточным,
если набор структурных линий содержит горизонтали. Чтобы в этом убедиться, нужно на треугольной сетке со слабыми ограничениями рассчитать изолинии с шагом изменения высоты в два раза меньшим, чем у исходных горизонталей (рис. 2). Из-за наличия в триангуляции горизонтальных треугольников (заштрихованных), все вершины которых лежат на одной изолинии, на исходных уровнях появятся паразитные линии, а на промежуточных уровнях линии окажутся чрезмерно спрямленными.
а б
Рис. 2. Изолинии на треугольной сетке со слабыми ограничениями (выделены горизонтальные участки модели рельефа): а - исходные, б - расчетные
Если набор исходных структурных линий содержит горизонтали, то интерполяционная поверхность должна не только полностью включать эти линии, но и удовлетворять следующим естественным ограничениям:
1. Точки изолинии не могут быть локальными минимумами или максимумами поверхности.
2. Изолиния не должна быть границей горизонтального участка поверхности, для выделения плоских участков рельефа нужно использовать структурные линии других типов (граница промышленной зоны, береговая линия озера).
3. Любая исходная отметка, высота которой совпадает с уровнем какой-либо изолинии, должна быть либо вершиной исходной горизонтали, либо являться вырожденной изолинией из одной точки.
4. Любой горизонтальный отрезок, соединяющий вершины изолиний одного уровня и целиком лежащий на поверхности, должен либо совпадать с отрезком изолинии (и, следовательно, ребром сетки), либо соединять конечные вершины одной или двух горизонталей.
5. Если некоторая линия начинается и заканчивается на изолиниях соседних уровней и 2В ,
< 2В , не пересекает структурных линий и целиком лежит на интерполяционной поверхности, то в любой ее внутренней точке М должно выполняться условие
< 2М < 2В .
В соответствии с данными требованиями введем сильные ограничения триангуляции:
на треугольной сетке выполняются слабые ограничения;
высотная отметка, которая не является узлом структурной линии, но имеет такую же высоту, как некоторая изолиния, считается вырожденной изолинией, содержащей одну конечную точку;
любое горизонтальное ребро пространственной триангуляции, соединяющее вершины изолиний, либо совпадает с отрезком изолинии, либо соединяет конечные вершины одной или двух изолиний.
Такие ограничения гарантируют отсутствие плоских горизонтальных участков на уровнях исходных горизонталей.
Построение треугольной сетки, удовлетворяющей сильным ограничениям, нужно начинать с триангуляции со слабыми ограничениями. Эта триангуляция может содержать горизонтальные ребра, на которых сильные ограничения не выполняются (для краткости будем называть такие ребра недопустимыми). Все треугольники, содержащие одно, два или три недопустимых ребра, необходимо перестроить, однако предварительно следует убедиться, что высоты исходных линий и точек определены корректно.
НАЧАЛЬНАЯ ПРОВЕРКА КОРРЕКТНОСТИ ЗАДАНИЯ ВЫСОТ
Ошибки в задании высот линий и точек характерны тем, что их практически невозможно обнаружить визуально или выделить все с помощью только автоматизированных проверок. Обработка высот должна включать несколько шагов и может потребовать корректировки исходных данных в интерактивном режиме. Начальным шагом, на котором выявляются самые очевидные ошибки, является проверка перепадов высот. Для ее проведения необходимо построить триангуляцию со слабыми ограничениями. Очевидно, что при правильном задании исходных данных перепад высот на любом ребре триангуляции не может превышать значения шага изменения высот для основных изолиний. Нарушение этого условия свидетельствует о неправильном задании высот или отсутствии некоторых горизонталей в исходном наборе линий. Но существуют и два исключения:
Если концы незамкнутых линий не выходят точно на границу триангуляции, то ребра на границе могут принадлежать очень вытянутым треугольникам и соединять точки, расположенные на большом расстоянии друг от друга. Необходимо либо вывести концы линий на границу, либо исключить из триангуляции граничные треугольники с недопустимым перепадом высот.
На участках с почти вертикальными обрывами большой высоты основные изолинии практически сливаются, поэтому некоторые из них могут быть пропущены. Такая ситуация не является ошибочной, более того, попытка восстановить все изолинии может привести к тому, что из-за совпадения их проекций будет невозможно построить триангуляцию со слабыми ограничениями.
При проверке сильных ограничений триангуляции потребуется перестраивать треугольники с недопустимыми горизонтальными ребрами. Для этого бывает необходимо определять знак приращения высоты во внутренних точках горизонтальных треугольников, что возможно только при корректном задании высот изолиний. Проверка корректности задания высот основана на простом правиле: точки изолинии не могут быть локальными минимумами или максимумами поверхно-
сти. Для триангуляции со слабыми ограничениями это можно сформулировать следующим образом: если AB
- отрезок изолинии (zA = zB), который является общим ребром треугольников ABC и BAD, причем zc ф Za и Zd Ф Za , тогда (zc - Za ) • (Zd - Za ) < 0 . Другими словами, если zc < zA, то zD > zA , и наоборот.
Пусть ABC - горизонтальный треугольник, т.е. zA = zB = zc . Для определения знака приращения высоты в его внутренних точках необходимо проверять окрестность ABC, переходя к смежным треугольникам строго по горизонтальным ребрам до тех пор, пока не будет найден негоризонтальный треугольник. Пусть это треугольник DEF, и zE = zF = zA , а zD ф zA . Если переход к DEF проходил только по недопустимым ребрам, то знак приращения относительно уровня zA во внутренних точках ABC такой же, как в точке D. При каждом пересечении отрезка изолинии знак приращения меняется на противоположный. Если высоты изолиний заданы корректно, то знак приращения не зависит от порядка проверки окрестности ABC. В этом можно убедиться, проверив все треугольники на границе горизонтального участка, в который входит ABC. Рис. 3 иллюстрирует данную проверку: знаки приращений для пар точек D,E и G,H должны быть разными.
E.
мыми ребрами, причем данные цепочки начинаются и заканчиваются треугольниками с одним недопустимым ребром.
Пусть некоторая цепочка начинается с треугольника ЛВС, а заканчивается треугольником иУЖ, причем ВС и УЖ - недопустимые ребра, т.е. zB = zC = zУ = zW , а точки В, С, У и Ж являются вершинами изолиний одного уровня. Предположим, что ломаная, проходящая из Л в и через центры недопустимых ребер, имеет длину Ь, а длина самого короткого недопустимого ребра цепочки равна I. Вычислим на ломаной новую точку М по следующим правилам:
Если zЛ Ф zи, то полагаем, что высота вдоль ломаной от Л до и изменяется линейно. Тогда М - это равноудаленная от Л и и точка ломаной (рис. 4, а), и
^ = ^л + %)/2 .
Если zЛ = zи Ф zB, то поверхность имеет седло-вую точку. Будем считать, что седловой точкой М является точка пересечения ломаной и самого короткого недопустимого ребра цепочки (рис. 4, б и в), причем Zв * Ь + z л * I Ь +1
Если zЛ = zи = zв , то цепочку образуют треугольники внутри замкнутой изолинии. Естественно предположить, что соответствующий участок рельефа также является почти горизонтальным. Будем считать, что наибольшее отклонение от горизонтальной плоскости z = zв достигается в точке М - центре самого длинного недопустимого ребра (рис. 4, г), а высота zM = zв ±Лz, где знак отклонения определяется так же, как и для треугольников с тремя недопустимыми ребрами.
U+»JV
Рис. 3. Проверка знака приращения высоты для центра треугольника ABC (выделены изолинии и треугольники с недопустимыми ребрами)
ПОСТРОЕНИЕ ТРИАНГУЛЯЦИИ С СИЛЬНЫМИ ОГРАНИЧЕНИЯМИ
Треугольная сетка, на которой выполняются слабые ограничения, может содержать треугольники с одним, двумя или тремя недопустимыми ребрами. Горизонтальный треугольник ABC с тремя недопустимыми ребрами на уровне zA интерполирует участок рельефа, который реально также является почти горизонтальным. Будем считать, что центр ABC (точка O) отклоняется от z A по вертикали на фиксированную величину Az. Для определения знака отклонения необходимо проверить треугольники в окрестности ABC (алгоритм проверки приведен выше). Задав в O высоту zA ±Az , можно заменить ABC на три треугольника ABO, BCO, CAO, содержащих по одному недопустимому ребру.
Перестроим все треугольники с тремя недопустимыми ребрами. После этого в триангуляции останутся цепочки смежных треугольников с двумя недопусти-
б
C A
г
Рис. 4. Проверка цепочек треугольников с недопустимыми ребрами и добавление новой точки М: а - между изолиниями разных уровней, б и в - на участках с седловой точкой, г - внутри замкнутой изолинии
Вычисленная точка М должна стать новой вершиной триангуляции, и это потребует соответствующего перестроения треугольной сетки. Для исключения возможного вырождения треугольников точку М следует добавлять, как в итеративном алгоритме триангуляции Делоне [2], но перестраиваться должны только недопустимые ребра. Полученная система треугольников уже не будет триангуляцией Делоне, зато в ней не будут нарушены слабые ограничения. После добавления точки М цепочка недопустимых треугольников разбивается на две более короткие (возможно, пустые) це-
а
почки, а общее число недопустимых ребер уменьшается. Процесс заканчивается, когда недопустимых ребер в триангуляции не останется.
Алгоритм проверки сильных ограничений триангуляции. Предполагается, что по исходному набору точек и линий построена триангуляция со слабыми ограничениями.
Шаг 1. Проверка всех треугольников сетки, выделение и отметка всех недопустимых ребер.
Шаг 2. Перестроение всех треугольников с тремя недопустимыми ребрами.
Шаг 3. Цикл, пока есть цепочки треугольников с недопустимыми ребрами:
- обработка очередной цепочки,
- вычисление новой точки M и включение ее в триангуляцию.
Конец алгоритма
Обработка очередной цепочки включает выделение начального треугольника и проход по недопустимым ребрам до конечного треугольника. При этом производится расчет длины ломаной вдоль цепочки, а также поиск самого длинного и самого короткого недопустимого ребра. Добавление точки M в триангуляцию является элементарным шагом, так как треугольник, в который попадает M, известен заранее. В целом приведенный алгоритм напоминает алгоритм быстрой сортировки Хоара [3]: включение точки M и разбиение цепочки на две более короткие - это полный аналог разделения сортируемого массива опорным элементом. Поэтому совпадут и оценки трудоемкости: если в цепочку входят m треугольников, то для ее полного перестроения потребуется в среднем O(m log m) операций. Общая трудоемкость будет вполне приемлемой даже в том случае, когда исходная триангуляция с n вершинами содержит O(n) треугольников с недопустимыми ребрами.
Для проверки корректности исходных данных и качества цифровой модели можно сравнить исходные и расчетные горизонтали. Для этого необходимо построить триангуляцию с сильными ограничениями и на полученной поверхности рассчитать изолинии исходных (основных) и промежуточных уровней. О возможных пропусках линий, а также ошибках при определении высот свидетельствуют:
- «лишние» участки основных расчетных изолиний, которым не соответствуют исходные линии;
- пересечения или касания расчетных изолиний, указывающие на наличие седловых точек поверхности;
- слишком спрямленные промежуточные расчетные линии, форма которых должна соответствовать форме двух соседних основных изолиний;
- изолинии очень малого размера, построенные вокруг высотных отметок.
При наличии существенных расхождений между исходными и расчетными горизонталями может потребоваться корректировка исходных линий и/или их высот и перестроение модели рельефа.
На рис. 5 приводятся исходные и расчетные изолинии участка местности, построенные на основе триангуляции с различными типами ограничений. Основные изолинии представлены сплошными, а дополнительные - штриховыми линиями, при этом исходные горизонтали выделены толстыми линиями. Для расчета линий строилась цифровая модель рельефа на основе триангуляции со слабыми (рис. 5, а) и сильными (рис. 5, б) ограничениями. В последнем случае также проводилось предварительное визуальное сглаживание исходных горизонталей.
а б
Рис. 5. Наборы исходных и расчетных изолиний участка рельефа: а - слабые, б - сильные ограничения триангуляции
Рис. 5, а показывает, что если при построении ЦМР учитываются только обычные (слабые) ограничения, то расчетные изолинии могут оказаться совершенно не похожими на исходные. И только учет сильных ограничений позволяет получить приемлемые результаты.
ЗАКЛЮЧЕНИЕ
Предложенный в статье новый тип сильных ограничений триангуляции позволяет более полно учитывать влияние горизонталей на форму рельефа. Эффективный алгоритм построения триангуляции с сильными ограничениями обеспечивает получение кондиционной цифровой модели рельефа в автоматическом режиме.
ЛИТЕРАТУРА
1. Скворцов А.В. Триангуляция Делоне и ее применение. Томск: Изд-во Том. ун-та, 2002. 128 с.
2. Препарата Ф., ШеймосМ. Вычислительная геометрия: Введение: Пер. с англ. М.: Мир, 1989. 478 с.
3. Вирт Н. Алгоритмы и структуры данных: Пер. с англ. М.: Мир, 1989. 360 с.
4. Костюк Ю.Л., Фукс А.Л. Предварительная обработка исходных данных для построения цифровой модели рельефа местности // Вестник ТГУ. 2003. № 280. С. 281-285.
Статья представлена кафедрой теоретических основ информатики факультета информатики Томского государственного университета, поступила в научную редакцию 15 мая 2003 г.