Космический мониторинг заснеженности территории речных бассейнов

Снежный покров отчётливо выделяется на снимках из космоса. Уже на первых снимках со спутника 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 (англ.) в оперативном режиме без участия оператора, а также реализует доступ к результатам обработки из сети Интернет. Система рассчитана на функционирование в автоматическом режиме в пунктах приема данных космического аппарата (КА) Terra и использование для бассейнов больших и средних рек (более 5 000–10 000 кв. км). Конечным результатом работы программы является распределение снежного покрова в относительном выражении по высотным зонам в речном бассейне, а также средней высотой кромки снега. Эти результаты используются для коррекции параметров гидролого-математических моделей для прогнозов уровней и расходов воды в реках [6]. Кроме того, данные сопровождаются обзорными снимками и картосхемами снежного покрова территории бассейнов. Система развёрнута и функционирует в Сибирском центре ФГБУ «НИЦ «Планета».

Описание методики

Выделение снежного покрова ведётся по алгоритму классификации, разработанному в NASA для глобального картографирования снежного покрова по данным с радиометра MODIS [7]. Данный алгоритм, основанный на дифференциальном нормализованном снежном индексе – 706 –

(NDSI, Normalized Differential Snow Index, англ.), использует специфическое свойство снега к поглощению излучения в диапазоне 1,6 мкм.

NDSI =

R4- R6

R4 + R6

Здесь R4 и R6 – альбедо в 4-м и 6-м каналах MODIS, соответственно. Эмпирическим путём по анализу одновременных снимков высокого разрешения установлено [7], что для пикселей, покрытых снегом более чем на 50 %, значение NDSI превышает значение 0,4.

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

NDVI =

R2 -R1

R2 +R1'

Здесь R1 и R2 - альбедо в 1-м и 2-м каналах MODIS соответственно. Лесной покров идентифицируется по превышению NDVI порогового значения 0,1. В этом случае значение порога по NDSI понижается по эмпирической табличной зависимости от NDVI [8] для возможности классификации снега при пониженном контрасте из-за экранирующего эффекта крон деревьев. Такая модификация алгоритма дает возможность обнаруживать от 50 до 98 % снега на лесистых территориях, в зависимости от типа и состава леса [9]. Анализ результатов применения алгоритма на территориях бассейнов сибирских рек за период с 2004 г. позволяет констатировать, что в условиях Сибири с её количеством снега он практически всегда обнаруживает снег на 100 % территории, так как для всех бассейнов существуют годы, имеющие плоский участок 100%-ной заснеженности в конце зимы – начале весны.

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

Полный алгоритм маскирования облачности весьма гибок и сложен одновременно, так как предназначен для применения на территории всего земного шара. Он содержит в себе 19 ветвей, применяемых сообразно характеру территории и времени суток, в которых выборочно используются несколько подходящих из 23 тестов и проверок. Учитывая высокие требования точности разделения снега и облачности для оценки заснеженности, из 19 ветвей признано обоснованным применение лишь одной, а именно облачность при наличии снега в светлое время суток. Данная ветвь алгоритма включает следующие тесты:

• –    тест в полосе поглощения углекислого газа на высокую облачность;

• –    тест в полосе поглощения воды на высокую облачность;

• –    тест разности температур на 11 и 3,75 мкм на туман и низкие облака;

• –    тест на высокую облачность в диапазоне 1,38 мкм;

• –    проверка на тонкие перистые облака, ближний ИК-диапазон;

• –    тест на тень от облачности;

• –    проверка по разности температур на 6,5, 11 мкм на чистое небо.

Результаты всех тестов и проверок в итоге сводятся к оценке вероятности наличия облачности в процентах. В дальнейшем при классификации снега при превышении установленного порога вероятности наличия облачности пиксели помечаются как облачные. Эксперименты с этим порогом ничего полезного не выявили, и порог был установлен в 0 %, что означает отказ от классификации снега при наличии малейшего подозрения на облачность.

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

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

Алгоритм, по сути, является попиксельным – обработка каждого пикселя ведётся независимо от других. В описываемой системе подсчёт сумм идёт в скользящем окне длиной в 16 суток до текущей даты плюс 16 суток после, а пороговое значение равно 3 случаям встречи одного класса. Каждый раз при обработке новых данных результаты прошлых дней (до 16 суток назад) пересчитывают с учётом новых данных и, возможно, получают новые значения в соответствии с решением по методу кумулятивных сумм. Помимо кумулятивных сумм для статистических

Рис. 1. Средняя и максимальная доля остаточной облачности при накоплении безоблачных данных в окне в 16 суток в бассейне Верхней Оби за октябрь 2004 – июнь 2006 г.

Fig. 1. Average and maximum share of the residual cloudiness cloudless with the accumulation of data in the window in 16 days in the Upper Ob basin during October 2004 – June 2006

оценок также за 16 суток накапливаются миноритарные композиты обоих классов, а также последнее значащее – безоблачное – значение.

Результаты статистического анализа доли остаточной облачности применяемой схемы композиции для бассейна Верховий Оби приведены на рис. 1, где показана средняя доля территории с остаточной облачностью за определённое число суток композиции, а также максимальная доля в виде Т-образных планок. В исследовании использовались данные с октября 2004 по июнь 2006 г. Характер полученного распределения близок к экспоненциальному. За 16 суток накапливается в среднем 99 % безоблачной территории, тогда как за одни сутки – только 63 %. В то же время в худшем случае даже за 16 суток накапливается только 65 % территории.

Результаты классификации снежного покрова и облачности, суточной и многосуточной композиции данных хранятся в растровом виде в равноплощадной конической проекции Альберса на каждый бассейн отдельно. В системе используются специально подготовленные покрытия бассейнов рек. Каждый бассейн должен оканчиваться створом, имеющим стационарный гидрологический пост, иначе невозможно построить прогностическую модель. Бассейн также может и начинаться со створа с гидрологическим постом. Каждый бассейн имеет собственные наиболее оптимальные параметры проекции. При их создании используются данные из различных источников:

• –    растровая цифровая модель рельефа ASTER GDEM2 [12];

• –    растровая гидрологическая модель водосборов HYDRO1K [13];

• –    гидрография векторной топографической карты масштаба 1:100 000.

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

Результаты подсчёта площадей заснеженности хранятся в виде распределений встречаемости классов по высотам, а высотные зоны модели для системы мониторинга заснеженности являются виртуальными. Площади заснеженности виртуальных высотных зон определяются «на лету» по текущим границам высотных зон районов и сохранённым в базе распределениям при запросе таблицы на сайте системы. Такой подход позволяет апостериорно менять границы высотных зон при оптимизации модели, не пересчитывая распределения.

Наряду с виртуальными высотными зонами в системе ландшафтно-гидрологические районы также виртуальны. Если при настройке модели по данным заснеженности будет выяснено, что два или более районов не отличаются друг от друга параметрами и нет никакой необходимости учитывать их отдельно, то в системе можно виртуально слить их в один. И они будут разными при подсчёте площадей, но пользователю будут показываться как один (с соответствующим пересчётом площадей, естественно) в системе мониторинга заснеженности. Данная возможность апостериорного слияния районов и изменения высотных зон без пересчёта всего многолетнего массива данных оказалась очень полезной и часто применяется на практике.

Описание реализации

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

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

На текущий момент в системе мониторинга заснеженности имеется 18 модельных бассейнов рек и водохранилищ Урала, Сибири и Дальнего Востока, входящих в бассейны рек Енисея, Оби, Лены, Амура и Волги, приведённых на рис. 2. Кроме модельных бассейнов в – 710 –

Рис. 2. Карта территорий 18 модельных бассейнов системы мониторинга заснеженности

Fig. 2. Map of the territories of 18 model basins monitoring system snow

Рис. 3. Многолетняя динамика относительной заснеженности бассейна Верхней Оби до Барнаула за 2005–2016 гг.: 1 – за 2016 г.; 2 – за 2015 г.; 3 – остальные годы; 4 – область наблюдавшихся значений

Fig. 3. Long-term dynamics of relative snow cover of the Upper Ob basin to Barnaul for 2005-2016.: 1 – for 2016; 2 – for 2015; 3 – other; 4 – area of observed values системе также есть три большие обзорные территории: в административных границах Сибирского ФО и европейской территории России, а также территория 107 бассейнов крупных (более 50 тыс. кв. км) рек России.

Все результаты обработки данных доступны через информационный сайт по сети Интернет. Сайт одновременно отображает данные одного выбранного бассейна, имеет сквозной календарь на весь период наблюдения (с 15 сентября 2004 г.), по которому легко двигаться в обе стороны с выбранным шагом до 30 суток, и позволяет составить любой набор информации для – 711 – отображения (снимки, картосхемы, таблицы). На рис. 3 приведён пример графика с информационного сайта – многолетняя динамика заснеженности бассейна Верхней Оби до Барнаула за 2005–2016 гг. Из-за большого числа кривых отдельно выделены только графики последних двух лет, а вся область наблюдаемых значений закрашена для более лёгкой общей оценки текущего года. По этому графику для текущего года можно легко качественно оценить снегозапасы на начало и сроки начала и конца периода снеготаяния, сравнить текущий год с прошлым и всем периодом наблюдений с 2005 г. Так, например, легко заметить, что вариация сроков начала схода снега в бассейне Верхней Оби превышает 1 месяц – с последней декады марта до начала мая. Также можно отметить, что весна 2016 г. была поздней, хотя её начало находилось среди среднемноголетних сроков.

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

Заключение

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

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

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