Построение зон фактического затопления местности при разливах рек на основе обработки радиолокационных космических снимков и анализа характерных точек рельефа Текст научной статьи по специальности «Строительство и архитектура»

УДК 777.555

СИСТЕМНЫМ АНАЛИЗ, УПРАВЛЕНИЕ И ОБРАБОТКА ИНФОРМАЦИИ

С. В. ПАВЛОВ, Ю.Н. КУНАКОВ

ПОСТРОЕНИЕ ЗОН ФАКТИЧЕСКОГО ЗАТОПЛЕНИЯ МЕСТНОСТИ ПРИ РАЗЛИВАХ РЕК НА ОСНОВЕ ОБРАБОТКИ РАДИОЛОКАЦИОННЫХ КОСМИЧЕСКИХ СНИМКОВ И АНАЛИЗА ХАРАКТЕРНЫХ ТОЧЕК РЕЛЬЕФА

Описывается метод построения на цифровой топографической карте зон фактического затопления территорий при разливах рек на основе радиолокационных космических снимков и анализа характерных точек рельефа местности. Приводятся результаты его применения и оценки. Зоны фактического затопления; радиолокационный космический снимок; характерные точки рельефа местности; географическая привязка космических снимков; цифровая топографическая карта

ВВЕДЕНИЕ

Одной из важных задач при осуществлении информационной поддержки мониторинга паводковой ситуации является построение на цифровой топографической карте зон фактического затопления местности при разливах рек, для выделения которых успешно применяются данные дистанционного зондирования Земли из космоса, особенно космические радиолокационные снимки, на качество которых облачность, сопровождающая паводок, оказывает минимальное влияние [1, 2]. Зоны фактического затопления, выделенные с привязанных по орбитальным данным радиолокационных космических снимков, полученных со спутника Radarsat-1 и являющимися одними из наиболее подходящих по своим параметрам для решения данной задачи [3], имеют погрешность совмещения с цифровой топографической картой масштаба 1:200000, достигающую 1000 м. Ввиду того, что радиолокационные космические снимки имеют большие размеры кадра (порядка 100^100 км), радиолокационная съемка ведется под углом к надиру и радиоволны более подвержены электромагнитным помехам, нужно «допривязать» радиолокационные космические снимки к векторной карте с высокой точностью, которая необходима для выделения зон затопления при мониторинге паводка, с помощью классических методов становится невозможным по следующим причинам. Во-первых, для точной привязки необходимо большое количество опорных точек (порядка 100 точек на снимок 100^100 км), что приводит к большим

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

Рис. 1. Космический снимок Radarsat, привязанный по орбитальным данным

Учитывая вышеперечисленные недостатки, авторами предложен метод построения зон фактического затопления местности при разливах рек на основе радиолокационных космических снимков и анализа характерных точек рельефа (далее для простоты будем называть его ПЗК), позволяющий исключить влияние описанных выше факторов и совмещать ин-

Контактная информация: (347) 273-77-33

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

1. ОПИСАНИЕ МЕТОДА

Метод ПЗК заключается в нахождении крутых берегов с помощью анализа характерных точек рельефа местности, создании множества опорных точек путем выбора на крутых берегах точек максимальной крутизны, выделении фрагментов космического снимка, содержащих по одному крутому берегу, привязки фрагментов, совмещая крутые берега на космическом снимке с теми же берегами на карте с использованием опорных точек, и объединении привязанных фрагментов с интерполяцией перекрывающихся областей. Далее подробнее рассмотрим этот метод, наглядно изображенный на рис. 2.

Рис. 2. Блок-схема алгоритма, основанного на разработанном методе

1. 1. Нахождение участков крутых берегов на цифровой топографической карте

Анализ характерных точек рельефа местности на цифровой топографической карте позволяет найти крутые берега рек, на которых выбираются опорные точки для привязки зон затопления с космического снимка. Полигональный водный объект можно представить как совокупность ломаных линий, ограничивающих этот объект.

С водными объектами в топологических отношениях находятся (непосредственно примыкают) линейные объекты - берега. Они представлены в виде ломаных линий

Сь = {({(-^,УЧ )ч},А*гЬ]).}у =1 /

= 1,4, д = 1,Ь}, (1)

где х 1 - координата долготы д-ой точки у-ой

ь,

линии, у 1 - координата широты д-ой точки у-

ь1

ой линии, Ліг = {аЦ,а1х2,..,а1;гпб} - совокупность атрибутивных характеристик каждой линии.

Ломаная линия состоит из точек, поэтому берег можно представить как множество точек

1 = 1, <і,

(2)

где Грбп1 - у-ая точка берега, хРпі 1 - координата

долготы 1-ой точки, у рп,

координата широты

у-ой точки, Ліг 1 - совокупность атрибутивных характеристик каждой точки.

Для удобства описания обозначим левый берег СЬ л и правый берег СЬ'П'. Принадлежность определяется относительно направления течения реки. Левый берег - берег, находящийся слева от направления течения реки, правый -соответственно справа (рис. 3).

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

Сьперп ={({(хь,уь)д},Лігь ),^

1 = 1, пГ«, д = 1^.

(3)

Среди всего множества точек пересечения

^ лив

можно найти координаты точек Срп , принад-

Си в

ь , находящихся на одной перпендикулярной прямой с соответствующими точками берегов С рП

Сри^' ={((хРп‘У, УРПУ ),Л1гРп‘У). } у =

Српперп ={((*Рп‘у,

у Рп‘у ),ЛігРпіі),} і = 1, (4)

Точки пересечения присутствуют как на левом, так и на правом берегу. Поэтому можно выделить следующие подмножества:

• Сб =

• ^Рпг

бл

і = 1, пРп1 - множество точек

пересечения левого берега со створами реки;

• Сб" =

к} і

= 1, пРп1 - множество точек

пересечения правого берега со створами реки.

Берега могут иметь как крутые участки, так и пологие. Для точек берега, находящихся на пересечениях со створами, определяется с помощью проверки правил принадлежность к определенному участку берега (крутому или по-

логому). Так как точность привязки не может быть выше пространственного разрешения космического снимка, а также необходимо отслеживать минимальное поднятие уровня воды на 1 метр, в правилах сравнивается пространственное разрешение космического снимка и затапливаемая площадь при изменении уровня воды на 1 метр

Если а < 90о и ^4к,к+1 = 0 и йки = 1 м, то 2\;

Я2: Если а > 90о и Рт > ШИкМ1 и йки = 1 м, то 2\;

Я3: Если а > 90о и Рт < й\ИкМХ и йки = 1 м, то 22;

где - имя /-го правила, а - угол между берегом и водной поверхностью, Рт - пространственное разрешение космического снимка, ^4к,к+1 - расстояние между точками пересечения

( Тш,

(5)

,и.ч и +.) створа Сьперп

и изолиниями вы-

^и.в.

сот С и

'и

С и

ближайшими к берегу,

разница между уровнями воды и5 и

и5+1, 21 - точка ГРп^ принадлежит крутому берегу, 22 - точка Грбп1 принадлежит пологому берегу.

Если при поднятии уровня воды на 1 м площадь затопления берега или расстояние между изолиниями меньше пространственного разрешения космического снимка, то будем считать, что точка принадлежит крутому берегу, и примем площадь его затопления равной нулю. Таким образом, опорные точки выбираются только на крутых берегах, чтобы сократить погрешность при совмещении зон затоплений с векторной картой.

В результате применения правил к множествам точек обоих берегов Срп1 и Срп1 получаем их подмножества: С р^- множество точек, принадлежащих крутым участкам левого берега, Срп - множество точек, принадлежащих

пологим участкам левого берега, Срп - множество точек, принадлежащих крутым участкам

правого берега, Ср^- множество точек, принадлежащих пологим участкам левого берега, для которых справедливы следующие выражения:

(6)

Таким образом определены крутые и пологие участки берегов.

Правый берег

Поперечный профиль реки

Перпендикулярная прямая (створ реки] к водному объекту Сперп

СЬ

С

л

Рис. 3. Схема реки с крутым и пологим берегом

1. 2. Выделение множества опорных точек

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

(7):

х'= /(х, у ); у'= g(х, у ), (7)

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

Наиболее простая модель функций / и g -линейная (8), позволяющая определить смещение, вращение, разницу масштабов в х- и у-направлениях и некоторые виды искажений:

х'= а1 + а 2 х + а3 у; У = а 4 + а5 х + а 6 у.

(8)

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

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

СРп; и СРп; точки с максимальной крутизной

б л б п

выделяются в подмножества СР1П:т“ и СРгк;т“ :

Српг = тІпСрп^} І = 1 Пперп , п (її г п -| _______

Српг = і =1, пперп.

(9)

Аналогично можно выделить подмножества

бл бп

точек с минимальной крутизной СР1т‘п и СР1™т :

(Ии;

перп ’

СрпгІп = тахКбп і }=1, п

(10)

(Ии;

перп

По приведенным ниже правилам определяются опорные точки:

г~Я г л г л ^л

Я,: Если Тр6;, е С& иТр6; е СРцгтах

гп гп гл

и Тр^ е С^, то / е СРппи

Гп гп гп р,п

Я-2- Если тРп,, е Срп и Т^пг, е СРпГ и ТЬ еСб;, то ТРИ; / е СРп:г;

г~л г л л ^л

Яз: ЕслиТр6п;/ е СЩ и Тр^ е С^Т

п

(11)

/-и гп гп Р,п

тт Тб а тт ТЬ а С к,тах

и грпіі Є Српі и грт; Є Српі

и (їЛ < (їм , то гр“ Є Срп1 ;

\ ,к+1 Ік ,і+1 ’ рп11 рп1 ’

Г Є С пр

Г>,л+ ^ Г>,л+

/~Л /~Л /~Л (^.Л

Я4-. Если Грбпіі є Српі и Грбпіі є Србпг

гп гП гП ап

т, Тб а С^х тл Тб а Г’бк,тах и ^рпг^ ^рпг и ^рпг^ ^рш

и (їк,>+1 > ^+l, то Грбп;,Є

где - имя і-го правила, і = 1, пперп , СрЩ - множество опорных точек привязки на векторной карте, Грбп1 - точка, принадлежащая левому берегу, Грбп1. - точка, принадлежащая правому берегу.

Таким образом, множество опорных точек привязки на векторной карте состоит из совокупности точек максимальной крутизны левого и правого берега:

Лп /^пр _^бк,тах . . ґ^бк,тах

^рпі “ ^рпі ^ ^рпі •

(12)

1. 3. Определение областей, выделяемых с космического снимка

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

бл

берег переходит в пологий и наоборот) СРп” и

Cpf1 , определяются на основе следующих правил:

,тЛ р.л г" рП

и . Тн лтттт T6extr _ 776к,ma^ T^extr _ T76",min

Rl: Еслиipn:? = lpn4 и-ipnt; = ^

^л бл 6" 6"

и Г extr = Г ",min и Г extr = Г K>max и JPn: і+і JPnt і+і и 1 Pnti+l Pnti+l

и dl% = dlhj, то Г* j є C^

тт T6 d Г*6 aver . и ^Pn:j є CPn: ,

п бп бл бл

Ї? • Th/*tttt T6extr _ T^^max T, TT6extr _ T76n,min r2: Если ^Pnt. = ipnt? и Jpn:? = ^

" б" бл бл

и Г extr = Г »>min и Г extr = Г к ’тах и 1 Pn: і+і JPn: i+i и Pnti+i Pnti+i

hj hj

лл

то Г є C aver

? то -'Pntj ^Pn:

и Г 6" є C 6"ver •

и Pn:j ^Pn: ?

бл бл п бп

/? • Th/*tttt 776extr _ T K,max 1Ж T7 6 extr _ T’°n,mm R3: Если ipnt? = Jpnt? и Jpn:? = ip

л л " "

и Г extr = Г »>mln и Г extr = Г »>mln

и JPn:і+і Pnti+i и Pnti+i Pnti+i

л л " "

и Г extr = Г K>max и Г extr = Г »>mln

и Pnti+2 Pnti+2 и Pnti+2 Pnti+2

(13)

Причем точки с индексами Трпк;,,тах и ТРп”-™п находятся между точками с индексами

л " л "

Г aver Г aver и Г aver Г aver

JPnt j Pn: j и JPnt j+^JPn: j+i ■

Области выделения фрагментов космического снимка описываются двумя парами координат (левый нижний угол и правый верхний угол):

A-0*1 Pol = {AOIPolt }= }

{( Xmin,l, ymin,l ), ( Xmax,l, Уmax,l ) }

Xmin,k = тІП (Ak (x)) k = 1 n6K ,

7, (іб)

min,k

= min(Ak (y)), k = 1,7

Xmax,k = max (Ak (x)) k = M

Ymax,k = max(Ak (y)\ k = 1, П6к .

В результате определены области, имеющие по одному крутому участку берега. Эти области используются для выделения фрагментов с космического снимка.

1. 4. Выделение фрагментов с космического снимка

Тп Т6 С reaver „ Т6 с. 6aver •

то ipn:j є CPn: и ipn:j є CPnt ,

6" 6" 6л 6л •

ї? • Тн/*тттт 776extr _ Т к,max т, T76extr _ T76n,min r4: Если ipn^ = Ipn:, и ipn:, = ipnt

1 Pnti

и Г6extr = Г »>min и Г6extr = Г »>mln

и J Pn: і+і JPn: і+і и 1 Pnti+i Pnti+i

a." 6Л г^л 6л

и Г extr = Г K>max и Г extr = Г »>mln

и 1 Pnti+2 1 PMi+, и 1 Pnti+2 1 Pnti+2

и dl'j = max dl' ,

/—л /—л /— Л r~ п

™ T6 d Г*6 aver TT T6 d (~*6 aver

то 7Pn: j є CPn: и ipn: j є CPnt

б л бп

где Я - имя /-го правила, Т^, Т^ - береговые точки-экстремумы, являющиеся точками максимальной либо минимальной крутизны,

1 < j < 1 + ^ 1 = 1, nextr , j = 1, Ппєрп .

Следовательно, справедливы выражения:

C 6л П C 6”er

Pn: i^Pn:

C бк П C aver і ■'-Pn:

бл aver

Pn:

л

C к u C

'-Pn: w ^

C6" u C6

'-TW ч~/ '—Г)

л

П Cp6l: =0,

"

П C =0

u C бл = C бл

^'-Pnt '-"Pnu

U C 6" = C 6"

w<^Pn: '-Pn:-

(14)

Таким образом, можно определить точки (15), которые должны попасть в области выделенных фрагментов космического снимка:

A = {Ak } =

л

Г aver

Pn: j ’

Г б aver

pn: j+1

Г б aver

JPn: j

л

Г б aver ,

Pn: j+i5 бл(") б"(л)

гГ,бк ,max ГП 6n,mln

1 Pn: k Г

Pn: k

(15)

б

Космический снимок можно представить как множество пикселей с парами координат, в которых хранятся значения спектральной яркости каждого спектрального канала:

Рг ={д/г }=_{( Хк, уі)і , К, 0^ _

i = 1,n6 ,k = 1,n,l = 1,m,s = 1,nk,t = 1,na.

(1l)

где пб - количество опорных точек (фрагментов

изображения), п - количество строк пикселей в изображении, т - количество столбцов пикселей в изображении, пк - количество спектральных каналов, па - количество значений спектральной яркости.

Множество фрагментов, которые выделяются с космического снимка (рис. 4), определяется пересечением множества точек космического снимка и областей выделения фрагментов:

п/г ={п/г}= п п{А01рок },/ = 1П; (18)

1. 5. Привязка выделенных фрагментов космического снимка

Для совмещения выделенных фрагментов космического снимка с векторной топографической картой используется полиномиальная трансформация (19), но с использованием только одной опорной точки для каждого фрагмента:

x = a 2 x, y = a6 y.

(19)

л

л

и

б

x

и dl л = max dl л

h

"

aver

Таким образом, использование только одной опорной точки при трансформировании фрагментов космического снимка позволяет избежать таких искажений, присущих аффинному преобразованию, как масштабирование по X и/или У, наклон по X и/или У, поворот, отражение (рис. 5).

Рис. 4. Мозаика выделенных фрагментов космического снимка и области выделения фрагментов АО1Ро1

1. 6. Объединение привязанных фрагментов космического снимка

При объединении привязанных фрагментов космического снимка нужно учитывать, что части фрагментов могут накладываться друг на друга и в этом случае области наложения требуют дополнительной обработки, которая заключается в определении фрагмента, с которого пиксели войдут в объединенное изображение. Выбираются пиксели фрагментов изображения с наименьшей яркостью, которая соответствует водной поверхности на космическом снимке:

дг.р™. = тт(Д& ■рг (а) п Д+рг (а)). (20)

а

Объединение привязанных фрагментов можно описать следующим образом:

прг = (п^ ии д+Г)/ (21)

(Дгрг п Д+рг), ( )

где Дгрг - привязанные фрагменты снимка, Прг -мозаика привязанных фрагментов.

Таким образом, формализован метод построения зон фактического затопления местности при разливах рек на основе радиолокационных космических снимков и анализа характерных точек рельефа. Разработанный алгоритм отражает все этапы построения зон фактического затопления местности.

2. ПРИМЕНЕНИЕ РАЗРАБОТАННОГО МЕТОДА

Разработанный метод был использован при построении зон фактического затопления местности на основе космических снимков Radarsat в период весеннего паводка 2008 года в окрестностях города Уфы (рис. 1, 4 и 5). Визуальный анализ результатов выделения зон затопления с обработанного космического снимка с помощью разработанного метода показывает эффективность применения данного метода совмещения зон затопления с цифровой топографической картой (рис. 5). Конвертирование зон затопления со снимков в векторный формат представления данных осуществляется с использованием известных методов дешифрирования изображений, сглаживания граничных областей между соседними зонами затоплений и создания векторных слоев с результатами обработки космических снимков [4]. На рис. 6 приведено сравнение результатов выделения зон затопления со снимков Radarsat, обработанных разработанным методом и классическим способом. Из рисунка видно, что увеличивается точность совмещения зон затопления, полученных с космического снимка, обработанного при использовании метода ПЗК.

3. ОЦЕНКА РАЗРАБОТАННОГО МЕТОДА

Для количественной оценки географической привязки космических снимков в геоинформа-ционных технологиях определяют точность совмещения всех опорных точек. Точность может быть измерена математически, путем сравнения известного местоположения точки в координатах карты и ее положения на трансформированном растре. Расстояние между этими двумя местоположениями точки называется остаточной ошибкой [5]. Общая среднеквадратическая ошибка (СКО) преобразования вычисляется путем извлечения квадратного корня из суммы квадратов погрешностей всех опорных точек и деления на их количество. Это значение показывает, насколько отличается результат трансформирования от идеального совмещения всех опорных точек. Но этот способ не позволяют оценить точность привязки, поскольку определяет только несовпадение известного местоположения опорных точек в координатах карты и их положения на трансформированном растре. При использовании метода ПЗК каждый фрагмент привязывается по одной опорной точке, в этом случае расстояние между парой точек всегда будет равно нулю и соответственно СКО будет тоже равна нулю.

При построении зон затопления важно оценить точность наложения выделенных зон затопления с космического снимка и участков крутых берегов рек на векторной топографической карте. Поэтому экспериментально было выяснено, что определять ошибку привязки необходимо с помощью оценки площади несоответствия между полигональным векторным слоем рек и зонами затопления, выделенными с космического снимка. Таким образом, площадь несоответствия между зонами затопления, полученными с космического снимка, обработанного с использованием метода ПЗК, и векторным слоем рек составляет 1990000 м , а между зонами затопления, полученными с космического снимка, привязанного классическим способом, с тем же векторным слоем рек 3310000 м , что говорит о повышении точности совмещения зон фактического затопления с цифровой топографической картой на 40%. Различие между са-

2

мими зонами затопления составляет 3410000 м (рис. 6). Эта количественная оценка подтверждает эффективность использования данного метода.

Рис. 5. Привязанные фрагменты с помощью одной опорной точки для каждого фрагмента

Рис. 6. Различие между зонами затопления, полученными с космического снимка, привязанного классическим способом (а), и зонами затопления, полученными с космического снимка, обработанного

с использованием метода ПЗК (б)

ЗАКЛЮЧЕНИЕ

Совместное использование данных дистанционного зондирования Земли из космоса и геоинформационных технологий при информационной поддержке мониторинга паводковой ситуации позволяет повысить качество построения зон фактического затопления местности и точность их совмещения с векторной топографической картой за счет разработанного метода обработки радиолокационных космических снимков с использованием анализа характерных точек рельефа местности. Предложенный в работе метод был использован при построении зон фактического затопления местности на основе космических снимков Radarsat в период весеннего паводка 2008 года в окрестностях города Уфы, благодаря чему точность совмещения зон затопления с векторной топографической картой повысилась на 40% по сравнению с результатами, полученными при использованием классических методов привязки космических снимков, что существенно при оценке затопленных и подтопленных территорий специалистами, осуществляющими мониторинг весеннего паводка.

СПИСОК ЛИТЕРАТУРЫ

1. Павлов С. В., Кунаков Ю. Н. Мониторинг

паводковой ситуации с использованием данных дистанционного зондирования // Геоинформацион-ные технологии в проектировании и создании корпоративных информационных систем: Межвуз.

науч. сб. Уфа: УГАТУ, 2007. С.176-181.

2. Рис У. Г. Основы дистанционного зондирования. М.: Техносфера, 2006. 336 с.

3. Кашкин В. Б., Сухинин А. И. Дистанционное зондирование Земли из космоса. Цифровая об-

работка изображений: Учеб. пособие. М.: Логос, 2001. 264 с.

4. Decision support system for actual spring tide flooding zones estimation using images from space with the help of geoinformation systems / S. V. Pavlov [and others] // Proc. of the 10th Intern. workshop on computer science and information technologies CSIT’2008. 2008. Vol. 1. P. 134-140.

5. Сайт неформального некоммерческого сообщества специалистов в области ГИС и ДЗЗ [Электронный ресурс] (http://gis-lab.info/).

ОБ АВТОРАХ

Павлов Сергей Владимирович,

проф., зав. каф. геоинф. систем. Дипл. математик по спец. «Вычислительная математика» (БГУ, 1977). Д-р техн. наук (УГАТУ, 1998). Иссл. в обл. обработки пространств. данных.

Кунаков Юрий Николаевич,

асп. Дипл. инженер по инф. системам и технологиям (УГАТУ 2006). Иссл. в обл. прикл. применения геоинформ. систем и методов обработки данных дистанционного зондирования в науке.