Основные геоморфометрические параметры: теория

Эта страница является черновиком статьи.

Расчет (Zevenbergen-Thorne, 1987) и интерпретация уклона, экспозиции, кривизн рельефа земной поверхности

Геморфометрический анализ растровых ЦМР базируется на двух исходных положениях. Первое основывается на математической формализации земной поверхности, а второе предусматривает расчет показателя в точке (пикселе) с учетом окружения.

Согласно первому положению, с математической точки зрения ЦМР является статистической поверхностью, которая характеризует пространственное распределение показателя высоты и может быть представлена функцией вида:

$z=f(x,y)$

(1)

где $z$ – значение высоты в точке с географическими координатами $(x,y)$, которое для лучшей аппроксимации рельефа может быть выражено более сложными функциями, например – полиномиальными (или многочленами). В таком случае многочлен 2-го порядка, используемый для аппроксимации земной поверхности, может иметь следующий вид:

$z=Ax^{2}+By^{2}+Cxy+Dx+Ey+F$

(2)

где $x$ и $y$ географические координаты точки, высоту $z$ которой необходимо определить, $A...F$ – коэффициенты уравнения аппроксимирующей поверхности 2-го порядка. Многочлены являются одними из наиболее простых и хорошо изученных функций в математике. Они характеризуются такими свойствами как непрерывность и сглаженность, благодаря чему их легко можно интегрировать и дифференцировать. Это открывает возможности использования математического анализа не только для более совершенного представления земной поверхности, но и для изучения ее свойств, например, на основе производных разных порядков.

Согласно второму положению, основной аналитической операцией в ГИС, которая используется для расчета большинства параметров на основе растровых ЦМР является анализ окружения. Он позволяет количественно описать связь между точкой (пикселем) и его ближайшим окружением, применяя для расчета локальное (чаще всего, размером 3×3 пиксела) скользящее окно (рис. 1).

рисунок 1 Рис. 1 Расчет большинства геоморфометрических параметров как правило производится на основе скользящего окна размером 3×3 пиксела (Geomorphometry…, 2008)

Окно двигается через все поверхность растра (в направлении от верхнего левого до нижнего правого угла) и последовательно применяет в каждой позиции одну и ту же математическую операцию (расчетную формулу) для ячеек основного растра. Таким образом, результат расчетов определяется формулой, которая используется для сравнения значений центральной ячейки с соседними. В результате получается новый растр, аналогичный по пространственному охвату исходной ЦМР, но с другим параметром.

В данной статье мы будем рассматривать особенности расчета основных геморфометрических параметров на примере алгоритма Zevenbergen-Thorne (Zevenbergen, Thorne 1987), который характеризуется расчетной эффективностью и высокой достоверностью результатов (Skidmore 1989, Jones 1998, Zhou, Liu 2004, Rodríguez, Suarez 2010). Кроме того, он реализован как в Открытых (SAGA), так и проприетарных ГИС (кривизны в ArcGIS, расширение для ArcGIS DEM Surface Tools от Jenness Enterprises).

Алгоритм Zevenbergen-Thorne использует модификацию (2) следующего вида:

$z=Ax^{2}y^{2}+Bx^{2}y+Cxy^{2}+Dx^{2}+Ey^{2}+Fxy+Gx+Hy+I$

(3)

где $A...I$ – коэффициенты аппроксимации, рассчитанные с помощью полиномов Лагранжа на основе 9 значений $z$ в ячейках окна 3×3. Геоморфометрические параметры получаются в результате дифференциациии (3) и решения соответствующих уравнений для центральной ячейки квадратной матрицы 3×3.

Основные геоморфометрические параметры, рассчитываемые на основе производных первого порядка

Фундаментальные геморфометрических параметры уклона и экспозиции взаимосвязаны, т.к. оба эти показателя характеризуют градиент поверхности, т.е. интенсивность изменения ее значений в пространстве, которая может быть выражена производной первого порядка. Как производная поверхности первого порядка, градиент характеризуется величиной (уклоном) и направлением (экспозицией).

Уклон поверхности (Slope)

Понятие

Уклон поверхности – угол наклона в точке пересечения между горизонтальной плоскостью и плоскостью касательной к земной поверхности; фиксирует интенсивность перепада высот (градиент) между двумя заданными точками (рис. 2)

рисунок 2 Рис. 2 XXX

Если земная поверхность представлена функцией $z=f(x,y)$, то уклон рассчитывается с учетом изменений значений $z$ в двух направлениях как ${\frac {\partial z}{\partial xy}}$:

$\alpha =\tan ^{-1}\left({\sqrt {\left({\frac {\partial z}{\partial x}}\right)^{2}+\left({\frac {\partial z}{\partial y}}\right)^{2}}}\right)$

(4)

где ${\frac {\partial z}{\partial x}}$ и ${\frac {\partial z}{\partial y}}$ - производные первого порядка, представляющие изменение значений абсолютной высоты $z$ с запада на восток ($x$) и с севера на юг ($y$).

Расчет по Zevenbergen-Thorne

Процедура определения уклона поверхности по алгоритму Zevenbergen-Thorne сводится к следующим шагам (рис. 3):

рисунок 3 Рис. 3 Определение уклона поверхности (Zevenbergen, Thorne 1987) XXX

1. определение уклона поверхности в направлении с востока на запад:

$G={\frac {\left(z_{6}+z_{4}\right)}{2l}}$

(5)

2. определение уклона поверхности в направлении с севера на юг:

$H={\frac {\left(z_{2}+z_{8}\right)}{2l}}$

(6)

где $z_{2...8}$ - высотные отметки в соответствующих ячейках растра, а $l$ - расстояние между индивидуальными элементами матрицы высот, другими словами – пространственное разрешение растра. При этом предусматривается, что единицы измерения абсолютной высоты и пространственного разрешения идентичны (как правило, метры);

3. определение интегрального значения уклона поверхности для центральной ячейки скользящего окна:

$\alpha =\tan ^{-1}\left({\sqrt {G^{2}+H^{2}}}\right)$

(7)

Рассмотрим алгоритм на примере Рис. 4:

рисунок 4 Рис. 4 Пример скользящего окна: ячейка для которой рассчитывается значение затенена, значения абсолютной высоты выражены в метрах, пространственное разрешение растра – 10 м

Пример расчета уклона поверхности

     
 Рассчитаем уклон по 5: 
$G={\frac {3-10}{2\times 10}}=-0.35$
Рассчитаем уклон по 6: 
$H={\frac {6-7}{2\times 10}}=-0.05$
Тогда уклон поверхности в центральной ячейке, рассчитанный по 7 составит:
$\alpha =\tan ^{-1}\left({\sqrt {\left(-0.35\right)^{2}+\left(-0.05\right)^{2}}}\right)=\tan ^{-1}\left({\sqrt {0.125}}\right)=\tan ^{-1}\left(0.3536\right)=19.47^{\circ }$

Интерпретация

Уклон поверхности фундаментальный морфометрический параметр, который закономерно связан со следующими процессами и характеристиками ландшафта:

• поверхностный сток и дренирование – чем более крутой склон, тем интенсивнее поверхностный сток и меньше инфильтрация влаги в почвенную толщу. Таким образом, уклон имеет принципиальное значение для режима увлажнения почв, особенно – верхних слоев;
• эрозия – интенсивность эрозии растет экспоненциально с увеличением уклона. Это объясняется тем, что с увеличением градиента кинетическая энергия осадков остается постоянной, но транспорт ускоряется в направлении подножья. В результате, кинетическая энергия стока превышает кинетическую энергию осадков, когда склон переходит отметку 8,5°, что и способствует проявлению эрозионных процессов;
• мощность почвенного профиля на склоне закономерно изменяется в соответствии с уклоном и относительной высотой. Как правило, почвенная толща меньше на возвышенных наклонных участках вследствие эрозионных процессов и гравитационного перемещения материала, и постепенно увеличивается в направлении пониженных участков с небольшим уклоном;
• количество солнечной энергии также зависит от уклона, поскольку он определяет угол падения солнечных лучей на земную поверхность. Увеличение уклона поверхности в направлении поступления солнечных лучей увеличивает угол их падения, а значит – количество энергии, которое получает поверхность. Это определяет микроклиматические особенности участка, в частности температуру, эвапотранспирацию и влажность верхних слоев почвы;
• особенности растительного покрова совокупно отражают все вышеперечисленные характеристики, поскольку они прямо или косвенно влияют на такие эдафические факторы как водный и температурный режим почвы, механический состав корнесодержащего слоя, содержание питательных элементов и т.д.

Простота расчета и информативность делают уклон поверхности наиболее употребимым показателем в моделировании процессов перераспределения поверхностного и внутрипочвенного стока, эрозии, определении эдафических условий, индикационном картографировании в физической географии и близких отраслях. Как правило, значения показателя измеряются в градусах (также это могут быть проценты или радианы) и колеблются в диапазоне от 0° (горизонтальная плоскость) до 90° (вертикальная плоскость). Часто используются следующие градации уклона поверхности:

Для равнинных территорий (Жучкова, Раковская, 2004)		Для горных территорий
		(Жучкова, Раковская, 2004)		(Міллер, 1996)
меньше 1°	плоские (субгоризонтальные) равнины	меньше 4°	плоские и почти плоские поверхности	меньше 3°	очень пологие склоны
1-3°	слабонаклонные равнины (очень пологие склоны)	4-10°	пологие склоны	3-6°	пологие склоны
3-5°	пологие склоны (наклонные равнины)	10-20°	покатые склоны	6-9°	слабопокатые склоны
5-7°	слабопокатые склоны	20-30°	склоны средней крутизны	9-12°	покатые склоны
7-10°	покатые склоны	30-45°	крутые склоны	12-15°	сильнопокатые склоны
10-15°	сильнопокатые склоны	45-60°	очень крутые склоны	15-30°	крутые склоны
15-20°	крутые склоны	больше 60°	скалистые (обрывистые) склоны	30-45°	очень крутые склоны
20-40°	очень крутые склоны			больше 45°	обрывистые склоны
больше 40°	обрывистые склоны

Экспозиция (Aspect)

Понятие

Экспозиция поверхности – угол по часовой стрелке между определенным направлением (как правило, на север) и проекцией уклона на горизонтальную плоскость; фиксирует направление (азимут) максимального уклона (градиента) земной поверхности (рис. 5).

рисунок 5 Рис. 5 ?Экспозиция с математической точки зрения?

Для земной поверхности представленной функцией $z=f(x,y)$ экспозиция рассчитывается как угол между двумя производными по формуле:

$\beta =-\tan ^{-1}\left({\frac {\partial z/\partial x}{\partial z/\partial y}}\right)$

(8)

где ${\frac {\partial z}{\partial x}}$ и ${\frac {\partial z}{\partial y}}$ - производные первого порядка, представляющие изменение значений абсолютной высоты $z$ с запада на восток ($x$) и с севера на юг ($y$).

Расчет по Zevenbergen-Thorne

Интегральное значение экспозиции поверхности в центральной ячейке скользящего окна по алгоритму Zevenbergen-Thorne определяется по формуле:

$\beta =-\tan ^{-1}\left({\frac {\pm H}{\pm G}}\right)$

(9)

В данном случае ± определяет в какой четверти находится β по отношению к выбранному направлению. Как правило, экспозиция отсчитывается по часовой стрелке от северного направления географического меридиана, потому (9) приобретает следующий вид:

$\beta =180^{\circ }-\tan ^{-1}\left({\frac {\pm H}{\pm G}}\right)+90^{\circ }\left({\frac {\pm G}{\left|G\right|}}\right)$

(10)

Пример расчета экспозиции поверхности

     
 Рассмотрим алгоритм на примере рис. 4. Уже известно, что по (5) G=-0.35, а по (6) H=-0.05. Тогда экспозиция в центральной ячейке, рассчитанная по (10) составит: 
$\beta =180^{\circ }-\tan ^{-1}\left({\frac {-0.05}{-0.35}}\right)+90^{\circ }\left({\frac {-0.35}{\left|0.35\right|}}\right)=180^{\circ }-\tan ^{-1}\left(0.1429\right)+90^{\circ }\left(-1\right)=180^{\circ }-8.1325^{\circ }-90^{\circ }=81.87^{\circ }$

Интерпретация

Функциональная интерпретация экспозиции может вестись в нескольких направлениях, поскольку она характеризует:

• основное направление линий тока, т.е. когда вода (или другой способный к перемещению материал) движется под действием силы тяжести вниз по склону, он делает это в направлении, которое определяется экспозицией. Эта зависимость положена в основу гидрологических алгоритмов моделирования поверхностного стока.

• ориентацию участка по отношению к потоку солнечных лучей, а значит и количество радиации, получаемой земной поверхностью – инсоляцию. Благодаря этому экспозиция существенно влияет на локальный климат (микроклимат) участка. Например, в северном полушарии склоны южной экспозиции прогреваются лучше, чем северной. Кроме того, вследствие большей эвапотранспирации южные склоны суше северных. Количество солнечной радиации непосредственно определяет интенсивность развития растений и их биологическую продуктивность. Такие закономерности иногда обуславливают довольно-таки существенные азональные и локальные отличия в почвенном и растительном покрове, сезонном течении функциональных процессов.

Отмеченные выше особенности распределения тепла дополнительно усложняются тем, что склоны восточной и западной ориентации формально имея одинаковую экспозицию (по 90°), также демонстрируют контрастность тепловых условий. Это объясняется тем, что на восточные склоны солнечные лучи попадают в утренние прохладные часы и тратятся на прогревание поверхности, а западне склоны освещаются во второй половине дня, когда поверхность уже прогрета. В результате, западне склоны несколько тепле и суше, чем восточные. М. Гродзинський (2013) отмечает, что при весьма значительной крутизне склонов (особенно в горных ландшафтах) склоны юго-западной экспозиции оказываются тепле и суше, чем южной экспозиции, а склоны северо-восточной экспозиции – холоднее и влажнее, чем северные. Еще раньше эту особенность отметил Р. Уиттекер (1980), который самым холодным считает не северное, а северо-восточное местоположение, а наиболее теплым – юго-восточное.

Типология инсоляционных местоположений проводится по сторонам горизонта (румбам)

северное	N	0-22.5°; 337.5-360°
северо-восточное	NE	22.5-67.5°
восточное	E	67.5-112.5°
юго-восточное	SE	112.5-157.5°
южное	S	157.5-202.5°
юго-западное	SW	202.5-247.5°
западное	W	247.5-292.5°
северо-западное	NW	292.5-337.5°

Определение соотношения по теплообеспеченности между этими типами может трактоваться по-разному с учетом того, как будут рассматриваться склоны восточной и западной экспозиций. Поэтому возможно несколько вариантов упорядочивания местоположений в инсоляционный ряд (табл. 1, рис. 4).

Варианты рядов тепло(влаго)обеспеченности по инсоляционной экспозиции (Гродзинський, 2013)

	Холодно => Тепло Влажно => Сухо
Ориентационный ряд	N → NE=NW → E=W → SE=SW → S
Традиционный "компасный" ряд	N → NE → NW → E → W → SE → SW → S
Ряд Уиттекера	NE → N → NW → E → W → SE → S → SW

рисунок 6 Рис. 6 Варианты типологии инсоляционных местоположений (Гродзинський, 2013): а – по ориентационному, б – по «компасному», в – по экологическому принципу (топографический ряд Уиттекера); 1-8 – типы местоположений (увеличение порядкового номера указывает на увеличение сухости местоположения). Стрелкой показано направление движения солнца в светлую пору дня.

Основные геоморфометрические параметры, рассчитываемые на основе производных второго порядка

Для рассмотренных выше производных функции $z=f(x,y)$ в свою очередь могут быть рассчитаны производные. И если первая производная описывала градиент поверхности (его величину и направление), то вторая фиксирует меру изменений этого градиента, т.е. является градиентом первой производной в заданном направлении. На производных второго порядка основывается система морфометрических кривизн, описывающих форму поверхности.

В общем, кривизну в некоторой точке поверхности можно описать как кривизну линии, образованную пересечением земной поверхности и плоскости определенной ориентации, которая проходит через заданную точку. Наиболее часто в геморфометрическом анализе используются горизонтальная (плановая), вертикальная (профильная) и (средняя, общая) кривизна. Рассмотрим особенности их расчета и физический смысл подробнее.

Горизонтальная (плановая) кривизна (Plane Curvature)

Понятие

Горизонтальная (плановая) кривизна – кривизна линии, образованной пересечением земной поверхности с плоскостью, перпендикулярной к направлению ориентации максимального градиента (экспозиции). Как производная второго порядка, горизонтальная кривизна описывает градиент экспозиции вдоль заданного контура (рис. 7).

рисунок 7 Рис. 7 Горизонтальная кривизна

Если земная поверхность представлена функцией $z=f(x,y)$, то плановая кривизна является функцией ее частичных производных:

$PLANC={\frac {\left({\frac {\partial z}{\partial y}}\right)^{2}{\frac {\partial ^{2}z}{\partial x^{2}}}-2{\frac {\partial z}{\partial x}}{\frac {\partial z}{\partial y}}{\frac {\partial ^{2}z}{\partial x\partial y}}+\left({\frac {\partial z}{\partial x}}\right)^{2}{\frac {\partial ^{2}z}{\partial y^{2}}}}{\sqrt {\left(\left({\frac {\partial z}{\partial x}}\right)^{2}+\left({\frac {\partial z}{\partial y}}\right)^{2}\right)^{3}}}}$

(11)

Расчет по Zevenbergen-Thorne (1987)

Для получения значения горизонтальной кривизны на основе алгоритма Zevenbergen-Thorne используют следующее выражение:

$PLANC={\frac {-2\left(DH^{2}+EG^{2}-FGH\right)}{\left(G^{2}+H^{2}\right)}}\times 100^{optionally}$

(12)

В данном случае G и Н рассчитываются по выражениям (5) и (6) соответственно, а другие частные производные как:

$D={\frac {\left[{\frac {\left(z_{4}+z_{6}\right)}{2}}-z_{5}\right]}{l^{2}}}$

(13)

$E={\frac {\left[{\frac {\left(z_{2}+z_{8}\right)}{2}}-z_{5}\right]}{l^{2}}}$

(14)

$F={\frac {\left(z_{3}-z_{1}+z_{7}-z_{9}\right)}{4l^{2}}}$

(15)

Пример расчета плановой кривизны

     
 Рассмотрим алгоритм на примере рис. 4. Уже известно, что по (5) G=-0.35, а по (6) H=-0.05. Рассчитаем другие частные производные:
по (13) $D={\frac {\left[{\frac {\left(10+3\right)}{2}}-4\right]}{10^{2}}}={\frac {2.5}{100}}=0.025$
по (14) $E={\frac {\left[{\frac {\left(6+7\right)}{2}}-4\right]}{10^{2}}}={\frac {2.5}{100}}=0.025$
по (15) $F={\frac {\left(9-4+8-1\right)}{4\times 10^{2}}}={\frac {12}{400}}=0.03$
Тогда горизонтальная кривизна для центральной ячейки, полученная по (12) составит:
$PLANC={\frac {-2\left(0.025\times \left(-0.05\right)^{2}+0.025\times \left(-0.35\right)^{2}-0.03\times \left(-0.35\right)\times \left(-0.05\right)\right)}{\left(-0.35\right)^{2}+\left(-0.05\right)^{2}}}={\frac {-2\left(0.025\times 0.0025+0.025\times 0.1225-0.000525\right)}{0.1225+0.0025}}={\frac {-0.0052}{0.125}}=-0.0416$

Интерпретация

Единицами измерения кривизны земной поверхности являются 1/м, но в результате расчетов получаются в основном очень маленькие значения, потому для облегчения интерпретации рекомендуется умножать их на 100, тогда окончательно единица измерения составит 1/100м. В таком случае кривизна характеризует изменение определенного градиента (для горизонтальной кривизны – экспозиции) на 100 м движения вдоль определенного направления. При интерпретации значений кривизны нужно обращать внимание как на величину самого значения (модуль), так и на его знак.

Модуль кривизны обратный радиусу кривой в данной точке, т.е. линии с плавным изгибом (и большим радиусом) имеют малую кривизну, и наоборот – большие значения кривизны соответствуют сильно изогнутым кривым с малым радиусом. Иными словами, чем больше значение кривизны (без учета знака), тем более выгнутой/ выпуклой является поверхность и наоборот.

рисунок 8 Рис. 8 Кривизна линий в заданной точке составляет 1/R , т.е. обратная радиусу R круга, вписанного в кривую в данной точке. В науках про землю согласовано, что кривизна положительная для выпуклых форм поверхности (R₂>0) и отрицательная – для вогнутых (R¹<0) (Geomorphometry…, 2008)

Знак, который присваивается значению кривизны, зависит от конкретной реализации расчетного алгоритма. Например в оригинальной работе Zevenbergen-Thorne (1987) выражение (12) не содержит знака «–», но реализация его в SAGA предусматривает умножение на –2. В результате, рассчитанная горизонтальная кривизна будет иметь положительные значения для выпуклых поверхностей и отрицательные для вогнутых, нулевые значения отвечают плоским в плане поверхностям (рис. 9).

рисунок 9 Рис. 9 Горизонтальная кривизна перпендикулярна направлению склона и влияет на конвергентность/ дивергентность поверхностного стока

Горизонтальная кривизна описывает так называемый первый механизм аккумуляции, который зависит от способности потока сворачиваться по мере движения по земной поверхности (рис. 10).

рисунок 10 Рис. 10 Первый механизм аккумуляции – реализуется за счет сближения линий тока в плане (Шарый, 2006)

Области с отрицательной плановой кривизной отвечают за вогнутые участки – зоны конвергенции, где происходит схождение линий тока. Области с положительной плановой кривизной характеризуют выгнутые участки – зоны дивергенции, где происходит расхождение линий тока. Благодаря этому плановая кривизна может быть использована для различия гребней, которым свойственен снос материала, и долин, которые этот материал аккумулируют.

Вертикальная (профильная) кривизна (Profile curvature)

Понятие

Вертикальная (профильная) кривизна – кривизна линии, образованной пересечением земной поверхности и вертикальной плоскости. Как производная второго порядка вертикальная кривизна описывает градиент уклона вдоль заданного контура (рис. 11).

рисунок 11 Рис. 11

Если земная поверхность представлена функцией $z=f(x,y)$, то профильная кривизна является функцией ее частных производных:

$PROFC={\frac {\left({\frac {\partial z}{\partial y}}\right)^{2}{\frac {\partial ^{2}z}{\partial x^{2}}}+2{\frac {\partial z}{\partial x}}{\frac {\partial z}{\partial y}}{\frac {\partial ^{2}z}{\partial x\partial y}}+\left({\frac {\partial z}{\partial x}}\right)^{2}{\frac {\partial ^{2}z}{\partial y^{2}}}}{\left(\left({\frac {\partial z}{\partial x}}\right)^{2}+\left({\frac {\partial z}{\partial y}}\right)^{2}\right){\sqrt {\left(1+\left({\frac {\partial z}{\partial x}}\right)^{2}+\left({\frac {\partial z}{\partial y}}\right)^{2}\right)^{3}}}}}$

(16)

Расчет по Zevenbergen, Thorne (1987)

Для получения значения вертикальной кривизны на основе алгоритма Zevenbergen-Thorne в SAGA используется следующее выражение:

$PLANC={\frac {-2\left(DG^{2}+EH^{2}+FGH\right)}{\left(G^{2}+H^{2}\right)}}\times 100^{optionally}$

(17)

По уже известным выражениям рассчитываются D (13), E (14), F (15), G (5) и H (6).

Пример расчета профильной кривизны

     
 Рассмотрим алгоритм на примере рис. 4. Уже известно, что по (5) G=-0.35, по (6) H=-0.05, по (13) D=0.025, по (14) E=0.025, по (15) F=0.03. Тогда горизонтальная кривизна для центральной ячейки, полученная по (17) составит:
$PROFC={\frac {-2\left(0.025\times \left(-0.35\right)^{2}+0.025\times \left(-0.05\right)^{2}+0.03\times \left(-0.35\right)\times \left(-0.05\right)\right)}{\left(-0.35\right)^{2}+\left(-0.05\right)^{2}}}={\frac {-2\left(0.025\times 0.1225+0.025\times 0.0025+0.000525\right)}{0.1225+0.0025}}={\frac {-0.0073}{0.125}}=-0.0584$