Оперативность в принятии природопользовательских и природоохранных решений при освоении и эксплуатации территорий требует наращивания имеющегося вычислительного и методического аппарата. Данные цифровой радарной интерферометрической модели рельефа SRTM, полученной совместными усилиями NASA и ESA, представляют собой массив высот в метрах с пространственным разрешением до 90 м, распространяемый в виде фрагментов мозаики шириной 5×5 градусов и доступный на всю территорию планеты до 60 С.Ш. и до 54 Ю.Ш. [1, 2]. Компьютерная обработка и анализ SRTM является источником разнородной информации о территории покрытия, которая может быть извлечена из параметров рельефа, в том числе: устойчивость подстилающих пород к размыву, наличие неоднородностей и разрывов в массиве грунта, гидрологический режим рек и другие сведения.
Несмотря на доступность данных ЦМР широкому кругу исследователей, анализ эрозионной сети и определение порядков ее притоков по Стралеру [3], представляет собой нетривиальную задачу, успешное выполнение которой осложняется наличием артефактов SRTM и тем, что помимо рельефа, в качестве поверхности, отражающей радиоволны, могут выступать растительный слой и местные предметы. Целью настоящей работы является разработка альтернативного алгоритма дешифрирования эрозионной сети на основе SRTM с применением вычислительных средств с открытым исходным кодом. Для достижения поставленной цели необходимо решение задач: формирования потока обработки исходных данных; подбор необходимого математического аппарата для выполнения каждого этапа анализа, а также пробная имплементация алгоритма с разработкой программы-прототипа, выполняющей дешифрирование эрозионной сети по спутниковой ЦМР и ее экспорт в ГИС-совместимый формат.
Материалы и методы исследования
Для корректной постановки задач исследования необходимо принимать во внимание, что при описании эрозионной сети, выделяемой на спутниковой основе, речь идёт не о буквальных «водотоках», а об их морфологическом выражении в цифровой модели рельефа – тальвегах, т.е. линиях, соединяющих наиболее пониженные участки эрозионной ложбины, вне зависимости от наличия или отсутствия в ней постоянного или временного потока воды.
В качестве существующих методов детектирования объектов эрозионной сети следует отметить: ручное дешифрирование; автоматизированное распознавание с помощью модуля Spatial Analysis ArcGIS [4–6], использование малоизвестного в отечественных исследованиях набора инструментов TopoToolbox проприетарной Mathworks Matlab [7], а также недавнюю разработку исследователей из Флоридского университета, использующую библиотеки Python Numpy (бесплатная программа, поддерживается сообществом) и ArcPy (лицензия ESRI) [8]. Основным недостатков перечисленных инструментов является их привязанность к проприетарным API Matlab и ArcGIS, ограничивающими их применением конечными пользователями. В связи с этим разрабатываемый алгоритм должен адекватно дешифрировать гидросеть с учетом вышеуказанных особенностей исходных данных, разбивать детектированную сеть на притоки и определять их степень ветвления; желательно также не привязывать его имплементацию к коммерческой среде выполнения.
Сведения о распределении высот, представленные в модели SRTM, позволяют судить о морфологии отражающей поверхности и характеристиках подстилающего субстрата, обусловивших наблюдаемую картину. Эрозионное расчленение рельефа, его густота, степень и характер разветвленности свидетельствуют о постоянной или сезонной обводнённости территории, а также о ее дренировании. Число Стралера позволяет установить иерархическую соподчинённость водотоков и порядок речных [3]. Этот параметр рассчитывается следующим образом: водотоки, не имеющие притоков, имеют первый порядок, те из них, в которые впадают исключительно водотоки первого порядка, имеют второй порядок, те, в которые впадает водоток второго порядка, имеют третий порядок и т.д.
Разрабатываемый алгоритм должен обладать необходимыми свойствами, включающими: имплементацию методик, описанных в открытых источниках, а также возможность выполнения проектируемых с ним приложений в свободных средствах программирования, средах и языках, что обеспечит повторяемость результатов. Адаптируемость алгоритма должна позволить применять его для анализа территорий, характеризующихся различными особенностями распределениями высот. Кроме того, значения ошибок определения высот изображения, являющиеся особенностью SRTM [9], не должны оказывать существенного влияния на детектирование водотоков.
Существующие и описанные в литературе методы распознавания водотоков в самом общем случае сводятся к определению источников («точек начала течения»), мест впадения («точек слияния линейных водотоков друг с другом или с площадным водным объектом»). Применяется нахождение производных и градиентов поверхности рельефа, при этом водоток детектируется как область смены направления градиента поверхности рельефа φ:
где i, j – масштабные коэффициенты.
Очевидная логичность данного метода при его практической реализации сталкивается с затруднением выделения границы направления градиента рельефа в «пойме», т.е. уплощённом днище ложбины эрозионного стока в связи с незначительным относительным перепадом высот, а также наличием помех для отражения: крон деревьев, кустарников и иных объектов. Этот факт снижает эффективность методов, основанных на выявлении направлений изменения высот модели, а также выявления границ градиента, а также делает привлекательными методы, основанные на анализе геометрии и морфологии радарной интерферометрической поверхности.
Предлагаемый алгоритм, использующий методы морфологического анализа, показан на схеме (рис. 1). Он включает два цикла, условно названных «последовательное затопление» и «поиск соседей».
Рис. 1. Блок-схема предлагаемого алгоритма выделения сети водотоков
В природе затопление речной долины приводит к формированию на ее участке водной поверхности, границы распространения которой повторяют очертания изолиний рельефа. При этом осевая линия этой поверхности, равноудалённая от изолиний, в начале затопления совпадает с руслом реки, а затем, по мере сокрытия долины под водой, смещается к оси речной долины. При затоплении долин ее притоков схожая картина наблюдается и в них. При этом происходит естественное сглаживание небольших неровностей рельефа. Предлагаемая методология использует анализ эрозионной сети, выбор подмножества диапазона высот радарного изображения, морфологические операции утончения, наращивания и скелетизации. Указанные функции входят в библиотеку Scikit-image, предназначенную для анализа изображений в Python [10]. Выбор и маркирование отдельных водотоков осуществлялись с помощью операции binary labels, позволяющей присваивать уникальный идентификатор кластеру бинарного (черно-белого) изображения.
Тестовая имплементация алгоритма осуществлялась в открытой среде научного программирования Spyder на языке Python3 с использованием, помимо функций ядра Python, а также упоминавшейся выше библиотеки Scikit-image, библиотеки Numpy (численные расчёты) и gdal (Geographic data abstraction library), которые обеспечивают ввод, обработку и вывод численных данных с пространственной привязкой.
Алгоритм включает 4 функциональных блока, работа которых не требует указания пользователем каких-либо параметров, влияющих на результативность анализа, за исключением указания пути к файлу исследуемого изображения (блок 1). Действия блока 2 выполняются циклически с шагом 1/10 перепада высот до достижения максимальной отметки и включают: «затопление», т.е. выбор и создание булевой матрицы для подмножества высот ЦМР, (2.1), морфологическое истончение skimage.thin (2.2), удаление части структурного рисунка, входящего в пределы предыдущих полей «затопления» (2.3), сохранение остаточного рисунка в матрице нераспределенных водотоков (2.4). По достижении максимальной высотной отметки осуществляется переход к блоку 3, минуя коррекцию массива нераспределенных водотоков с закрытием однопиксельных разрывов фрагментов водотоков, полученных при последовательном затоплении и очисткой структурного рисунка от ошибок дешифрирование («висящих водотоков»). Блок 3 включает поиск точек бифуркации рек и «истоков» водотоков первого порядка с помощью «подсчета соседей» (3.1), морфологическое наращивание (дилатансию) точек бифуркации с помощью дискового «структурного элемента» (3.2), разъединение дешифрированного рисунка в точках бифуркации (3.3), маркирование фрагментов гидросети (3.4) и выбор тех ее фрагментов, которым принадлежат точки начала, в качестве водотоков извлекаемого текущего порядка в матрицу порядков водотоков (МПВ) (3.5). После изъятия водотоков текущего порядка матрица нераспределенных водотоков нуждается в очистке от утолщений в точках бифуркации. Для этого применяется функция 3D скелетизации skimage.skeletonize_3d (3.6). После этого текущий порядок водотоков увеличивается на 1 и происходит возврат к шагу 3.1. Цикл суперблока 3 завершается при отсутствии точек «истока», доступных для определения. Нераспределенному массиву водотоков присваивается последний достигнутый в процессе итераций порядок (число Стралера). В завершении реализующая алгоритм программа формирует и записывает в файл изображение порядка водотоков в формате с геопривязкой (GeoTiff).
Результаты исследования и их обсуждение
Применение описанного алгоритма позволило дешифрировать эрозионную сеть с использованием цифровой модели рельефа SRTM для тестового региона. Оценка качества работы модели производилась визуальным сопоставлением полученной крупномасштабной модели гидросети и фактических данных (рис. 2).
Рис. 2. ЦМР и выявленные с помощью прилагаемого алгоритма водотоки на тестовой площади (а), диаграмма распределения длин водотоков по их порядкам (числу Стралера) (б)
Результирующая матрица порядков водотоков позволяет охарактеризовать соответствие длин выделенных водотоков и их порядков. Идентифицированные точки бифуркации и «истоки» показаны в виде дополнительного растрового слоя (рис. 2, а). Они позволяют нагляднее представить работу алгоритма. Сопоставление результатов дешифрирования с фактическими данными позволяет говорить о его достаточной достоверности: точность выделения водотоков достигает 83 % процентов (в сопоставлении с глазомерным дешифрированием водотоков на том же участке). Вместе с тем используемый алгоритм не требует применения ресурсоемких коммерческих настольных продуктов, таких как ArcGIS. Интеграция алгоритма и библиотеки GDAL позволит получать в результате работы программы файл с географической привязкой, который можно включить в ГИС-проект.
Для тестового района (пересеченная горно-таежная местность, перепад высот до 997 м) было проведено сопоставление рассматриваемого алгоритма с ближайшим аналогом, дополнением Matlab TopoToolbox [7], результат работы которого показан на рис. 3. Необходимо отметить, что обоими алгоритмами, не требующими пользовательских действий при анализе, были допущены неточности в распознавании водотоков. Так, модуль TopoToolbox, несмотря на высокую точность распознавания водотоков 1-го порядка, допустил разрывы в рисунке гидросети (рис. 3, а, цифры на схеме 1 и 2), что привело к снижению максимального распознанного им порядка (числа Стралера) для эрозионной сети (распознан 4-й, в действительности 5-й порядок). Морфологический алгоритм, предлагаемый нами, установил максимальный распознанный порядок верно, однако им неверно был оконтурен один из притоков (рис. 3, б, цифра на схеме 1), а также объединены в действительности разобщенные элементы гидросети (рис. 3, б, цифра на схеме 2).
Рис. 3. Сопоставление результатов работы Matlab TopoToolbox (a) и рассматриваемого морфологического алгоритма (б)
Несмотря на то, что оба алгоритма вполне справились с поставленной задачей, для интерпретации результата дешифрирования важно верно определить максимальный порядок водотоков, что не вполне удалось сделать близкому аналогу. Выявленные несовершенства в работе предлагаемого нами алгоритма морфологического дешифрирования послужат основой для дальнейшего развития с целью создания на его основе прикладного программного обеспечения.
Заключение
Рассмотрены особенности и преимущества морфологического анализа изображений для выделения рисунка эрозионной сети и сегментации порядков ее водотоков.
Проведенное сопоставление результатов работы предлагаемого алгоритма с ближайшим аналогом, модулем Matlab TopoToolbox [7], выявило недостатки в работе обеих программ, однако разрабатываемый алгоритм корректнее установил максимальный порядок водотоков для территории. С помощью описанного морфологического алгоритма, при его дальнейшем развитии, может быть создано или модернизировано программное обеспечение для географических, геологических и природоохранных целей на основе свободно распространяемых данных ЦМР SRTM. Разработанный алгоритм, опробованный в виде консольного приложения на языке Python, может быть реализован самостоятельно или в виде дополнения настольной ГИС.
Библиографическая ссылка
Шевырёв С.Л. АЛГОРИТМ МОРФОЛОГИЧЕСКОГО АНАЛИЗА И ОПРЕДЕЛЕНИЯ ПОРЯДКОВ ВОДОТОКОВ НА ОСНОВЕ ЦИФРОВОЙ МОДЕЛИ РЕЛЬЕФА ДЛЯ РАЗРАБОТКИ ОТКРЫТЫХ ГЕОИНФОРМАЦИОННЫХ СИСТЕМ // Успехи современного естествознания. – 2018. – № 12-2. – С. 390-395;URL: https://natural-sciences.ru/ru/article/view?id=37026 (дата обращения: 19.04.2024).