Информация о форме поверхности Земли востребована при решении многих геоэкологических задач и традиционно представляется в виде карт (горизонталей, изогипс, изобат и т.п.). В настоящее время основным способом представления формы земной поверхности является цифровая модель рельефа (ЦМР). ЦМР может использоваться в качестве геопространственной основы при прогнозировании опасных процессов и явлений, для оценки скорости протекания эрозионных процессов, выделения площадей аккумуляции осадков, загрязняющих веществ и пр. В свободном доступе в сети Интернет находятся следующие глобальные ЦМР: SRTM, SRTM-30, ASTER GDEM v.2, GTOPO30, ACE2, ETOPO2 и GMTED2010. Авторами произведен детальный анализ представленных моделей по таким параметрам, как покрытие, разрешение, точность, наличие артефактов и однородность материала. Установлено, что для территории Архангельской области наиболее подходящей является ASTER GDEM v.2.
ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) – усовершенствованный спутниковый радиометр теплового излучения и отражения, установленный на борту спутника NASA Terra и запущенный в декабре 1999 г. Продукт ASTER GDEM (Global Digital Elevation Model) разработан совместно METI (Ministry of Economy, Trade, and Industry) и NASA (National Aeronautics and Space Administration). Сенсор ASTER имеет возможность стереоскопической съемки вдоль полосы пролета с помощью двух телескопов, снимающих в надир и назад в ближнем инфракрасном диапазоне с отношением база-высота (base-to-height ratio) – 0,6. Пространственное разрешение в плане – 15 м. Одна сцена ASTER в видимом или ближнем инфракрасном диапазоне имеет размер 4100 на 4200 элементов, что соответствует 60x60 км на поверхности Земли (ASTER Global DEM…, 2009). Однако в исходном виде модель содержит большое количество артефактов, пустот и иных ошибок и требует значительной корректировки и подготовки для проведения на ее основе дальнейших исследований.
Подготовка ЦМР Архангельской области включала в себя следующие этапы:
1) извлечение и объединение данных;
2) замена аномальных значений высотных отметок и «пустот»;
3) устранение шероховатостей, шума и ошибок, возникающих при наложении снимков;
4) заполнение впадин.
Все операции производились в свободном программном продукте SAGA GIS.
Архангельская область покрыта 91 сценой ASTER GDEM v2. Сцены ASTER GDEM v.2 распространяются в формате GeoTIFF. Программный продукт SAGA GIS имеет в своей структуре модуль GDAL: Import Raster, использующий открытую библиотеку преобразования растровых и векторных форматов GDAL (Geospatial Data Abstraction Library), позволяющую работать с большим количеством различных форматов, в том числе с GeoTIFF. Сцены были получены с сервера Геологической службы США и загружены в SAGA GIS с помощью библиотеки GDAL.
Следующий шаг – извлечение только тех данных, которые располагаются в границах области исследования. Граница Архангельской области в векторном формате доступна в наборе данных Vmap0. Набор данных Vmap0 (Vector Map Level 0) разработан Национальным агентством геопространственной разведки США (NGA, National Geospatial Agency) на основе набора данных DCW (Digital Chart of the World) в 1990-х гг. Векторные данные в наборе оцифрованы с карт масштаба 1:1000000 и снабжены атрибутами и текстовой информацией. Набор состоит из 10 тематических разделов по 50–70 карт: политические границы, оценка качества данных, высота, гидрография, население, транспорт, растительность и т.д. (Neteler, 2005).
Слой bnd-political-boundary-a из библиотеки Европы и Северной Азии состоит из набора полигонов, представляющих политические границы различных административно-территориальных единиц, в том числе и границы Архангельской области. Каждый полигон имеет совокупность атрибутов, в которых представлена различная информация о каждой административно-территориальной единице. Таким образом, воспользовавшись методом выборки по атрибутам (ARKHANGEL'SKAYA OBLAST'), из всего набора Vmap0 получили полигон Архангельской области, с помощью которого извлечены данные ASTER GDEM в границах области исследования (рис. 1).
Рис. 1. Граница Архангельской области в наборе данных Vmap0
Библиотека GDAL определяет значение отсутствующих данных в наборе ASTER GDEM равным -32767. Однако некоторые ошибочные данные в ASTER GDEM были удалены разработчиками и заполнены значением -9999. Для дальнейшей работы значения пустот были приведены к единому виду (-32767) с помощью модуля Grid Calculator. Данный модуль из существующих растров позволяет конструировать новый, в ячейках которого значения вычисляются из значений соответствующих ячеек входных растров по математической функции. Переменные функции – это значения входных растров, обозначенных либо в алфавитном порядке, либо с помощью буквы g и порядкового номера входного растра. Кроме того, SAGA GIS предоставляет набор математических операторов, математических функций и функций сравнения и ветвления. Таким образом, для обозначения пустот в наборе данных ASTER GDEM v.2 в модуле Grid Calculator была построена следующая формула:
Для заполнения недостающих данных в ASTER GDEM v.2 были использованы данные глобальной ЦМР GMTED2010 (рис. 2).
Рис. 2. Заполнение пустот ASTER GDEM v.2 (сверху) данными GMTED2010 (снизу)
Архангельская область покрывается одной сценой GMTED2010. Таким образом, указав в качестве параметров модуля Clip Grid With Polygon сцену GMTED2010 и полигон границы области из набора Vmap0, извлекли данные GMTED2010 в пределах Архангельской области.
ASTER GDEM v.2 имеет размер ячейки, равный 1 угловой секунде (~30 м.), размер ячейки GMTED2010 равен 7,5 угловым секундам (~225 м.). Для дальнейшей композиции этих двух моделей размер ячейки был приведен к единому виду, т.е. искусственно уменьшен размер ячейки GMTED2010 до 1 угловой секунды. С использованием модуля Patching отсутствующие данные в ASTER GDEM v.2 были заполнены значениями GMTED2010. Далее проводился статистический анализ для обнаружения аномальных значений.
По данным ASTER GDEM v.2 на территории области диапазон высот варьирует от –393 до 6255 м. Понятно, что граничные значения данного диапазона являются погрешностями и не представляют высоту естественного рельефа. Уточнить границы диапазона значений высоты можно с помощью методов обнаружения выбросов (артефактов, аномальных значений) в совокупности наблюдений. В теории статистики выбросом называется наблюдение, значение которого численно удалено от остальных значений в совокупности наблюдений. В целом наличие выбросов негативно влияет на анализ и получение корректных данных (Seo, 2006; Olewuezi, 2011).
Существует два типа методов обнаружения выбросов: формальные и неформальные тесты. Формальные тесты основываются на предположении о законе распределения и определяют выброс как отклонение от предполагаемого распределения. Неформальные тесты генерируют интервал показателей (или некий критерий), значения за пределами которого рассматриваются как выбросы. При этом неформальные тесты не позволяют сделать предположения о законе распределения (Olewuezi, 2011). Так как невозможно сделать корректное предположение о законе распределения высоты рельефа земной поверхности, следует воспользоваться одним из неформальных методов: методом среднеквадратического отклонения (3SD); методом Z-Score; модифицированным методом Z-Score (Modified Z-Score); методом абсолютного медианного отклонения (3MADe); правилом медианы (Median Rule).
Для определения аномальных значений в ЦМР области был написан сценарий на языке программирования Python, который позволяет рассчитывать вероятные диапазоны значений (и отбрасывать выбросы) по всем из пяти рассмотренных методов. В качестве входящего параметра сценарий принимает статистические данные растрового слоя (модуль Zonal Grid Statistics SAGA GIS), выводит в упорядоченном виде таблицу, в одном столбце которой находятся значения растрового слоя, а в другом – количество их вхождений.
Неформальные методы дают приблизительно одинаковый результат на относительно небольших выборках. Однако модифицированный метод Z-Score, метод абсолютного медианного отклонения и правило медианы, т.е. методы, которые не основываются на среднем арифметическом и стандартном отклонении, более устойчивы к наличию экстремальных выбросов в больших наборах данных, поэтому лучше использовать их. Поскольку метод абсолютного медианного отклонения работает быстрее на больших наборах данных, то именно он был применен для обнаружения выбросов в подготавливаемой ЦМР.
Для устранения шероховатостей, шума и ошибок наложения снимков в SAGA GIS использовался модуль DTM Filter. Работа модуля основывается на предположении о том, что большая разность высот между двумя соседними ячейками вряд ли является крутым склоном, скорее ошибкой в данных. Вероятность того, что ячейка с большими значениями высот не является ошибочной, увеличивается пропорционально расстоянию между двумя рассматриваемыми ячейками, поэтому допустимая разность высот рассчитывается с помощью функции от этого расстояния. Рассматриваемая ячейка считается верной, если в пределах плавающего окна не встретится другая ячейка, разность высоты которой с высотой рассматриваемой ячейки не будет превышать допустимую разность (Vosselman, 2000).
Одним из показателей точности ASTER GDEM v.2 является размер стэка, т.е. количество стереопар, использованных для расчета значения высоты в конкретном пикселе. Размер стэка представлен в файле оценки качества (QA-файл) для каждой сцены ASTER GDEM v.2. В тех областях, где размер стэка мал, наиболее вероятны ошибки и неточности в данных. Также вероятны ошибки на линейных или криволинейных границах между областями с разным размером стэка. Такие ошибки принимают форму прямых линий, бугров, ям и т.д. С помощью модуля DTM Filter были устранены погрешности в данных на границах между областями с разным размером стэка и в областях, где размер стэка мал.
После выполнения всех процессов, связанных с устранением ошибок, образованные на месте неверных значений пустоты были также заменены данными GMTED2010.
Для сглаживания ЦМР применен фильтр в модуле Simple Filter, основанный на среднеквадратическом отклонении. Работа фильтра заключается в том, что вокруг каждой ячейки модели строится окно определенного размера и для полученного окна высчитывается среднеквадратическое отклонение. Далее всем ячейкам окна, значение которых превышает двойное среднеквадратическое отклонение (по модулю), задается новое значение, равное среднеквадратическому отклонению окна.
На последнем этапе подготовки модели необходимо произвести ее гидрологическую коррекцию. В большинстве ЦМР содержится множество понижений рельефа, выраженных в виде площадей без стока и часто называемых впадинами (Zhu et al., 2013). Определение и удаление поверхностных впадин – это важный этап подготовки цифровой модели рельефа для дальнейшего гидрологического моделирования и анализа на ее основе. Впадиной называют локальный минимум, не имеющий стока вниз по склону ни к одной из смежных ячеек цифровой модели. Впадина может состоять из одной или из группы смежных ячеек с одинаковым значением высоты, которое меньше значения окружающих ячеек. В качестве впадин могут выступать как элементы естественного рельефа, так и ложные артефакты. Ложные впадины – это недостатки цифровой модели, их наличие обусловлено ошибками входных данных, ошибками интерполяции в процессе генерации цифровой модели, округления интерполированных значений, усреднения значений высоты (Wang, Liu, 2006).
Существует достаточное количество автоматических методов устранения впадин из ЦМР. Одни методы основаны на простом последовательном увеличении значений высоты в ячейках впадины до тех пор, пока не станет возможным сток вниз по склону с тем ограничением, что сток не может вернуться обратно в ячейку впадины. Другие методы предлагают простое размытие значений цифровой модели рельефа для уменьшения количества и размера впадин (Garbrecht, Martz, 2000). Эффективные методы устранения впадин разработаны авторами L. Wang и H. Liu (Wang, Liu, 2006), и O. Planchon и F. Darboux (Planchon, Darboux, 2002). Оба этих метода реализованы в SAGA GIS в виде соответствующих модулей Fill Sinks (Planchon/Darboux, 2002) и Fill Sinks (Wang/Liu).
Метод (Planchon, Darboux, 2002) заключается в том, что вместо постепенного заполнения впадин в первую очередь рельеф заполняется «толстым слоем воды», а затем удаляется избыток. В основе метода (Wang, Liu, 2006) лежат такое понятие, как «высота заполнения», и алгоритм поиска направления потоков с наименьшими затратами. Данный метод способен за один проход одновременно определить как направления потоков, так и пространственное разделение водоразделов (рис. 3). Для заполнения впадин в ЦМР Архангельской области был использован метод Wang/Liu как наиболее быстро выполняющий процесс обработки больших наборов данных по сравнению с другими методами.
Рис. 3. Иллюстрация работы метода L. Wang и H. Liu (2006):
а) Профиль рельефа с выходом потока в точке A;
б) высота заполнения в ячейках между точками A и C равна исходной высоте для каждой ячейки;
в) высота заполнения между точками C и I равна исходной высоте точки C;
г) высота заполнения в ячейках между точками I и H равна исходной высоте
В результате проведенной коррекции получена подготовленная для дальнейших исследований ЦМР Архангельской области. Результат сравнения исходной и подготовленной модели приведен на рисунке 4.
Рис. 4. Исходная (слева) модель ASTER GDEM v.2 и подготовленная ЦМР Архангельской области
Исследование выполнено при финансовой поддержке ФАНО России в рамках проекта № 0410-2014-0024 «Разработка комплексной физико-геоэкологической количественной модели взаимодействия (литосфера, гидросфера, биосфера, атмосфера и, частично, ионосфера) в районах тектонических узлов севера Русской плиты и оценка их влияния на окружающую среду»; программы УрО РАН в рамках темы 0410-2015-0016 «Отражение проявлений кимберлитового магматизма и зон глубинного нефтегазообразования в современном геодинамическом режиме Арктического сегмента земной коры».
Рецензенты:
Беленович Т.Я., д.г.-м.н., старший научный сотрудник лаборатории глубинного геологического строения и динамики литосферы Института экологических проблем Севера УрО РАН, г. Архангельск;
Малов А.И., д.г.-м.н., заведующий отделом по Наукам о Земле Института экологических проблем Севера УрО РАН, г. Архангельск.