The automated geographical information system for monitoring of snow cover of river watersheds based on Terra and Suomi NPP satellite data is developed for providing data for hydrological mathematical models of river runoff forecasts. System is built on standard algorithms of snow and cloud classification and include new multiday data composition algorithm based on method of cumulative sums for detecting abrupt changes in a casual process, which makes possible to achieve 99% cloud free territory accumulation and unbiased assessment of moment of snow melting.

Keywords: snow cover, image classification, river watershed

Ромасько В.Ю.(1), Бураков Д.А.(2) Для обеспечения гидролого-математических моделей прогноза стока рек создана автоматическая геоинформационная система мониторинга заснеженности речных бассейнов по данным КА Terra и Suomi NPP. Система основана на стандартных алгоритмах классификации снега и облачности, а также оригинальном алгоритме композиции данных за несколько суток на основе метода кумулятивных сумм для обнаружения разладки, с помощью которого достигается накопление 99% безоблачной территории и несмещённая оценка момента стаивания.

Ключевые слова: снежный покров, классификация изображений, речной бассейн.

Введение. Одним из важных направлений гидрологических исследований является прогнозирование уровней и расходов воды в створах бассейнов рек и водохранилищ с оборудованными гидропостами. Для сибирских рек, текущих в условиях длительной зимы и накопления снежного покрова, характерен режим преимущественно снегового и смешанного питания с максимумом годового стока в конце весны – начала лета [1].

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

Снежный покров отчётливо выделяется на снимках из космоса. Уже на первых снимках со спутника TIROS-1, последовавших после запуска в апреле 1960 г был зарегистрирован снежный покров [ 2 ]. Регулярное картографирование снежного покрова ведётся NOAA (National Oceanic and Atmospheric Administration, англ., США) c 1966 г [ 3 ]. В рамках программы EOS (Earth Observing System, англ., США) с 1998 г выпускаются ежедневные и 8-дневные карты снежного и ледового покровов и температуры льда морей в NSIDC (National Snow and Ice Data Center, англ., США) в Университете Колорадо [ 4 ].

Однако картографирование снежного покрова само по себе не даёт и не может дать никаких оценок запасов воды в снеге, ибо переход от площади к толщине и плотности снегового покрова исключительно по снимку невозможен. Тем не менее, в конце 1990-х гг. авторами была предпринята успешная попытка применить площади заснеженности, полученные по спутниковым данным, для коррекции параметров гидролого-математической модели речного стока. Суть методики, описанной в работе [5], сводится к моделированию в течение периода снеготаяния, вместе с другими параметрами, и площадей заснеженности высотных зон бассейна, которые можно сравнить с их оценкой по спутниковым данным. Сравнение модельной и фактической заснеженности в середине сезона снеготаяния может выявить расхождение, вызванное ошибками начальных снегозапасов, или коэффициентов стаивания, или других параметров модели, которое можно использовать для коррекции этих параметров модели (путём подбора), тем самым уточняя прогнозы. Ошибки в оценке начальных снегозапасов корректируются оперативно в течение сезона, а ошибки параметров модели снеготаяния корректируются единовременно – при подготовке модели по многолетнему массиву данных о заснеженности. В результате совершенствования автоматизированной системы оперативной оценки площади снегового покрытия речных бассейнов, разработанной в Красноярском НИЦ и Институте леса СО РАН в 1995-2003 гг. для вышеописанных целей [5], создана новая, полностью автоматическая «Система мониторинга заснеженности», которая обеспечивает полный цикл обработки данных радиометра MODIS (MODerate resolution Infrared Sensor, англ., США) со спутник Terra (англ.) и радиометра VIIRS (Visible InfraRed Scanner, англ., США) со спутника Suomi NPP в оперативном режиме без участия оператора, а также реализует доступ к результатам обработки из сети интернет. Система рассчитана на функционирование в автоматическом режиме в пунктах приема данных этих космических аппаратов (КА) и использованию для бассейнов больших и средних рек (более 5 000—10 000 кв. км). Конечным результатом работы программы является распределение снежного покрова в относительном выражении по высотным зонам в пределах речного бассейна, а также средней высотой кромки снега. Эти результаты используются для коррекции параметров гидролого-математических моделей для прогнозов уровней и расходов воды в реках [6]. Кроме того, данные сопровождаются обзорными снимками и картосхемами снежного покрова территории бассейнов. Система развёрнута и функционирует в Сибирском центре ФГБУ «НИЦ «Планета».

Описание методики. Выделение снежного покрова ведётся с использованием алгоритма классификации, разработанного в NASA для глобального картографирования снежного покрова по данным с радиометра MODIS [ 7 ] и аналогичного алгоритма для данных радиометра VIIRS [ 8 ]. Оба алгоритма основаны на дифференциальном нормализованном снежном индексе (NDSI) и используют специфическое свойство снега к поглощению излучения в диапазоне 1,6 мкм.

0,5 + 1,6 Здесь R0,5 и R1,6 — альбедо в диапазоне 0,5 и 1,6 мкм, соответственно (для VIIRS используется альбедо в 0,6 мкм). Эмпирическим путём, по анализу одновременных снимков высокого разрешения установлено [ 7 ], что для пикселей, покрытых снегом более, чем на 50%, значение NDSI превышает значение 0,4.

Применение NDSI для выделения снега наталкивается на трудности при наличии экранирующих помех, частично заслоняющих снежный покров, например, ветвей деревьев в лесу. Для коррекции алгоритма и возможности надёжного выделения снега на лесистых территориях в работе [ 9 ] предложено понижать пороговое значение NDSI на лесных территориях. Для этого дополнительно рассчитывается значение дифференциального нормализованного растительного индекса (NDVI):

0,5 − 1,6 1,0 − 0,6 1,0 + 0,6 Здесь R0,6 и R1,0 — альбедо в диапазоне 0,6 и 1,0 мкм, соответственно. Лесной покров идентифицируется по превышению NDVI порогового значения 0,1. В этом случае значение порога по NDSI понижается по эмпирической табличной зависимости от NDVI для MODIS [ 9 ] или проверяется на попадание в диапазон значений, рассчитываемых по двум полиномиальным зависимостям от NDVI — для VIIRS [ 8 ]. Такая модификация алгоритма позволяет обнаруживать от 50 до 98% снега на лесистых территориях, в зависимости от типа и состава леса [ 10 ]. Анализ результатов применения алгоритма на территориях бассейнов сибирских рек за период с 2004 г позволяет констатировать, что в условиях Сибири он практически каждый год обнаруживает снег на 100% территории всех модельных бассейнов.

Основным мешающим фактором для классификации снега является облачность. Из процесса классификации снега необходимо с максимальной точностью исключить районы, экранированные облачностью, для этого используется стандартные алгоритмы маскирования облачности [ 11, 12 ]. Полные версии алгоритмов маскирования облачности весьма гибкие и сложные одновременно, так как предназначен для применения на территории всего Земного шара. Учитывая высокие требования точности разделения снега и облачности для оценки заснеженности в условиях неявного наличия снега, из всех ветвей признано обоснованным применение лишь одной, а именно: облачность при наличии снега в светлое время суток. Результаты всех тестов и проверок в итоге сводятся к оценке вероятности наличия облачности. Данное значение, отличное от нуля, влечёт за собой отказ от классификации снега.

Наилучшая точность определения относительной заснеженности достигается на безоблачных снимках, что очевидно. Однако таких снимков за сезон снеготаяния попадается не много, и тем меньше, чем больше бассейн. Иногда для бассейна за всю весну с трудом можно выбрать один снимок. В существующих системах [ 4 ] для исключения влияния облачности используется миноритарная композиция данных за 8 суток в отношении снега — хотя бы один случай встречи снега за период даёт в результате снег. Однако для применения в гидрологических прогнозах, с одной стороны, шаг в 8 суток является грубым, и, с другой стороны, важна несмещённость оценки, тогда как миноритарный композит смещён в сторону переоценки наличия снега. Поэтому для данной системы был разработан алгоритм мониторинга динамики заснеженности на основе метода обнаружения разладки, известного как метод кумулятивных сумм [13].

Разладкой случайного процесса называется существенное изменение его параметров [13]. Метод кумулятивных сумм позволяет определить момент изменения среднего значения величины по серии наблюдений. Оценка момента является несмещённой, но требует знания наблюдений из серии после момента изменения, а также априорного знания среднего значения до и после изменения. В данном случае случайный процесс — это классификация снежного покрова при наличии помех в виде облачности, а существенное изменение параметров — это смена результата классификации со «снег» на «не снег» или обратно. Применительно к оценке заснеженности метод кумулятивных сумм заключается в последовательном подсчёте фактов встречи классов до достижения некоторого порога. Факт встречи одного класса увеличивает счётчик одного класса и, одновременно, сбрасывает счётчик другого. Факты встречи облачности, как и факты отсутствия данных, никак не учитываются и счёта не меняют. В момент достижения порога делается заключение о том, что на момент отрыва суммы этого класса от нуля процесс уже имел тот класс, сумма которого и превысила порог. То есть решение принимается «задним числом».

Алгоритм, по сути, является попиксельным — обработка каждого пикселя ведётся независимо от других. В описываемой системе подсчёт сумм идёт в скользящем окне длиной в 16 суток до текущей даты плюс 16 суток после, а пороговое значение равно 3 случаям встречи одного класса. Каждый раз при обработке новых данных результаты прошлых дней (до 16 суток назад) пересчитываются с учётом новых данных и, возможно, получают новые значения в соответствие с решением по методу кумулятивных сумм. Статистический анализ доли остаточной облачности применяемой схемы композиции показал близкий к экспоненциальному характер распределения. За 16 суток накапливается в среднем 99% безоблачной территории, тогда как за одни сутки только 63%, однако, в худшем случае даже за 16 суток накапливается только 65% территории.

Результаты классификации снежного покрова и облачности, суточной и многосуточной композиции данных хранятся в растровом виде в равноплощадной азимутальной проекции Ламберта на каждый бассейн отдельно. Подсчёт площадей заснеженности ведётся отдельно для каждого ландшафтно-гидрологического района каждого бассейна. Число районов и высотных зон для бассейна определяется из потребностей гидролого-математической модели на этапе её подготовки и в системе мониторинга заснеженности не ограничено. Возможности применяемых на сегодняшний день гидролого-математических моделей ограничены 7 районами и 10 высотными зонами для каждого района.

Описание реализации. Геоинформационная система мониторинга заснеженности реализована на основе архитектуры «клиент-сервер», причём клиент является веб-приложением, работающем в обычном обозревателе интернета. Серверная часть реализована как приложение-служба в среде операционной системы Microsoft Windows, работающая в связке с системой управления базами банных (СУБД) PostgreSQL и веб-сервером Apache. Хранение данных в системе реализовано на основе смешанного файлового и реляционного подхода. На долю реляционной СУБД отводится хранение параметров бассейнов и индексов файлового хранилища, а все порции каждого вида данных хранятся в файловом хранилище. Всего на текущий момент в системе рассчитывается 44 различных вида информации, начиная с обзорных повитковых изображений и заканчивая таблицами малооблачных снимков и картами многолетней статистики встречаемости классов.

В реализации серверной части используются алгоритмы с крупномодульной параллельностью по данным на основе парадигмы «портфеля задач» [ 14 ]. Отдельной задачей в «портфеле» — единицей обработки данных — является вычисление одного продукта для одного бассейна. Исполнение независимых друг от друга задач для одного или разных бассейнов производится параллельно, число параллельных вычислителей задаётся в конфигурации и обычно соответствует числу логических ядер процессора в системе. По опыту эксплуатации развёрнутой системы на обработку 3–4 витков данных спутника Terra в настоящее время тратится примерно 12–16 ядро-часов в сутки, что при 24 ядрах двух процессоров E5620 с тактовой частотой 2,4 ГГц занимает меньше 1 часа в сутки.

На текущий момент в системе мониторинга заснеженности имеется 18 модельных бассейнов рек и водохранилищ Урала, Сибири и Дальнего Востока, входящих в бассейны рек Енисея, Оби, Лены, Амура и Волги. Кроме модельных бассейнов в системе также есть 3 больших обзорных территории: в административных границах Сибирского ФО и европейской территории России, а также территория 107 бассейнов крупных (более 50 тыс. кв. км) рек России.

Результаты мониторинга в виде относительной заснеженности районов и высотных зон модельных бассейнов используются для коррекции параметров моделей в процессе весеннего снеготаяния, а именно — коррекции начальных снегозапасов по районам и высотным зонам [6]. При возникновении расхождений в расчётной и фактической наблюдаемой заснеженности в высотных зонах производится коррекция начальных снегозапасов, и расчёт повторяется. Данная итеративная процедура повторяется до достижения наименьших различий между моделью и наблюдениями. Коррекция начальных снегозапасов по этой методике в оперативной практике прогнозов всегда приводит улучшению критерия качества прогнозов, что служит наглядным подтверждением достоверности и данных о заснеженности и работоспособности всей методики в целом.

Заключение. Стандартный алгоритм классификации снежного покрова для данных MODIS и VIIRS с коррекцией порогов для лесистой местности способен выявлять 100% заснеженной территории для бассейнов сибирских рек.

Предложенный оригинальный попиксельный алгоритм композиции данных с обнаружением разладки на основе метода кумулятивных сумм вместе со стандартным алгоритмом маскирования облачности для данных MODIS и VIIRS способен обеспечить накопление 99% безоблачной территории за 16 суток, при этом опаздывая, в среднем, на 3 дня, и при этом обладает свойством несмещённости оценки.

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