Введение

Современные задачи обработки данных дистанционного зондирования Земли характеризуются экспоненциальным ростом объёмов пространственной информации. Согласно обзору Ma et al. [1], архивы данных NASA превышают 7,5 петабайта, ежедневно пополняясь на 4 ТБ, а глобальный объём архивных данных спутниковых наблюдений, по оценкам OGC, может превысить один ЭБ [1]. Эффективная обработка таких массивов требует специализированного программного обеспечения, способного загружать и обрабатывать данные над конкретными территориями, в том числе и речными бассейнами.

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

Для решения данной проблемы целесообразно использовать пакетную обработку и преобразование пространственных данных, которое в последние годы удобно реализовать с использованием языка программирования R.

Язык программирования R, изначально разработанный как среда для статистических вычислений, за последние два десятилетия превратился в мощный инструмент для работы с пространственными данными. Как отмечается в обзоре Hesselbarth et al. [2], проведённый опрос сообщества показал, что 90,8% исследователей в области ландшафтной экологии оценивают R как «очень полезный» или «полезный» для решения своих задач [2]. Slater et al. [3] подчёркивают, что R играет центральную роль в гидрологических науках благодаря демократизации вычислительных методов, поддержке воспроизводимости и доступу к широкому спектру статистических инструментов [3]. Ключевым преимуществом R является его открытая архитектура и возможность расширения через пакеты. Совокупность этих возможностей делает среду RStudio рекомендованным инструментом как для научных исследований, так и для практического мониторинга и управления речными бассейнами [3-5]. Bivand [6] отмечает, что более 900 пакетов в CRAN так или иначе используют классы пространственных данных, причём количество пакетов, зависящих от пакета «sp», продолжает расти [6]. Naveena et al. [7] отмечают, что рост использования R в науке о воде отражается в количестве новых пакетов, посвященных проблемам водоснабжения, и это решает задачи, связанные с осадками, засухой, эвапотранспирацией, анализом речного стока и т. д. К примеру, Kaya et al. [8] с помощью программы RStudio анализируют качество воды у плотины Мамасын (Турция). В исследовании Kobal et al. [9] демонстрируются возможности R на примере пространственного анализа лесного массива в Снежнике (Южная Словения), где показаны возможности импорта, преобразования и экспорта данных в различные форматы ГИС (геоинформационные системы), а также возможности геостатистики, пространственного моделирования и пространственной визуализации.

Цель исследования – разработать и программно реализовать автоматизированный конвейер обработки геопространственных данных на языке R для балльной геоэкологической оценки территории речного бассейна, выполнить вычислительный эксперимент на девяти бассейнах Юго-Западного Прикаспия и оценить влияние метода классификации и числа классов на интегральную карту. Практическая цель – сократить долю ручных операций при массовой подготовке согласованных растровых слоёв и обеспечить воспроизводимый пересчёт интегральной карты при смене классификационной схемы. Методика ориентирована на тематическое картографирование по бассейновому принципу; архитектура конвейера допускает адаптацию к иным полигональным объектам при сохранении схемы именования файлов.

Материалы и методы исследования

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

Геоэкологическая оценка речного бассейна представляет собой комплексный подход к изучению и анализу состояния природной и антропогенной среды в конкретном речном бассейне. Геоэкологическая оценка речного бассейна является комплексной и позволяет выявить уязвимые точки, прогнозировать возможные негативные изменения и разрабатывать стратегии управления и использования ресурсов с учетом принципов и целей устойчивого развития. Геоэкологическая оценка включает в себя оценку различных критериев и их последующее совместное воздействие [10-12]. Геоэкологической оценке территорий речных бассейнов посвящено незначительное количество работ. Практически все они выполнены по различным методикам и трудно сопоставимы между собой. Например, Карпенко Н. П. [13] исследовал бассейн реки Медвенка в Московской области с использованием балльной оценки, которая включала в себя 10 различных показателей; Сысоева П. С. [14] при оценке геоэкологической ситуации в бассейне реки Днепр при балльной оценке использует 6 показателей; Белов Н. С. в работе [15] при оценке геоэкологической ситуации в речных бассейнах Калининградской области использовал балльный подход и коэффициенты, рассчитанные по 16 показателям.

Проведение геоэкологической оценки в настоящей работе реализуется для девяти речных бассейнов Юго-Западного Прикаспия. В качестве критериев выбраны 34 фактора, суммарное воздействие которых и представляет собой результат балльной геоэкологической оценки территории речного бассейна. Методика получения данных подробно описана в [16, с. 46]. В то же время в [16] проведение геоэкологической оценки выполнялось в полуавтоматическом режиме, для территорий речных бассейнов с малой площадью и требовало больших затрат времени и ресурсов.

Геоэкологическая оценка в данной работе понимается как балльная интегральная характеристика территории водосбора по набору критериев природной и антропогенной среды. Настоящее исследование реализует исполнение геоэкологической оценки речного бассейна как воспроизводимый вычислительный конвейер. В [16] расчёт выполнялся не как полностью автономная пакетная процедура, и это требовало полного участия оператора на стадиях подготовки и классификации слоёв. Предложенная методика полностью автономна и производит вычисления на основе исходных растровых данных критериев оценивания (34 показателя) и границы водосборного бассейна. Настоящая статья не вводит новую предметную шкалу, а задаёт вычислительную реализацию: приведение слоёв к единой сетке, маскирование, перекодировку в баллы и суммирование.

Вычислительный эксперимент. Вычислительный эксперимент состоял в пакетном прогоне четырёхмодульного конвейера на полном наборе исходных данных девяти речных бассейнов Юго-Западного Прикаспия (Сунжа, Сулак, Уллучай, Карачай, Атачай, Хераз, Горган, Самур и Шураозень) и в сопоставлении интегральных растров, полученных при разных классификационных схемах. Предметная интерпретация критериев не варьировалась: менялись только вычислительные правила перекодировки и последующего суммирования. Имена растров имеют вид XX_YY_…, где XX – номер критерия, YY – номер бассейна.

Эталонной сеткой бассейна принят растр критерия 26 (26_YY_…), и к его разрешению и экстенту после перепроецирования в CRS шейпа приводятся все остальные слои того же бассейна. Перепроецирование выполняется методом ближайшего соседа. После приведения к сетке каждый растр маскируется полигоном бассейна; значения вне границы получают NA и в расчёт границ классов не входят. Входной набор содержит 34 растровых критерия на бассейн.

Пусть на территории бассейна заданы растровые критерии x1,…,xm, m=34. Обозначим множество включённых критериев через S.

Для метода M ∈ {equal, quantile} и числа классов k ∈ {4,5,6,7} слой-критерий i ∈ S переводится в баллы ci(p). Интегральный индекс в ячейке p

I(M, k, p) = c1(p) + c2(p) + … + cn(p),

где суммируются только критерии из множества S, n – число таких слоёв-критериев, ci(p) – балл i-го критерия в ячейке p. Если в ячейке значение критерия отсутствует (NA), это слагаемое пропускается. Веса включённых слоёв-критериев равны единице. Равные веса позволяют сравнивать классификационные варианты, не смешивая эффект метода разбиения с экспертным взвешиванием. В таблице 1 представлен список критериев оценивания.

Итого 8 классификационных схем на каждый включённый в обычную процедуру слой. Полный факториальный план по бассейнам имеет 72 интегральных растра. Методы fisher и kmeans в программном модуле предусмотрены, но в данной работе не представлены в результатах (в описываемый прогон не входили и в сопоставление карт не включались).

Разработанный программный конвейер реализован на языке R (версия 4.5.1) и состоит из четырёх последовательных этапов, каждый из которых отвечает за свой определённый этап обработки геопространственных данных.

Таблица 1

Перечень критериев оценивания

Выбор R обусловлен наличием высокопроизводительных пакетов для пространственных вычислений («terra»), работы с векторными данными («sf»), статистической классификации («classInt») и параллельной обработки («parallel»). Для вычислений используется высокоскоростной персональный компьютер (AMD Ryzen 9 9950X3D OEM; GIGABYTE B850 GAMING X WIFI6E; 128 ГБ оперативной памяти; твердотельный накопитель Solid State Drive (SSD) 4 TБ). Для пакета terra задавалось terraOptions(memfrac = 0.85, todisk = FALSE), то есть процедуре был доступен контролируемый объём оперативной памяти, а не весь физический объём памяти.

1. Исходными данными являются каталог с шейп-файлами бассейнов (9 речных бассейнов) и каталог с растрами (34 показателя для каждого речного бассейна). Модуль выполняет сканирование файловой системы, извлечение идентификаторов бассейнов из имён файлов с использованием регулярных выражений и построение ассоциативных словарей. Для каждого бассейна, для которого найдены и шейп, и эталонный растр, выполняется перепроецирование с помощью функции project() пакета «terra». Используется метод ближайшего соседа.

2. На втором этапе выполняется рекурсивный обход корневого каталога для выявления всех растровых файлов. Для каждого файла, удовлетворяющего формату именования, производится: загрузка растра в объект SpatRaster; извлечение соответствующего шейпа из словаря; проверка и при необходимости преобразование системы координат шейпа к CRS растра; применение операции маскирования mask(r, vect) с заменой значений вне полигона границы бассейна на NA; сохранение результата в формате GeoTIFF со сжатием LZW.

3. Третий этап заключается в проведении переклассификации растровых карт критериев и является центральным вычислительным звеном разработанного конвейера. Для каждого замаскированного растра извлекаются немаскированные значения, для каждой пары (M,k) вычисляются границы классов, строится матрица перекодировки и выполняется переклассификация. Полученные классифицированные растры сохраняются в иерархическую файловую структуру исходного вида. Обработка файлов происходит в параллельном режиме с использованием всех ядер процессора и оперативной памяти, за счет чего достигается значительное ускорение вычислений. По завершении обработки формируется сводная статистика, включающая число успешно классифицированных растров, количество ошибок и причины пропусков, что обеспечивает прозрачность и контроль качества на каждом шаге.

4. На 4-м этапе выполняется суммирование всех классифицированных растров внутри каждой конечной папки (комбинация «бассейн – метод – число классов»). Для каждой группы производится загрузка всех растров в коллекцию sds(); проверка идентичности размеров и экстентов; вычисление поэлементной суммы с игнорированием пропусков; сохранение результата с именем, содержащим все идентификационные параметры. В результате создается растр со значениями геоэкологической оценки территории речного бассейна.

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

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

В результате исследования проведена геоэкологическая оценка для девяти речных бассейнов Юго-Западного Прикаспия (Сунжа, Сулак, Шураозень, Уллучай, Самур, Карачай, Атачай, Хераз, Горган). На рисунке представлены результаты геоэкологической оценки для бассейна реки Самур (карты построены в проекции WGS-84).

Разработанный программный конвейер хорошо справляется с массовой обработкой растровых геоданных. Его главные плюсы заключаются в том, что он разделён на независимые этапы, поэтому отдельные части скриптов можно менять или заменять, не ломая всю систему обработки данных. Как видно из рисунка, бассейн реки Самур характеризуется неблагоприятным геоэкологическим состоянием в центральной своей части в районе крупных населенных пунктов Касумкент, Рутул, Ахты и др. В целом использование метода классификации quantile показывает более низкие значения неблагоприятной геоэкологической оценки, чем использование метода equal (табл. 2).

Таблица 2

Доля площади бассейна реки Самур (%), занятая каждым из классов геоэкологической оценки, выполненной по методам equal и quantile

Геоэкологическая оценка бассейна реки Самур, полученная по результатам различных методик классификации и числа классов: а) равные интервалы, 4 класса; б) равные интервалы, 5 классов; в) равные интервалы, 6 классов; г) равные интервалы, 7 классов; д) квантиль, 4 класса; е) квантиль, 5 классов; ж) квантиль, 6 классов; з) квантиль, 7 классов. Примечание: составлено авторами по результатам данного исследования

Расчет площадей выполнен с использованием перепроецирования растровых данных в метрическую систему координат (UTM 38N) и округления до одного знака после запятой, в связи с чем получена итоговая сумма, соответствующая площади бассейна.

Валидация и ограничения исследования. Способы валидации интегрального показателя в этой постановке ограничены устройством самого индекса. I (M, k, p) есть сумма баллов при равных весах и особых правилах. Для такой величины нет эталонного числового поля «истинной» геоэкологической оценки, с которым можно было бы сравнить полученную карту или долю неверно классифицированных ячеек. Поэтому классические критерии погрешности, устойчивости и сходимости, принятые для итерационных численных методов, здесь неприменимы: алгоритм конечен, не решает уравнение относительно неизвестного поля и не уточняет приближение от шага к шагу.

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

Во-первых, конвейер зависит от имен исходных файлов. Обрабатываются только растры с именем вида XX_YY_… и только те бассейны, для которых есть и шейп-файл с границами. Нарушение шаблона даёт пропуск файла, а не автоматическое сопоставление по другому признаку. В связи с этим необходимо тщательно готовить исходные данные и выполнять их проверку. Приведение разнородных исходных материалов к этому соглашению остаётся ручной операцией.

Во-вторых, веса всех входящих в сумму слоёв равны единице. Экспертное или статистическое взвешивание не проводилось, это является перспективой будущих исследований.

В-третьих, в опубликованном эксперименте использованы только методы equal и quantile. В коде конвейера предусмотрена реализация других методов классификации. Это также является перспективным направлением анализа.

В-четвертых, границы интервалов классификации вычисляются по выборке значений в пределах отдельного бассейна. Вследствие этого одному и тому же номеру класса в разных бассейнах соответствуют различные диапазоны физических величин, и интегральные оценки сопоставимы внутри бассейна, но несопоставимы между бассейнами. Для внутрибассейнового зонирования, на которое ориентирована предложенная методика, такая схема корректна и предпочтительна, поскольку обеспечивает максимальную дифференциацию в пределах анализируемой территории. Для межбассейновых сопоставлений требуется переход к единым границам интервалов, вычисляемым по объединённой выборке значений всех рассматриваемых бассейнов; такой режим работы конвейера составляет направление дальнейшей работы. Два критерия обрабатываются вне общей схемы классификации в силу своей структуры: один задан в бинарной шкале, другой изначально представлен в порядковой шкале и переклассификации не требует. Расширение множества агрегируемых критериев требует предварительного приведения их к общей порядковой шкале с согласованным направлением. Исходные критерии имеют различное пространственное разрешение и приводятся к единой сетке методом ближайшего соседа. Такой выбор сохраняет исходное множество значений и статистическую структуру растра, необходимую для корректного определения границ интервалов, но вносит геометрическое смещение до половины размера ячейки.

В-пятых, оценка является статической. Частные критерии получены из источников и характеризуют состояние территории на некоторый осреднённый момент времени. Динамика геоэкологического состояния, сезонная изменчивость и тренды в рамках предложенной модели не рассматриваются, при этом архитектура конвейера допускает построение разновременных серий оценок при наличии соответствующих данных, и это является перспективой дальнейших исследований.

Хотя язык программирования R изначально предназначался для статистических вычислений, он проложил себе путь в область ГИС благодаря разработке специализированных пакетов расширения. Язык программирования R предлагает широкий спектр функций на всех уровнях ГИС: сбор данных, обработка данных, графическое представление и количественный анализ. В данной статье R представлен как альтернатива с открытым исходным кодом существующему коммерческому программному обеспечению ГИС. Анализ более 60 тыс. статей [4] в 30 ведущих экологических журналах показал, что доля публикаций, использующих R для анализа данных, выросла линейно: с 11,4% в 2008 г. до 58% в 2017 г. Язык программирования R особенно хорошо зарекомендовал себя, когда необходимы передовые количественные методы обработки пространственных данных (например, пространственное моделирование) [9].

Заключение

В работе реализован автоматизированный конвейер балльной геоэкологической оценки территории речного бассейна на языке R в среде RStudio. Предложена формализованная модель интегральной балльной геоэкологической оценки территории речного бассейна. В отличие от словесных описаний балльного подхода, принятых в практике геоэкологических исследований, модель задаёт в явном виде оператор классификации частных критериев, правило приведения шкал к единому направлению для критериев и формирования итоговой карты, у которой рост значения соответствует ухудшению состояния территории. Разработан алгоритм многовариантной классификации разнородных растровых данных и реализующий его комплекс программ на языке R, состоящий из четырёх функционально независимых модулей: приведения растров к эталонной сетке бассейна, маскирования по векторным границам водосбора, многовариантной классификации и агрегирования с построением итоговой карты. Проведён вычислительный эксперимент на данных по девяти речным бассейнам Юго-Западного Прикаспия. Полученные растры пригодны для бассейнового тематического картографирования. Дальнейшее развитие связано с введением обоснованных весов критериев, расширением набора классификаторов и интеграцией выходных слоёв в геопортал.