CC BY
124
МЕТОДЫ ГЕОГРАФИЧЕСКИХ ИССЛЕДОВАНИЙ
УДК 004.942; 551.43:004.9 С.М. Кошель1, А.Л. Энтин2
ВЫЧИСЛЕНИЕ ПЛОЩАДИ ВОДОСБОРА ПО ЦИФРОВЫМ МОДЕЛЯМ РЕЛЬЕФА НА ОСНОВЕ ПОСТРОЕНИЯ ЛИНИЙ ТОКА
Гидрологический анализ рельефа на основе цифровых моделей - важная область морфометри-ческого анализа рельефа средствами ГИС. Такие параметры, как площадь водосбора и объем стока, активно применяются в географических исследованиях (напрямую или опосредованно). Их вычисление опирается на процедуру определения направления и пропорции перемещения водных масс по регулярной сетке цифровых моделей рельефа (ЦМР). К настоящему времени разработано достаточно большое количество алгоритмов, реализующих разные подходы к этой процедуре, и их применение приводит к различающимся (иногда на несколько порядков) результатам расчета площади водосбора. Все эти алгоритмы реализуют концепцию дискретного стока «из ячейки в ячейку», что, во-первых, не вполне соответствует реальному перемещению водных масс по повнрхности, и, во-вторых, затрудняет верификацию результатов расчета: ячейки регулярной сетки не всегда можно однозначно соотнести с элементами и формами рельефа.
Предлагается принципиально другой подход к расчету площади водосбора и объема стока, основанный на построении линий тока на непрерывной поверхности, восстановленной методом интерполяции по значениям высот в дискретных узлах регулярной сетки ЦМР.
Ключевые слова: цифровая модель рельефа (ЦМР), геоинформационная система (ГИС), алгоритм, гидрологический анализ рельефа, площадь водосбора, объем стока, построение линий тока.
Введение. Площадь водосбора и объем поверхностного стока - морфометрические показатели, рассчитываемые по цифровым моделям рельефа (ЦМР), - широко используются в геоинформационном анализе местности. Их принято относить к блоку гидрологических параметров, где они являются одними из базовых понятий. На их основе, например, может быть построена сеть потенциальных водотоков, или выполнен расчет различных производных показателей (таких, как компоненты универсального уравнения смыва почв - USLE), необходимых для оценки и моделирования процесса эрозии почв.
В современных программных средствах ГИС вычисление площади водосбора и объема стока опирается на расчет направлений поверхностного стока на основе сеточной ЦМР, трактуемой в качестве матрицы квадратных (реже прямоугольных) ячеек с заданной высотой. Моделируемая водная масса перемещается из ячейки с большей высотой в соседнюю ячейку с меньшей высотой. Направления стока для данной ячейки указывают, куда и каким образом происходит перемещение. Возможны следующие варианты: перемещение потока из ячейки строго в одну соседнюю ячейку или распределение потока между несколькими соседними ячейками в некоторой пропорции. К настоящему времени предложено и реализовано множество алгоритмов расчета направлений стока, реализующих различные
варианты и пропорции распределения потока между соседними ячейками.
Площади водосбора и объемы стока, рассчитанные на базе различных алгоритмов, могут отличаться на несколько порядков. При этом до сих пор не предложено надежного способа оценки достоверности расчета направлений стока и пропорций распределения потока на реальных поверхностях. Соответственно, нет возможности выбрать среди «сеточных» алгоритмов расчета направлений стока такой алгоритм, который дает наиболее достоверный результат.
Возможным способом решения описанной проблемы может являться изменение концепции расчета показателей таким образом, чтобы перемещение потоков не было привязано к регулярной сетке. В настоящей работе представлен алгоритм расчета базовых гидрологических параметров, основанный на построении линий тока на непрерывной поверхности. Сама поверхность при этом восстанавливается путем интерполяции по значениям высот в дискретных узлах регулярной сетки ЦМР.
Материалы и методы исследований. Площадь водосбора (англ. Catchment Area, Contributing Area) - морфометрический показатель, выражающий общую площадь (в плане) территории, с которой осуществляется потенциальный сток в данную ячейку регулярной сетки. Этот показатель отража-
1 Московский государственный университет имени М.В. Ломоносова, географический факультет, кафедра картографии и геоинформатики, вед. науч. сотр., канд. геогр. наук; e-mail: [email protected]
2 Московский государственный университет имени М.В. Ломоносова, географический факультет, кафедра картографии и геоинформатики, аспирант; e-mail: [email protected]
ет только геометрические формы (свойства) поверхности рельефа. Процедуры расчета площади водосбора входят в блоки гидрологического анализа рельефа практически всех программных средств ГИС. Значение показателя вычисляется для каждой ячейки сетки цифровой модели рельефа. В итоге получается матрица значений, аналогичная матрице высот исходной ЦМР. Используя рекурсивный подход, идею расчета этого показателя можно выразить так: площадь водосбора ячейки (ij) равна сумме площадей водосбора соседних ячеек, из которых осуществляется сток в (i j), площади всех ячеек при этом считаются равными. Площадь самой ячейки (i,j) может включаться или не включаться в рассчитываемую площадь водосбора - в модулях анализа рельефа современных ГИС реализованы оба варианта. Для ячеек, находящихся на водораздельных позициях или локальных максимумах, площадь водосбора равна площади ячейки (или нулю), по мере понижения высот значения этой величины быстро растут и в тальвегах достигают уже порядка 104-105 площади одной ячейки ЦМР. Максимальное значение в пределах модели зависит от особенностей строения речной сети и может превышать значения в вершинах тальвегов еще на несколько порядков.
Объем стока (англ. Flow Accumulation) рассчитывается по аналогичной схеме. Разница заключается в том, что исходное значение рассчитываемой величины в каждой ячейке представляет не ее площадь, а слой стока (количество осадков, участвующее в формировании поверхностного стока). Это исходное значение может быть различным для разных ячеек и подается на вход алгоритмам расчета показателя в виде матрицы той же размерности, что и матрица высот.
На практике «истинные», то есть соответствующие определению, значения показателя получаются только для ячеек внутри водосборных бассейнов, целиком входящих в экстент ЦМР. Для остальных ячеек значения получаются заниженными, из-за отсутствия в ЦМР части водосборного бассейна, территория которого могла бы внести свой вклад в значение рассчитываемого показателя.
Практически все разработанные к настоящему времени алгоритмы используют для расчета площади водосбора (или объема стока) приведенную выше схему. Различия заключаются в том, какие ячейки считаются соседними и какая процедура расчета направлений стока используется для распределения воды из текущей ячейки в соседние.
Существующие алгоритмы расчета направлений стока. Алгоритм Deterministic Eight-Neighbor (D8) был впервые применен в ГИС-ана-лизе в середине 80-х годов [O'Callaghan, Mark, 1984]. Поток из рассматриваемой ячейки целиком направляется в ту из восьми соседних, которая имеет, во-первых, меньшую высоту и, во-вторых, наибольший уклон линии, соединяющей центр текущей ячейки с центром соседней. Такой подход (с переносом потока только в одну из восьми соседних ячеек) позволяет построить сеть ячеек, связанных между
собой потоками воды (в работе [O'Callaghan, Mark, 1984] эта сеть называется drainage graph; в более поздних работах встречаются другие варианты, например flow network), и далее выделить сеть потенциальных водотоков. Алгоритм D8 наиболее прост в реализации, и в силу этого присутствует во всех программных средствах ГИС, имеющих модули гидрологического анализа рельефа. Однако он имеет ряд недостатков. Один из них связан с неопределенностью, возникающей в ситуациях, когда максимальный уклон соответствует более чем одной соседней ячейке. Существует модификация алгоритма D8, позволяющая кодировать такие ситуации и затем выбирать единственное направление на основе специальной перекодировочной таблицы [Greenlee, 1987]. Другим недостатком является распределение потоков строго под углами, кратными 45°, что приводит к формированию нереалистичной сети потоков, особенно на достаточно ровных наклонных поверхностях.
Алгоритм Random Eight-Neighbor (p8, Rho8) был создан как попытка решить описанную выше проблему нереалистичной сети потока, которая возникает при использовании алгоритма D8 [Fairfield, Leymarie, 1991]. По замыслу разработчиков, если реальное направление стока поверхности имеет азимут, например, 30°, то можно распределить поток таким образом, чтобы сток части ячеек был направлен на северо-запад (азимут 45°), а часть - на север (азимут 0°). Алгоритм p8 вычисляет уклоны между соседними ячейками таким же образом, как и D8, но затем величины этих уклонов корректируются - умножаются на некоторую случайную величину, разную для различных ячеек. Наибольшие значения скорректированных уклонов в большинстве случаев будут соответствовать направлению наибольшего уклона, как и в D8; в тех же случаях, когда результат будет отличаться, и происходит устранение нежелательного эффекта «дискретности». Однако, поскольку в расчете участвуют случайные величины, обработка одной и той же ЦМР по алгоритму p8 может приводить к разным результатам в зависимости от используемого датчика случайных чисел. В настоящее время этот алгоритм реализован только в SAGA GIS.
Алгоритмы D8 и p8 направляют поток из ячейки строго в одну из восьми соседних ячеек. Более сложные алгоритмы могут распределять поток между несколькими соседними ячейками. Deterministic Infinity (D<x>) также рассчитывает для ячейки единственное значение азимута направления стока, но, в отличие от предыдущих вариантов, это направление не ограничено условием кратности 45° [Tarboton, 1997]. Для расчета азимута направления стока рассматриваются не пары «центральная ячейка - соседняя ячейка», а треугольники, образованные центральной ячейкой и двумя ее соседями (которые, в свою очередь, являются соседями друг друга). Всего рассматривается восемь таких треугольников; каждому треугольнику соответствует единственная наклонная плоскость, Для каждой из этих восьми
плоскостей вычисляется вектор, противоположный градиенту, с началом в середине центральной ячейки. В качестве направления стока выбирается направление того из построенных векторов, который имеет наибольшую длину и при этом находится в пределах угла, образованного сторонами соответствующего треугольника.
Метод Ли (Lea's method; известен также под названиями Kinematic Routing Algorithm или способ «катящегося шарика») представляет поверхность ячейки как наклонную плоскость, положение которой определяется с учетом высот соседних ячеек [Lea, 1992]. Этой плоскости соответствует определенное направление вектора градиента. Азимут противоположного направления выбирается в качестве направления стока в данной ячейке. В алгоритме предполагается, что поток начинается из центра произвольной ячейки, далее течет по азимуту, рассчитанному для этой ячейки, до ее границы. Переходя границу, поток попадает в новую ячейку и поворачивает так, чтобы соответствовать новому азимуту. Движение элементарной водной массы, таким образом, уподобляется шарику, катящемуся по поверхности ЦМР. Суммирование всех потоков, проходящих через некоторую ячейку, позволяет вычислить площадь водосбора или объем стока для этой ячейки.
На основе метода Ли (KRA) был разработан алгоритм DEMON [Costa-Cabral, Burges, 1994]. Направления в нем рассчитываются так же, как в алгоритме KRA, разница же заключается в интерпретации потока. В алгоритме KRA водная масса считается сконцентрированной в центре ячейки и перемещающейся далее как единое целое («катящийся шарик»), DEMON использует введенное авторами понятие «полоса потока» (flow tube). Эта полоса ограничена двумя линиями тока, и характеристики рассчитываются на основе ее ширины.
Группа алгоритмов, распределяющих поток из ячейки между всеми соседними нижележащими ячейками, носит название Multiple Flow Direction (MFD). Различия между этими алгоритмами заключаются в способах вычисления пропорций распределения потока между соседними ячейками. В наиболее ранних вариантах такого подхода [Freeman, 1991; Quinn et al., 1991] доля потока, приходящегося на данного соседа, определяется формулой, содержащей эмпирические константы. В более поздних работах [Seibert, McGlynn, 2007; Qin et al., 2007] предлагаются специальные процедуры для расчета пропорции, в которых может учитываться даже кривизна поверхности. В целом алгоритмы группы MFD позволяют распределять потоки более достоверно, жертвуя при этом возможностью строить связную сеть ячеек.
Более подробно, с графическими иллюстрациями и указанием формул, изложенные выше алгоритмы рассматриваются в обзорной статье [Кошель, Энтин, 2016]. Там же можно ознакомиться с методами, применяемыми для предварительной обработки ЦМР с целью устранения локальных понижений
и других возможных препятствий корректной работе алгоритмов гидрологического анализа рельефа.
Помимо разных алгоритмов расчета направлений стока, применяются также разные подходы при вычислении собственно площади водосбора или объема стока. Эти подходы отличаются порядком просмотра и обработки ячеек модели. При просмотре «сверху вниз» (Top-Down Processing) ячейки сортируются по убыванию высоты и обрабатываются по очереди от самой высокой к самой низкой; взаимное расположение ячеек не учитывается. Такой способ может комбинироваться с любыми алгоритмами расчета направлений стока, кроме KRA и алгоритма DEMON. Другой вариант - рекурсивная обработка ячеек «от устья к истоку» [Olaya, 2004; Jasiewicz, Metz, 2011]: ячейки просматриваются в порядке возрастания высоты, для каждой нижележащей ячейки в просмотр включаются все вышележащие, вплоть до водоразделов или максимумов (то есть до ячеек, для которых не существует соседей с большей высотой). Присвоение значений производится рекурсивно в процессе просмотра. Для KRA и алгоритма DEMON применяется специальный способ «трассировки потока» (Flow Tracing). Обработка начинается с произвольной ячейки, от нее строится линия тока (для алгоритма Ли) или ложе потока (DEMON) вниз по склону до границы ЦМР. В процессе построения производится изменение значений всех ячеек, по которым проходит линия тока. Одна и та же ячейка, если через нее проходят несколько линий тока, обрабатывается несколько раз. Построение линий тока и вычисление площади водосбора осуществляется параллельно.
Различные реализации алгоритмов расчета площади водосбора (разный порядок просмотра ячеек), применяемые с одними и теми же алгоритмами расчета направлений стока генерируют одинаковые результаты. Однако распределения показателей, рассчитанных по одной и той же модели, но с применением разных алгоритмов расчета направлений стока, могут различаться очень значительно. Вопрос о большей или меньшей достоверности получаемой с помощью различных алгоритмов картины распределения показателя в настоящее время решается следующим образом: показатель рассчитывается по модели, представляющей некоторую простую поверхность (наклонная плоскость, эллипсоид, конус), для которой возможно теоретически предсказать распределение его значений. Чем ближе рассчитанное распределение оказалось к ожидаемому, тем более достоверным считается алгоритм [Qin et al., 2007; Seibert, McGlynn, 2007]. Этот подход позволяет оценить достоверность через количественные характеристики (например, средне-квадратическое отклонение), однако не вполне ясно, можно ли переносить выводы, сделанные для относительно простых модельных поверхностей, на сложные формы реального рельефа. Иными словами, низкие отклонения результатов на модельных поверхностях еще не являются достаточным доказа-
тельством большей достоверности того или иного алгоритма.
Построение линий тока. Перечисленные проблемы могут отчасти быть решены, если отказаться от концепции дискретного стока «из ячейки в ячейку» и моделировать распространение потоков на основании того факта, что водные массы под действием гравитации перемещаются по непрерывной поверхности в направлении, противоположном направлению вектора градиента в каждой точке. Выражаясь формальным языком, имея модель поверхности рельефа в виде непрерывной дифференцируемой функции Дх,у), линию тока в параметрическом представлении (ху(0), выходящую из точки х у0 вниз по склону, можно определить как решение задачи Коши для системы дифференциальных уравнений первого порядка при значениях параметра ^0 [Кошель, 2004]:
х'(0 =-~Т (х, у)
^ ; х(0) = Х0, у(0) = У0. (1)
x'(t) = ~ (х, у) Ох
С практической точки зрения, чтобы избежать слишком маленьких расстояний между последовательными точками решения, правую часть системы желательно пронормировать, то есть разделить на длину вектора градиента. Нормировка выполняется в том случае, если длина градиента не меньше некоторого заданного порога, иначе правой части присваиваются нулевые значения. В итоге правая часть системы (1) будет представлять собой компоненты единичного вектора, противоположного вектору градиента. Конечный результат численного интегрирования этой системы определяется двумя факторами: способом вычисления правой части в произвольной точке исследуемой области (при том, что нам известны значения высот только в узлах регулярной сетки) и выбранным методом численного интегрирования. Восстановление непрерывных частных производных на основе высот в узлах сеточной модели рельефа можно выполнять разными способами, мы использовали билинейную интерполяцию по значениям в четырех ближайших узлах сетки. При выборе способа интегрирования были опробованы разные алгоритмы, начиная от простейшего метода Эйлера первого порядка и заканчивая методом Рунге-Кутты четвертого-пятого порядка [Каханер, Моулер, Нэш, 1998; Форсайт, Малькольм,
Моулер, 1980]. Иллюстрация работы метода Эйлера с постоянным шагом представлена на рис. 1. В первую очередь, задается шаг интегрирования (шаг, с которым увеличивается параметр 0. Для построения плавной линии тока желательно, чтобы величина этого шага была меньше половины размера ячейки ЦМР. Далее, в начальной точке вычисляются компоненты вектора градиента (частные производные модельной функции) и в противоположном направлении откладывается отрезок длины, равной шагу интегрирования. Конец этого отрезка принимается за новую начальную точку. Процесс продолжается до тех пор, пока ломаная линия не достигнет границы ЦМР, либо когда значение правой части не станет равным нулю (попадание в точку локального минимума или на горизонтальную площадку). К сожалению, этот простейший подход оказался непригодным для наших целей, поскольку зачастую у построенных таким образом ломаных наблюдаются осцилляции и даже зацикливание в районе некоторых тальвегов и на участках с близкими к нулю значениями модуля градиента. В результате экспериментов с использованием ЦМР на участки с различными типами рельефа для решения системы (1) в итоге был выбран неявный метод Эйлера. Он оказался оптимальным с точки зрения соотношения получаемого качества результата и затрат на вычисления.
Расчет площади водосбора на базе построенных линий тока. Для расчета площади водосбора и объема стока выполняется построение линий тока из центра каждой ячейки ЦМР. С линией тока ассоциируется значение рассчитываемой величины (площади водосбора или объема стока), соответствующее той ячейке, из которой она строится.
Когда линия построена, анализируется ее положение относительно ячеек регулярной сетки. Предлагается два варианта такого анализа, назовем их условно «простой» и «точный». В «простом» способе используется только положение вершин построенной ломаной; к значениям рассчитываемой величины в тех ячейках, куда попала хотя бы одна вершина, добавляется значение, ассоциированное с данной линией. Если в ячейку не попала ни одна вершина, значение в ней не увеличивается, даже если часть линии фактически проходит по этой ячейке. Такой подход значительно снижает вычислительные затраты, и при этом вполне правомерен, поскольку шаг, с которым строится ломаная, представляющая линию тока, меньше, иногда значительно, размера
! a df\ \ \ б I \ \ в \ | \ г
I Эх | || | it | || |
j£?.......: : : V:. : : X: :: X;
Г""7-----i i----------i : \ - ; \
I I II I II I II I
I I II I II I II I
I______I______I I______I______I Г______I______I I______|_
Рис. 1. Построение линии тока из центра ячейки. Объяснение а - г - в тексте Fig. 1. Flowline tracing from the center of a cell. (For a - г see the text)
ячейки сетки. «Точный» вариант основан на подсчете длин сегментов линии тока. Для каждой ячейки, через которую проходит линия, определяется длина той части ломаной, которая находится в пределах ячейки. Значение рассчитываемой величины, ассоциированное с линией тока, затем умножается на коэффициент Касс, зависящий от длины сегмента ломаной и добавляется к значению, уже имеющемуся в данной ячейке. Коэффициент К вычисляется по формуле (2):
(
K„„ = min
1,
L,„
\
L
(2)
'0 у
Здесь L - длина сегмента линии в пределах рассматриваемой ячейки, L0 - длина условного единичного сегмента линии. Значение L0 устанавливается из следующих соображений: коэффициент Kacc должен быть равен единице в том случае, если линии тока параллельны одной из осей координат. В таком случае длина сегмента линии в пределах ячейки равна шагу сетки вдоль данной оси. Следовательно, для квадратной сетки значение L0 принимается равным шагу сетки. Для прямоугольных сеток с различными шагами вдоль разных осей предлагается устанавливать в качестве значения L0 полусумму шагов сетки. Коэффициент Kacc ограничивается сверху единицей, чтобы вносимая прибавка не превосходила значение, ассоциированное с линией тока. Процедура анализа положения линии тока относительно ячеек сетки ЦМР проиллюстрирована на рис. 2.
Алгоритм, названный нами FLBA (FlowLine-Based Algorithm), реализован как самостоятельная программа на языке Fortran (стандарт Fortran 90). Исходными данными для программы являются сеточная цифровая модель рельефа и параметры ме-
Рис. 2. Расчет гидрологического показателя на базе построения линии тока - «простой» способ (вверху) и «точный» способ (внизу). Более темный цвет соответствует большему значению показателя
Fig. 2. Computation of a hydrological parameter using the flowline tracing: «simple» way (above) and «accurate» way (below). Darker color stands for higher parameter value
тода, на выходе программа создает цифровую модель площади водосбора, со значениями показателя, рассчитанными в узлах регулярной сетки той же размерности, что и исходная.
Результаты исследований и их обсуждение. Тестирование разработанного алгоритма проводилось на ЦМР, представляющей несколько небольших эрозионных форм рельефа в северной части гор Западного Саяна. Абсолютная высота в пределах рассматриваемого участка изменяется от 700 до 900 м. Цифровая модель рельефа построена по данным воздушного лазерного сканирования путем интерполяции методом естественного соседа (Natural Neighbor) [Sibson, 1981] с получением значений высот на регулярной сетке с шагом 1 м. Модель была обработана при помощи алгоритма заполнения локальных понижений Ванга-Лю [Wang, Liu, 2006]; замкнутые локальные понижения были устранены, на их месте сформированы наклонные поверхности с углом наклона 0,1°. Другие способы предварительной обработки (сглаживающие фильтры и т. п.) не применялись.
Работа алгоритма оценивалась визуально на примере результата вычисления площади водосбора. Для сопоставления с другими методами этот же параметр был рассчитан в пакете SAGA GIS с использованием алгоритмов D8, MFD (наиболее распространенные и часто используемые), MFD-md, KRA (идеологически близкий к разработанному FLBA).
Сравнение результатов, полученных двумя вариантами алгоритма FLBA («простым» и «точным») показало, что существенные отличия между ними наблюдаются только вблизи тальвегов, где ячейки задеваются «по касательной» большим количеством линий тока (рис. 3, д, е). При простом варианте расчета одни ячейки «недополучают» площадь водосбора, поскольку ни одна вершина ломаной не попадает в них, а другие, наоборот, «получают больше, чем должны», если общая протяженность линии тока в них невелика. На участках, удаленных от тальвегов, существенных отличий не наблюдается.
Сравнение результатов, полученных разными алгоритмами, показывает, что картина распределения показателя, формируемая FLBA (с любыми параметрами), отличается довольно низким «разбросом» потока (рис. 3). Области высоких значений площади водосбора вблизи тальвегов формируют линии, ширина которых равна размеру ячейки ЦМР. Эти области, таким образом, гораздо более узкие, чем при использовании алгоритмов MFD и MFD-md (рис. 3, б, в) и сопоставимы с результатом расчета по D8 (рис. 1, а). Однако предлагаемый алгоритм не демонстрирует характерного для D8 неестественного расположения потоков под углами, кратными 45°. Похожая картина наблюдается также при использовании идейно близкого алгоритма KRA (рис. 3, г), но из-за особенностей построения линий тока в нем часто возникает ситуация, когда на линии тальвега оказывается несколько ячеек с низкими значениями площади водосбора (образно
Рис. 3. Результаты расчета площади водосбора с использованием разных алгоритмов определения направлений тока: D8 (а), MFD (б), MFD-md (в), KRA (г), предложенного алгоритма (FLBA) - «простой» расчет (д), «точный» расчет (е)
Fig. 3. Catchment area computation results using different algorithms for identification of flowline direction: D8 (а), MFD (б), MFD-md (в), KRA (г), suggested FLBA - «simple» (д), FLBA, «accurate» (е)
эта ситуация похожа на остров посреди реки). Результаты расчета по разработанному алгоритму лишены этого недостатка, что позволяет их использовать в качестве основы для выделения сети потенциальных водотоков (сейчас для этих целей используется исключительно алгоритм D8).
На рис. 4 приведены поперечные разрезы через характерные участки рельефа, позволяющие более детально и наглядно проиллюстрировать особенности работы перечисленных выше алгоритмов. Профиль I (линия 1а-1б) пересекает склон водосборного понижения параллельно горизонталям, профиль II (Па-Пб) пересекает тальвег наиболее крупной эрозионной формы в верхнем течении, профиль III (Ша-Шб) пересекает тот же тальвег в нижнем течении. Плановое положение линий профилей показано на левой части рис. 4; сами профили приведены на правой части рис. 4. Видно, что для профиля I, ближайшего к водоразделу, для результатов расчета всеми алгоритмами характерны невысокие значения площади водосбора (порядка 100 м2). При этом алгоритмы MFD и MFD-md распределяют сток практически равномерно по склону, KRA и обе разновидности FLBA, на фоне общей равномерности, концентрируют сток вдоль несколь-
ких слабо выраженных эрозионных борозд. Распределение же, получаемое по алгоритму D8, отличается высокой неравномерностью и связано с его ограничениями, описанными выше, - именно так на профиле выглядит формируемая этим алгоритмом сеть параллельных водотоков. На профиле II можно наглядно видеть ширину полосы относительно высоких значений водосборной площади вблизи тальвега. Для алгоритмов группы MFD ширина этой полосы составляет около 3-5 м (напомним, что размер ячейки ЦМР равен 1 м), алгоритмы KRA, FLBA-s, FLBA-a концентрируют сток в единственной ячейке (самой низкой в поперечном профиле). Результат, показываемый алгоритмом D8 (три отдельных невысоких пика водосборной площади вблизи тальвега), вновь связан с особенностями определения направлений стока, но не с рельефом местности. Наконец, профиль III демонстрирует дивергенцию потока при переходе от ложбины к конусу выноса. Для алгоритмов группы MFD «полоса потока» становится еще шире (7-10 м). Алгоритм KRA рассчитывает два максимума вдоль линии профиля (та самая ситуация «острова посреди реки»); больший из этих максимумов соответствует «настоящему» тальвегу, меньший появляется
1а 16 Па Пб Ша 1116
О 10 20 30 40 50 0 10 20 30 40 SO 0 10 20 30 40 50
Расстояние, м
Рис. 4. Поперечные профили характерных элементов и форм рельефа: плановое положение (слева), значения площади водосбора, рассчитанной различными алгоритмами (D8, MFD, MFD-md, KRA, FLBA-s, FLBA-a) вдоль линий профилей I, П, Ш (справа)
Fig. 4. Profile lines across typical elements and landforms: location of profile lines (left), catchment area values along I, II and III profile lines (right). Catchment areas were computed using D8, MFD, MFD-md, KRA, FLBA-s, and FLBA-a algorithms
вследствие особенностей расчета в алгоритме КЯА и соответствует «псевдотальвегу», появление которого связано скорее с особенностями реализации алгоритма, чем с рельефом местности. Отметим также, что положение меньшего из отмеченных локальных максимумов полностью совпадает с положением максимума, рассчитанного алгоритмом D8. Алгоритм FLBA по-прежнему, как и на профиле II выше по склону, концентрирует сток в одной-двух ячейках. Эту особенность FLBA следует учитывать при выборе разрешения ЦМР для анализа - если размер ячейки будет заметно меньше ширины днища ложбины, поток будет концентрироваться в наиболее низкой ее части, что приведет к недооценке водосборной площади в пределах относительно высоких участков ложбины.
Таким образом, предложенный алгоритм сочетает достоинства разработанных к настоящему времени методов, и при этом лишен их недостатков. К отрицательным моментам в работе алгоритма можно отнести, пожалуй, только его более низкую вычислительную эффективность (работает медленнее остальных). Однако это нельзя назвать существенным недостатком, особенно учитывая возможный потенциал для его дальнейшего развития. В частности, предложенный подход позволяет легко интегрировать в алгоритм расчета площади водосбора возможность прохождения линий тока через локальные понижения, что позволит использовать его и для моделей, не подготовленных предварительно с помощью специальных процедур. Для этого достаточно внести изменения только в способ вычисления градиента в уравнении (1), все остальные этапы алгоритма остаются неиз-
менными. Помимо того, промежуточная стадия алгоритма - построение линий тока - представляет самостоятельный интерес, особенно при реализации в интерактивном режиме.
Выводы:
- представлен алгоритм расчета площади водосбора и объема стока, базирующийся на построении линий тока, исходящих из узлов регулярной прямоугольной сетки в пределах исходной ЦМР, и последующем анализе прохождения построенных линий через ячейки результирующей регулярной сетки. Важной отличительной особенностью алгоритма является отказ от дискретного подхода к распределению стока из «ячейки в ячейку» ЦМР. Вместо этого используется непрерывная модельная функция рельефа, интерполирующая данные в узлах регулярной сетки исходной ЦМР, на основе которой строятся линии тока как решение некоторой системы дифференциальных уравнений;
- реализованный алгоритм формирует компактный поток, подобный результату расчета по алгоритму D8. Это позволяет, при достаточном качестве ЦМР, учитывать влияние форм рельефа, размер которых сопоставим с шагом сетки модели, а также использовать результаты расчета в качестве исходных данных для построения сети потенциальных водотоков. Однако, в отличие от D8, разработанный алгоритм не создает нереалистичные сети потоков, текущих под углами, кратными 45°. Это сближает его с подходами группы MFD. Таким образом, предложенный алгоритм сочетает достоинства разработанных к настоящему времени методов расчета площади водосбора, и при этом лишен их недостатков.
Благодарности. Исследование выполнено при поддержке гранта Президента Российской Федерации для молодых российских ученых - кандидатов наук МК-4829.2016.5.
СПИСОК ЛИТЕРАТУРЫ
Каханер Д., Моулер К., Нэш С. Численные методы и математическое обеспечение: Пер. с англ. М.: Мир, 1998. 575 с.
Кошель С.М. Теоретическое обоснование структуры и функций блока моделирования рельефа в ГИС. Дис. ... канд. геогр. н. М., 2004. 105 с.
Кошель С.М., Энтин А.Л. Современные методы расчета распределения поверхностного стока по цифровым моделям рельефа // Геоморфология: Современные методы и технологии цифрового моделирования рельефа в науках о Земле. М.: Медиа-ПРЕСС, 2016. Вып. 6. С. 24-34.
Форсайт Д., Малькольм М., Моулер К. Математические методы машинных вычислений. Пер. с англ. М.: Мир, 1980. 280 с.
Costa-Cabral M.C., Burges S.J. Digital elevation model networks (DEMON): A model of flow over hillslopes for computation of contributing and dispersal areas // Water resources research. 1994. V. 30(6). P. 1681-1692.
Fairfield J., Leymarie P. Drainage networks from grid digital elevation model // Water Resources Research. 1991. V 27(5). P. 709717.
Freeman Т. Calculating catchment area with divergent flow based on a regular grid // Computers and Geosciences. 1991. V. 17. P. 413-422.
Greenlee D.D. Raster and vector processing for scanned linework // Photogrammetric Engineering and Remote Sensing. 1987. V. 53. P. 1383-1387.
Jasiewicz J., MetzM. A new GRASS GIS toolkit for Hortonian analysis of drainage networks // Computers and Geosciences. 2011. V. 37(8). P. 1162-1173.
Lea N.L. An aspect-driven kinematic routing algorithm // Overland Flow: Hydraulics and Erosion Mechanics / Eds.: A.J. Parsons, A.D. Abrahams. New York, 1992.
O 'Callaghan J.F., Mark D.M. The extraction of drainage networks from digital elevation data // Computer vision, graphics, and image processing. 1984. V. 28(3). P. 323-344.
Olaya V. HidrologHa Computacional y Modelos Digitales de Terreno: Teoría, práctica y filosofea de una nueva forma de análisis hidrolygico [Электронный ресурс] // Режим доступа: http:// www.gabrielortiz.com/descargas/Hidrologia_Computacional_ MDT_SIG.pdf. Дата обращения: 24.09.2015.
Qin C., Zhu A.-X., Pei T., Li B., Zhou C., YangL. An adaptive approach to selecting a flow partition exponent for a multiple flow direction algorithm // International Journal of Geographical Information Science. 2007. V. 21(4). P. 443-458.
Quinn P., Beven K., Chevallier P., Planchón O. The prediction of hillslope flow paths for distributed hydrological modelling // Hydrological Processes. 1991. V. 5(5). P. 59-79.
Seibert J., McGlynn B.L. A new triangular multiple flow direction algorithm for computing upslope areas from gridded digital elevation models // Water Resources Research. 2007. V. 43(4).
Sibson R. A brief description of natural neighbour interpolation // Interpreting multivariate data. 1981. V. 21. P. 21-36.
Tarboton D.G. A new method for the determination of flow directions and upslope areas in grid digital elevation models // Water Resources Research. 1997. V. 33(2). P. 309-319.
Wang L., Liu H. An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling // International Journal of Geographical Information Science. 2006. V. 20(2). P. 193-213.
Поступила в редакцию 04.07.2016 Принята к публикации 09.02.2017
S.M. Koshel s A.L. Entin2
CATCHMENT AREA DERIVATION FROM GRIDDED DIGITAL ELEVATION MODELS USING THE FLOWLINE-TRACING APPROACH
Hydrological terrain analysis on the basis of the digital elevation models is an important part of GIS-based geomorphometry. Hydrological parameters, such as the catchment area or flow accumulation, are frequently used in geographical studies, both directly and indirectly. Their computation is based on the determining of water flow through the regular grid of digital elevation model (DEM). There are many algorithms developed for this procedure, and their application gives extremely different computation results for the catchment area. All these algorithms apply a discrete, «cell-to-cell» flow conception. This, first, does not adequately represent real water movement, and, second, makes it difficult to verify results, because cells of a regular grid do not usually «suit» to the landforms.
The paper presents a new approach for catchment area and flow accumulation computation. It is based on flowline tracing on a continuous surface. The surface could be obtained through interpolation of discrete elevation values from a regular grid of DEM.
Keywords: digital elevation model (DEM), geographical information system (GIS), algorithm, hydrological terrain analysis, catchment area, flow accumulation, flowline tracing.
Acknowledgements. The study was financially supported by the President of Russia grant for young PhD scientists (MK-4829.2016.5).
1 Lomonosov Moscow State University, Faculty of Geography, Department of Cartography and Geoinformatics, Leading Scientific Researcher, PhD. in Geography; e-mail: [email protected]
2 Lomonosov Moscow State University, Faculty of Geography, Department of Cartography and Geoinformatics, postgraduate student; e-mail: [email protected]
REFERENCES
Costa-Cabral M.C., Burges S.J. Digital elevation model networks (DEMON): A model of flow over hillslopes for computation of contributing and dispersal areas // Water resources research. 1994. Vol. 30(6). P. 1681-1692.
Fairfield J., Leymarie P. Drainage networks from grid digital elevation model // Water Resources Research. 1991. Vol. 27(5), P. 709-717.
Forsythe G.E., Malcolm M.A., Moler C.B. Matematicheskiye metody mashinnykh vychisleniy (per. s angl.) [Computer Methods for Mathematical Computations (translated from English)]. Moscow: Mir, 1980. 280 p. (in Russian).
Freeman T. Calculating catchment area with divergent flow based on a regular grid // Computers and Geosciences. 1991. Vol. 17. P. 413-422.
Greenlee D.D. Raster and vector processing for scanned linework // Photogrammetric Engineering and Remote Sensing. 1987. Vol. 53. P. 1383-1387.
Jasiewicz J., MetzM. A new GRASS GIS toolkit for Hortonian analysis of drainage networks // Computers and Geosciences. 2011. Vol. 37(8). P. 1162-1173.
Kahaner D., Moler C., Nash S. Chislennye metody i matematicheskoe obespecheniye (per. s angl.) [Numerical Methods and Software (translated from English)]. Moscow: Mir. 1998. 575 p. (in Russian).
Koshel S.M. Teoreticheskoye obosnovaniye struktury i funktsij bloka modelirovaniya reljefa v GIS [Theoretical basis of structure and functions of relief modeling block in GIS applications]. PhD Thesis. Moscow, 2004. 105 p. (in Russian).
Koshel S.M., Entin A.L. Sovremennye metody raschyota raspredeleniya poverhnostnogo stoka po tsifrovym modelyam reljefa [Contemporary methods of DEM-based surface flow distribution modelling]. Geomorfologiya: sovremennye metody i
tehnologii tsifrovogo modelirovaniya reljefa v naukah o Zemle. Issue 6. Moscow: Media-PRESS, 2016. P. 24-34 (in Russian).
Lea N.L. An aspect-driven kinematic routing algorithm // Overland Flow: Hydraulics and Erosion Mechanics / Eds. A.J. Parsons, A.D. Abrahams. New York, 1992.
O 'Callaghan J.F., Mark D.M. The extraction of drainage networks from digital elevation data // Computer vision, graphics and image processing. 1984. Vol. 28(3). P. 323-344.
Olaya V. HidrologHa Computacional y Modelos Digitales de Terreno: Teorna, pr6ctica y filosofHa de una nueva forma de an6lisis hidrolygico [Electronic resource] // URL: http://www.gabrielortiz. com/descargas/Hidrologia_Computacional_ MDT_SIG.pdf. (Accessed 24.09.2015).
Qin C., Zhu A.-X., Pei T., Li B., Zhou C., YangL. An adaptive approach to selecting a flow partition exponent for a multiple flow direction algorithm // International J. Geographical Information Science. 2007. Vol. 21(4). P. 443-458.
Quinn P., Beven K., ChevallierP., Planchon O. The prediction of hillslope flow paths for distributed hydrological modelling // Hydrological Processes. 1991. Vol. 5(5). P. 59-79.
Seibert J., McGlynn B.L. A new triangular multiple flow direction algorithm for computing upslope areas from gridded digital elevation models // Water Resources Research. 2007. Vol. 43(4).
Sibson R. A brief description of natural neighbour interpolation // Interpreting multivariate data. 1981. Vol. 21. P. 21-36.
Tarboton D.G. A new method for the determination of flow directions and upslope areas in grid digital elevation models // Water Resources Research. 1997. Vol. 33(2). P. 309-319.
Wang L., Liu H. An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling // International Journal of Geographical Information Science. 2006. Vol. 20(2). P. 193-213.
Received 04.07.2016 Accepted 09.02.2017
