CC BY
50
УДК 004.942:556.507
С.Г. Яковченко, И.С. Постнова, В.А. Жоров ИВЭП СО РАН, Барнаул
АЛГОРИТМЫ ПОСТРОЕНИЯ СОГЛАСОВАННОГО ПРОДОЛЬНОГО ПРОФИЛЯ РУСЛА В ГИС ДЛЯ ПРОВЕДЕНИЯ ГИДРОЛОГИЧЕСКИХ РАСЧЕТОВ
Введение
Продольный профиль русла (ППР), под которым в контексте данной работы понимается совокупность высот водной поверхности (ВП) в местоположении речной сети (РС), необходим для проведения многих гидрологических расчетов. Картографическая информация для прямого определения ППР включает в себя урезы воды и пересечения РС с горизонталями, а также данные по геометрическому положению РС [3, 4]. Необходимым условием для ППР, как части модели рельефа местности является гидрологическая согласованность, выражающаяся в понижении высоты в направлении течения реки и корректном сопряжении ППР в местах слияния рек [3]. Как будет показано далее, задача автоматизированного определения высот ППР может быть решена в стандартных ГИС линейки ESRI путем использовании встроенного в них инструментария картографической алгебры.
Способы построения Ш1Р по рельефу в окрестности РС
Наиболее простым с точки зрения реализации представляется способ расчета ППР по гидрологически корректной цифровой модели рельефа (ЦМР) [3, 4]. В модуль 3D Analyst ArcView заложена принципиальная возможность присвоения высот РС по триангуляционной (TIN) или растровой (GRID) ЦМР с сохранением результата в форме 3D файла. Растровая гидрологически корректная ЦМР может быть создана в модуле TopoGrid ГИС ArcInfo, если включена опция принудительного согласования ЦМР с РС [3, 4, 7, 8]. Необходимо отметить, что РС, представляемая на карте как область, в этом подходе аппроксимируется линейной РС, причем все дуги РС должны быть сомкнуты и иметь корректные направления оцифровки. Высоты ППР, рассчитанные по растровой ЦМР, из-за особенностей алгоритмов интерполяции могут иметь большие ошибки на территории днищ долин, где практически отсутствуют горизонтали.
Использование в качестве источника высот для РС TIN, созданной на основе горизонталей, наталкивается на трудность следующего характера. Для такой TIN вдоль профиля реки характерны резкие перепады высот при переходе от триангуляционных площадок, вершины (nodes) которых являются точками одной горизонтали, к площадкам, у которых вершины принадлежат другой горизонтали [1]. Петлеобразная форма горизонталей вблизи пересечения с рекой еще более усугубляет этот эффект. В [3, 5] предложено простое решение этой проблемы: использовать для присвоения высот точкам на РС TIN, созданный на основе только урезов воды и точек горизонталей, лежащих на пересечении с РС (далее [TlNh]). Такой [TlNh] будет иметь корректное изменение высоты вдоль РС (линейное между отдельными
горизонталями и урезами воды) практически всюду, где есть хорошая обеспеченность горизонталями. Область, в которой [TlNh] можно использовать для присвоения высот РС, определяется простым условием [TlNh]<[TlN] (т.е. высоты, рассчитываемые по [TlNh], должны быть меньше, чем по [TIN], рассчитываемому по горизонталям и отметкам высот). Оба рассмотренных способа требуют наличия информации по рельефу вне РС.
Метод расчета высот ППР по отметкам точек пересечения с горизонталями и урезов воды
В ряде случаев возникает необходимость расчета ППР только по гидрографии и отметкам точек пересечения с горизонталями и урезов воды без использования информации по рельефу вне РС. Такая ситуация возникает, например, при использовании ППР для подготовки гидрологически корректных ЦМР по разномасштабным данным, необходимых для автоматизированного деления поверхности водосбора на расчетные участки и вычисления усредненных морфометрических характеристик этих участков для оценки стока [2, 3]. Основная проблема расчета ППР без использования данных по рельефу вне РС заключается в том, что средние уклоны рек имеют тенденцию уменьшаться с ростом порядка (принимается схема порядков притоков по Хортону-Штралеру), вследствие чего плотность точек пересечения РС с горизонталями для притока, в общем, всегда больше, чем для главного русла (см. рис. 1). При любом интерполяционном процессе такая разная плотность точек на притоках может приводить к ошибкам в расчете высот ВП на реке впадения.
Рис. 1. Пример построения продольного профиля русла (р. Самыш, бассейн р.
Бия, Алтае-Саянская горная система)
Для решения этой проблемы предлагается подход, в основе которого положена идея раздельной интерполяция высот для рек разного порядка:
а) В расчете высот реки-притока используется информация о высоте реки впадения,
б) В расчете высоты реки впадения не используется информация о высотах РС притока, и тем самым при интерполяции учитываются только исходные данные по основному потоку.
Далее рассматривается программная реализация указанного подхода в ГИС ArcView средствами внутреннего языка Avenue.
Исходными данными для построения ППР являются: РС (линия) составляющая единую неразрывную сеть, урезы (точки), пересечения горизонталей с РС (точки), пересечения РС с границей области (точки) с оценочными величинами высот (эти данные необходимы для однозначной интерполяции высот на границах области). Последние три покрытия далее рассматриваются как единое покрытие [Rpoint].
Для расчета порядка водотоков используется промежуточная модель рельефа, вычисляемая как сумма интерполированного по [Rpoint] рельефа (далее называемого "базис эрозии") и добавки, быстро растущей (по сравнению с падением высоты вдоль русла реки) при удалении от РС. В качестве добавки используется умноженное на большую величину поле расстояний от РС. С целью присвоения порядков РС используется грид аллокаций вдоль РС, рассчитываемый из [Rpoint] (функцией CostDistance), значением которого для каждой точки РС является значение высоты ближайшей вдоль РС к данной точке точки [Rpoint]. Предлагается следующая последовательность действий:
1. Расчет грида исходных источников рельефа, расчет гридов аллокаций и расстояния до РС.
2. Расчет на основе грида аллокаций (после преобразования его в точки) грида базиса эрозии с использованием сплайн интерполяции.
3. Расчет на основе гридов базиса эрозии и расстояния до РС промежуточного грида.
4. Заливка промежуточного грида и расчет направлений потоков на нем.
5. Расчет на основе направлений потоков грида порядков РС.
6. Выделение диапазона порядков и построение пополняемого далее грида текущих источников рельефа (лежащих на дренаже, принадлежащем выбранному диапазону порядков), используемых для интерполяции высот РС (вначале пустой).
7. Цикл по порядкам i, начиная от самого высшего:
а) Запросом выделяется область, занимаемая порядками > i;
б) Рассчитывается грид базовых данных по рельефу путем "подтягивания" грида исходных источников рельефа (функция SnapPourPoint) к дренажу на величину размера ячейки;
в) В гриде базовых данных по рельефу выделяется область, занимаемая порядками > i, и присоединяется к гриду текущих источников рельефа;
г) Грид текущих источников рельефа сплайн-интерполируется, его значения сравниваются с гридом аллокаций и, если разность по абсолютной величине превышает шаг горизонталей h, значение ячейки принимается равным гриду аллокаций ±И (знак выбирается как знак разности). По результатам расчетов лучше всего подходит реализованная в ArcView сплайн интерполяция по типу упругой мембраны, когда в критерии минимизации учитывается первая производная с весом w [6];
д) Из результирующего грида вырезается область, занимаемая порядками > i, и сохраняется в грид текущих источников рельефа.
8. Грид текущих источников рельефа (после преобразования его в точки) используется как основа для создания скорректированного грида базиса эрозии.
9. Рассчитывается корректированный грид на основе гридов полученного базиса эрозии и расстояния до РС.
10. Корректированный грид подвергается заливке (Fill). Вычисляется разность между ним и гридом текущих источников рельефа, используемая для оценки корректности исходных данных.
11. Область, занимаемая РС залитого скорректированного грида, конвертируется в покрытие точек, используемое как ППР.
Данный алгоритм был апробирован для территории бассейна оз. Телецкое. Ошибка расчета, оцениваемая как разность между значениями высот в интерполированном DEM и [Rpoint], по порядку величины не превышала уклона реки, умноженного на размер ячейки (см. табл. 1).
Таблица 1. Характеристики расчетных ППР, сгенерированных по данным
разных масштабов для участка вблизи оз. Телецкое
Масшта б Число точек Шаг горизонтал ей, h м Стандартное отклонение высоты в точках пересечений РС с горизонталями, м/ доля нулевых уклонов %
w=0.1 w=5 w=10 w=10 (c расчетом высот вершин дуг)
1:50К 3159 20 0.411.65 0.12/1.28 0.05/1.44 0.04/1.00
1:100К 2444 20 0.72/2.76 025/2.07 0.12/2.08 0.05/1.85
1:200К 965 40 1.66/2.94 0.10/2.4 0.14/2.72 0.15/1.76 (0.11/1.61 при w=5)
Для повышения точности ППР на участках истоков РС величины высоты РС в наивысшей точке отдельной дуги каждого порядка могут рассчитываться из условия сохранения среднего уклона данного потока для
чего предложен ряд дополнительных ГИС операций в цикле между стадиями 7в и 7г:
1. Выделяется область, занимаемая данным порядком РС.
2. Строится последовательность уникальных непрерывных дуг, в которые разбивается область (функция StreamLink).
3. Для каждой дуги методом пространственной статистики рассчитывается максимальное Hup и минимальное Hdown значение высот текущих источников рельефа. Применением функции FlowLength рассчитываются длина каждой дуги L, длина отрезка от наинизшего источника рельефа Hdown до устья дуги, обозначаемая Ldown, длина отрезка до устья дуги от наивысшего источника рельефа Hup, обозначаемая Lup.
4. Для каждой дуги максимальная высота на дуге Hmax оценивается формулой:
нгш, = нtop + min( К(Htop -Hdown)■(/.- Ltop)/(Ltop - Ldown)
5. Для каждой дуги вычисляется наивысшая точка дуги по признаку равенства нулю аккумуляции для неё и ей присваивается значение максимальной высоты Нтах, рассчитанное ранее.
6. Полученные точки добавляются в покрытие текущих источников рельефа.
Расчет вершин дуг указанными операциями существенно понижает количество плоских участков РС (оцениваемое долей нулевых значений уклона, см. табл. 1), делая распределение уклонов более равномерным. Как видно из результатов (см. табл. 1), при увеличении масштаба и уменьшении числа точек оптимальное значение w уменьшается, а применение алгоритма предварительной оценки значений высоты в вершинах дуг не приводит к уменьшению отклонения высот. Таким образом, алгоритм оценки высот на концах дуг крупномасштабных карт имеет смысл использовать только для данных, снятых с карт масштаба более крупного, чем 1:100К. Оптимальное значение параметра w =10 для таких карт и w =5 для средне- и мелкомасштабных.
Выводы
Предложенный алгоритм создания ППР не использует данные по рельефу вне гидрографической сети, что позволяет применять его при создании согласованной информации по рельефу на гидрографии крупного масштаба для уточнения рельефа мелкого масштаба и последующем построении гидрологически согласованных ЦМР в рамках настольных ГИС, без привлечения специальных программ.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Яковченко, С.Г. Об одном способе получения цифровой модели рельефа / С.Г. Яковченко, Н.М. Ковалевская // Материалы междунар. конф. "ГИС для устойчивого развития территорий", июнь 1999 г., Якутск ч. 2. - Якутск, 1999. - С. 72-77.
2. Яковченко, С.Г. Подготовка смешанных ЦМР с использованием данных по рельефу русла и речной сети для расчета морфометрических характеристик водосборов /
С.Г. Яковченко, В.А. Жоров, И.С. Постнова // Сб. материалов науч. конгр., ГЕО-СИБИРЬ-2005, Т. 4, Геоинформатика. - Новосибирск, 2005. - С. 83-88.
3. Яковченко, С.Г. Создание и использование цифровых моделей рельефа в гидрологических и геоморфологических исследованиях / С.Г. Яковченко, В.А. Жоров, И.С. Постнова. - Кемерово: Изд-во ИУУ СО РАН, 2004. - 92 с.
4. Яковченко, С.Г. Технология создания гидрологически корректных моделей рельефа/ С.Г. Яковченко, В.А. Жоров, И.С. Постнова // Материалы междунар. конф. "ГИС для устойчивого развития территорий", 28 мая - 1 июня 2002 г., Санкт-Петербург. - СПб., 2002. - С. 137-142.
5. "SNIPCALC" - геоинформационное приложение для автоматизированного расчета параметров водосборов / С.Г. Яковченко, В.А. Жоров, И.С. Постнова, О.В. Ловцкая, Е.К. Воробьев // Материалы междунар. конф. "ГИС для устойчивого развития территорий", 2003. - С. 216-221.
6. ArcView GIS. Руководство пользователя, ESRI, 1996, с. 377.
7. Arc/Info Help.
8. Hutchinson M.F., A new procedure for gridding elevation and stream line data with automatic removal of spurious pits, Journal of Hydrology, 106 (1989) pp. 211-232.
© С.Г. Яковченко, И.С. Постнова, В.А. Жоров, 2007
