ИССЛЕДОВАНИЕ ЗЕМЛИ ИЗ КОСМОСА, 2008, № 6, с. 43-51
МЕТОДЫ И СРЕДСТВА ОБРАБОТКИ И ИНТЕРПРЕТАЦИИ КОСМИЧЕСКОЙ ИНФОРМАЦИИ
УДК 528:551.4
РАСЧЕТ ПРОИЗВОДЯЩЕЙ ФУНКЦИИ ВЫСОТЫ ДЛЯ ВЫДЕЛЕНИЯ СТРУКТУРНЫХ ЛИНИЙ РЕЛЬЕФА ПО СПУТНИКОВЫМ ДАННЫМ
И ТОПОГРАФИЧЕСКИМ КАРТАМ © 2008 г. И. В. Флоринский
Институт математических проблем биологии РАН, Пущино Московской обл.
Тел.: (4967) 73-24-08, e-mail: iflor@mail.ru Поступила в редакцию 13.03.2008 г.
Выделение линий экстремальной кривизны земной поверхности — структурных линий рельефа — является важным этапом при решении многих прикладных и фундаментальных задач геоморфологии, геологии и смежных дисциплин. Структурные линии могут быть аналитически описаны производящей функцией высоты, которая зависит от первых, вторых и третьих частных производных высоты. Представлен вывод формул для расчета третьих частных производных по цифровой модели высот, заданной на квадратной сетке. Приводится пример выделения линий экстремальной кривизны рельефа с помощью расчета цифровой модели производящей функции.
ВВЕДЕНИЕ
Выделение линий экстремальной кривизны земной поверхности — структурных линий рельефа (водоразделов, тальвегов, бровок и подножий) — является важным этапом при решении многих прикладных и фундаментальных задач геоморфологии, геологии и смежных дисциплин. Для автоматизированного выделения структурных линий по цифровым моделям высот (ЦМВ) предложено большое число алгоритмов [1—7]. В отличие от
методов расчета локальных морфометрических характеристик, основанных на подходах дифференциальной геометрии [8, 9], выделение тальвеговой и водораздельной сети, как правило, основано на логических процедурах анализа ЦМВ, таких как пошаговое определение направления линии тока и пр.
Вместе с тем структурные линии рельефа могут быть аналитически описаны с помощью производящей функции высоты (Т):
T = 1 {qa - 3pq2b + 3p2qc -p3d +
J(p2 + q )3(i + p2 + q2)
2 2 2 2 2 2 1 (1) + ( q2r - 2pqs + p21) [pq(t - r) + s(p2 - q2) ] ( 2 + 3p 2 + 3q ) I
(p2 + q2)(i + p2 + q2) j
Если высота г задана в виде г = /(х, у), где х, у — декартовы координаты, то
a=
d3z дх3
d=
д3 z д j3
b=
д3 z дх2д y
r=
д2 z
д x
2
t=
д2 z д y2
д-z
дх
p = -т, q = -г.
д-z д y
c=
д3z
дхд y
2
s=
дх дy,
(2)
Теоретически обосновано, что линии водоразделов и бровок являются множеством точек Т = 0 при положительных значениях горизонтальной кривизны земной поверхности (кь), а линии тальвегов и подножий — множеством точек Т = 0 при къ < 0 [10].
Значения первых и вторых частных производных высоты г, ж, р и q (2) могут быть рассчитаны по ЦМВ, заданной на квадратной сетке, с помощью методов аппроксимации частных производных конечными разностями. Благодаря высокой точности и устойчивости к шуму наибольшее рас-
16
21
12
17
22
13
18
23
14
19
24
10
15
20
25
и
(0, —и, ^18), (и, —и, ^19), (2м, -и, £20), (—2м, —2м, £21), (—и, —2м, £22), (0, —2м, £23), (и, —2м, £24) и (2м, —2м, £25) известны их декартовы координаты и высоты.
Приблизим полином третьего порядка [14]
1 3,1,3,1,2 ,1 21 21.2, z = - ах + - ау + - Ьх у + - сху + - гх + - 1у + 6 6 2 2 2 2 (3)
+ зху + рх + ду + и
к скользящему окну 5 х 5 с помощью метода наименьших квадратов [15]. Записав полином (3) для 25 точек скользящего окна, получаем систему 25 линейных нормальных уравнений, которая может быть представлена в виде
Рис. 1. Скользящее окно 5 х 5. 1, ..., 25 — номера точек, м — шаг сетки.
а = рр,
(4)
у
3
4
8
х
и
пространение получил метод Эванса [11, 12]. В соответствии с этим методом [13], полином второго порядка приближается методом наименьших квадратов к точкам скользящего окна (матрице высот) размера 3 х 3 с шагом и. Для точек окна известны декартовы координаты и высоты. Значения г, ж, р и # определяются для центральной точки окна. Перемещая скользящее окно 3 х 3 вдоль ЦМВ, можно рассчитать значения г, ж, р и # для всех точек ЦМВ, кроме крайних строк и столбцов.
Насколько нам известно, до сих пор не опубликованы формулы для вычисления третьих частных производных высоты а, й, Ь и с (2) по ЦМВ. Кроме того, отсутствуют работы, где показана практическая возможность использования производящей функции (1) для выделения структурных линий рельефа по ЦМВ. В данной статье дается вывод формул третьих частных производных высоты для их расчета по ЦМВ, заданной на квадратной сетке. Приводится пример выделения линий экстремальной кривизны рельефа по ЦМВ с помощью расчета цифровой модели производящей функции.
ВЫВОД ФОРМУЛ
Рассмотрим скользящее окно размера 5 х 5 (рис. 1). Для точек этого окна (—2м, 2м, £1), (—и, 2м, £2), (0, 2м, £3), (и, 2м, £4), (2м, 2м, £5), (—2м, и, £6), (—и, и, £7), (0, и, £8), (и, и, £9), (2м, и, £10), (—2м, 0, £П), (—и, 0, £12), (0, 0, £13), (и, 0, £14), (2м, 0, £15), (—2м, —и, £16), (—и, —и,
где а — матрица значений высот в точках скользящего окна размера 25 х 1:
а=
¿1 ¿2
¿25
(5)
в — матрица неизвестных коэффициентов ряда (3) размера 10 х 1:
в =
(6)
и Р — матрица известных коэффициентов системы уравнений (4) размера 25 х 10:
4 3 — Ц 3 43 - Ц 3 3 4Ц 3 -4 Ц 2Ц2 2 Ц 2 -4Ц -2 Ц 2Ц
1 3 --Ц 6 43 - Ц 3 3 Ц -2 ц3 1 2 -ц 2 2 Ц -2ц2 -Ц 2Ц
0 43 - Ц 3 0 0 0 2 Ц 0 0 2Ц
13 - Ц 6 43 - Ц 3 3 Ц 2 ц3 12 -ц 2 2 Ц 2ц2 Ц 2Ц
43 - Ц 3 43 - Ц 3 3 4Ц 3 4ц 2Ц2 2 Ц 2 4ц 2Ц 2Ц
43 --Ц 3 13 - Ц 6 2Ц3 3 -Ц 2Ц2 12 -ц 2 -2Ц2 -2 ц Ц
13 --Ц 6 13 - Ц 6 13 - Ц 2 13 --Ц 2 12 -ц 2 12 -ц 2 2 -Ц -Ц Ц
0 13 - Ц 6 0 0 0 12 -ц 2 0 0 Ц
13 - Ц 6 13 - Ц 6 13 - Ц 2 13 -Ц 2 12 -ц 2 12 -ц 2 2 Ц Ц Ц
43 - Ц 3 13 - Ц 6 2Ц3 3 Ц 2Ц2 12 -ц 2 2Ц2 2Ц Ц
43 --Ц 3 0 0 0 2Ц2 0 0 -2 ц 0
13 --Ц 6 0 0 0 12 -ц 2 0 0 -Ц 0
0 0 0 0 0 0 0 0 0
13 - Ц 6 0 0 0 12 -ц 2 0 0 Ц 0
43 - Ц 3 0 0 0 2Ц2 0 0 2 ц 0
43 --Ц 3 13 --Ц 6 -2Ц3 3 -Ц 2Ц2 12 -ц 2 2Ц2 -2 ц -Ц
13 --Ц 6 13 --Ц 6 13 — Ц 2 13 --Ц 2 12 -ц 2 12 -ц 2 2 Ц -Ц -Ц
0 13 --Ц 6 0 0 0 12 -ц 2 0 0 -Ц
13 - Ц 6 13 --Ц 6 13 — Ц 2 13 -Ц 2 12 -ц 2 12 -ц 2 2 -Ц Ц -Ц
43 - Ц 3 13 --Ц 6 -2Ц3 3 Ц 2Ц2 12 -ц 2 -2Ц2 2Ц -Ц
43 --Ц 3 43 --Ц 3 3 -4Ц 3 -4 ц 2Ц2 2 Ц 2 4Ц -2 ц -2 ц
13 --Ц 6 43 --Ц 3 3 -Ц -2 ц3 12 -ц 2 2 Ц 2Ц2 -Ц -2 ц
0 43 --Ц 3 0 0 0 2 Ц 0 0 -2 ц
13 - Ц 6 43 --Ц 3 3 -Ц 2 ц3 12 -ц 2 2 Ц -2Ц2 Ц -2 ц
43 - Ц 3 43 --Ц 3 3 -4Ц 3 4ц 2Ц2 2 Ц 2 -4ц 2Ц -2 ц
(7)
ь
Чтобы найти неизвестные коэффициенты полинома (3), необходимо решить уравнение
в = ()-У а, (8)
где Ьт — транспонированная матрица Ь, а (ЬТЬ)-1 — матрица, обратная матрице ЬТЬ. С помощью матричных операций находим матрицу (ЬТЬ)-1ЬТ размера 10 х 25:
т -1 т (FTF) FT
5w 1
5w3 l
1
lOw
__1_
5w l
3
70w 35w'
__1___1
3 „ - 3
70 w l
35 w 1
35 w l
lOOw 31 2lO w
31 2 10 w
H
l75
35w' 2
35 w2
__l_
35 w1 1
50w2 1
84w 11 lO5 w
175
1 1 O 1 1 1 1 O
lOw3 5w3 5w3 10 w3 lOw3 5 w3
1 1 1 1 1 1 1 1
lOw3 lOw3 lOw3 lOw3 10 w3 5 w3 5 w3 5w3
2 1 2 1 2 1 1 1
35w3 35 w3 35w3 35w3 35 w3 35 w3 70w3 35 w3
2 1 O 1 2 1 1 O
35w3 35 w3 35w3 35 w3 35 w3 70 w3
2 1 2 1 2 2 1 2
35w2 35 w2 35w2 35w2 35w2 35 w2 35w2 35 w2
2 2 2 2 2 1 1 1
35w2 35 w2 35w2 35w2 35w2 35 w2 35w2 35 w2
1 1 O 1 1 1 1 O
25w2 5Ow2 5Ow2 25w2 5Ow2 100 w2
31 ll O 11 3l 1 31 O
42O w 105 w 1 05w 4 2O w 8 4w 210 w
3l 1 17 1 3l 11 31 17
42O w 84w 420 w 84w 4 2O w 105 w 210 w 105 w
13 2 1 2 1 3 2 17 22
175 l75 25 l75 175 l75 175 175
1 1 O 1 1 1 1 O
lOw3 5w 5w lOw3 10 w3 5w
O O O O O 1 5w3 1 5 w3 1 5w3
O O O O O 1 35w3 1 70 w3 1 35 w3
2 1 O __l_ __2_ 1 1 O
Q - 3 35 w 35w3 35w3 35 w3 35w3 70 w3
2 1 2 1 2 2 1 2
35 w2 35w2 35 w2 35w2 35 w2 35w2 35 w2 35 w2
2 2 2 2 2 1 1 1
2 35 w 35w2 35 w2 35w2 35 w2 35w2 35 w2 35 w2
O O O O O 1 5Ow2 1 100 w2 O
l7 17 O 17 17 1 31 O
4 2O w 1 05w 105 w 420 w 84w 210 w
O O O O O 11 105 w 31 210 w 17 105 w
1 22 27 22 1 2 17 22
(9)
25
l75 l75 l75
25
l75 l75
l75
1 1 1 1 0 1 1
С 3 5 у 10 у3 10у3 С 3 С 3 10у3
1 1 1 1 1 1 1
5у3 - 3 5у 10у3 10 у3 10у3 10у3 10 у3
1 1 2 1 2 1 2
70 у3 35у3 35 у3 35 у3 35у3 35у3 35 у3
__1_ __1_ __2_ __1_ 0 1 2
70у3 35у3 35 у3 35у3 35у3 35 у3
1 2 2 1 2 1 2
35у2 35 у2 35 у2 35 у2 35 у2 35 у2 35 у2
1 1 2 2 2 2 2
35у2 35 у2 35 у2 35 у2 35у2 35 у2 35 у2
1 1 1 1 0 1 1
100 у2 50 у2 25 у2 50 у2 50у2 25 у2
31 1 31 11 0 11 31
210 у 84 у 420 у 105 у 105 у 4 20 у
31 11 31 1 17 1 31
210 у 105 у 420 у 84 у 4-20 у 8 4у 420 у
17 2 1 3 2 1 2 1 3
175 175 175 175 25 175 175
Решая уравнение (8), получаем следующие выражения для частных производных высоты а, й, Ь, с, г, ж, р и
а = ¿5 + ¿1 0 + ¿15 + ¿20 + ¿25 - ¿1 - ¿6 - ¿1 1 - ¿1 6 - ¿2 1 +
10у
(10)
+ 2 (¿2 + ¿7 + ¿ 12 + ¿17 + ¿22 - ¿4 - ¿9 - ¿14 - ^ 19 - ¿24 )
10 у
^ _ ¿1 + ¿2 + ¿3 + ¿4 + ¿5 - ¿2 1 - ¿22 - ¿2 3 - ¿24 - ¿25 +
+
10 у
2 (¿ 16 + ¿17 + ¿ 18 + ¿19 + ¿20 - ¿6 - ¿7 - ¿8 - ¿9 - ¿10 )
10у3
(11)
ь = ¿1 7 + ¿ 19 - ¿7 - ¿9 + 4 (¿ 1 + ¿5 + ¿23 - ¿3 - ¿21 - ¿25) +
70у3 (12) + 2 (..¿6 + ¿ 10 + ¿1 8 + ¿22 + ¿24 - ¿2 - ¿4 - ¿8 - ¿16 - ¿20)
70 у3
с _ ¿7 + ¿17 - ¿9 - ¿ 19 + 4 (¿5 + ¿ 11 + ¿25 - ¿ 1 - ¿15 - ¿21 ) + С = 3 +
70у (13)
2(¿4 + ¿6 + ¿12 + ¿16 + ¿24 - ¿2 - ¿10 - ¿14 - ¿20 - ¿22)
1
;[ 2( ¿1 + ¿2 + ¿3 + ¿4 + ¿5 + ¿21 + ¿22 + ¿23 +
г =
35 у'
+ ¿24 + ¿25 ) - 2(¿11 + ¿12 + ¿13 + ¿14 + ¿15) " (15) - ¿6 —¿7 — ¿8 - ¿9 - ¿10 - ¿16 - ¿17 - ¿18 - ¿19 - ¿20 ] ,
5 =
1
: [ ¿9 + ¿17 - ¿7 - ¿19 + 4 (¿5 + ¿21 - ¿1 - ¿25) +
100у (16)
+ 2 (¿4 + ¿10 + ¿16 + ¿22 - ¿2 - ¿6 - ¿20 - ¿24) ] ,
1
420у
{ 44(¿4 + ¿24 - ¿2 - ¿22) + 31 ¡¿1 + ¿21 - ¿5 -
-¿25 + 2(¿9 + ¿19- ¿7- ¿17)] + 17 [¿15- ¿11 + (17) + 4(¿14 - ¿12) ] + 5(¿10 + ¿20 - ¿6 - ¿16) } ,
1
420у
{ 44(¿6 + ¿10 - ¿16 - ¿20 ) + 31 [¿21 + ¿25 - ¿1 -
-¿5 + 2(¿7 + ¿9- ¿17- ¿19)] + 17 [¿3- ¿23 + (18) + 4(¿8 - ¿18) ] + 5 (¿2 + ¿4 - ¿22 - ¿24) } •
+
Г =
1
35 у
70у3
■-,[ 2 (¿1 + ¿5 + ¿6 + ¿10 + ¿11 + ¿15 + ¿16 + ¿20 +
Как и в методе Эванса, в нашем случае полином третьего порядка (3) приближается к значениям высот в точках скользящего окна 5 х 5, а не проходит строго через них. Это ведет к локальному подавлению шума и может оптимизировать + ¿21 + ¿25) - 2 (¿3 + ¿8 + ¿13 + ¿18 + ¿23) - (14) расчет частных производных, так как они весьма
чувствительны к высокочастотной
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.