UAV-based monitoring of the thermal structure of heterogeneous landscapes

Full Text

ВВЕДЕНИЕ

Температура подстилающей поверхности является важнейшим параметром окружающей среды, определяющим широкий спектр природных процессов, в том числе режим взаимодействия атмосферы с подстилающей поверхность [Li et al., 2013b]. Неоднородность распределения температуры поверхности в ряде ландшафтов связана с неоднородностью влажности (например, в болотах), отражательной способности, наклонов рельефа и отражает таким образом мозаичность распределения потоков явного и скрытого тепла. Эта неоднородность прослеживается на широком диапазоне пространственных масштабов, причем при размерах <1 км ее отдельные элементы являются подсеточными для всех современных и перспективных математических моделей погоды и климата. Существует два подхода параметризации неоднородности поверхности в прогностических моделях: метод агрегирования по параметрам и метод агрегирования по потокам (мозаичный или плиточный подход) [Avissar, Pielke, 1989; Koster, Suarez, 1992; Bartlett et al., 2002]. В гидродинамических моделях обычно используется комбинация двух методов, когда внутри ячейки сетки производится осреднение по параметрам для десятков типов растительности с целью оценки параметров одного–двух “средних” типов растительности (см., например, [Part IV: Physical Processes, 2020]), а агрегирование по потокам применяется к наиболее контрастным типам внутри ячейки: растительность, обнаженная почва, снежный покров, города, и т. д. Все предложенные варианты мозаичного метода [Arola, 1999; Molod et al., 2003, 2004; Li et al., 2013a; de Vrese et al., 2016] подразумевают наличие развитого и горизонтально статистически однородного приземного слоя над каждым типом поверхности. Но при размерах неоднородностей, сопоставимых с характерной толщиной приземного слоя при данных метеорологических условиях (менее десятков метров) условие горизонтальной однородности не выполняется, а внутренние пограничные слои могут быть сильно ниже характерных высот расположения приборов или уровней гидродинамической модели. Методы агрегирования для таких неоднородностей практически не разрабатывались. При разработке и верификации методов агрегирования для случая сильной мелкомасштабной неоднородности очень полезно использовать фактические данные о распределении температуры поверхности высокого разрешением. Это касается как поверхностей естественных равнинных ландшафтов, так и антропогенных геометрически сложных поверхностей. Детальные данные о термической неоднородности городских ландшафтов полезны для верификации и калибровки параметризаций взаимодействия атмосферы с городской средой [Тарасова и др., 2023].

Задача совершенствования методов описания термически неоднородных поверхностей в моделях атмосферы накладывает ряд требований к данным наблюдений о температуре поверхности. Пространственное разрешение и точность таких данных должны быть достаточным для идентификации основных элементов подсеточной термической неоднородности. Это означает, что дискретность таких данных по пространству должна быть минимум на порядок меньше шага сетки современных и перспективных моделей прогноза погоды и климата, который для региональных моделей уже сейчас может составлять менее 1 км. Кроме того, пространственное разрешение таких данных должно соответствовать характерному масштабу основных структурных элементов, формирующих неоднородные поверхности. Например, в случае городской поверхности, включающей здания, дороги и элементы растительности, необходимы данные с дискретностью порядка первых метров. Временная дискретность температурных данных должна разрешать суточный ход, играющий ключевую роль в процессах взаимодействия атмосферы и поверхности. В качестве минимально-достаточной дискретности по времени можно принять 3 ч – дискретность стандартных наблюдений на метеорологических станциях. При этом погрешность измерений должна быть мала по сравнению с дисперсией температуры для рассматриваемых видов подсеточной неоднородности, а также по сравнению с характерными значениями разности температуры между поверхностью и приземным слоем атмосферы, определяющей турбулентные потоки тепла. Получение соответствующих таким требованиям данных о температуре поверхности остается нетривиальной задачей.

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

Дистанционное термическое зондирование земной поверхности возможно в СВЧ- и ИК-диапазонах [Эткин и Шарков, 1976]. В ИК-диапазоне излучение формируется в очень тонком поверхностном слое и коэффициент излучения близок к единице, что позволяет достаточно точно определять поверхностную температуру. Но возможности ИК-измерений с искусственных спутников Земли существенно ограничивает влияние облачности и атмосферного аэрозоля. При этом разрешение спутниковых СВЧ радиометров, влияние атмосферы на результаты измерений которых минимизировано, не позволяет отслеживать мелкомасштабные контрасты температуры. И при СВЧ-измерениях проблемой, осложняющей их интерпретацию, является сильное влияние на наблюдаемый спектр яркостных температур вариаций коэффициента излучения, которые определяются, как правило, вариациями увлажненности грунта [Шутко, 1986].

В настоящее время в научных исследованиях широко используются данные ИК космической съемки среднего пространственного разрешения (≈1 км²), например спектрометров MODIS и Sentinel. Восстановление температуры поверхности по их данным считается достаточно надежным (±1 °C), что подтверждено сравнением с наземными измерениями [Coll et al., 2016; Wan, 2014]. Но репрезентативность таких данных ограничена для неоднородных поверхностей, поскольку они не дают представления о температуре их отдельных элементов [Yu, Ma, 2015]. Тепловые космические снимки высокого разрешения (≈100 м), например системы TIRS и ETM на спутниках Landsat, более репрезентативны для неоднородных территорий, в частности для городских ландшафтов [Weng, 2009; Ho et al., 2014; Варенцов и др., 2021]. Но даже снимки высокого разрешения не позволяют оценить термические контрасты типичных элементов неоднородного ландшафта, например отдельных зданий [Garcia-Santos et al., 2019]. Кроме того, временное разрешение спутниковых данных недостаточно для многих задач в области метеорологии и физики атмосферы. Даже в условиях безоблачной погоды спутниковые данные среднего разрешения обычно доступны 2 раза в сутки (для одной съемочной системы), а данные высокого разрешения – 1 раз в несколько дней, что не вызволяет использовать их для анализа процессов взаимодействия атмосферы с подстилающей поверхностью на внутристрочном масштабе.

Новые возможности для мониторинга, детализированного картирования и анализа термической неоднородности подстилающей поверхности открывает стремительно развитие беспилотных летательных аппаратов (БПЛА). В настоящее время исследования с использованием БПЛА набирают все большую популярность в науках о Земле, например для построения цифровых моделей рельефа, мониторинга состояния растительности [Sizov et al., 2022], измерения альбедо [Canisius et al., 2019; Ryan et al., 2017] и других задач. В области метеорологии и физики атмосферы БПЛА, оборудованные контактными датчиками, применяются для измерения метеорологических величин в пограничном слое атмосферы, в первую очередь их вертикальных профилей [Chilson et al., 2019; Kral et al., 2021; Segales et al., 2020; Varentsov et al., 2021; 2023; Репина и др., 2020]. Также развиваются методы использования БПЛА для исследования атмосферной турбулентности. При этом применяются как БПЛА самолетного типа, оборудованные многоканальными приемниками воздушных давлений [Чечин и др., 2021; Rautenberg et al., 2019], так и БПЛА мультироторного типа (квадрокоптеры), для которых возможно восстановление компонент скорости ветра и их пульсаций по данным измерений углов наклона, доступных в телеметрии БПЛА [Шелехов и др., 2021; Shelekhov et al., 2022; 2023].

Активно развивается и направление исследований, связанное со съемкой с БПЛА в инфракрасном диапазоне. Оборудованные тепловизорами БПЛА применяются для тепловой съемки различных объектов, включая геотермальные зоны, элементы флоры и фауны [Курамагомедов и др., 2015]. Тепловизионная съемка с помощью БПЛА представляет собой перспективный метод для наблюдения за водными объектами, в том числе за выходами подземных вод [Abolt et al., 2018], сельскохозяйственным орошением [Bellvert et al., 2014].

Достаточно узкие углы обзора современных тепловизоров и невысокое разрешение получаемых изображений (десятые доли мегапикселя) огранивают площадь территории, попадающий на один снимок тепловизора. Увеличение детальности получаемых данных и их пространственного охвата возможно при использовании методов аэрофотограмметрической съемки с последующим объединением множества отдельных снимков в ортомозаику. Такой подход является стандартным инструментом аэрофотосъемки в видимом диапазоне, но для тепловой съемки его применение сопряжено с большим количеством сложностей и недостаточно отлажено. В качестве примеров его применения в задачах тепловой съемки можно привести работы по картографированию температуры поверхности геотермального поля в Новой Зеландии [Nishar et al., 2016], ледника в Гималаях [Kraaijenbrink et al., 2018], урбанизированной территории [Garcia-Santos et al., 2019, Feng et al., 2019], прибрежной территории озера Севан [Medvedev et al., 2020], поля солнечных панелей в Южной Корее [Lee, Park, 2019]. К сожалению, во всех этих работах не представлены технические детали по части выполнения тепловой съемки и построения ортомозаик, что препятствует воспроизведению предлагаемых авторами методов. Кроме того, проблема термической неоднородности поверхности в этих работах рассматривается с технической и географической точек зрения, но не в контексте проблем взаимодействия атмосферы с поверхностью; при этом статистические характеристики термической неоднородности детально не анализируются.

Данная работа направлена на то, чтобы дать представление о технических нюансах задачи теплового картографирования подстилающей поверхности с использованием БПЛА, о разработанной авторским коллективом методике измерений и обработке данных, продемонстрировать результаты апробации этих методик для различных неоднородных ландшафтов, и сравнить статистические характеристики термической неоднородности этих ландшафтов.

РАЙОНЫ ИССЛЕДОВАНИЯ

Опыт мониторинга термической структуры неоднородных ландшафтов рассмотрен в статье на примере четырех полигонов, различающихся по типу земельного покрова, климатическим условиям, степени антропогенной нагрузки, типам и масштабам неоднородности подстилающей поверхности. Первый полигон представляет собой группу песчаных барханов в аридной природной зоне, расположенную в окрестностях поселка Нарын-Худук, Республика Калмыкия (45.43 N, 46.47 E, рис. 1а). Неоднородность свойств поверхности для него определяется главным образом микро- и мезорельефом (чередованием песчаных гряд с перепадов высот в несколько метров). Второй полигон представляет собой участок верхового болота в природной зоне таежных лесов на территории научно-исследовательского стационара Мухрино в Ямало-Ненецком Автономном округе (60.89 N, 68.68 E, рис. 1б). Неоднородность свойств поверхности здесь определяется микрорельефом и тесно связанными с ним условиями увлажнения в пределах характерного для верхового болота грядово-мочажинного комплекса [Kupriianova et al., 2023]. Третий полигон охватывает район субарктического города Надым в Ямало-Ненецком автономном округе (65.53 N, 72.51 E, рис. 1в) и представляет собой пример урбанизированной территории, где неоднородность поверхности определяется чередованием различных видов жилых и промышленных зданий, дорог, элементов городского озеленения, включая участок парка им. Е. К. Козлова [Варенцов и др., 2022; Sizov et al., 2022]. Четвертый полигон расположен на территории поселка Баренцбург на арктическом архипелаге Шпицберген в Норвегии (78.06 N, 14.21 E, рис. 1г). Здесь неоднородность свойств поверхности обуславливается большим набором природных и антропогенных факторов, включая расположение поселка на склоне в долине фьорда, наличие формам рельефа меньших масштабов, различные виды построек и другие антропогенно-измененные виды поверхности.

 

Рис. 1. Спутниковые снимки четырех полигонов, для которых выполнялась термическая съемка: песчаные барханы в районе пос. Нарын-Худук (а), верховое болото в районе исследовательской станции Мухрино (б), центральная часть г. Надым (в), пос. Баренцбург (г). Зеленым цветом показаны маршруты полета квадрокоптера над данными полигонами во время съемки

 

Измерения в Калмыкии выполнялись в июле 2021 г., в г. Надым – в августе 2021 г., в поселке Баренцбург – в сентябре 2021 г., в Мухрино – в июне 2022 г. Каждый из полигонов представляет собой близкую к прямоугольнику фигуру, длинна и ширина которой составляют несколько сот метров, а площадь порядка 10–20 га. Измерения для каждого из полигонов выполнялись в теплое время года в ходе полевых измерительных кампаний. На каждом из полигонов выполнялись серии измерений различной временной протяженности, в настоящем исследовании рассмотрены наиболее показательные примеры из этих серий. Для первого полигона доступны результаты измерений только для дневных условий, для остальных полигонов – для контрастных дневных и ночных условий.

МЕТОДИКА ИЗМЕРЕНИЙ И ОБРАБОТКИ ДАННЫХ

Нами использован измерительный комплекс на базе квадрокоптера DJI Mavic 2 Zoom с установленным тепловизором Flir TAU2R (рис. 2). Для установки тепловизора на квадрокоптер использован подвес производства голландской компании DroneExpert (https://dronexpert.nl/en/), стабилизирующий положение тепловизора относительно земли. Для подачи питания для тепловизора и его подвеса на корпус квадрокоптера выведен дополнительный разъем (монтаж выполнен сотрудниками компании DroneExpert). Используемая модификация тепловизора Flir TAU2R имеет размер матрицы 336 × 256 пикселей, размер пикселя 17 мкм, фокусное расстояние 9 мм, угол обзора 35°. Тепловизор работает в спектральном диапазоне 7.5–13.5 мкм.

 

Рис. 2. Фотография тепловизора FLIR Tau 2 (а) и измерительного комплекса на базе квадрокоптера DJI Mavic 2 Zoom c подвесом компании Drone Experts (б). Фото с сайта https://dronexpert.nl/en/

 

Для получения детализированных данных о пространственном распределении яркостной температуры поверхности использован метод построения ортомозаик, широко применяемый в современной геодезии и картографии. Он заключается в получении набора снимков в ходе пролета квадрокоптера над полигоном по заранее спланированному маршруту, обеспечивающему заданный уровень перекрытия снимков, и последующее объединение (сшивку) снимков в специальном программном обеспечении (ПО). Для планирования маршрута использовано ПО DJI GS Pro. Параметры полетного задания подбирались таким образом, чтобы обеспечить съемку полигона с достаточным уровнем перекрытия (не менее 50%) на одном аккумуляторе. Для корректного определения уровня перекрытия, при планировании полетного задания параметры камеры явно прописывались в ПО GS Pro. Далее квадрокоптер выполнял полет по спланированному маршруту в автоматическом режиме, без участия оператора в управлении. Управление тепловизором в нашем случае не синхронизировано с управлением квадрокоптером, поэтому перед взлетом на тепловизоре включалась серийная съёмка с максимально возможной частотой 1 раз в секунду. Для рассматриваемых полигонов съемка выполнялась с высоты от 100 до 160 м, зависящей от площади полигона. Это обеспечивало пространственное разрешение получаемых снимков в диапазоне от 15 до 30 см. Время выполнения съемки полигона во всех случаях составляло порядка 10–12 мин.

Принцип действия тепловизора основан на преобразовании инфракрасного теплового излучения в электрический сигнал, который усиливается и затем воспроизводится на ЖК-экране в виде цветной картины распределения температуры [Афонин и др., 2000; Госсорг, 1988]. Современные тепловизоры в своей работе используют матричные фотоприемники (матрицы ПЗС-датчиков), дающие сразу полное изображение температурного поля объекта. При этом следует учитывать, что тепловизионная камера измеряет не температурное поле исследуемого объекта, а яркостную температуру в ИК области, которая связана с тепловым потоком, излучаемым телом по закону Стефана-Больцмана. Реальная температура вычисляется по формуле: $T=T_{яр}/\sqrt[4]{{\epsilon }_T}$, где ε_T – интегральный коэффициент излучательной способности или черноты, который зависит от теплофизических свойств поверхности, ее геометрии, угла визирования, шероховатости и прочих сложноопределимых параметров. При этом считается, что изучаемая поверхность является “серой”, то есть коэффициент черноты постоянен во всем диапазоне волн. При тепловизионной сьемке коэффициент излучения задается пользователем в меню тепловизора, и тот факт, что различные поверхности, попадающие в поле зрения прибора, могут обладать различным коэффициентом черноты, не учитывается.

Результатом съемки тепловизионной камерой является серия снимков, каждый из которых сохраняется в отдельный графический файл в формате JPEG. Сохраненное изображение – визуализация поля автоматически рассчитанной яркостной температуры на основе предустановленной цветовой шкалы. Для получения точного значения яркостной температуры в каждом пикселе необходимо обратиться к дополнительной информации, записанной внутри того же графического файла в формате EXIF (Exchangeable Image File Format). Чтение EXIF-данных осуществлялось при помощи ПО EXIFTool и Matlab, после чего в Matlab рассчитывалась яркостная температура для каждого пикселя выбранного снимка по следующей формуле:

$T=\frac{B}{\ln \left(\frac{R_1}{R_2\left(S+O\right)+F}\right)}-273.15$, (1)

где T – яркостная температура (°C), B, R₁, R₂, O, F – параметры, характеризующие спектральные свойства излучения, S – измеренное камерой излучение, скорректированное на излучательную способность поверхности:

$S=\frac{S_0-S_{RAT}\left(1-\epsilon \right)}{\epsilon }$, (2)

где S₀ – измеренное камерой излучение, S_RAT – отраженное от объекта съемки излучение других объектов и атмосферы (RAT – Reflected Apparent Temperature), ε – коэффициент излучательной способности (коэффициент черноты) поверхности объекта съемки.

Указанные в формулах параметры для каждого снимка считываются из его EXIF-данных: S₀ и S_RAT (в EXIF-данных записана температура RAT, из которой можно получить S_RAT, подставив RAT в формулу 1 вместо T) – для каждого пикселя, B, R₁, R₂, O, F, ε – для всего снимка. При всех описываемых далее съемках коэффициент излучательной способности ε принимался равным 1.

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

Параллельно производится геопривязка инфракрасных снимков с помощью данных из полетных записей (flightlog) квадрокоптера. По времени съемки кадра определяются пространственные координаты дрона в этот момент времени: широта, долгота, высота над поверхностью земли и углы вращения. Между настройками времени квадрокоптера и инфракрасной камеры может возникать расхождение в несколько секунд, поправка на это вводится отдельно для каждой серии снимков. После определения координат проводится отсеивание лишних снимков серии – кадров, снятых во время взлета, движения к точке начала маршрута, возвращения от точки конца маршрута, посадки. В итоге остаются только снимки, сделанные дроном во время съемки полигона, и сохраняется таблица (CSV-файл), в которой для каждого снимка маршрута указываются пространственные координаты съемки, необходимые для дальнейшего объединения снимков.

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

1. Резкое (обычно на 1–3 °C) изменение средней по снимку яркостной температуры между последовательными снимками, происходящие в произвольные моменты времени (обычно несколько раз за съемку). С учетом того, что съемка велась с частотой 1 с и значительным перекрытием, можно с уверенностью говорить о том, что такие изменения связаны не с изменением температуры поверхности или внешних условий, а с особенностями работы прибора.
2. Изменение регистрируемой тепловизором яркостной температуры при изменении ориентации квадрокоптера относительно сторон света или при изменении его наклона в солнечную погоду. Это приводит к появлению систематических различий между измерениями на соседних полосах сканирования (на которых направление движения квадрокоптера различается на 180°). Подобные проблемы ранее отмечались и в других работах [Kelly et al., 2019; Malbéteau et al., 2021]. Они объясняются отсутствием у используемых моделей тепловизоров систем стабилизации внутренней температуры и чувствительностью результатов измерений к изменениям интенсивности солнечного нагрева и аспирации (что обычно происходит при изменении направления движения квадрокоптера).

Для устранения подобных артефактов предложен алгоритм коррекции, состоящий из нескольких этапов (рис. 3). На всех этапах анализируются и корректируются значения средней температуры каждого снимка T_i. На первом этапе осуществляется коррекция резких скачкообразных изменений температуры соседних снимков. Для этого вначале рассчитываются разности значений средних температур соседних снимков ΔT_i = T_i – T_i–1, и характеристики статистического распределения их значений ΔT, и на основе критерия трех сигм определялись выбросы в ряду значений ΔT:

$\left|\Delta T_i\right|>3\sigma \left(\Delta T\right)$. (3)

 

Рис. 3. Схема алгоритма коррекции. Цифры в круглых скобках советуют номеру уравнения в тексте. Пунктиром выделены этапы опциональные коррекции, необходимость которых определяется экспертным образом для каждой съемки

 

В случае выполнения такого критерия температура данного и последующего снимков корректируется следующим образом:

$T_k^{new}=T_k-∆T_i, где k\ge i$. (4)

Операция (4) далее повторятся для всех выбросов, в результате получается временной ряд скорректированных значений T^C0. За счет накопленной суммы значений ΔT_i возможно смещение среднего значения температуры временного ряда, поэтому далее необходима коррекция смещения среднего:

$T_i^{C1}=T_i^{C0}-\left(\overline{T^{C0}}-\overline{T}\right)$, (5)

где вертикальная черта обозначает осреднение по всему временному ряду. Опционально после этого также возможно удаление из ряда данных линейного тренда. Работу первого шага алгоритма коррекции иллюстрирует пример на рис. 4а.

 

Рис. 4. Пример применения алгоритма коррекции к данным съемки поверхности верхового болота в районе станции Мухрино вечером 17 июня 2022 г.: (а) временной ряд исходных значений температуры (цветные кривые на графике, цвет соответствует полосе сканирования на рис. 6) и значений, скорректированных на этапе алгоритма № 1 (сверху), временной ряд разности температуры между соседними точками (снизу); (б) идентифицированные полосы сканирования, пунктирные линии обозначают переходы между полосами; (в) зависимость температуры от положения вдоль полос сканирования и результаты ее аппроксимации локально-линейной регрессией на этапе алгоритма № 3

 

Следующий шаг коррекции направлен на устранения систематических различий между соседними полосами сканирования съемки исходя из предположения о малых различиях их температуры ввиду значительного перекрытия. Ранее в работе [Malbéteau et al., 2021] предложена идея корректировать различия средних температур соседних полос, мы же допускаем возможность непостоянства коррекции вдоль полосы. Полосы сканирования идентифицируются на основе данных о направлении движения дрона (рис. 4б). Коррекция для i-й точки задана следующим образом:

$T_i^{C2}=T_i^{C1}-\frac{1}{2w+1}\sum _{j=i-w}^{i+w}{T}_j^{C1}-\frac{1}{2}\left({T_{prev}^{C1}}_j-{T_{next}^{C1}}_j\right)$, (6)

где $T_{prev}^{C1}$ и $T_{next}^{C1}$ – средние значения температуры для n точек, входящих соответственно в предшествующую и следующую полосы сканирования, ближайших к i-ой точке (в текущем исследовании задается n = 3), оператор суммирования обеспечивает сглаживание бегущим средним c параметром w = 2 и шириной окна 2w + 1 = 5.

Третий шаг коррекции заключается в исключении систематической зависимости регистрируемой температуры от положения квадрокоптера относительно полосы сканирования:

$T_i^{C3}=T_i^{C2}-\left(f\left(x_i\right)-\overline{T^{C2}}\right)$, (7)

где f – функция, аппроксимирующая зависимость температуры от расстояния вдоль полосы сканирования полета x_i (отсчитываемого для нечетных полос от их начала, для четных – от конца), горизонтальная черта обозначает осреднение по всему ряду. В текущем исследовании в качестве функции f использована локально-линейная регрессия (модель LOWESS) с относительной шириной окна, равной 0.5. Пример такой аппроксимации показан на рис. 4в. В случае, если полученная аппроксимация хорошо описывает зависимость температуры от положения вдоль полосы, она применяется для коррекции температуры. Опыт показал, что данный шаг коррекции актуален при условии выполнения съемки в дневное время в солнечную погоду.

В результате работы алгоритма коррекции рассчитываются значения поправок для каждого снимка $c_i=T_i^C-T_i$, где T^C – ряд скорректированных значений после второго или третьего этапов коррекции. Далее эти поправки прибавляются к растрам яркостной температуры соответствующих снимков. Подразумевается, что в пределах отдельно взятого снимка значение поправки не меняется. Скорректированные растры сохраняются в виде одноканальных изображений в формате TIFF.

Построение карты температуры поверхности всего полигона производится объединением скорректированных снимков с использованием программного обеспечения Agisoft Metashape. Вначале в программу импортируются одноканальные TIFF-изображения со значениями скорректированной температуры поверхности, полученные ранее, а также CSV-файл с геопривязкой каждого изображения. В настройках программы указываются параметры инфракрасной камеры (абсолютный размер пикселя в мм, фокусное расстояние в мм) и производится трансформация канала со значениями температуры в изображение для визуализации результатов сшивания снимков. Далее производится выравнивание снимков с учетом их геопривязки – программа уточняет позиции снимков друг относительно друга, ищет на них общие точки и составляет общее разреженное трехмерное поле точек на основе всех снимков. Далее строится плотное облако точек, а на его основе трехмерная полигональная модель поверхности, необходимая для построения ортомозаики – изображения поверхности (в нашем случае температуры поверхности), сшитого из множества отдельных снимков и спроектированного на плоскую горизонтальную поверхность. Результаты промежуточных этапов сшивания снимков и построения ортомозаики представлены на рис. 5. Полученная ортомозаика в виде геопривязанного поля значений температуры поверхности экспортируется в одноканальный GeoTIFF-файл, который далее может быть, например, импортирован в какую-либо ГИС для уточнения привязки и дальнейшей обработки или визуализирован множеством различных программ.

 

Рис. 5. Результаты различных этапов сшивания изображений температуры поверхности и построения ортомозаики в Agisoft Metashape: геопривязка изображений (а), разреженное облако точек (б), плотное облако точек (в), ортомозаика (г)

 

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

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

На всех рассматриваемых полигонах массивы исходных данных содержали перечисленные выше артефакты, но их устранение оказывалось возможным средствами предложенного алгоритма коррекции. Его применение позволило устранить необъяснимые с физической точки зрения артефакты, такие как резкие скачки температуры, полосы в получаемых ортомозаиках, а также наличие и крупномасштабного градиента (плавное изменение температуры от одного края полигона к другому), что демонтирует пример для полигона на станции Мухрино (рис. 6). При этом применение алгоритма коррекции ожидаемо уменьшает характеристики статистической изменчивости температуры в пределах полигона. Для показанного на рис. 6 примера по мере применения шагов алгоритма коррекции стандартное отклонение температуры уменьшается с 3.2 до 1.4 °C, а разность между 1-м и 99-м перцентилями – с 16.1 до 7.2 °C.

 

Рис. 6. Тепловые ортомозаики, построенные по исходным данным тепловой схемки (а) и по данным после применения первого (б) и второго (в) и третьего (г) шагов алгоритма коррекции для болотного полигона на территории исследовательской станции Мухрино вечером 17 июня 2022 г. (17:30). Снизу от каждой ортомозаики приведено среднее значение (mean), стандартное отклонение (std) и разность 99-го и 1-го перцентилей (IQR)

 

Безусловно, возникает вопрос об объективной проверке адекватности алгоритма коррекции сравнением скорректированных результатов измерений с альтернативными способами оценки неоднородности температуры поверхности. На текущий момент такая задача не решена ввиду сложности и неоднозначности сравнения различных подходов измерения температуры поверхности. В перспективе для ее решения планируется привлечь данные измерений ручных ИК-радиометров, контактных термодатчиков и спутниковые тепловые снимки высокого разрешения.

Выполненные измерения позволяют проанализировать и сравнить закономерности пространственной изменчивости яркостной температуры поверхности для различных ландшафтов и в различное время суток. Для аридного полигона основной вклад в термическую неоднородность вносит чередование склонов дюн различной экспозиции. Утром (08:25 по местному времени) контрасты между освещенными (восточными) и теневыми (западными) склонами дюн достигают 10 °C (рис. 7а). Масштаб элементов неоднородности определяется размером дюн и составляет в нашем случае десятки метров. За пределами зоны барханов структура термической неоднородности резко меняется – там она определяется отдельными элементами аридной растительности и более мелкими формами рельефа. Результаты съемки в тот же день на пике дневного прогрева демонстрируют схожие закономерности (рис. 7б). В этом случае хорошо заметны различия в площади теплых и холодных элементов термической неоднородности: наиболее высокие температуры локализованы на вершинах дюн, склоны и ложбины холоднее. Статистические характеристики изменчивости температуры в пределах полигона утром и днем близки: стандартное отклонение температуры σ(T) составляет 2.0 и 2.3 °C утром и днем соответственно, разности между 99-м и 1-м перцентилями ΔP_99–1(T) составляют 10 и 11.3 °C (табл. 1).

 

Рис. 7. Тепловые ортомозаики для природных полигонов: для зоны барханов в районе поселка Нарын-Худук (Калмыкия) по данным съемки утром (а) и днем (б) 22 июля 2021 г., для верхового болота на территории исследовательской станции Мухрино по данным съемки ночью (в) и днем (г) 16–17 июня 2022 г. Снизу от каждой ортомозаики приведено среднее значение (mean), стандартное отклонение (std) и разность 99-го и 1-го перцентилей (IQR). (а) Барханы, утро (22.07.2021, 08:25), (б) Барханы, день (22.07.2021, 14:15), (в) Верховое болото, ночь (16.06.2022, 23:30), (в) Верховое болото, день (17.06.2022, 11:25)

 

Таблица 1. Статистические характеристики изменчивости яркостной температуры поверхности в пределах рассматриваемых полигонов по данным тепловой съемки с квадрокоптера

Полигон	Дата и время суток	Среднее значение, °C	Стандартное отклонение, °C	Разность 99-го и 1-го перцентилей, °C
Барханы (Нарын-Худук)	Утро (рис. 7а)	38.0	2.3	10.3
	День (рис. 7б)	55.5	2.9	11.8
Верховое болото (Мухрино)	Ночь (рис. 7в)	11.9	2.1	9.0
	День (рис. 7г)	26.8	2.1	9.1
	Вечер (рис. 6г)	25.1	1.7	7.2
Городская среда (Надым)	Ночь (рис. 8а)	9.1	2.6	12.4
	День (рис. 8б)	23.7	7.2	30.7
Городская среда (Баренцбург)	Ночь (рис. 8в)	0.5	1.7	7.2
	День (рис. 8г)	8.0	3.3	14.6

 

Для болотного полигона термическая неоднородность поверхности выражена не менее ярко, чем для барханов, при ином механизме ее формирования. Как и для барханов, ключевую роль в этом играет рельеф, но не через экспозицию склонов, а через увлажненность. Температура поверхности болота влияет на рост наземной растительности (мхи, кустарники и др.), формирование микроклимата и является ключевым фактором, контролирующим многие биотические и абиотические процессы, протекающие в почвах (торфяных, минеральных): разложение и минерализацию органического вещества почв, эмиссию парниковых газов, выделение растворенного органического углерода [Головацкая и др., 2008; Молчанов, 2015]. Режим увлажнения отдельных элементов болотного ландшафта напрямую зависит от его высоты над уровнем воды [Киселев и др., 2019]. Поэтому деятельный слой возвышенных элементов микроландшафта – болотных гряд – оказывается значительно суше и имеет меньшую теплоемкость по сравнению с более низкими участками – мочажинами и топями. Следовательно, в условиях ночного радиационного выхолаживания болотные гряды остывают значительно мочажин и топей, а днем, напротив, медленнее нагреваются. По данным выполненных измерений контрасты яркостной температуры между грядами и мочажинами ночью (23:30 по местному времени) достигают 5–10 °C, при этом наиболее теплыми оказываются участки открытой воды в топях (рис. 7в). Днем (11:25) наблюдаются сравнимые по величине контрасты противоположного знака, при этом минимумы температуры характерны для открытой воды, максимумы – на южной стороне болотных гряд, наименее затененной деревьями и кустарниками (рис. 7г). Статистические характеристики термической неоднородности полигона в ночных и дневных условиях для рассматриваемых примеров почти совпадают: σ(T) составляет 2.1 °C, ΔP_99–1(T) – около 9 °C. В переходное время суток термические контрасты между различными типами ландшафтов уменьшаются, а вместе с ними ослабевает термическая однородность территории в целом, что показывает пример для вечерних условий на рис. 6г (для этого примера σ(T) составляет 1.7 °C, ΔP_99–1(T) – 7.6 °C).

Результаты съемок для г. Надым и пос. Баренцбург позволяют получить представление о термической неоднородности антропогенно-измененных городских ландшафтов (рис. 8). Полигон в Надыме особенно показателен, поскольку позволяет сравнить городской парк и район с плотной застройкой, представленной жилыми, административными и промышленными зданиями. Парк и застроенная территория сильно различаются по характеру термической неоднородности. Ночью территория парка достаточно однородна, прослеживаются лишь слабые различия между более теплыми кронами деревьев и более холодной поверхностью. Застроенная территория более неоднородна. Во-первых, искусственные поверхности (стены и крыши зданий, заасфальтированные дороги и тротуары) в целом теплее естественных вследствие их способности первых более эффективно поглощать тепло в течение дня. При этом хорошо прослеживается эффект “тепловой след” вокруг зданий. Он обусловлен переизлучением длинноволновой радиации между стенами зданий и поверхностью земли и играет важную роль в формировании характерных особенностей городского микроклимата. Максимальная температура поверхности дорог наблюдается на узких улицах в северной части полигона, формирующих практически идеальные “городские каньоны” [Nunez, Oke, 1977]. Отдельные здания характеризуются аномально низкой яркостной температурой крыш, что, видимо, связано с низкой излучательной способностью их материалов. Например, металлические крыши могут иметь чрезвычайно низкую излучательную способность порядка 0.25 [Oke et al., 2017; Stewart et al., 2014]. В дневное время суток ключевым фактором термической неоднородности становится затенение поверхности зданиями и деревьями. Контрасты яркостной температуры между затененными и незатененными участками одной и той же поверхности могут достигать 20 °C. Из-за совместного влияния перечисленных факторов городской полигон в Надыме выделяется на фоне остальных полигонов по статистическим характеристикам изменчивости температуры. Для дневных условий эти характеристики превосходят соответствующие значения для природных полигонов более чем в 3 раза: σ(T) и ΔP_99–1(T) составляют 7.2 и 30.7 °C соответственно (табл. 1). Ночью неоднородность яркостной температуры в г. Надым также максимальна среди всех полигонов, но различия не столь велики: σ(T) и ΔP_99–1(T) составляют 2.6 и 12.4 °C. Полигон в Надыме выделяется не только на фоне других полигонов из данной работы, но и на фоне литературных данных. Например, для урбанизированной территории в Майорке разность максимальной и минимальной яркостной температуры по данным тепловой съемки с БПЛА не превышала 18 °C [Garcia-Santos et al., 2019], что меньше разности 99-го и 1-го перцентилей в Надыме.

 

Рис. 8. Тепловые ортомозаики для полигонов с антропогенно-измененной поверхностью: для центральной части г. Надым (ЯНАО) по данным съемки ночью (а) и днем (б) 11 августа 2021 г., для территории пос. Баренцбург (архипелаг Шпицберген, Норвегия) по данным съемки ночью (в) и днем (г) 10 сентября 2021 г. Снизу или сбоку от каждой ортомозаики приведено среднее значение (mean), стандартное отклонение (std) и разность 99-го и 1-го перцентилей (IQR). (а) г. Надым, ночь (11.08.2021, 01:00), (б) г. Надым, день (11.08.2021, 14:15), (в) пос. Баренцбург, ночь (10.09.2021, 23:45), (г) пос. Баренцбург, день (10.09.2021, 14:10)

 

Полигон в пос. Баренцбург сочетает виды термической неоднородности поверхности, обусловленные природными и антропогенными факторами. Ночью (23:45) наибольший вклад в наблюдаемую неоднородность вносит контраст между сушей и более теплым морем, а также рассмотренные ранее на примере Надыма факторы влияния застройки, включая “тепловые следы” вокруг зданий и различия излучательной способности искусственных материалов. Также прослеживается влияние высоты над уровнем моря. Участки естественной поверхности на удалении от зданий в западной части полигона на 2–3 °C холоднее, чем на склоне долины в западной части полигона, приподнятой над уровнем моря более чем на 100 м. Это ожидаемо, поскольку во время выполнения съемки наблюдалась ясная безветренная погода, способствующая радиационному выхолаживанию и формированию приземной инверсии температуры. В дневных условиях (14:10) одним из важных факторов термической неоднородности поверхности становится экспозиция склонов. Несмотря на разнообразие факторов, усиливающих термическую неоднородность, ее статистические характеристики уступают надымскому полигону (табл. 1), что может быть связано как с более слабой инсоляцией, так и значительно меньшей озелененностью ландшафта. Скудная растительность здесь представлена преимущественно мхами, лишайникам, карликовыми кустарничками, отдельными видами злаковых трав [Мосеев и др., 2014], при этом значительные территории в поселке и его окрестностях полностью лишены растительности.

ЗАКЛЮЧЕНИЕ И РЕКОМЕНДАЦИИ

В работе представлена разработанная авторами методика измерения температуры неоднородной подстилающей поверхности с использованием квадрокоптера и обработки получаемых данных, а также результаты ее апробации для различных неоднородных ландшафтов. Методика предполагает использование квадрокоптера DJI Mavic 2 Zoom, оборудованного тепловизором FLIR Tau 2 R, и включает следующие шаги: планирование маршрута полета над выбранным полигоном в ПО DJI GS Pro, выполнение измерений, препроцессинг тепловых снимков и коррекцию аппаратных погрешностей, построение ортомозаик в ПО Agisoft Metashape. Методика такого типа не нова, однако, в отличие от предшествующих работ, нами дано подробное описание всех ее этапов, что облегчает их воспроизведение другими исследователями. Проведенная апробация предлагаемой методики в различных условиях позволяет рекомендовать ее для построения тепловых карт поверхности, в том числе для микроклиматических исследований.

Предложен оригинальный алгоритм коррекции аппаратных погрешностей, включая резкие случайные флуктуации регистрируемой тепловизором яркостной температуры, а также различий между полосами сканирования, возникающими из-за чувствительности регистрируемых значений к инсоляции и аспирации тепловизора. Подобные артефакты возникали на всех полигонах, вне зависимости от типа подстилающей поверхности и условий съемки. Тем не менее, предложенный алгоритм коррекции во всех случаях позволил значительно снизить их влияние на результирующие характеристики термической неоднородности. Снизить зависимость результатов от используемых методов коррекции в дальнейшем можно путем выполнения рекомендаций из работы [Kelly et al., 2019], в том числе путем увеличения времени на стабилизацию работы тепловизора перед полетом и установки дополнительной защиты тепловизора от солнечного нагрева и ветра.

Важным ограничением технического характера, выявленным по результатам апробации, стали случаи присутствия в пределах съемочных полигонов крупных термически-однородных или слабоконтрастных участков, соизмеримых по размеру с пространственным охватом отдельных снимков тепловизора. Такими поверхностями могут быть водные объекты, ледники, участки открытой почвы в условиях отсутствия солнечного освещения. Для таких поверхностей сложно идентифицировать контрольные точки, необходимые для совмещения отдельных изображений. В результате совмещение изображений выполняется с ошибками (как для территории парка в г. Надым ночью, рис. 8а) или оказывается вовсе невозможным (тогда для соответствующий области ортомозаики задается идентификатор отсутствия данных или используется интерполяция, как для водной поверхности на рис. 8в, г). Решение подобной проблемы возможно путем размещения на местности искусственных контрольных точек, например алюминиевых пластин [Nishar et al., 2016].

Для четырех рассматриваемых полигонов с помощью предлагаемой методики получены данные о термической неоднородности поверхности разрешением 15–30 см (при протяженности полигонов 300–500 м по каждой из сторон). На основе этих данных впервые получены однотипные взаимно сравнимые метрики термической неоднородности поверхности для ландшафтов четырех различных типов (песчаные барханы, верховое болото, городская застройка в субарктическом городе, сочетание естественных и антропогенных ландшафтов в Арктике) в дневное и ночное время суток. Для всех рассмотренных ландшафтов выявлена высокая степень термической неоднородности. Установлено, что тепловые контрасты могут достигать первых десятков °C на площади порядка 10–20 га как на фоне дневного прогрева, так и ночного выхолаживания поверхности. Дневные значения метрики термической неоднородности, как правило, превосходят ночные, однако в отдельных случаях могут иметь близкие с ними значения (например, для верхового болота). Среди рассмотренных ландшафтов, максимальная термическая неоднородность наблюдалась в озелененной городской среде субарктического г. Надым, причем в дневное время суток ее метрики почти в 3 раза превосходят аналогичные метрики для других полигонов.

Полученные результаты в очередной раз демонстрируют важность учета подсеточной термической неоднородности в параметризациях взаимодействия атмосферы с поверхности, используемых в численных моделях прогноза погоды и климата. Во всех рассмотренных примерах пространственный масштаб отдельных элементов неоднородности не превышает первые десятки метров, то есть речь идет о заведомо подсеточных эффектах для современных и перспективных моделей атмосферы. В рамках моделей деятельного слоя суши к настоящему времени разработаны параметризации различных типов поверхности, включая городскую застройку и болота. Однако не решенным остается фундаментальный вопрос о методах агрегирования потоков с контрастных типов подстилающей поверхности для рассмотренных конфигураций ее термической неоднородности, когда пространственный масштаб ее отдельных элементов не превышает первые десятки метров. Разработанная методика мониторинга и получаемые на ее основе термические карты, в сочетании с данными пульсационных измерений потоков тепла и влаги, в перспективе могут позволить более аккуратно оценить применимость существующих методов агрегирования потоков для различных видов неоднородных ландшафтов [Lee et al., 2017]. Кроме того, термические карты могут быть использованы для задания реалистичных граничных условий в вихреразрешающих моделях атмосферы, позволяющих явно учесть влияние такой неоднородности на структуру пограничного слоя атмосферы [Глазунов, Степаненко, 2015; Варенцов и др., 2022], что также будет способствовать решению задачи проверки и улучшения методов агрегирования потоков для неоднородных ландшафтов.

Предложенная методика позволяет получать температурные данные, соответствующие большей части требований по точности и разрешению, выдвигаемых задачей совершенствования параметризаций для моделей погоды и климата сформулированных во Введении. Под вопросом остается присутствие в данных измерений систематических, постоянных по пространству погрешностей, которые могут исказить градиенты температуры в приземном слое воздуха. В дальнейших исследованиях вопрос присутствия таких погрешностей будет исследован путем сравнения результатов измерений с БПЛА с данными эталонных контактных и радиометрических измерений. Однако даже в случае наличия таких погрешностей они могут быть скорректированы при наличие эталонных данных в отдельных точках в пределах полигона. Отдельным достоинством рассматриваемой методики является возможность получения экспериментальных данных о термической неоднородности поверхности с достаточно высоким временным разрешением, ограниченным лишь размером полигона, возможностями оборудования и операторов. У авторов уже имеется опыт выполнения подобных съемок раз в 3 ч в течение нескольких суток, представляется возможным выполнение съемок и раз в 1 ч. Результаты таких измерений позволяют анализировать суточный ход характеристик термической неоднородности и учитывать его при разработке параметризацией подсеточных процессов для моделей прогноза погоды и климата. Это выгодно отличает подобные измерения на фоне данных дистанционного зондирования Земли из космоса.

БЛАГОДАРНОСТИ

Разработка методики измерений и алгоритма обработки данных выполнена при поддержке Российского научного фонда, проект № 21–17–00249. Экспедиционные работы на архипелаге Шпицберген выполнены в рамках темы госзадания 1022032200002–2–1.5.10. Экспедиционные работы на станции Мухрино выполнены при поддержке гранта Правительства Тюменской области в соответствии с программой Западно-Сибирского межрегионального научно-образовательного центра мирового уровня в рамках национального проекта “Наука”. Анализ данных измерений на станции Мухрино выполнен при поддержке Российского научного фонда, проект № 22–47–04408. Авторы благодарны Д. Г. Чечину (ИФА РАН) за содействие в проведении измерений на станции Мухрино; О. Г. Чхетинни (ИФА РАН) за содействие в проведении измерений в республике Калмыкия; П. И. Константинову (МГУ имени М. В. Ломоносова) и А. В. Соромотину (ТюмГУ) за содействие в проведении измерений в г. Надым.

About the authors

M. I. Varentsov

Lomonosov Moscow State University; Obukhov Institute of Atmospheric Physics of the Russian Academy of Sciences

Author for correspondence.
Email: mikhail.varentsov@srcc.msu.ru
Russian Federation, Moscow; Moscow

A. I. Varentsov

Lomonosov Moscow State University; Obukhov Institute of Atmospheric Physics of the Russian Academy of Sciences

Email: mikhail.varentsov@srcc.msu.ru
Russian Federation, Moscow; Moscow

I. A. Repina

Lomonosov Moscow State University; Obukhov Institute of Atmospheric Physics of the Russian Academy of Sciences; Yugra State University

Email: mikhail.varentsov@srcc.msu.ru
Russian Federation, Moscow; Moscow; Khanty-Mansiysk

A. Yu. Artamonov

Obukhov Institute of Atmospheric Physics of the Russian Academy of Sciences; Yugra State University

Email: mikhail.varentsov@srcc.msu.ru
Russian Federation, Moscow; Khanty-Mansiysk

I. D. Drozd

Lomonosov Moscow State University; Obukhov Institute of Atmospheric Physics of the Russian Academy of Sciences

Email: mikhail.varentsov@srcc.msu.ru
Russian Federation, Moscow; Moscow

A. E. Mamonotov

Obukhov Institute of Atmospheric Physics of the Russian Academy of Sciences

Email: mikhail.varentsov@srcc.msu.ru
Moscow

V. M. Stepanenko

Lomonosov Moscow State University; Obukhov Institute of Atmospheric Physics of the Russian Academy of Sciences; Yugra State University

Email: mikhail.varentsov@srcc.msu.ru
Russian Federation, Moscow; Moscow; Khanty-Mansiysk

References

Supplementary files