Scientific journal
Advances in current natural sciences
ISSN 1681-7494
"Перечень" ВАК
ИФ РИНЦ = 0,775

ALGORITHM OF MORPHOLOGICAL ANALYSIS AND IDENTIFICATION OF STREAM ORDER IN DIGITAL ELEVATION MODEL FOR DEVELOPMENT OF OPEN-SOURCE GEOINFORMATIC SYSTEMS

Shevyrev S.L. 1, 2
1 Far Eastern Geological Institute
2 Far Eastern Federal University
Nowadays, the freely available satellite digital relief model (DRM) for the most part of the globe makes this source of remote sensing data affordable for the great variety of applied activities and the research. Methods of remote detecting of erosion network elements and establishing of their hierarchy are applicable for assessment of the erosion dissection grade, neotectonic and landscape reconstructions as well as for resource management and various mapping. Despite of rapid development of information technologies and image analysis recently, extraction of erosion network and dividing it into the generations (orders) remains challenging task, which is complicated by artifacts in digital relief model and shortcomings in current detection algorithms or by their limitations. Errors in absolute elevations and leaving decision making onto the user responsibility also decreases reliability of available methods. Given paper considers simple and efficient algorithm of erosion network detecting and computing of Strahler number with morphological threshold-less processing in open-source scientific computing environments on a sample of Python. Output result containing spatial reference could be used for its stacking in geographic information systems (QGIS, ArcGIS etc.).
morphological image processing
geoinformatics
digital relief model
remote sensing of the Earth
Strahler numbe

Оперативность в принятии природопользовательских и природоохранных решений при освоении и эксплуатации территорий требует наращивания имеющегося вычислительного и методического аппарата. Данные цифровой радарной интерферометрической модели рельефа 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], не должны оказывать существенного влияния на детектирование водотоков.

Существующие и описанные в литературе методы распознавания водотоков в самом общем случае сводятся к определению источников («точек начала течения»), мест впадения («точек слияния линейных водотоков друг с другом или с площадным водным объектом»). Применяется нахождение производных и градиентов поверхности рельефа, при этом водоток детектируется как область смены направления градиента поверхности рельефа φ:

hevir01.wmf

где i, j – масштабные коэффициенты.

Очевидная логичность данного метода при его практической реализации сталкивается с затруднением выделения границы направления градиента рельефа в «пойме», т.е. уплощённом днище ложбины эрозионного стока в связи с незначительным относительным перепадом высот, а также наличием помех для отражения: крон деревьев, кустарников и иных объектов. Этот факт снижает эффективность методов, основанных на выявлении направлений изменения высот модели, а также выявления границ градиента, а также делает привлекательными методы, основанные на анализе геометрии и морфологии радарной интерферометрической поверхности.

Предлагаемый алгоритм, использующий методы морфологического анализа, показан на схеме (рис. 1). Он включает два цикла, условно названных «последовательное затопление» и «поиск соседей».

hevir1.tif

Рис. 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).

hevir2.tif

Рис. 2. ЦМР и выявленные с помощью прилагаемого алгоритма водотоки на тестовой площади (а), диаграмма распределения длин водотоков по их порядкам (числу Стралера) (б)

Результирующая матрица порядков водотоков позволяет охарактеризовать соответствие длин выделенных водотоков и их порядков. Идентифицированные точки бифуркации и «истоки» показаны в виде дополнительного растрового слоя (рис. 2, а). Они позволяют нагляднее представить работу алгоритма. Сопоставление результатов дешифрирования с фактическими данными позволяет говорить о его достаточной достоверности: точность выделения водотоков достигает 83 % процентов (в сопоставлении с глазомерным дешифрированием водотоков на том же участке). Вместе с тем используемый алгоритм не требует применения ресурсоемких коммерческих настольных продуктов, таких как ArcGIS. Интеграция алгоритма и библиотеки GDAL позволит получать в результате работы программы файл с географической привязкой, который можно включить в ГИС-проект.

Для тестового района (пересеченная горно-таежная местность, перепад высот до 997 м) было проведено сопоставление рассматриваемого алгоритма с ближайшим аналогом, дополнением Matlab TopoToolbox [7], результат работы которого показан на рис. 3. Необходимо отметить, что обоими алгоритмами, не требующими пользовательских действий при анализе, были допущены неточности в распознавании водотоков. Так, модуль TopoToolbox, несмотря на высокую точность распознавания водотоков 1-го порядка, допустил разрывы в рисунке гидросети (рис. 3, а, цифры на схеме 1 и 2), что привело к снижению максимального распознанного им порядка (числа Стралера) для эрозионной сети (распознан 4-й, в действительности 5-й порядок). Морфологический алгоритм, предлагаемый нами, установил максимальный распознанный порядок верно, однако им неверно был оконтурен один из притоков (рис. 3, б, цифра на схеме 1), а также объединены в действительности разобщенные элементы гидросети (рис. 3, б, цифра на схеме 2).

hevir3.tif

Рис. 3. Сопоставление результатов работы Matlab TopoToolbox (a) и рассматриваемого морфологического алгоритма (б)

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

Заключение

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

Проведенное сопоставление результатов работы предлагаемого алгоритма с ближайшим аналогом, модулем Matlab TopoToolbox [7], выявило недостатки в работе обеих программ, однако разрабатываемый алгоритм корректнее установил максимальный порядок водотоков для территории. С помощью описанного морфологического алгоритма, при его дальнейшем развитии, может быть создано или модернизировано программное обеспечение для географических, геологических и природоохранных целей на основе свободно распространяемых данных ЦМР SRTM. Разработанный алгоритм, опробованный в виде консольного приложения на языке Python, может быть реализован самостоятельно или в виде дополнения настольной ГИС.