Алгоритм вычисления границ области затопления для речной сети с моделированием распространения воды по растровому представлению рельефа

Доклады БГУИР

2016 № 5

УДК 004.94

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

А.А. ВОЛЧЕК, Д О. ПЕТРОВ, Д А. КОСТЮК

Брестский государственный технический университет Московская, 267, Брест, 224017, Беларусь

Поступила в редакцию 24 мая 2016

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

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

Введение

Проблема прогноза затопления территории в речных поймах в результате сезонных паводков, либо из-за возникновения аварийных ситуаций, является важной и решается на основе различных технологий. Существует несколько подходов к вычислению области затопления: геометрический, гидродинамический и основанный на моделировании распространения воды непосредственно по цифровой модели рельефа местности (ЦМР), представленной регулярной матрицей высотных отметок [1].

Геометрический подход предполагает создание трехмерной модели водной поверхности и последующее ее пересечение с ЦМР для определения контура границы области затопления. Наиболее известным методом расчета границ затопления территории, использующим геометрический подход, является метод створов. В методе створов трехмерная модель водной поверхности строится в виде сеточной интерполяционной поверхности между последовательностью поперечных сечений речной долины при измеренном либо вычисленном уровне подъема воды в каждом из поперечных профилей [2, 3]. К недостаткам геометрического подхода можно отнести следующие три: чрезмерное упрощение гидрологических и гидродинамических процессов; ограничение расчета области шириной поперечных сечений речной долины; нетривиальность выбора и расположения данных сечений (особенно при необходимости вычисления затопления речной системы). К достоинствам геометрического подхода относятся нетребовательность к вычислительным ресурсам, удовлетворительное качество прогнозирования области затопления при наличии густой сети гидрологических постов и при недоступности гидродинамических характеристик речной долины.

Методы вычисления области затопления, основанные на гидродинамическом подходе, построены на прямых расчетах динамики поверхностных вод (численное интегрирование уравнений Сен-Венана) [4]. Чаще всего используется пространственно-двумерная гидродинамическая модель, позволяющая вычислить высоту уровня воды и вектор скорости ее движения в узлах сети конечных элементов, покрывающих поверхность речной поймы. Гидродинамический подход обеспечивает более точное решение задач определения области затопления за счет значительной вычислительной сложности и необходимости проведения частых гидрологических изысканий для обеспечения адекватности моделирования.

Подход к вычислению области затопления, основанный на моделировании распространения воды непосредственно по ЦМР в большинстве случаев комбинирует различным образом методы геометрических и гидродинамических подходов с включением элементов теории клеточных автоматов. Методология указанного подхода в настоящее время активно развивается [5-9]. В настоящей работе авторами предлагается алгоритм моделирования паводка, специально ориентированный на речную сеть и находящийся в рамках третьего подхода.

Особенности алгоритма

Для задач моделирования затопления распределенной речной сети был разработан алгоритм, использующий математический аппарат клеточных автоматов, и принимающий в качестве входных данных объекты трех типов (рис. 1): ЦМР в виде матрицы Я [х,у],

хранящей высоту точек рельефа над уровнем моря, списка элементов растра а £ А, лежащих на осевых линиях всех водотоков речной системы, а также точек измерения уровня воды вдоль течения реки ^ е О (элемент gi хранит координаты соответствующей точки растра гх £ Я, а также высоту уровня воды И, измеренную либо полученную в результате прогноза).

Рис. 1 Представление данных в алгоритме расчета затопления речной сети

Построению области затопления предшествует подготовительный этап, состоящий из следующих шагов.

Шаг 1. Создается список Ж, содержащий водотоки речной системы, отсортированные

по возрастанию их порядков: каждый элемент wi еЖпредставляет собой список Ак,

элементы которого принадлежат осевой линии 7-го водотока. Далее выполняются следующие действия:

- для каждого водотока элементы соответствующего ему списка а е А сортируются согласно очередности прохождения ячеек по осевой линии в направлении от истока к устью (включая точку впадения в водоток более низкого порядка). Аналогично сортируются элементы ^ е О;

- список О дополняется новыми элементами, каждому из которых соответствует не менее двух элементов а: а £ £ A2w1 ^ w2 (т. е. элементы, соответствующие точкам пересечения водотоков).

Шаг 2. Для каждого элемента гх £ Я вычисляется кратчайшее Евклидово расстояние

до ближайшего элемента щ, и вычисленное значение сохраняется как элемент е е Е .

Шаг 3. Для вычисления зоны затопления необходимо:

- создать матрицу V, совпадающую по размеру с Я, но проинициализированную значением «-1» во всех ячейках.

Шаг 4. Для каждого водотока wi еЖ выполняется следующая последовательность действий:

- значения уровня воды из ^ е О присваиваются элементам V , координаты которых совпадают с координатами, хранящимися в ^ ;

- для таких соседних элементов ^ и , которым соответствуют элементы ат е А и а е Ак (т. е. точек измерения уровня воды, принадлежащих одному общему водотоку Ак ) находятся соответствующие им по положению на ЦМР элементы Я(xm,ym) и Я(xn,уп) , и интерполированные значения уровня воды присваиваются элементам матрицы V(хк,ук),

таким, что ак е Ам>,к = ш,п. Для водотоков Ак с индексом w>1 (т. е. являющихся притоками) интерполяция между заключительной парой элементов списка О не проводится, а последняя ячейка списка О заимствует значение из водотока А-1.

Шаг 5. Выполняется сравнение гх е Я с соответствующими значениями V : в случае, если не соблюдается неравенство Уху<~Уху, то значение V следует уменьшить таким

образом, чтобы оно выполнялось (это позволит нивелировать возможные погрешности ЦМР).

Шаг 6. Для определения зоны затопления, начиная от ячеек, представляющих собой осевую линию и входящих в матрицу V, необходимо произвести циклическое распространение величины вычисленной высоты уровня воды в соседние ячейки, имеющие положительное значение Евклидова расстояния от осевых линий водотоков речной системы.

Рассмотрим расчет распространения фронта воды более подробно.

1. Значения ячеек Vx > 0 копируются в множество соседних ячеек , , удовлетворяющих следующим условиям:

- находится в окрестности фон Неймана первого порядка относительно V ;

- е < ех,у, (Евклидово расстояние для ячейки больше такового для V , см. рис. 2);

- ^у < 0 (вычисление высоты воды не проводилось ранее);

- Vxy > гх, , (вычисленная высота воды превышает высоту рельефа).

Рис. 2. Схема распространения фронта воды между ячейками ЦМР по градиенту кратчайшего Евклидова расстоянию от осевой линии реки

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

полученный результат, попадая при этом в пределы зоны затопления и становясь в свою очередь источником распространения фронта водной поверхности;

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

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

Результаты и их обсуждение

Первичная проверка работоспособности алгоритма была выполнена на основе данных мониторинга паводка в ЕС, случившегося после нескольких дней проливных дождей в конце мая - начале июня 2013 года. Затопления и разрушения затронули преимущественно восток и север Германии, Чехию и Австрию, Польшу, Венгрию. Паводок распространился вниз по течению Эльбы, Дуная и в бассейне их притоков, что привело к затоплению берегов этих рек.

Благодаря наличию свободного доступа к массиву наблюдений за изменением уровня воды на гидропостах речной сети Германии, для моделирования было выбран участок р. Эльба длиной 181 км, расположенный между гидропостами Шёна, Дрезден, Майсен, Риза, Торгау, Прецш-Маукен, по состоянию на 03.06.2013. По окончании моделирования выполнялся контроль совпадения вычисленной границы с линией уреза воды, фактически зарегистрированной в трех районах (см. рис. 3). В качестве модели рельефа местности использованы данные ЦМР SRTM v4.1 (Shuttle Radar Topographic Mission) с горизонтальным разрешением до 90 м и абсолютной ошибкой определения высоты до 16 м.

Границы наблюдавшихся областей затопления были переведены из проективной системы координат EPSG:32633 (универсальная поперечная проекция Меркатора, зона 33N) в географическую систему координат на основе эллипсоида WGS-84 и растеризованы согласно пространственному разрешению элементов регулярной широтно-долготной решетки ЦМР SRTM (3x3 угловые секунды).

Рис. 3. Участок моделирования паводка на р. Эльба с расположением гидропостов; прямоугольниками отмечены районы проверки качества моделирования в окрестностях городов

Торгау (1), Риза (2) и Дрезден (3).

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

„ А п В

корректно рассчитанных областей затопления применен критерий Ъ =- [9], а для

А

А п В - В \ А

оценки адекватности расчета областей затопления - критерий г7 =- [10].

А п В + А \ В + В \ А

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

Сравнение площадей (в пикселах) рассчитанной (модель) и фактической (снимок) областей затопления для трех районов бассейна р. Эльба

Номер района А В А п В А \ В В \ А

1 972 492 448 524 44 0,46 0,4

2 1523 1124 810 314 713 0,72 0,05

3 629 842 571 271 58 0,68 0,57

Рис. 4. Сравнение наблюдавшихся и вычисленных зон затопления с отметкой расположения ближайшего гидропоста: а - Торгау, б - Риза, в - Дрезден

Согласно критерию Ъ , предложенный геометрический алгоритм расчета области

затопления показывает удовлетворительные результаты - в особенности с учетом того, что использовалась ЦМР с абсолютной ошибкой определения высоты элементов рельефа местности в 16 м [11]. Поскольку гидропосты, выполнявшие измерение высоты подъема

уровня воды, располагаются в пределах областей затопления со средней длиной до 10 км, следует предположить, что наибольший вклад в колебания величины критерия F вносят именно погрешности ЦМР.

Заключение

Разработан алгоритм для моделирования сценария разлива воды, пригодный к применению на разветвленной речной сети, отличающийся от существующих вариантов ориентацией на использование растрового представления данных и автоматическим объединением зоны затопления основного русла и его притоков в общую картину. Несмотря на крайне ограниченное разрешение использованной публично-доступной ЦМР, моделирование показывает приемлемые результаты, а именно, демонстрирует совпадение площади расчетной зоны затопления с площадью наблюдаемой зоны более чем на 62 %. Можно также отметить применимость использованной ЦМР в пилотных проектах; однако ее использование для продуктивного решения задачи прогнозирования паводков не оправдано из-за достаточно большой ошибки определения высоты [11]. Ожидается, что использование модели рельефа с вертикальным разрешением в десятки сантиметров [12] позволит существенно минимизировать погрешности расчета.

ALGORITHM FOR THE FLOOD INUNDATION ZONE BONDARIES FOR THE RIVER NETWORK BASED ON WATER SPREAD OVER THE RASTER TERRAIN MODEL

A.A. VOLCHAK, DO. PETROV, DA. KOSTIUK

Abstract

The calculation algorithm is presented for the river network inundation based on the modeling of the water front spread over the raster terrain model. A specific of its implementation is reviewed related to the data representation model, which allows consistent calculation for the river tributaries. The algorithm of primary approbation carried out on EU flooding data is discussed.

Keywords: flooding of the territory, digital elevation model, geometric calculation.

Список литературы

1. Митакович С.А. Моделирование затопления территорий и ArcGIS // Esri CIS [Электронный ресурс]. - Режим доступа: http://esri-cis.ru/blogs/? page=post&blog=arcgis&post_id= modelirovanie-zatopleniya-rek-i-arcgis. - Дата доступа : 20.04.2016.

2. Серебряков С.В., Гущин А.Н., КоршуновМ.Е. и др. // Геопрофи. 2005. № 5. С. 53-55.

3. Серебряков С.В., Гусев В.В., Зраенко Ю.Д. // Вестн. СГУГиТ. 2010. № 1 (12). С. 139-144.

4. Кобелев И.А., Писарев А.В., Храпов С. С. и др. // Матер. Междунар. конф «ИнтерКарто-ИнтерГИС 17». Белокуриха, Денпасар, 14-19 декабря 2011 г. С. 190-197.

5. Bates P.D., De Roo A.P.J. // Journal of hydrology. 2000. Vol. 236 (1). P. 54-77.

6. HorrittM.S., Bates P.D. // Hydrological processes. 2001. Vol. 15(5). P. 825-842.

7. Krupka M., Wallis S., Pender G. et. al. // Transport phenomena in hydraulics, Publications of the Institute of Geophysics, Polish Academy of Sciences. 2007. Vol. E-7 (401). P. 129-135.

8. KrupkaM., Wallis S., Pender G. et. al. // Proceedings of the congress-international association for hydraulic research. 2007. Vol. 32, № 1. P. 28.

9. Борщ С.В., Самсонов Т.Е., Симонов Ю.А. и др. // Труды Гидрометеорологического научно-исследовательского центра Российской Федерации. 2013. № 349. С. 47-62.

10. Horritt M.S. // Journal of Hydrology. 2006. Vol. 326, № 1. P. 153-165.

11. Карионов Ю.И. Оценка точности матрицы высот SRTM // Геопрофи. 2010. № 1. С. 48-51.

12. Sanders R., Shaw F., MacKay H. et. al. // Hydrology and Earth System Sciences Discussions. 2005. Vol. 9, № 4. P. 449-456.