Применение спутниковой радарной интерферометрии для изучения и моделирования полей смещений на склонах вулканов полуострова Камчатка

Глава 1 посвящена первому защищаемому положению: Технология обработки РСА
снимков, эффективность которой в условиях полуострова Камчатка (низкая когерентность
радарных снимков, горный рельеф, существенное влияние атмосферы) обоснована путем
массового тестирования на снимках C и L диапазона. Снимки, выполненные спутниковыми
радарами с синтезированной апертурой (РСА), позволяют оценить смещения земной
поверхности за период съемки. Эти результаты дают возможность исследовать вулканические и
сейсмические процессы. В п. 1.1 приводятся краткие сведения из теории радарной
интерферометрии. В основе интерферометрических методов лежит технология обработки пары
снимков, так называемая дифференциальная интерферометрия DInSAR. Ключевым моментом
является то, что повторная съёмка должна проводиться с относительно близких орбит в случае
одного спутника или в случае группировки спутников. Отраженные сигналы, полученные в
двух точках последовательной съёмки от одного и того же отражающего объекта, приходят с
разной фазой. В общем случае, после вычитания фазовой компоненты, содержащей
информацию о топографии (a priori известная цифровая модель рельефа – ЦМР), остаётся
разность фаз двух отражённых от земной поверхности сигналов, принятых в близко
расположенных точках съёмки. По разности фаз можно определить смещение за время между
съёмками. Важный этап интерферометрической обработки состоит в удалении фазовых
компонент, связанных с шумами, влиянием атмосферы, и неточным знанием орбиты и ЦМР.
На применении технологии DInSAR основаны более сложные технологии оценки
смещений по серии радарных снимков, например, метод малых базовых линий SBAS (Small
BAseline Subsets) – наблюдение динамики смещений во времени по серии снимков. При работе
со множеством радарных изображений, можно устранять некоррелируемые во времени фазовые
сигналы, которые в основном связаны с влиянием атмосферы, определять коррекции высоты к
ЦМР. Применение радарной съёмки для изучения и мониторинга вулканических районов
Камчатки осложняется рядом обстоятельств. Это горный рельеф, наличие мощного снежного
покрова зимой и густой растительности летом, сложные погодные условия и т.д.
Для Камчатки в основном мы используем снимки спутников Sentinel-1A (C-диапазон),
находящиеся в открытом доступе и снимки L-диапазона спутников ALOS-1 и ALOS-2 (радары
PALSAR-1, PALSAR-2 соответственно), полученные в рамках научного проекта ER2A2N075,
поддержанного Японским аэрокосмическим агентством JAXA. Исходя из характеристик
радаров и накопленного опыта работы в различных районах Камчатки, можно сделать
следующие выводы об использовании радарных снимков и применяемых к ним методов
обработки. Для определения смещений единичных событий и/или там, где невозможно
получить серию из достаточного количества снимков, рекомендуется применять технологию
DInSAR на основе пары радарных изображений дециметрового спутника ALOS-1 или ALOS-2.
В ряде случаев удается получить интерферограммы по парам снимков спутника Sentinel-1, но
это сильно зависит от благоприятного совпадения условий, влияющих на когерентность
снимков или на качество сигнала, связанного с атмосферными эффектами. Примером работы в
благоприятных условиях являются, например, наши исследования землетрясения под озером
Хубсугул [Тимошкина и др., 2022]. В большинстве случаев предпочтительнее использовать
подход SBAS для наблюдения динамики смещений во времени и получения более точного
результата. Технология SBAS применяется, если серия составляет более 12 снимков, что
возможно при использовании радарных изображений спутника Sentinel-1, работающего в
сантиметровом диапазоне длин волн.
Процесс обработки РСА снимков по существу состоит в выделении полезного
деформационного сигнала на фоне различных помех. Для сложных, с точки зрения
интерферометрии, районов, приходится проводить большие серии пробных расчетов для
выработки наилучшей последовательности действий, позволяющих эффективно устранять
помехи и региональный набег фазы, проводить развертку фазы в условиях низкой
когерентности и т.д. В п. 1.4 главы 1 приведено описание технологии оценки полей смещений
методом SBAS, которая была разработана и обоснована численными расчетами в процессе
исследований и рекомендуется для применения в северных и/или горных районах.
Для обработки серии снимков в пакете SARscape сложность выбора параметров и
настроек алгоритмов заключается в том, что при расчете серии интерферограмм, их развёртке,
последующем выравнивании и уточнении результатов, все параметры задаются одновременно
для всего стека данных. В связи с этим рекомендуется предварительно по некоторым парам
отдельно выполнить дифференциальную интерферометрическую обработку для того, чтобы
подобрать наиболее подходящие общие параметры.
При выборе области исследования важно, чтобы ее границы выходили за пределы
исследуемого горного массива, что актуально в условиях Камчатского региона. В дальнейшем
это скажется при определении местоположения контрольных точек (Ground Control Points).
В зависимости от целей применения метода SBAS необходимо учитывать следующие
ограничения на временную и пространственную базовые линии. В случае построения ЦМР
интерферограммы с очень малыми значениями перпендикулярной базовой линии (< 30 м) оказываются практически бесполезны из-за их высокой чувствительности к фазовому шуму и атмосферным эффектам. Для расчёта смещений, напротив, значение пространственной базовой линии должно быть как можно меньше. Поскольку для Камчатки редко удается собрать большое количество снимков с хорошей когерентностью, максимальная временная базовая линия (т.е. временной интервал между снимками) должна покрывать весь период съёмки, чтобы получить максимальное количество интерферограмм. Как правило, серия состоит из 12–14 снимков, с интервалом выполнения съёмки спутником Sentinel-1A 12 дней, следовательно, пороговое значение временной базовой линии должно превышать 168 дней, чтобы гарантированно учесть все возможные пары. Важным моментом на этапе построения интерферограмм является выполнение операции осреднения по пространству (multilooking). Целесообразно увеличивать размер ячейки разрешения как минимум до одинакового размера вдоль и поперек орбиты. Коэффициенты осреднения задают больше 1, чтобы получить интерферограммы с более высоким соотношением сигнал/шум. Осреднение также значительно уменьшает время работы алгоритмов. При расчёте интерферограмм рекомендуется использовать опцию "Coregistration with DEM" – корегистрация (пространственное совмещение) снимков с учетом ЦМР. Основная цель опции состоит в том, чтобы улучшить корегистрацию между основным снимком (master) и его парой (slave). Это важно, например, в областях с крутым горным рельефом и при обработке данных, полученных в высоких широтах. Также необходимо знание опорной ЦМР. Главные характеристики ЦМР это точность задания высот и пространственное разрешение. Предпочтительно использовать ЦМР с разрешением 1 угловая секунда: SRTM1 и AW3D30. Для обоснованного выбора ЦМР необходимо проводить сравнительный анализ имеющихся моделей на предмет наличия пустот, как в SRTM1 в некоторых горных районах (например, на вулканах Толбачик, Безымянный и некоторых др.) или наличия аретефактов с разницей в высотах между SRTM1 и AW3D30 до 1000 метров, например, в районе вулкана Мутновский. Так как рассчитанные интерферограммы содержат значения, свёрнутые по модулю 2π (относительные значения фазы), возникает нетривиальная задача восстановления полного цикла периодов (абсолютного значения фазы), называемая разверткой фазы. При проведении операции развёртки фазы рекомендуется включать в рассмотрение только те пиксели, в которых значения когерентности выше заданного порога. Тем самым исключаются пиксели, сильно подверженные воздействию помех и, следовательно, не несущие полезной информации о фазе. Эти пиксели ухудшают результаты развёртки. При фазовой развёртке обычно хорошо работает метод MCF (Minimum Cost Flow), но в сложных условиях, лучше всего выбрать метод Delaunay MCF. Этот метод обеспечивает большую надёжность в областях с низкой когерентностью и преимущественно не допускает ошибок фазовой развёртки (по сравнению с MCF). Здесь следует уточнить, что ошибкой фазовой развертки считается наличие скачков в фазовом поле более чем на π радиан в соседних пикселях. В ряде случаев на интерферограммах появляется региональная компонента, называемая набег фазы. Он может возникнуть из-за неточного задания орбиты или атмосферных эффектов. Устранение набега фазы возможно путем задания контрольных точек, т.е. точек, в которых задаются значения развернутой фазы, например, нулевые значения в областях, где смещения поверхности не ожидаются. Для расстановки контрольных точек с целью уточнения развёрнутых интерферограмм и выравнивания фазового уклона (если он имеется) существуют общие правила. Главное, рекомендуется выбирать ровные, стабильные зоны, с хорошей когерентностью, подальше от крутого рельефа и возможных областей существенных смещений. Однако в реальной ситуации на Камчатке зачастую наилучшую когерентность имеют скальные выступы, горные хребты, лавовые поля, т.е. как раз те области, которых надо избегать при выборе местоположения контрольных точек. Хорошие результаты при расчёте смещений на вулкане Большая Удина [Сенюков и др., 2020] были нами получены расстановкой более ста контрольных точек по периферии области исследования (не затрагивая склоны горных массивов) и без учета когерентности. После вычисления развёрнутой фазы создаётся последовательность карт смещений во времени и формируются временные ряды фазовых смещений для каждого пикселя. Далее начинается работа с временными рядами. При расчёте поправок к заданной опорной ЦМР (коррекция высот) используется следующий факт. Если имеется ошибка в задании ЦМР в данном пикселе, то на парной интерферограмме она породит фазовый сдвиг, величина которого будет прямо пропорциональна величине пространственной базовой линии двух снимков, по которым эта интерферограмма построена. Поэтому поправка к ЦМР определяется путем оценки линейной корреляции фазовых смещений и величины базовой линии снимков. Одновременно может быть задана модель зависимости смещений от времени, например, линейная. Атмосферные помехи и топографические артефакты, как правило, эффективно удаляются с помощью встроенных в пакет SARscape атмосферных и пространственных фильтров на последующих шагах обработки. Чем меньше размер окна атмосферного фильтра, тем фильтрация «жёстче». Здесь важно оценить эффект сглаживания, внесённого фильтром. Это можно сделать путём сравнения полученных полей смещений до и после атмосферной коррекции, что позволяет понять, насколько сглаживание атмосферным фильтром повлияло на детали, имеющие отношение к деформациям поверхности. Если области деформаций имеют небольшой размер, не рекомендуется чрезмерно уменьшать размер окна, чтобы не отнести полезный сигнал к атмосферной помехе. Точность расчётов контролируется получаемыми на соответствующих этапах картами точности оценок скоростей смещений, точности определения высот рельефа и среднеквадратическим отклонением временного ряда от линейного тренда. Как правило, хорошие результаты получаются при величине коррекции высоты менее 5 м, оценке точности временного ряда менее 5 мм/год для областей с когерентностью выше 0.3. Далее, в п. 1.5 приводится технология интерпретации полей смещений. Описывается подход к комплексной интерпретации различных геологических и геофизических данных, основанный на построении и использовании при решении обратных задач численной математической модели исследуемого геодинамического процесса. Выводы по главе 1. Для обработки спутниковых радарных снимков хорошо зарекомендовал себя L-диапазон (спутник ALOS-1/2). Ограничения к применению снимков ALOS: у длинноволнового диапазона хуже радиометрическое разрешение, и большой интервал съёмки территории Камчатки (1-2 раза в год). Снимки L-диапазона спутника ALOS-1/2 подходят для применения метода парных интерферограмм (DInSAR). Для исследования динамики смещений во времени (применение метода SBAS) используются радарные изображения спутника Sentinel-1, C-диапазона. Отработанная технология применения РСА- интерферометрии позволила определить сценарий обработки снимков: диапазон параметров и методов для применения на территории полуострова Камчатка. Были рассмотрены почти все заметные события вулканической и сейсмической активности, происходившие на полуострове Камчатка с 2006 года (год запуска ALOS-1) по настоящее время (2021). За этот период удалось получить поля смещений и/или временные ряды смещений на вулканах Плоский Толбачик, Корякский и Шивелуч. Глава 2. Моделирование процесса термического остывания поверхности лавового поля, сформированного в результате Трещинного Толбачинского Извержения 2012 – 2013 гг. Предложена новая численная математическая модель процесса термического остывания поверхности лавы. В начале главы приводится подробный обзор литературы по теме изучения деформаций земной поверхности в вулканических районах. Многие исследования посвящены моделированию термического остывания лавового потока и последующему оседанию поверхности лавы. Модели для разных вулканических районов, используемые в работах Shaw et al. (1977); Keszthelyi и Denlinger (1996); Neri (1998); Patrick et al. (2004); Chaussard (2016); Wittmann et al. (2017); Carrara et al. (2019) отличаются своей сложностью и диапазоном учитываемых процессов и термических параметров. Но в этих работах не рассматривалась скорость формирования мощности лавового слоя и не учитывались конкретные параметры лав Толбачика. Трещинное Толбачинское извержение началось 27 ноября 2012 года и продолжалось почти 10 месяцев, до 15 сентября 2013 года. В результате извержения сформировалось крупное лавовое поле площадью 45.8 км2 с суммарным объёмом лавы около 0.6 км3. Оценка смещений лавовой поверхности проводилась методами радарной интерферометрии по серии снимков 2017–2019 гг. По результатам обработки методом SBAS, максимальные скорости субвертикальных смещений поверхности лавового слоя зарегистрированы в районе конуса Клешня и имеют значение 285 мм/год для 2017 года, 249 мм/год для 2018 года и 261 мм/год для 2019 года. Скорость оседания поверхности в центральных и нижних частях лавовых потоков не превышает 80 мм/год. В целом Ленинградское поле оседает немного быстрее, чем Толудское (рисунок 2.1). При моделировании была использована карта мощности лавового поля, полученная путем вычисления разности двух моделей рельефа, построенных до и после извержения по данным ArcticDEM [Dai, Howat, 2017]. Рисунок 2.1 – Карта средних скоростей смещений поверхности лавы извержения 2012–13 по спутниковым данным (цветная шкала, мм/год), совмещённая с картой мощности лавы, построенной по данным, предоставленным Dai и Howat (2017) (чёрно-белая шкала, м). На рисунке 2.2 приведено сравнение скорости оседания, полученной по данным радарной интерферометрии с мощностью лавы. По плотному облаку точек в верхней части этого графика прослеживается явный тренд, из которого следует, что скорость оседания лавы пропорциональна её мощности и достигает значений ~100 мм/год для толщины потока около 60–70 м. На рисунке выделяется отдельная группа точек (ниже красной линии) с “аномально" высокими значениями скорости оседания поверхности лавы (до 260 мм/год). Для объяснения такого распределения скоростей построена новая термическая модель остывания лавового потока и выполнено сравнение результатов расчетов с результатами, полученными по спутниковым данным. Рисунок 2.2 – Сравнение скорости оседания, полученной по данным радарнойинтерферометриис мощностью лавового слоя. Красная линия обозначает границу между "нормальной"и"аномальной" скоростями оседания лавы. Для моделирования остывания лавового поля был применён подход иерархии моделей. Сначала рассматривалась простая модель (Модель I) для оценки величины термического погружения поверхности лавы с постоянными коэффициентами. По результатам расчетов по Модели I получено, что оценки термического погружения поверхности лавы в процессе ее остывания в целом близки к величинам оседания на значительной части лавового поля, но в ряде областей его центральной части реальные величины погружения в несколько раз превосходят термические оценки. Для выяснения причин высоких скоростей погружения требуются более детальные расчеты. Анализ результатов позволил выполнить постановку задачи для более детального моделирования. В модель были внесены изменения: уточнены параметры начальных условий, учтена скрытая теплота плавления, зависимости от температуры физических свойств среды, рассмотрена скорость формирования мощности лавового слоя (Модель II). Постановка задачи (Модель II). По физическим свойствам моделируемая среда разделена на слой лавы и подстилающий слой породы. Лавовый поток будем рассматривать как изначально однородную среду, состоящую из магматического расплава и некоторого начального объёма кристаллов. Поскольку горизонтальные размеры потока существенно больше вертикальных, пренебрегаем горизонтальными изменениями температуры и решаем одномерную задачу вдоль вертикальной координаты z. Для дискретизации области моделирования введём сетку элементов с шагом ∆zi (t), i = 1, ... , N; где N – количество элементов сетки, определяющее размер модели. Ось Oz направлена вниз, а начало координат совпадает с подошвой слоя лавы. Поскольку фазовый переход (образование кристаллов из расплава) сопровождается скрытым выделением тепла во время охлаждения лавы, задача сводится к решению одномерного уравнения теплопроводности с массопереносом за счет усадки и объемным источником тепла: (2.1) с начальными и граничными условиями: при при (2.2) гдеи W(z, t) – температура и скорость оседания лавы на глубине z в момент времени t; – количество скрытой теплоты, выделяемой в единице объёма в процессе фазового перехода (кристаллизации); Xfaz(z,T) - доля твёрдой фазы на глубине z; – удельная теплота кристаллизации;– удельная теплоёмкость; – температуропроводность. Начальная мощность лавового слоя равна Hini. Глубина нижней границы модели, Hmax , задаётся достаточно большой, чтобы считать, что изменением температуры ниже этой границы можно пренебречь;– начальная температура расплава,– градиент температуры в грунте к моменту начала извержения. Плотность лавы зависит от ее состава, температуры, наличия пор и количества закристаллизованного (кристаллы) и незакристаллизованного (стекло/расплав) материала. Введём долю затвердевшей части лавы, равную равновесной доле кристаллов, рассчитанной для кристаллизации наиболее примитивной лавы ТТИ-50 в программе Petrolog 3 [Danyushevsky, Plechov, 2011]. Пусть Xcr, Xgl и Xpor – объемные доли кристаллов, стекла и пор (газа) в затвердевшей части лавы, которые будем считать не зависящими от времени и температуры: Xcr + Xgl + Xpor = 1. Расплавленная часть состоит из магматического расплава и газа в порах, причем пористость расплава также равна Xpor. Речь идет о средних по толщине потока величинах, хотя в верхней части (корке) и массивной основной части доли пор и кристаллов будут различаться. Параметры Xgl и Xpor будем считать свободными, и подбирать в процессе сопоставления с данными об оседании поверхности лавы. Плотность расплавленной лавы () равна: где– плотность расплава,– плотность материала, заполняющего поры. Затвердевшая часть состоит из кристаллов с плотностью, стекла с плотностью и газа в порах (. Тогда плотность затвердевшей части –, равна: , и плотность лавы, содержащей расплав и затвердевшую часть (обе с долей газа), равна: . Подставляя сюда плотности лавы и затвердевшей части, получим: (2.3) Теперь в уравнении (2.1) выделенное в процессе кристаллизации тепло будет равно: поскольку часть расплава при затвердевании переходит в стекло и выделяемой при этом энергией можно пренебречь. При моделировании плотность кристалловрассчитана из равновесной плотности кристаллического ансамбля (Pl, Ol, Cpx) с учётом их составов (An, Fo and Mg#). Эти свойства были получены по расчетам в программе Petrolog 3 (рисунок 2.3А и Б). Плотность стекла близка к плотности расплава и меняется только за счет охлаждения: , где– коэффициент термического расширения стекла, равный 5.8·10-7 1/град, плотность расплава считается постоянной=2.6428 г/см3 и соответствует плотности при=1100 оC (рисунок 2.3A), плотность газа в порахтоже можно считать постоянной 0.001 г/см3 , поскольку ее вклад пренебрежимо мал. АБ Рисунок 2.3 – Зависимости плотности (г/см ) магмы (А) и содержания кристаллов (доли единицы) в расплаве (Б) от температуры, определенные по лабораторным данным с помощью программы Petrolog 3 [Danyushevsky, Plechov, 2011]. Высота всех элементов сеткипересчитывается на каждом временном шаге в соответствие с изменением температуры. Это позволяет рассчитать скорость смещения W(z,t) в уравнении (2.1) с использованием закона сохранения массы . Следовательно: .(2.4) Плотность кристаллов и коэффициент термического расширениязависят от кристаллического состава. Здесь используется значение, усредненное по минералам, составляющим Толбачикскую лаву. Ниже солидуса плотность рассчитывается как функция температуры: ,(2.5) с параметрами= 505оС;- равно равновесному значению плотности 2.9752 г/см3 при температуре. Зависимость (2.5) соответствует наклону кривой на рисунке 2.3А. В модели также учтён процесс постепенного формирования лавового потока. Предполагается, что слой лавы мощностью H(t) формируется в течение первого года (год извержения) с постоянной скоростью dh/dt, где dh – заданное приращение в метрах за время dt так, чтоПри расчете H(t) в первый год, в модель сверху добавлялся слой ячеек заданной мощности dh через равные интервалы времени dt до тех пор, пока не достигалась заданная конечная мощность лавового слоя Hfin. При этом считалось, что в момент времени t0 некоторый начальный слой Hini уже существует. При проведении численного моделирования было рассмотрено несколько сценариев формирования лавового слоя. При быстром формировании лавовый слой наращивался со скоростью 6 метров каждые 2 дня, при медленном – каждые 20 дней. Такие значения обеспечивают временные интервалы накопления лавы, не выходящие за пределы хронологии извержения. Так же было выполнено моделирование без учета наращивания слоя, когда считалось, что лавовый слой в начальный момент времени t=0 уже имел заданную мощность Hfin. Результаты Модели I показали, что оценки термического оседания поверхности лавы подтверждают фактические данные о смещении поверхности в областях лавового поля с малой мощностью. В центральных же частях лавовых полей с мощностями лавы > 30 м реальные
величины погружения в несколько раз превосходят значения, полученные по термической
модели (таблица 2.1).

Таблица 2.1 – Сравнение расчетных данных по Модели I для максимальной мощности слоя лавы
100 м с интерферометрическими оценками

Результаты термического
моделированияВертикальное
оседание по
ГодПогружениеданным
Уменьшение
дневнойинтерферометрии,
мощности лавового
поверхности (∆S,м
потока (∆H, м)
м)

20170.0570.0830.285
20180.0510.0740.249
20190.0460.0670.261
Модель II позволила выполнить детальное сравнение интерферометрических оценок с
результатами моделирования. В начале были рассчитаны оседания для медленного
формирования лавового потока при dh/dt = 6 м за 20 дней при различном процентном
содержании стекол и пористости в магматическом веществе и получена оценка влияния
пористости и содержания стекла в лаве на ее оседание. Расчеты показывают, что увеличение
количества стекла в конечном составе лавы, так же, как и увеличение пористости, приводит к
снижению скорости оседания поверхности лавы в течение 7-го года после извержения. На
рисунке 2.4 показаны результаты моделирования с параметрами, которые хорошо
соответствуют “нормальному” тренду оседания лавы (зеленые круги).

Рисунок 2.4 – Результаты
моделированияс
параметрами,наилучшим
образом соответствующие
“нормальным”скоростям
оседаниялавы(зелёные
круги).Жёлтыекруги
соответствуют области с
“аномальными” скоростями.
Результатырассмотрены
для сценария “медленного”
формирования слоя (dh/dt =
6м/20дн).

Выводы по главе 2: 1) Оценки оседания поверхности лавовых потоков, основанные на
данных РСА интерферометрии, показали наличие двух различных типов характера оседания:
“нормальный” и “аномальный”. 2) Более чем на 90% поверхности лавового поля зависимость
оседания от толщины лавы соответствует “нормальной” скорости оседания и может быть
хорошо объяснена моделью, основанной на термическом охлаждении и уплотнении отложений
лавы. 3) Наилучшее соответствие данным РСА было достигнуто для медленно
формирующегося слоя лавы (dh/dt = 6 м в течение 20 дней) с пористостью между 10-30% и
содержанием стекла между 30-50%. 4) Возможные причины появления области с “аномально”
высокими скоростями оседания лавовой поверхности в окрестности кратера Набоко и конуса
Клешня состоят в особых условиях формирования лавового потока. Например, эта часть потока
могла образоваться очень быстро после раскрытия трещины, в отличие от более отдаленных
частей лавовых полей, формировавшихся более медленно. Это соответствует расчётам,
основанным на сценариях “быстрого” или “мгновенного” образования слоя лавы, которые
показывают значительно более высокие скорости оседания для слоев толщиной более 25 м
(рисунок 2.5). 5) Значительное количество “аномальных” точек в части потока с относительно
небольшой мощностью не могут быть объяснены исключительно моделью, основанной на
термическом уплотнении. Следует рассмотреть другие процессы, например, необходимо
учитывать образование во время извержения в верхней части потока прямо под кратером
Набоко системы лавовых труб и каналов. Эти трубы быстро опустошались во время
извержения, и лава стекала во внешнюю часть поля, что привело к образованию множества
полостей. Области аномально быстрого оседания, вероятно, располагаются над такими
“лавовыми пещерами”.
Рисунок2.5–Сравнение
результатов,полученныхпо
спутниковым оценкам (зелёные
круги – “нормальная” скорость,
жёлтые круги – “аномальная”
скорость)срезультатами
моделирования (цветные линии) при
10% пористости и 20% содержания
стекла,которыенаилучшим
образомсоответствуют
“аномальному”оседанию.На
вертикальнойосипоказана
скорость оседания поверхности в
2019 г. (в мм/год). Красная линия – результаты для моментального формирования слоя
(Hini=Hfin, dh/dt = 0). Зелёная линия – “быстро” сформированный лавовый слой (dh/dt = 6 м в
течение 2 дн). Синяя линия – “медленно” сформированный лавовый слой (dh/dt = 6 м в течение
20 дн).

В главе 3 представлена модель формирования дайки (т.е. трещины при поступлении
магматического расплава). Эта модель была построена для проверки гипотезы о внедрении
дайки во время извержения вулкана Корякский 2008–2009 гг. Корякский стратовулкан
расположен в 30 км от г. Петропавловск-Камчатский. Извержение 2008–2009 гг. вызвало
дискуссию по вопросу: с чем связана активизация вулкана и есть ли свидетельства подъема
магмы в вулканическую постройку? Решение вопроса важно не только для понимания
динамики вулканических процессов, но и для прогноза дальнейшего развития событий и их
опасности для населения, инфраструктуры и авиации. По снимкам спутника ALOS-1,
выполненных 16.08.2007 и 06.10.2009, интервал съемки которых полностью покрывает период
извержения, выявлены смещения поверхности вулкана Корякский, с помощью
дифференциальной интерферометрии (DInSAR). В пределах снимка обнаруживается всего одна
область поднятий, на северо-западном склоне вулкана Корякский, совпадающая с областью
извержения 2008–2009 гг. Значения смещений в направлении на спутник у подножия
составляют 9–15 см и увеличиваются к вершине до 20–22 см в проекции на линию обзора
радарного луча (LOS). В предположении чисто вертикальных смещений, у вершины смещения
превосходят 28 см. По полученной карте смещений суммарный объём поднятий на северо-
западном склоне равен 1.3·106 м3. В работе [Гордеев, Дрознин, 2010] сделаны расчёты тепловой
мощности вулкана за период извержения. Энергия извержения составила 109 МДж. Такое
количество выделившейся энергии может соответствовать кристаллизации магмы объёмом ~106
м3. В работе [Кирюхин и др., 2017], выполнено моделирование инжекции магмы в трещину и
так же получена величина раскрытия трещины объёмом 1.2·106 м3.
Для интерпретации поля смещений выполнено математическое моделирование с
использованием решения Окада [Okada, 1985] о деформации поверхности упругого
полупространства в результате смещений по расположенной в нём дислокации. Проблема
применения этого решения в данной задаче состоит в том, что оно получено для дислокации,
расположенной в упругом полупространстве с горизонтальной свободной поверхностью. В
пределах области смещений высота рельефа меняется от 1500 м до 3300 м, поэтому
пренебрежение реальным рельефом может привести к ошибкам. Для этого был осуществлён ряд
преобразований системы координат так, чтобы аппроксимирующая рельеф плоскость совпадала
со свободной поверхностью упругого полупространства, а отклонения остаточного рельефа
были малы по сравнению с глубиной трещины. После решения обратной задачи по
определению параметров площадки разрыва и смещений на ней, полученные результаты
переводятся обратно в локальную систему координат и рассчитывается смещение в
направлении на спутник с использованием азимута орбиты и угла наклона зондирующего луча
(для восходящей орбиты спутника ALOS-1 они равны:= 8.16 и= 38.69 градусов
соответственно), так что смещения в направлении LOS равны:

Для аппроксимации площадки разрыва мы ограничились одним прямоугольным
элементом, что обеспечивает численную устойчивость обратной задачи (корректность по
Адамару). Прямоугольный элемент характеризуется десятью параметрами. Это три координаты
середины нижней грани прямоугольника (x, y, d); его размеры по падению и простиранию (W,
L); углы падения и простирания (); три компоненты вектора смещений (u1, u2, u3). При
выборе параметров прямоугольника мы ориентировались на результаты анализа распределения
сейсмических событий. В частности, из работы [Кирюхин и др., 2017] следует, что верхний
кластер землетрясений располагается под вулканом Корякский на уровне от 0.5 км над уровнем
моря и выше. Размер дайки, для которой выполнялось моделирование в этой работе был 2·2 км2,
угол падения 60о.
Исходя из размеров области смещений, нижняя грань прямоугольного элемента задана
равной d = 0.5 км над уровнем моря, размер по падению W = 2.4 км, по простиранию L = 1.0 км.
Угол падения варьировался в пределах 45–80о. Угол простирания, после поворота площадки
разрыва на 43.030 составляет = 133.03o. Координаты середины нижней грани прямоугольника
и его простирание легко выбрать путем смещения друг относительно друга карт расчетного и
измеренного полей смещений. В результате координаты центра дайки выбраны равными:
широта 53.328о, долгота 158.696о.
Следуя [Кирюхин и др., 2017] смещения вдоль простирания дайки (u1) заданы равными
нулю (u1=0). В этом случае необходимо найти две компоненты: раздвиг (u3) и смещение по
падению (u2), путем решения системы линейных уравнений из [Okada, 1985].
Модель дайки, аппроксимированная одним прямоугольным элементом, хорошо
приближает реальное поле смещений. Наилучшее согласование с полем смещений на спутник
достигается при углах падения трещины от 45 до 60 о. Во всех моделях сбросовая компонента
смещений составляет первые сантиметры, т.е. в пределах точности равна нулю. Расширение при
угле падения 45о составляет 82 см, при 60о – 71 см, при 80о – 64 см. Следовательно, объём
внедрившегося материала составляет 2.0·106, 1.7·106 и 1.5·106 м3, что согласуется с
приведенными выше данными других авторов.
Выводы по главе 3:
1.Смещения поверхности вулкана Корякский, полученные методами РСА
интерферометрии, превосходят 25 см и не могут быть объяснены слоем пепловых отложений,
сформированным в период извержения 2008–2009 гг. Склоновые процессы и таяние ледника
должны были привести к смещениям противоположного знака. Следовательно, наиболее
вероятной причиной смещений следует признать внедрение магматического материала в
постройку вулкана. На это указывают также данные анализа сейсмических каталогов и
тепловизионных исследований.
2.Модель дайки с глубиной нижней кромки 0.5 км над уровнем моря, шириной по
простиранию 1.0 км, по падению 2.4 км, с углом падения от 45 до 60о хорошо соответствует
смещениям, определенным по РСА интерферометрии. Объем трещины согласуется с оценками
работы [Кирюхин и др., 2017], описывающей моделирование инжекции магмы в трещину, а так
же с оценками в работе [Гордеев, Дрознин, 2010], где были выполнены расчёты энергии
извержения. Все результаты соответствуют объёму магмы по порядку величины равной 10 6 м3.
В нашей модели глубина трещины может быть увеличена на 1 км при соответствующем
уменьшении ее геометрических размеров.
3.Полученные результаты поддерживают гипотезу о том, что активизация вулкана
Корякский была связана с подъемом вулканического материала и его внедрением, в том числе, в
вулканическую постройку самого вулкана. Поэтому происходящие под вулканом процессы
могут создавать опасность для расположенных в его окрестности населенных пунктов и
объектов инфраструктуры и требуют непрерывного мониторинга.

Глава 4 посвящена исследованию пирокластического потока вулкана Шивелуч,
образовавшегося во время извержения 29.08.2019 г. С начала века вулкан извергается почти
постоянно, а 29 августа 2019 г. произошло мощное извержение взрывного типа, в результате
которого вниз по склону на 12 км распространился пирокластический поток, отложившийся на
площади 7.2 км2, объёмом 87·106 м3 по оценкам [Shevchenko et al., 2015].
По серии снимков спутника Sentinel-1А за 2020 и 2021 гг. методом малых базовых линий
(SBAS) была выявлена область с большими оседаниями, совпадающая с областью
пирокластического потока на юго-восточном склоне. В пересчёте в вертикальные смещения
максимальные значения составили 385 мм и 257 мм за 2020 и 2021 год соответственно. По
разности цифровых моделей рельефа, построенных по сериям снимков Sentinel-1 2018 года и
2020 года, была оценена мощность слоя пирокластических отложений. Их максимальная
мощность достигает 30 метров со средним значением 6.2 м. Из рисунка 4.1 четко видно, что
скорости оседаний поверхности пирокластического потока коррелируют с его толщиной
(коэффициент корреляции равен -0.69) и их можно считать зависимыми величинами. При этом
дисперсия относительно линейного тренда довольно существенная. Для объяснения этого
явления были рассмотрены влияние термического остывания и влияние компакции на усадку
пород потока, и проведено сравнение результатов моделирования с данными, полученными по
РСА интерферометрии.
Рисунок 4.1 – Зависимость скорости оседания
в мм/год (вертикальная ось) от мощности
пирокластическогопотокавметрах
(горизонтальная ось). Прямой линией показана
линейная регрессия.

Для математического моделирования и интерпретация полученных смещений была
поставлена следующая задача. Рассмотрим термическое остывание и уплотнение горячего
пирокластического материала, расположенного на слое грунта. Поток достаточно однородный
по горизонтали, его горизонтальные размеры существенно больше его толщины. Поэтому
будем рассматривать одномерную задачу, вдоль вертикальной оси Oz, которая направлена вниз.
Разобьём моделируемое пространство по вертикали на N элементов с шагом дискретизации
, i = 1, …, N, который зависит от времени . Для решения задачи термического остывания
будем решать одномерное уравнение теплопроводности со следующими начальными и
граничными условиями:
(4.1)
при
(4.2)
при

при
(4.3)
при

гдеи W(z, t) – температура и скорость оседания материала на глубине z в момент
времени t; H(t) – мощность пирокластических отложений; Hmax – глубина нижней границы
модели (задается достаточно большой, чтобы считать, что изменением температуры ниже этой
границы за 2 года можно пренебречь);– начальная температура материала;– градиент
температуры в грунте к моменту начала извержения принятый равным 25 град/км;
– температуропроводность; – теплопроводность;– плотность;– удельная
теплоёмкость.
Пусть коэффициент пористости асимптотически стремится к величине k0 при t ∞, где t –
время. Тогда– средняя начальная пористость в момент времени t = 0, где–
суммарное изменение пористости; k(t) – пористость в момент времени t > 0; (1 – ) и (1 – k(t)) –
коэффициенты скелетной фракции в начальный и текущий моменты времени соответственно;
– плотность породы в начальный момент времени t = 0. Определим деформацию ячейки
сетки за счёт уплотнения среды при изменении пористости. Из закона сохранения массы
скелета для каждого элемента сетки имеем:
,
откуда получаем изменение размера ячейки сетки по вертикали в момент времени t:
.(4.4)

Пусть средняя по разрезу пористость экспоненциально зависит от времени как:

,(4.5)
а плотность скелета зависит от температуры:
,(4.6)

где в (4.9) – (4.10) a – временной масштаб;– коэффициент термического расширения.
Обозначим термическую компоненту оседания, тогда из формулы (4.8) получаем:
(4.7)

Если k(t) определяет среднее изменение пористости по вcей толщине пирокластического
потока, то можно сначала рассчитать термическую компоненту оседаний, а затем
определить суммарное оседание по формуле (4.7).
Для нижележащего слоя грунта деформация ячеек по z учитывает только температурное
расширение/сжатие. Деформации, вызванные литостатическим давлением, при малой мощности
слоя отложений (<30 м) составляют пренебрежимо малую величину порядка 10 -5, поэтому в расчётах не учитываются. Поскольку продукты извержения вулкана Шивелуч имеют достаточно постоянный состав [Горбач, Портнягин, 2011], использованы значения физических свойств изверженной породы (начальная температура, пористость, плотность, теплоёмкость), известные по прошлым извержениям. Выбранные параметры модели приведены в таблице 4.1. Суммарное изменение пористостибыло определено методом подбора. Таблица 4.1 – Начальные параметры термомеханической модели Физическая величинаЗначение , г/см32.4 850 (горячий поток) T, о C 400 (холодный поток) , Вт/(мК)1.1 103 (лава) , Дж/(кг·К) 1.4·103 (грунт) 3.4·10-5 k00.5 ,м0.1 , дни0.5 При заданных выше термофизических параметрах среды величина погружения поверхности составляет в среднем за первый год 123 мм, а за второй год 52 мм (175 мм накопленной усадки за два года). Величина оседаний по результатам интерферометрической обработки составляет 385 и 257 мм за первый и второй год (суммарно 642 мм), что, соответственно по годам примерно в 3 и в 5 раз больше, чем получено по термической модели. Следовательно, одним термическим оседанием невозможно объяснить наблюдаемое погружение, выявленное по данным РСА интерферометрии. Добавим к термической усадке деформации, связанные с уплотнением вещества за счёт изменения пористости со временем. Наилучшее согласование расчётных и фактических данных для 30 метрового слоя получается при изменении пористости по закону (4.7) с параметрами a = 5 лет,= 0.023. Пористость при этом снизилась от = 0.523 до= 0.515 в 2021 г. В случае если лава при извержении остыла более чем в 2 раза, и средняя температура вещества в момент отложения пирокластического потока составила 400 0С, термическое оседание за первый год оказывается меньше измеренного более чем в 6.6 раз, а за второй год – более, чем в 10.7 раз. Для наилучшего совпадения результатов расчётов с фактическими результатами необходимо увеличить параметрдо 0.027, временной масштаб остается при этом таким же a = 5 лет. В этом случае пористость изменяется от 0.527 в начальный момент до 0.518 через 2 года. Поскольку скорости оседания слоев толщиной от 5 до 30 метров изменяются в широком диапазоне (рисунок 4.1), следует допустить влияние других, в частности эрозионных процессов. На рисунке 4.2А приведена средняя амплитуда отраженного радарного сигнала по снимкам спутника Sentinel-1А за период исследования (2020–2021 гг.). Видно, что склоны изрезаны множеством потоков, которые размывают пирокластические отложения. Следы этих потоков трассируются в тех местах, где на карте мощности видно резкое уменьшение толщины отложений (рисунок 4.2Б). АБ Рисунок 4.2 – А – Амплитудное изображение, среднее по снимкам Sentinel-1А за период 2020– 2021гг. Б – Карта мощности слоя пирокластических отложений извержения 29 августа 2019 года. Чёрный пунктир – следы поверхностных потоков, оцифрованные по амплитудному изображению. Изолинии – скорости оседания поверхности, полученные по РСА снимкам за 2020 год. Выводы по главе 4: 1. Снимки спутника Sentinel-1А позволили оценить мощность пирокластического потока и скорость оседания его поверхности. Карты скоростей смещений поверхности вулкана, построенные по данным РСА-интерферометрии по сериям снимков спутника Sentinel-1А за 2020 и 2021 годы в периоды с мая по октябрь, показывают оседания в области пирокластического потока с максимальной скоростью в 2020 году равной 385 мм/год, в 2021 году равной 257 мм/год. Мощность пирокластических отложений, оцененная по тем же снимкам, достигает 30 м. 2. Для исследования процесса оседания построена термомеханическая модель, в которой учтена компакция отложений за счёт изменения во времени их пористости и плотности. Сопоставление расчётных и реальных данных показало, что для объяснения зависимости скорости оседаний поверхности потока от мощности пород, достаточно предположить небольшое изменение пористости с 2019 по 2021 гг., которое в зависимости от начальной температуры потока составило от 1.5 до 1.7%. Возможно, что уменьшение пористости зависит не только от незначительного для данных толщин литостатического давления, а так же от изменения влажности породы. 3. Разброс зависимости скорость оседаний – мощность потока объясняется процессами размыва пирокластических отложений. Следы многочисленных потоков хорошо видны на спутниковых снимках. Их положение согласуется с картой мощности потока, построенной по данным спутниковой РСА интерферометрии. ЗАКЛЮЧЕНИЕ Полученные в диссертационной работе результаты являются основой для дальнейших исследований ИФЗ РАН в области применения спутниковой радарной интерферометрии, в том числе в районе Курило-Камчатской зоны субдукции. Основные результаты диссертационной работы: Разработана и обоснована путем массового тестирования на снимках C и L диапазона технология обработки РСА снимков, адаптированная к условиям Камчатки. Сформулированы методические рекомендации по применению спутниковых технологий в горных районах и в условиях северных широт. Определены величины смещений на лавовом поле Толбачик за 2019 гг, уточнены оценки смещений за 2017 и 2018 годы. Для интерпретации поля смещений построена новая математическая модель остывания лавового потока. Получены оценки различных физических свойств лавы (содержания пор, стёкол), а также скорости формирования мощности потока при моделировании остывания Толбачинского лавового поля. Оценки оседания поверхности лавовых потоков по данным РСА интерферометрии выявили наличие небольшой области с "аномальной" скоростью оседания поверхности. Термическое охлаждение и уплотнение хорошо объясняет зависимость оседания от толщины лавы более чем 90% поверхности лавового поля, которая соответствует “нормальной” скорости оседания. Показано, что наилучшее соответствие данным РСА было достигнуто для медленно формирующегося слоя лавы (dh/dt = 6 м в течение 20 дней) с пористостью между 10–30% и содержанием стекла между 30–50%. Учёт в модели скорости формирования толщины потока позволило объяснить образования "аномальной" области погружения в окрестности кратера Набоко и конуса Клешня: расчёты, основанные на сценариях “быстрого” или “мгновенного” образования слоя лавы, показывают значительно более высокие скорости оседания для слоев толщиной более 25 м. Значительное количество “аномальных” точек не могут быть объяснены исключительно моделью, основанной на термическом уплотнении. Следует рассмотреть другие процессы, например, необходимо учитывать образование во время извержения в верхней части потока прямо под кратером Набоко системы лавовых труб и каналов. Результаты смещений, полученных по РСА-интерферометрии, на склоне вулкана Корякский не могут быть объяснены слоем пепловых отложений, сформированным в период извержения 2008–2009 гг. или склоновыми процессами, что даёт основания предположить внедрение магматического вещества в постройку вулкана. Модель трещины с глубиной нижней кромки 0.5 км над уровнем моря, шириной по простиранию 1.0 км, по падению 2.4 км, с углом падения от 45 до 60о хорошо соответствует смещениям, определенным по РСА интерферометрии. Объём трещины согласуется с оценками работы [Кирюхин и др., 2017], описывающей моделирование инжекции магмы в трещину, а также с оценками работы [Гордеев, Дрознин, 2010], где были выполнены расчёты энергии извержения. Все результаты соответствуют объёму магмы по порядку величины равной 106 м3. Модель внедрения дайки в постройку Корякского вулкана подтвердила гипотезу подъёма магмы к поверхности во время извержения 2008–2009, следовательно, происходящие под вулканом процессы могут создавать опасность для расположенных в его окрестности населенных пунктов и объектов инфраструктуры и требуют непрерывного мониторинга, в том числе с применением спутниковых технологий. По снимкам спутника Sentinel-1А впервые определена скорость оседаний слоя пирокластических отложений на вулкане Шивелуч, сформировавшихся во время извержения 29 августа 2019 года. В 2020 году скорость оседаний достигала 385 мм/год, в 2021 году – 257 мм/год. Оценена мощность пирокластических отложений, оцененная по тем же снимкам, достигает 30 м. Для исследования процесса оседания построена термомеханическая модель, в которой учтена компакция отложений за счёт изменения во времени их пористости и плотности. Сопоставление расчётных и реальных данных показало, что для объяснения зависимости скорости оседаний поверхности потока от мощности пород, достаточно предположить небольшое изменение пористости с 2019 по 2021 гг., которое в зависимости от начальной температуры потока составило от 1.5 до 1.7%. Разброс зависимости скорость оседаний – мощность потока объясняется процессами размыва пирокластических отложений. Следы многочисленных потоков хорошо видны на спутниковых снимках. Их положение согласуется с картой мощности потока, построенной по данным спутниковой РСА интерферометрии.