Вулканология и сейсмология, 2022, № 4, стр. 47-66

Экспериментальное изучение и моделирование вулканических структур с использованием активных вибросейсмических методов

Б. М. Глинский a*, В. В. Ковалевский a, М. С. Хайретдинов a, А. Г. Фатьянов a, В. Н. Мартынов a, Д. А. Караваев a**, А. Ф. Сапетина a, А. Л. Собисевич b, Л. Е. Собисевич b, Л. П. Брагинская a, А. П. Григорюк a

a Институт вычислительной математики и математической геофизики СО РАН
630090 Новосибирск, просп. Акад. Лаврентьева, 6, Россия

b Институт физики Земли им. О.Ю. Шмидта РАН
123242 Москва, ул. Большая Грузинская, 10, стр. 1, Россия

* E-mail: gbm@sscc.ru
** E-mail: kda@opg.sscc.ru

Поступила в редакцию 16.11.2021
После доработки 15.02.2022
Принята к публикации 08.04.2022

Полный текст (PDF)

Аннотация

В данной статье представлен обзор работ авторов по экспериментальному изучению и математическому моделированию сейсмического поля в вулканических структурах с использованием вибраторов в качестве источников возбуждения упругих колебаний. Обобщены результаты экспериментальных исследований грязевых вулканов, проведенных ИВМиМГ СО РАН, ИФЗ РАН и КубГУ в Таманской грязевулканической провинции с помощью вибраторов. Проведено математическое моделирование в неоднородных геофизических средах для уточнения информации о структуре исследуемого объекта, а также об отличительных свойствах сейсмического поля. Разработан математический подход к моделированию вибропросвечивания грязевого вулкана произвольной геометрии с учетом глубинных разломов у вулкана, перекрывающихся слоев и т. п. На основе численных методов решения системы уравнений теории упругости разработаны параллельные алгоритмы, программные пакеты и проведены численные эксперименты на высокопроизводительных вычислительных системах. Приведены результаты расчетов сейсмического поля очаговой зоны грязевого вулкана Шуго. В данной статье представлены разработанные 3D и 2D геофизические модели и результаты моделирования сейсмического поля грязевого вулкана Карабетова гора и магматического вулкана Эльбрус. Показано, что предложенный подход с использованием активных вибросейсмических методов может успешно применяться на практике для уточнения особенностей сейсмического поля, глубинной структуры геофизических моделей и изучения влияния геометрии магматического очага и наличия выходных каналов на данные, получаемые системой наблюдения на свободной поверхности. Проведенные исследования доказывают возможность использования вибросейсмических источников с высокой точностью периодического излучения для исследования вулканических структур и активного мониторинга вулканической активности.

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

ВВЕДЕНИЕ

Сейсмология вулканов как область науки достигла значительного прогресса в изучении и прогнозировании вулканических извержений. Этому способствовало детальное изучение многих извержений и строение различных типов вулканических структур, предпринятое научными группами различных стран. Тем не менее, несмотря на увеличение количества и точности исследований, надежного метода прогнозирования извержений вулканов на основе мониторинга не существует. Например, группа ученых из различных стран в своей работе [Brenguier et al., 2011] сообщает, что, несмотря на значительные усилия, точное прогнозирование извержений и их интенсивности оказалось затруднительным. Следовательно, существует постоянная потребность в новых методах наблюдений для получения информации о продолжающихся вулканических процессах. Данная работа посвящена применению вибросейсмического исследования вулканов для разработки метода мониторинга вулканических структур.

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

Многие опубликованные работы касаются пассивной сейсмики, когда источниками возбуждения являются процессы, происходящие в каналах и магматических очагах, например, [Fehler, Aki, 1978; Chouet, 1985; Aki, Ferrazzini, 2000]. В работе [Aki et al., 1977] отмечается полезность кодовых волн для прогнозирования извержения вулкана, что было впервые признано в работе [Fehler et al., 1988] для горы Сент-Хеленс и убедительно продемонстрировано [Londono et al., 1988] для вулкана Невадо-дель-Руис. Параметр кода Q, используемый в этих исследованиях, однако, имеет высокую изменчивость [Aster et al., 1996].

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

Особое место в изучении вулканических структур занимает сейсмическая томография, когда в большинстве случаев в качестве источников используются землетрясения. Его применение к вулканической среде позволило получить изображения с высоким разрешением вулканических структур, таких как Этна [Cardaci et al., 1993], Редаут [Benz et al., 1996], Килауэа [Dawson et al., 1999], вулкан Ключевской [Koulakov et al., 2013] и других. В сейсмической томографии трехмерная сейсмическая структура определяется временами прихода продольных и поперечных волн. Следует отметить, что этот метод применим только в районах с длительными наблюдениями с использованием плотных регулярных сетей и равномерной непрерывной сейсмической активности. Землетрясения – это неизвестный и непредсказуемый источник, что ограничивает разрешающую способность вычислительных моделей и затрудняет их использование для целей мониторинга.

Эти проблемы частично исчезают при проведении эксперимента с активным сейсмическим источником на основе искусственных сейсмических сигналов, зарегистрированных плотной сетью сейсмических станций, работающих в течение относительно короткого времени.

Для вулканических регионов, расположенных близко к морю, выстрелы из пневматических пушек являются активными источниками. Например, этот метод использовался в экспериментах по изучению таких вулканических зон, как Везувий [Zollo et al., 2002], Асама [Aoki et al., 2009], Десепшн [Zandomeneghi et al., 2009], Монтсеррат [Paulatto et al., 2010] и Тенерифе [Garcia-Yeguas et al., 2012]. На суше в качестве источников используются искусственные взрывы, для которых известны параметры времени и мощности. Они успешно применяются для сейсморазведки, например, [Landru et al., 1999], а также для изучения некоторых вулканов (вулкан Асама [Aoki et al., 2009], вулкан Авачинский [Koulakov et al., 2014]). Однако взрывы не являются повторяющимися источниками и, кроме того, они слишком дороги и малоприменимы в заселенных вулканических областях, что делает их непригодными для целей мониторинга.

Альтернативный подход – использование специальных мощных прецизионных вибраторов, которые используются для вибросейсмического мониторинга Земли. Например, в начале 21 века в Японии была разработана технология ACROSS (Accurately-Controlled Routinely-Operated Signal System – точно контролируемая регулярно работающая сигнальная система) [Kumazawa et al., 2000; Kunitomo, Kumazawa, 2004]. Использование этой технологии позволило отслеживать изменение физических свойств среды с течением времени из-за изменения напряженного состояния среды, миграции жидкости в потенциальной очаговой зоне и т.п. Эта технология успешно использовалась в нескольких задачах EOR (Enhanced Oil Recovery – повышение нефтеотдачи) и CCS (Carbon Capture and Storage – улавливание и хранение углерода) [Kasahara et al., 2010а, 2010b, 2011а, 2011b, 2013]. Двухмерное численное моделирование было выполнено для обнаружения изменяющихся во времени областей, в частности, фокальной области вдоль границы опускающейся плиты, с использованием нового метода визуализации, который предполагает непрерывный активный мониторинг по технологии ACROSS и сейсмической группой [Kasahara et al., 2015].

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

ЭКСПЕРИМЕНТАЛЬНОЕ ИССЛЕДОВАНИЕ ГРЯЗЕВЫХ ВУЛКАНОВ ТАМАНСКОЙ ГРЯЗЕВУЛКАНИЧЕСКОЙ ПРОВИНЦИИ

Керченско-Таманский район отличается активной грязевулканической деятельностью и является вторым в СНГ после Азербайджана по количеству грязевых вулканов. О вулканах этого типа говорится, например, в работах [Тверитинова и др., 2015; Собисевич и др., 2015; Овсюченко и др., 2017; Канарейкин и др., 2019; Преснов и др., 2020; Mazzini, Etiope, 2017]. Грязевой вулканизм, обусловлен флюидодинамическими процессами в земной коре, в результате которых на поверхность в большом количестве извергаются перетертые осадочные породы (сопочная брекчия), вода и газ.

Серия экспериментальных исследований с вибросейсмическим зондированием этих природных явлений была проведена Институтом вычислительной математики и математической геофизики СО РАН, Институтом физики Земли РАН и Кубанским государственным университетом. В частности, впервые в практике геофизических исследований грязевых вулканов эти эксперименты проводились на грязевых вулканах Тамани: Шуго, Ахтанизовский, Карабетова гора [Глинский и др., 2007, 2008, 2010б; Ковалевский и др., 2012]. В ходе экспериментов методический план предусматривал изучение особенностей волновых процессов в окрестностях вулканического сооружения и оценку их вклада в параметры регистрируемых сейсмических полей; изучение условий трансформации волновых полей и резонансных явлений на геологических структурах вулкана; построение скоростного разреза геологической среды в районе вулканического сооружения грязевого вулкана.

Грязевой вулкан Шуго

Кратер вулкана Шуго стоит особняком среди грязевых вулканов Тамани и Северо-Западного Кавказа. Он имеет вид громадной чаши, в центре которой возвышается вулканическая постройка из светло-серой брекчии. Сама чаша напоминает собой кальдеру проседания, что отличает вулкан Шуго от других грязевых вулканов Тамани и Северо-Западного Кавказа. В этой связи он представляется наиболее интересным не только с геологической, но и с сейсмической точки зрения, что и явилось определяющим при выборе объекта для опробования предлагаемой технологии мониторинга сложно-построенных геологических объектов. Такое исследование было проведено вибросейсмическим методом ИФЗ РАН, ИВМиМГ СО РАН и КубГУ впервые в практике геофизических зондирований грязевых вулканов [Глинский и др., 2007].

При сейсмическом зондировании этого вулкана методика предусматривала изучение условий трансформации волновых полей в геологических структурах со съемкой продольных линий на траверсах “вибратор–сейсмическая станция–грязевой вулкан” и “вибратор–вулкан–сейсмическая станция”. В качестве сейсмического источника использовался вибратор СВ-10/100. В качестве зондирующих сигналов вибратора использовались сигналы развертки частотного диапазона 10–64 Гц с линейной разверткой частоты. Схема зондирования представлена на рис. 1.

Рис. 1.

Расположение вибратора и датчиков при вибрационном зондировании вулкана Шуго: точки зондирования 62, 63, 64, 55, треугольниками обозначены линии расположения линий сейсмоприемников.

В методическом плане работы выполнялись путем прохождения продольных профилей на траверсах “вибратор–регистрирующая сейсмостанция–вулкан” и “вибратор–вулкан–регистрирующая сейсмостанция”. В качестве вибратора использовался сейсморазведочный вибратор типа СВ-10/100. В состав сейсмостанции входили 12 групп сейсмических датчиков. Группы располагались на линейном профиле с шагом 30 м. В составе каждой группы было по 5 датчиков типа СВ-5. Методика проведения экспериментальных работ предусматривала проведение сеансов зондирования в 4-х пространственно разнесенных пунктах, условно обозначаемых как пункты 62, 63, 64, 55, с расстояниями между соседними пунктами 550, 1175 и 2115 м соответственно.

Помимо переноса вибратора, сейсмические группы были перегруппированы по указанным маршрутам в пределах базы “источник–приемник” 500–5500 м. При этом каждое расположение сейсмических групп соответствовало зондированиям в отмеченных точках. Количество повторных сеансов зондирования в каждой точке установки вибратора варьировалось в пределах 5–7 раз. Этого количества было достаточно для повышения помехоустойчивости вибрационных сейсмограмм при наличии повышенного сейсмического шума в районе вулкана. Расчет вибрационных сейсмограмм (аналогов взрывных сейсмограмм) проводился по алгоритму

(1)
$\begin{gathered} \bar {r}(m) = \frac{1}{L}\sum\limits_{i = 1}^L {\sum\limits_{n = 0}^{N - 1} {x_{n}^{{(t)}}{{S}_{{n - m}}}} } , \\ m = 0,K,M - 1;\,\,\,\,i = 1,K,L, \\ \end{gathered} $
где M – число дискретных отсчетов вибрационной сейсмограммы, L – число усреднений, N – множество дискретных отсчетов регистрируемого входного сигнала ${{x}_{n}} = x({{t}_{n}}),$ ${{S}_{n}} = S({{t}_{n}})$ – опорный сигнал с линейной частотной модуляцией (ЛЧМ) вида $S(t) = \alpha (t)\cos (2\pi {{f}_{0}}t + {{\beta {{t}^{2}}} \mathord{\left/ {\vphantom {{\beta {{t}^{2}}} 2}} \right. \kern-0em} 2})$, параметрами которого являются $\alpha (t)$ – огибающая, ${{f}_{0}}$ – начальная частота развертки, $\beta $ – скорость развертки по частоте, равная $\beta = {{({{f}_{{\max }}} - {{f}_{0}})} \mathord{\left/ {\vphantom {{({{f}_{{\max }}} - {{f}_{0}})} T}} \right. \kern-0em} T}{\kern 1pt} {\text{,}}$ где ${{f}_{{\max }}}$ – максимальная частота, T – длительность развертки. При проведении эксперимента за основу были выбраны следующие величины: ${{f}_{0}}$ = 10 Гц, ${{f}_{{\max }}}$ = 64 Гц, T = 60 с.

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

Важной характеристикой вибрационных сейсмограмм являются спектрально-временные функции (СВФ), позволяющие характеризовать спектральную структуру волнового сейсмического поля в зоне вулкана как функцию времени и частоты. Вычисление СВФ производится с помощью уравнения (2):

(2)
${{F}_{r}}({{\omega }_{k}},{{t}_{l}}) = {{F}_{r}}(k,l) = \sum\limits_{m = 0}^{{{M}_{l}} - 1} {{{r}_{{m + l}}}\exp ({{ - j2\pi km} \mathord{\left/ {\vphantom {{ - j2\pi km} {{{M}_{1}}}}} \right. \kern-0em} {{{M}_{1}}}})} ,$
где rm + l – текущие дискретные отсчеты вибрационных сейсмограмм, $m = 0.1,...,$ ${{M}_{1}} - 1,$ $l = 0,...,{{M}_{2}};$ ${{M}_{1}}{{M}_{2}} = M.$

С помощью уравнения (2) вычислены СВФ в зонах до вулкана и за вулканом. Вид этих функций представлен на рис. 2. Виден ряд частот, связанных с излучением вибратора (10–64 Гц).

Рис. 2.

Спектрально-временные функции коррелограмм, полученных при R = 3290 м (вулкан находится по одну сторону от вибратора и сейсмической станции) (а), R = 3380 м (вулкан находится между вибратором и сейсмической станцией) [Глинский и др., 2008, рис. 7, 8] (б).

Отметим некоторые наиболее характерные особенности сейсмических полей, полученные при вибросейсмическом зондировании тела грязевого вулкана Шуго.

1. Наблюдается существенное усложнение структуры волнового поля за вулканом, обусловленное прохождением сейсмических колебаний через сложно построенную геологическую среду вулканической постройки в сравнении с картиной волновых процессов, регистрируемых на сейсмограммах до вулкана. Во временной области это приводит к появлению на вибрационных сейсмограммах вслед за головными волнами хаотически распределенных вторичных волн, которые отсутствуют в структуре волнового поля до вулкана, т.е. когда источник и приемник находятся по одну сторону вулкана. Следствием этого является почти пятикратное уширение длительности отклика среды за вулканом по отношению к длительности до вулкана при незначительной разнице в расстояниях “источник–приемник” (3380 м по отношению 3290 м). Вследствие этого при практически одинаковых расстояниях “вибратор–сейсмостанция”, соответствующих обоим случаем зондирования, во втором наблюдается увеличение длительности преобладающей по амплитуде части сейсмограмм до 1.5–2.0 c вместо исходных 0.3 с во втором.

2. Отмеченное явление сопровождается значительным расширением спектров откликов среды за вулканом в точках профиля, что было доказано на основе сопоставления двух множеств спектрально-временных функций видов (2), полученных “до вулкана” и “за вулканом” на сравнимых расстояниях “источник–приемник”. В частности, это следует из сопоставления СВФ вибрационных сейсмограмм, представленных на рис. 2. Сравнительный анализ обоих спектров откликов достаточно убедительно показывает вклад флюидо-магматических и трещиноватых структур вулканической постройки в процесс обогащения спектра колебаний, прошедших через них. При этом наблюдается почти 10-ти кратное расширение спектров колебаний. Подобные эффекты могут быть объяснены трансформацией спектров вибросейсмических колебаний, обусловленной воздействием двух взаимосвязанных факторов: нарастающей нелинейностью и многомасштабностью геологической среды. Влияние первого из факторов по отношению к трещиноватым средам обосновано в работах А.В. Николаева и других авторов [Проблемы …, 1987]. Влияние второго фактора проявляется в нелинейности динамики микронеоднородных сред, для которых масштаб неоднородностей мал по сравнению с длиной волны поля возмущений, а число неоднородностей на длине волны велико [Володин, 2003]. Оба из указанных факторов характерны для трещиноватых зон вулканов.

3. Редуцированные годографы головных сейсмических волн (скорость редукции V = 4550 м/c), полученные по данным вибрационного зондирования, имеют хорошо выраженные линейные зависимости “время–расстояние” в диапазонах 300–5000 м.

Годографу, полученному на профиле “вибратор–приемник–вулкан”, соответствует кажущаяся средняя скорость 1740 м/с, которая с точностью 0.5% повторяется на разных участках зондирования. Второму годографу, охватывающему профили за вулканом и соответствующему схеме расстановки “вибратор–вулкан–сейсмостанция” (т.е. вулкан оказывается между источником и приемником) с точностью до 3% соответствует кажущаяся скорость 3400 м/с. Оба годографа получены на основе использования вибрационных сейсмограмм с хорошо выраженными вступлениями головных волн.

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

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

Данные вибросейсмических зондирований вулкана Шуго позволяют предположить трехслойную горизонтально-стратифицированную модель Земли с параметрами слоев, представленными в табл. 1.

Таблица 1.  

Трехслойная модель Земли в районе Шуго на основе экспериментов

Объект Глубина раздела, м Скорость преломления, м/с
Слой 1 0 1130
Слой 2 70 1740
Слой 3 632 3400

Более подробно проведенные эксперименты и полученные результаты описаны в работах [Глинский и др., 2007, 2008].

Грязевой вулкан Карабетова гора

На рис. 3 показан геологический разрез вулкана Карабетова гора по данным, представленным в работе [Шнюков и др., 2006]. Строение этого вулкана характеризуется слоистой антиклинальной структурой, сложенной различными отложениями.

Рис. 3.

Геологический разрез вулкана Карабетова гора [Шнюков и др., 2006].

В экспериментах с активным вибросейсмическим зондированием грязевого вулкана Карабетова гора (рис. 4), по [Собисевич и др., 2008] исследовалась структура сейсмического поля вибраторов при прохождении флюидосодержащих зон грязевого вулкана с использованием профильной регистрации волнового поля и методов обработки сейсмических данных.

Рис. 4.

Геологическая схема района Карабетова гора и общая схема профилей системы наблюдений при вибросейсмических исследованиях грязевого вулкана. 1 – оси антиклинальных складок и их номера; 2 – достоверно установленные разломы; 3 –разломы, идентифицированные по структурным и геоморфологическим особенностям; 4 – брекчия кургана грязевого вулкана; 5 – активные грифоны; 6 – неактивные грифоны; 7 –сальза; 8 – крупные конусовидные грифоны высотой несколько метров; 9 – очаг эксплозивного извержения 6 мая 2001 г.; 10 – точки вибросейсмического излучения вибратора СВ-10/100; 11 – пикеты профиля вибросейсмической зондирования системами RefTek (показан каждый пятый датчик) [Горбатиков и др., 2008, рис. 4; Ковалевский и др., 2012, рис. 2].

В качестве источника сейсмических сигналов использовался вибратор СВ-10/100, излучающий свип-сигналы в диапазоне частот 10‒64 Гц с линейной разверткой частоты.

Пять точек излучения сигнала вибратором располагались на линии длиной 2 км с интервалом 500 м (точки T1–T5, см. рис. 4). Поперечный профиль регистрации (профиль I, см. рис. 4) длиной 3.4 км, пересекающий вулкан, был отработан при указанном выше расположении точек излучения. Профилирование выполнялось в два этапа: с использованием 37 и 31 регистраторов RefTek с интервалом 50 м на двух последовательных участках профиля. Чтобы излучать сигналы от вибратора в противоположном направлении, точка возбуждения сигнала была помещена в конце профиля (точка T8) в дополнение к точкам T1–T5.

Регистрация вибросейсмических сигналов проводилась регистрирующими комплексами RefTek-125A и РОСА.

Для определения структуры и скоростных характеристик среды вне зоны вулкана волновое поле от вибратора регистрировалось на профиле II длиной 1–8 км с помощью 37 самописцев RefTek с шагом 50 м. Сигналы от вибратора излучались в 3 точках профиля II (точки Т8, Т10, Т11) с последующей обработкой по методике сейсморазведки по отраженным волнам.

Было выполнено множество идентичных сеансов излучения с последующей записью файлов волновых форм длиной 3300 с. Синхронизация по времени излучающей и регистрирующей систем осуществлялась с помощью GPS-приемников. Одним из сейсмоприемников RefTek осуществлялась регистрация сигнала от вибратора СВ-10/100 (генератора сигналов). Запись формы излучаемого вибратором сигнала велась с целью его использования в качестве опорного при вычислении корреляционных (вибрационных) сейсмограмм. Такие сейсмограммы были получены с помощью стандартной процедуры корреляционной свертки регистрируемых и излучаемых сигналов от вибратора, что позволило преобразовать виброзапись в импульсную форму.

На рис. 5 представлены корреляционные сейсмограммы, полученные на профилях I от источников в точках T1 и T8. На оси времени отраженная волна наблюдается в моменты времени 430–450 мс. Это время пробега хорошо аппроксимируется теоретической кривой времени пробега отраженной волны со скоростью 2200–2700 м/с.

Рис. 5.

Корреляционные сейсмограммы, полученные на профилях I от источников в точках T1 и T8.

Увеличение эффективных скоростей в направлении вулкана по сравнению с параллельным профилем можно интерпретировать как поднятие первой отражающей границы на глубине 350 м. Это связано с выклиниванием на поверхность пласта, отражение от кровли которого наблюдается и хорошо согласуется с геологическим строением вулкана Карабетова гора (см. рис. 3).

На рис. 6 представлен вариант временного разреза ОГТ (общей глубинной точки) по профилю I. Кратность суммирования равна единице, поскольку точки излучения (Т1 и Т8) располагались на концах профиля, а топография не позволяла проводить стандартные наблюдения методом ОГТ с движением вибратора по всему профилю. Временной разрез на рис. 6 показывает совпадение фаз, связанных с подъемом границ слоя к центральной части сейсмического разреза. Границы и толщины этих слоев соответствуют значениям, полученным из анализа данных профиля II (см. рис. 4). Предлагаемая в геологической модели (см. рис. 3) слоистая куполообразная структура грязевого вулкана не просматривается в центральной части сейсмического разреза. Более подробно проведенные эксперименты и полученные результаты описаны в работе [Ковалевский и др., 2012].

Рис. 6.

Пример временного разреза по методу ОГТ вдоль составного профиля I [Ковалевский и др., 2012, рис. 4].

В целом анализ сейсмограмм на профилях регистрации сигналов от различных точек излучения вибратора позволил уточнить основные особенности вулканического строения Карабетовой горы: наличие выклинивающих слоев, центрального флюидозаполненного канала вулкана, проявления резонансных эффектов сейсмических волн в выходном канале вулкана. Подробная схема натурных экспериментов и результаты обработки данных представлены в работах [Глинский и др., 2010б, 2011].

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

СЕЙСМИЧЕСКОЕ МОДЕЛИРОВАНИЕ ГЕТЕРОГЕННЫХ СРЕД НА ОСНОВЕ УРАВНЕНИЙ ТЕОРИИ УПРУГОСТИ

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

Численное моделирование распространения сейсмических волн в неоднородных упругих средах проводится на основе решения замкнутой системы уравнений теории упругости с соответствующими начальными и граничными условиями. В данной статье задача решается в терминах скоростей смещения и напряжений. В случае использования смещений мы столкнемся с необходимостью аппроксимировать производные второго порядка, включая смешанные производные. Реализация условий на границах раздела сред и граничных условий для жидко-упругой среды будет нетривиальной задачей, и потребуются дополнительные методы. Для формулировки в “скоростях-напряжениях” задача существенно упрощается, поскольку мы имеем дело с производными первого порядка, а также с набором простых средств между скоростями перемещений и напряжений (закон Гука с уравнениями, дифференцированными относительно времени). В этом случае условия на границах среды внутри области моделирования получаются естественным образом.

Моделируемой областью считается изотропная трехмерная неоднородная упругая среда сложной геометрии. Эта область имеет форму параллелепипеда с характерными линейными размерами в прямоугольной системе координат, где ось Oz направлена вертикально. При этом одна из граней моделируемой области представляет собой свободную поверхность (плоскость z = 0), на которой выполняются нулевые граничные условия. Геометрия свободной поверхности в данной работе не рассматривается. Это связано с тем, что основные изученные эффекты точно связаны с глубинным строением вулканов.

Основные уравнения рассматриваемой системы в векторном виде можно представить следующим образом:

(3)
$\rho \frac{{\partial u}}{{\partial t}} = A\sigma + F(x,y,z,t),$
(4)
$\frac{{\partial \sigma }}{{\partial t}} = Bu,$
$\begin{gathered} A = \left[ {\begin{array}{*{20}{c}} {\partial {\text{/}}\partial x} \\ 0 \\ 0 \end{array}\begin{array}{*{20}{c}} 0 \\ {\partial {\text{/}}\partial y} \\ 0 \end{array}\begin{array}{*{20}{c}} 0 \\ 0 \\ {\partial {\text{/}}\partial z} \end{array}\begin{array}{*{20}{c}} {\partial {\text{/}}\partial y} \\ {\partial {\text{/}}\partial x} \\ 0 \end{array}\begin{array}{*{20}{c}} {\partial {\text{/}}\partial z} \\ 0 \\ {\partial {\text{/}}\partial x} \end{array}\begin{array}{*{20}{c}} 0 \\ {\partial {\text{/}}\partial z} \\ {\partial {\text{/}}\partial y} \end{array}} \right], \\ B = \left[ {\begin{array}{*{20}{c}} {(\lambda + 2\mu )\partial {\text{/}}\partial x}&{\lambda \partial {\text{/}}\partial y}&{\lambda \partial {\text{/}}\partial z} \\ {\lambda \partial {\text{/}}\partial x}&{(\lambda + 2\mu )\partial {\text{/}}\partial y}&{\lambda \partial {\text{/}}\partial z} \\ {\lambda \partial {\text{/}}\partial x}&{\lambda \partial {\text{/}}\partial y}&{(\lambda + 2\mu )\partial {\text{/}}\partial z} \\ {\mu \partial {\text{/}}\partial y}&{\mu \partial {\text{/}}\partial x}&0 \\ {\mu \partial {\text{/}}\partial z}&0&{\mu \partial {\text{/}}\partial x} \\ 0&{\mu \partial {\text{/}}\partial z}&{\mu \partial {\text{/}}\partial y} \end{array}} \right], \\ \end{gathered} $
где $u = {{(u,{v},w)}^{T}}$ – вектор скорости смещения, $\sigma = {{({{\sigma }_{{xx}}},{{\sigma }_{{yy}}},{{\sigma }_{{zz}}},{{\sigma }_{{xy}}},{{\sigma }_{{xz}}},{{\sigma }_{{yz}}})}^{T}}$ – тензор напряжений, t – время, $\rho (x,y,z)$ – плотность, $\lambda (x,y,z),\mu (x,y,z)$ – коэффициенты Ламе. Начальные и граничные условия на свободной поверхности:

(5)
$\begin{gathered} {{\left. {{{\sigma }_{{xx}}}} \right|}_{{t = 0}}} = 0,\,\,\,\,{{\left. {{{\sigma }_{{yy}}}} \right|}_{{t = 0}}} = 0,\,\,\,\,{{\left. {{{\sigma }_{{zz}}}} \right|}_{{t = 0}}} = 0, \\ {{\left. {{{\sigma }_{{xz}}}} \right|}_{{t = 0}}} = 0,\,\,\,\,{{\left. {{{\sigma }_{{yz}}}} \right|}_{{t = 0}}} = 0,\,\,\,\,{{\left. {{{\sigma }_{{xy}}}} \right|}_{{t = 0}}} = 0, \\ {{\left. u \right|}_{{t = 0}}} = 0,\,\,\,\,{{\left. v \right|}_{{t = 0}}} = 0,\,\,\,\,{{\left. w \right|}_{{t = 0}}} = 0, \\ \end{gathered} $
(6)
${{\left. {{{\sigma }_{{xz}}}} \right|}_{{z = 0}}} = 0,\,\,\,\,{{\left. {{{\sigma }_{{yz}}}} \right|}_{{z = 0}}} = 0,\,\,\,\,{{\left. {{{\sigma }_{{zz}}}} \right|}_{{z = 0}}} = 0.$

Компоненты тензора напряжений связаны с компонентами тензора деформаций известными соотношениями, использующими соотношения Вольтерра, которые включают неупругие эффекты [Фатьянов, 1990, 1991].

Для численного решения задачи (3)–(6) применяются хорошо зарекомендовавшие себя разностные схемы на сдвинутых сетках. В данной работе используются схемы второго порядка приближения по времени и второго и четвертого порядков приближения по пространству [Virieux, 1986; Levander, 1988; Kristek et al., 2002] для численного решения трехмерных задач. Расчет коэффициентов сетки в виде коэффициентов Ламе осуществляется на основе интегральных законов сохранения. Однако, чтобы исключить отражения упругих волн от границ расчетной области во внутренней ее части, мы используем специальные приемы для реализации неотражающих границ. При численном решении уравнений каждая из границ расчетной области (параллелепипеда), кроме свободной поверхности, окружена поглощающим слоем. Внутри волновое поле рассчитывается с помощью первичных конечно-разностных уравнений. В так называемых поглощающих слоях расчет производится по формулам с параметрами демпфирования. Описание этого подхода и выбор значений параметров демпфирования для расчетов в соответствующих поглощающих слоях приведены в работе [Komatitsch, Martin, 2007].

Для расчета 2D-задачи используется аналогичная система уравнений и комплекс методов, но с редуцированной y-координатой. В качестве модели среды для 2D-задач обычно рассматриваются различные плоскости сечения полных 3D моделей.

На основе этих методов разработан пакет параллельных алгоритмов и программ для численного решения рассматриваемой задачи [Karavaev et al., 2015; Sapetina, 2016; Глинский и др., 2015]. Разработанный программный комплекс адаптирован для кластеров, оснащенных мощными многоядерными процессорами и графическими ускорителями (GPU). Основными технологиями, используемыми для реализации параллельных вычислений, являются MPI (интерфейс передачи сообщений) и CUDA (программно-аппаратная архитектура с использованием графических процессоров). Основные модули программного пакета состоят из следующих блоков:

(i) решатель трехмерной задачи динамической упругости для кластера с классической массивной параллельной архитектурой и гибридного кластера с GPU (графическим процессором);

(ii) решатель для двумерной задачи динамической упругости с реализацией для одного вычислительного узла гибридного кластера с графическими процессорами;

(iii) модули для построения 2D и 3D массивов значений коэффициентов сетки для среды, типичной для вулканов;

(iv) модули для записи результатов моделирования (снимки сейсмического поля и теоретические сейсмограммы) для задаваемой системы наблюдения.

Для распараллеливания на гибридном кластере мы используем усеченную двумерную декомпозицию расчетной области, масштабируемость которой в одном из направлений ограничена количеством GPU, присутствующих в вычислительном узле. Чтобы минимизировать время обмена, данные передаются между узлами с использованием соответствующих неблокирующих асинхронных функций MPI и использования функции асинхронного копирования CUDA для обмена между GPU. Для эффективного использования графических ускорителей мы оптимизируем работу с различными типами памяти, включая общую и постоянную память. Также мы изучили зависимость производительности программного обеспечения от размерности и размера блока нитей GPU и нашли оптимальные параметры. Разработанный пакет программ имеет высокую масштабируемость около 80% при имитационном моделировании исполнения программы на 106 ядер и эффективно использует архитектуру многоядерных процессоров и графических ускорителей [Glinskiy et al., 2017].

Эксперименты показали, что расчет полномасштабной 3D-модели в узлах гибридного кластера ССКЦ ИВМиМГ СО РАН занимает несколько часов, а время расчета 2D-модели составляет всего около получаса на одном GPU. В этом случае 2D-моделирование можно использовать в качестве корректировки для последующего 3D-моделирования. Таким образом, можно идентифицировать основные особенности сейсмического поля на двухмерном срезе и скорректировать модель для более точных трехмерных расчетов. Это позволяет быстро выбирать кинематические характеристики исследуемой среды.

Предложенные подход, методы и программы были успешно применены для моделирования распространения сейсмических волн в неоднородных средах, в частности, для геофизических моделей грязевых и магматических вулканов [Глинский и др., 2010а; Мартынов и др., 2018; Glinskii et al., 2015; Sapetina et al., 2021].

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВИБРОЗОНДИРОВАНИЯ ГРЯЗЕВЫХ ВУЛКАНОВ: ШУГО И КАРАБЕТОВА ГОРА

Численное моделирование грязевого вулкана Шуго

Мы дополнили экспериментальные исследования в Таманской грязевулканической провинции моделированием волновых полей для произвольных расположений вибратора, систем регистрации и других параметров. Как отмечалось в работе [Лаверов и др., 1997], можно достаточно надежно предсказать разлом в ядре синклинали, являющейся питающим каналом вулкана Шуго. Именно в этой синклинали находится грязевулканический очаг вулкана, а его “магматический” очаг находится на большой глубине. Мы стремились выяснить основные особенности волнового поля путем моделирования, проведенного в несколько этапов. В работе [Фатьянов, Мирошников, 1996] разработан метод расчета функции Грина для многомерно-неоднородных моделей сред. Он позволяет моделировать волновые поля для произвольной блоковой геометрии сред. Сначала для сравнения мы использовали простую модель нормального разлома без вулканического очага [Фатьянов, Мирошников, 1996].

Следующим шагом было включение вулканической камеры. На рис. 7а приведена сейсмограмма с источником, расположенным над центром очаговой зоны вулкана Шуго. База наблюдения составляет 3 км. По вертикали приведено время в секундах. Расчет осуществлен до времени 1.2 с. На рис. 7б показан фрагмент сейсмического профиля для трех источников, расположенных на дневной поверхности.

Рис. 7.

Моделирование волнового поля: сейсмограмма с источником над центром зоны очага (а) и фрагмент сейсмического профиля при наличии вулканической зоны (б).

Проведенные численные эксперименты показали, что математическое моделирование может быть эффективным инструментом для уточнения глубинного строения среды. Так, например, уточнена модель среды, полученная по данным вибросейсмических зондирований [Глинский и др., 2008].

Численное моделирование грязевого вулкана Карабетова гора

Для изучения и анализа зарегистрированного волнового поля при вибросейсмическом зондировании грязевого вулкана Карабетова гора была разработана трехмерная структурная модель верхней части вулкана [Глинский и др., 2010а]. Для этой модели было выполнено моделирование распространения упругих волн. Отметим, что для грязевых вулканов высота рельефа дневной поверхности незначительна по сравнению с глубиной каналов извержения, поэтому мы можем пренебречь реальной топографией в наших исследованиях. В результате серии численных экспериментов впервые была получена 3D-модель вулкана (значения упругих параметров), которая является наилучшей при сравнении с результатами геофизических экспериментов (табл. 2).

Таблица 2.  

Значения параметров геофизической модели вулкана Карабетова гора

Объект ${{V}_{p}}$, км/с ${{V}_{s}}$, км/с $\rho $, г/см3
Слой 1 1.7 0.68 2.55
Слой 2 2.2 1.1 2.9
Слой 3 2.5 1.4 2.7
Слой 4 2.6 1.25 2.8
Слой 5 2.5 1.25 3.0
Цилиндр и конус 1.3 0.81 1.2

Значения упругих параметров получены из следующих источников [Шнюков и др., 2006]. Полученная модель вулкана Карабетова гора включает пять слоев с цилиндрическим и коническим включением, которые соответствуют антиклинальной структуре, центральному заполненному флюидом канала вулкана и области его выхода на поверхность. Для моделирования был выбран источник типа “центр давления” с доминирующей частотой 25 Гц, расположенный вблизи свободной поверхности. Система наблюдения представлена профилем А–А1 (рис. 8). Геофоны размещаются на свободной поверхности. Такое размещение наблюдательной системы соответствует схеме проведенного натурного полевого эксперимента. Моделирование проводилось на кластере НКС-30Т ССКЦ.

Рис. 8.

Трехмерная структурная модель верхней части разреза грязевого вулкана Карабетова гора [Ковалевский и др., 2012, рис. 8].

Снимки рассчитанного сейсмического поля (результаты сейсмического моделирования) в разные моменты времени (рис. 9) показывают процесс прохождения прямой волны через территорию вулкана. На представленных снимках можно видеть процесс резонансного возбуждения волн в канале вулкана (сечение по профилю A–A1 вертикально ориентированных цилиндрического и конического включений в центральной части снимков) и их медленное затухание после выхода прямой и отраженной волн из этой зоны.

Рис. 9.

Снимки Z-компоненты сейсмического поля в разные моменты времени для профиля A–A1 модели (см. рис. 8). Линиями отмечены границы слоев антиклинальной структуры вулкана и центрального канала [Мартынов и др., 2018, рис. 2].

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

Рис. 10.

Теоретическая сейсмограмма от источника после обработки и сборки редуцированных экспериментальных сейсмограмм от источников T1, T6, T8 (а) и после обработки и преобразования Гильберта (б). Линии отмечают кривые времени пробега волн, отраженных от наклонных слоев [Ковалевский и др., 2012, рис. 6].

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ВИБРОЗОНДИРОВАНИЯ МАГМАТИЧЕСКОГО ВУЛКАНА ЭЛЬБРУС

Рассмотрим возможность применения технологии активного виброзондирования грязевых вулканов для изучения и мониторинга магматических вулканов на примере вулкана Эльбрус.

Двуглавый вулкан Эльбрус расположен на северной стороне Большого Кавказа. Рядом находятся густонаселенные районы Северного Кавказа, прилегающие территории юга России и северной Грузии. Период формирования вулкана составляет 20 млн лет. Последнее извержение произошло в 50 г н. э. Наличие обширного снежного и ледяного покрова делает вулкан Эльбрус еще более опасным, поскольку в случае будущих извержений любой силы и типа к вулканической опасности добавиться катастрофическая опасность в виде селей и крупномасштабных наводнений, которым может предшествовать водная волна высотой до нескольких десятков метров [Лаверов и др., 2005].

Существует возможность проводить исследования вулкана Эльбрус с помощью системы наблюдений, установленной в туннеле Баксанской нейтринной обсерватории (БНО) Института ядерных исследований РАН. Координаты входа в туннель 43°16.338′ N 42°40.878′ E. Высота над уровнем моря 1740 м. Расстояние до вершины Эльбруса 21.9 км. Азимут туннеля составляет 150°37′. Тоннель находится в поселке Нейтрино под горой Андырчи и имеет длину более 4 км [Собисевич, 2012]. Следует обратить внимание, что расположение датчиков в штольне горы сведет к минимуму промышленные помехи.

В штольне БНО ИЯИ РАН в совместных экспериментах ИВМиМГ СО РАН и ИФЗ РАН впервые была развернута и опробована в режиме непрерывной работы линейная сейсмическая группа с апертурой 2.5 км, состоящая из шести трехкомпонентных сейсмоприемников СК-1П с автономными цифровыми регистраторами “Байкал” [Ковалевский, 2013]. Результаты изучения записей слабых сейсмических событий показали, что по уровню микросейсмического шума и регистрируемых сигналов, а также по их корреляционным свойствам сейсмическая группа в туннеле БНО ИЯИ РАН имеет характеристики в короткопериодном диапазоне сопоставимые с характеристиками современных сейсмических групп, действующих в рамках Международной Системы Мониторинга [Беляшова, Михайлова, 2007; Ковалевский и др., 2014]. Это позволяет эффективно использовать такую группу для мониторинга микросейсмической активности в рамках активного мониторинга района вулкана Эльбрус с целью определения областей активизации сейсмических процессов, связанных с геодинамикой магматического очага вулкана, а также для проведения сейсмоэмиссионной томографии активных районов вулкана Эльбрус.

Рассмотрим возможность активного мониторинга вулкана Эльбрус с использованием вибросейсмических источников с высокой точностью повторного излучения при помощи численных экспериментов. В частности, мы рассмотрим возможность использования сейсмической группы, расположенной в тоннеле БНО ИЯИ РАН, для целей мониторинга и рассмотрим несколько возможных геофизических моделей строения вулкана Эльбрус.

По данным [Milyukov et al., 2008; Gurbanov et al., 2004; Авдулов, 1962, 1993; Богатиков и др., 2004; Нечаев, 1999] построены приближенные геолого-геофизические модели внутреннего строения вулкана Эльбрус с различной формой магматических очагов (рис. 11).

Рис. 11.

Приближенная геофизическая модель строения вулкана Эльбрус с эллиптическим (справа) и конусовидным (слева) магматическими очагами в слоистой среде и схема вибросейсмического мониторинга [Sapetina et al., 2021, Fig. 1; Глинский и др., 2015, рис. 1].

Модели представляют собой многослойные среды с магматическими включениями. Вулканическая структура залегает на гранитном блоке (слой +I); эффузивные породы образуют вулканическую постройку (слой +II); ниже нулевой отметки мы можем выделить еще пять слоев. Параметры слоев взяты из монографии [Собисевич, 2012] и представлены в табл. 3.

Таблица 3.  

Параметры слоистой среды геофизической модели вулкана Эльбрус

Объект Интервал глубин, км ${{V}_{p}}$, км/с ${{V}_{s}}$, км/с $\rho $, г/см3
Слой +II +2…+5 2.85 1.65 2.4
Слой +I 0…+2 3.1 1.79 2.66
Слой I 0–1 3.2 1.82 2.7
Слой II 1–6 5.9 3.42 2.85
Слой III 6–11 6.22 3.59 2.62
Слой IV 11–15 5.82 3.37 2.7
Слой V 15–22 5.97 3.45 2.75

Первая модель магматического очага основана на модели вулканов центрального типа (см. рис. 11, справа). В ней очаг определяется как эллипсоид с горизонтальной и вертикальной осями 9 км и 6 км. Другая модель вулканического очага предложена в [Богатиков и др., 2004]. По мнению авторов, очаг имеет форму конуса с направленной вверх вершиной (см. рис. 11, слева), нижняя граница корового очага находится примерно на 8–15 км ниже уровня моря. Нижние острые углы конуса неточно представляют форму камеры с точки зрения физики расплава пород, поэтому они сглаживаются вписанными окружностями. Во включениях магмы примем $\rho = 2.1~\,\,{{\text{г}} \mathord{\left/ {\vphantom {{\text{г}} {{\text{с}}{{{\text{м}}}^{{\text{3}}}}}}} \right. \kern-0em} {{\text{с}}{{{\text{м}}}^{{\text{3}}}}}},$ ${{V}_{p}} = 2.2\,\,{{{\text{км}}} \mathord{\left/ {\vphantom {{{\text{км}}} {\text{с}}}} \right. \kern-0em} {\text{с}}},$ диаметр верхнего канала – 130 м, нижнего – 260 м.

Отметим, что проведенные по этим моделям 2D- и 3D-расчеты не учитывают геометрию свободной поверхности, которая из-за своей статичности и линейности волнового поля не будет иметь достаточного влияния на мониторинг изменений волнового поля, связанных с изменениями внутренней структуры вулкана.

Мы предполагаем, что в неактивном состоянии верхний канал закупорен и сливается со слоями, которые его включают. Основываясь на этом допущении, выполнены трехмерные расчеты для модели с эллиптической камерой и двумя состояниями верхнего канала: закупоренным (состояние покоя) и заполненным магмой (момент извержения) [Glinskii et al., 2015]. Система возбуждения представляет собой точечный источник типа “центр давления” с частотой 8 Гц, расположенный вблизи свободной поверхности в левой части расчетной области в одной плоскости с осью вращения каналов и магматического очага.

Каждый расчет проводился на 11 узлах гибридного кластера, оснащенного графическими ускорителями, всего было использовано 16 896 графических ядер. Чтобы смоделировать процесс на 12 000 шагов по времени на такой вычислительной системе потребовалось в среднем чуть больше полутора часов. Результаты расчетов представлены на рис. 12 в виде теоретических сейсмограмм, снятых с профиля на поверхности, проходящего через ось вращения магматических включений. Наличие значительного различия между теоретическими сейсмограммами (а) и (б), которое наглядно представлено в виде их разницы (в), позволяет сделать вывод о возможности фиксации таких изменений при вибросейсмическом мониторинге магматических вулканов в случае значительных изменений их внутренней структуры (например, подъема магмы по каналу).

Рис. 12.

Теоретические сейсмограммы по компоненте u, полученные в ходе 3D расчетов для модели с эллиптической камерой и двумя состояниями верхнего канала: а – без верхнего канала (спящий вулкан), б – с верхним каналом, заполненным магмой (извержение), в – разница сейсмограмм (а) и (б) [Glinskii et al., 2015, Fig. 7].

Для исследования информативности сигналов, поступающих в систему наблюдения в туннеле БНО, проводилось 2D моделирование для эллипсоидальных и конических магматических очагов с закупоренным верхним каналом [Sapetina et al., 2021].

Чтобы использовать для активного мониторинга сейсмическую группу, расположенную в тоннеле БНО, рассмотрим сейсмический источник, расположенный на противоположной стороне вулкана (рис. 13).

Рис. 13.

Карта Приэльбрусья и схема вибросейсмического мониторинга вулкана Эльбрус для численных экспериментов.

Наличие подходящих дорог на картах предполагает возможность размещения вибратора в районе долины Битик-Тюбю на расстоянии около 9 км от пика Эльбруса и 30 км от системы наблюдения в тоннеле БНО. В этом случае в качестве источника может быть использован портативный вибратор ЦВ-40, разработанный в СО РАН. Он развивает колебательные силы 40–50 т и имеет рабочий диапазон 5–15 Гц [Еманов и др., 1986]. При этом потенциальное расположение вибратора оказывается примерно на линии вершины Эльбруса и тоннеля БНО. Таким образом, рассмотрим для 2D моделирования сечение описанных выше 3D моделей по линии “гора Андырчи–источник” параллельно оси Oz.

Расчеты проводились на узле с графическим ускорителем NVIDIA K40. Размер расчетной сетки 7200 × 4800 узлов. Один расчет для 56 000 временных шагов занимает в среднем около 25 мин.

Рассчитаны волновые поля для модели с эллиптической камерой (a) и модели с конической камерой (б). Для анализа информации о включениях, поступающих на приемники, мы устраняем влияние слоев на волновое поле. Для этого вычтем результаты расчета для слоистой среды без включений из расчетного поля для случаев (a) и (б). Результаты этой манипуляции показаны на рис. 14 (снимки волнового поля) и на рис. 15 (теоретические сейсмограммы).

Рис. 14.

Разница между снимками волнового поля по компоненте u для слоистой среды с магматическими включениями и без них в разные моменты времени: a – эллипсоидальная камера, б – конусовидная камера. Серая линия в правом верхнем углу указывает расположение линии сейсмоприемников. Кругами отмечена зона выхода волн из магматического очага (слева) и приход этих же волн в район системы наблюдения (справа) [Sapetina et al., 2021, Fig. 3].

Рис. 15.

Различие теоретических сейсмограмм по компоненте u для слоистой среды с магматическими включениями и без них на линии сейсмоприемников, соответствующей туннелю БНО: а – эллипсоидальная камера, б – конусообразная камера. Вертикальная ось представлена в секундах [Sapetina et al., 2021, Fig. 4].

Сравнивая их, мы можем проанализировать кинематические характеристики волнового поля с точки зрения влияния магматических включений. На снимках волнового поля линии указывают границы слоев и включений. Серой линией в правом верхнем углу отмечено положение профиля приемников в туннеле БНО. Из снимков видно, что волновое поле имеет сложную структуру со значительным влиянием формы камеры на форму проходящих через нее волн. В случае (б) камера генерирует множественные отражения внутрь. Это приводит к застреванию поля в ней, и волны, отраженные от нижней части камеры и выходящие за ее границу, достигают поверхности, не доходя до линии приемников. На снимках поля а1 и б1 отмечена зона выхода волн из магматического очага, а на снимках поля а2 и б2 отмечается приход этих же волн в область системы наблюдения (см. рис. 14). На сейсмограммах (см. рис. 15) также можно наблюдать приход этих волн в соответствующие моменты времени. Также следует отметить, что в корпусе эллиптической камеры (см. рис. 15а) на сейсмограммах мы наблюдаем приход волн, отраженных от нижней части эллипса, а в корпусе камеры конической формы (см. рис. 15б) большое количество волн, приходящих в приемники, отражаются от верхних точек дифракции камеры и распространяются внутри слоев + I и I.

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

ДИСКУССИЯ

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

В течение нескольких лет научные сотрудники СО РАН разрабатывают методы и приборы для вибросейсмического зондирования Земли. Созданы мощные источники систем возбуждения и регистрации упругих колебаний. Также были проведены многочисленные полевые эксперименты. Наиболее полно результаты этих исследований представлены в монографиях [Алексеев и др., 2004; Alekseev et al., 2010] и в работе [Seleznev et al., 2004]. Следует особо отметить, что эти источники являются экологически чистыми и, в отличие от взрывных источников, не наносят вреда природе [Glinskii et al., 2000].

Создание мощных вибросейсмических источников позволило реализовать более широкие возможности глубинного зондирования Земли. Достигнутый эффект позволяет говорить о волновом эквиваленте мощных источников с длительным накоплением слабых колебательных сигналов сравнимом со средними землетрясениями.

Также можно сказать, что данный метод сейсмологии приближен к методике сейсморазведки по разрешающей способности изучения строения земной коры и мониторинга очаговых и вулканических зон [Алексеев и др., 2004; Alekseev et al., 1995, 2010].

Например, отметим высокую амплитудно-фазовую идентичность сигналов от мощных вибраторов [Юшин и др., 1981], что позволило регистрировать лунно-солнечные вариации [Глинский и др., 1999], при этом были изучены монохроматические сигналы и сигналы с частотной разверткой (свип сигналы). Воспроизводимость вибросейсмического воздействия можно оценить как 10–3 с [Соловьев и др., 2016] и для монохроматических сигналов 10–5 с [Kovalevsky et al., 2006]. Это было определено в специальных многодневных экспериментах в стационарном сейсмологическом павильоне “Ключи” на расстоянии около 50 км от вибратора, где была измерена разница во времени повторных сеансов опорных P- и S-волн от поверхности кристаллических пород.

Из этого можно сделать вывод, что чувствительность активного мониторинга с использованием мощных вибраторов достаточна для мониторинга вулканических структур. Это также косвенно подтверждается рассмотренными выше численными экспериментами для вулкана Эльбрус. Так, например, численные расчеты для варианта с подъемом магмы показывают наличие значительного различия между сейсмограммами при изменениях в верхнем канале вулкана (см. рис. 12).

Представленные в работе результаты реальных и численных экспериментов показывают актуальность и возможность вибросейсмического мониторинга вулканических структур с использованием переносных и стационарных источников упругих волн. Таким образом, регулярное вибросейсмическое зондирование позволит по результатам обработки корреляционных сейсмограмм выявлять в динамике существенные изменения напряженно-деформированного состояния исследуемого объекта, которые могут привести к катастрофическим явлениям. В свое время ИФЗ РАН и ИВМиМГ СО РАН планировали поставить работы по мониторингу вулкана Эльбрус с применением 40-тонного передвижного вибратора ЦВ-40, разработанного в АСОМСЭ СО РАН, но, к сожалению, ограниченность финансовых ресурсов не позволила поставить эти работы.

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

ЗАКЛЮЧЕНИЕ

В работе приведен обзор и анализ ранее выполненных экспериментальных исследований волновых полей методом вибрационного просвечивания Земли грязевых вулканов в Таманской грязевулканической провинции: Шуго, Ахтанизовский, Карабетова гора и их математических моделей. Благодаря высокой повторяемости актов просвечивания зон вулканов накоплен значительный объем вибрационных сейсмограмм, позволяющий установить общие для вулканов характерные закономерности распространения волновых полей. Основные из них следующие.

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

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

Экспериментально выявлен эффект пространственного перераспределения грязепроявлений на грифонах вулкана в период вибросейсмического воздействия на среду.

Данные вибросейсмических зондирований и математическое моделирование, основанное на расчете функции Грина для среды, в которой расположен вулкан Шуго, позволяют предположить наличие трехслойной горизонтально-стратифицированной среды, в которой расположен вулкан.

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

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

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

Проведенные экспериментальные и модельные исследования показали, что вибраторы могут быть успешно применены для мониторинга вулканической среды в динамике.

Список литературы

  1. Алексеев А.С., Геза Н.И., Глинский Б.М. и др. Активная сейсмология с мощными вибрационными источниками. Новосибирск: ИВМиМГ СО РАН, 2004. 387 с.

  2. Авдулов М.В. О геологической природе гравитационной аномалии Эльбруса // Изв. АН СССР. Сер. геол. 1962. № 9. С. 67–74.

  3. Авдулов М.В. О геологической природе Эльбрусского гравитационного минимума // Вестник МГУ. Сер. 4. Геология. 1993. № 3. С. 32–39.

  4. Беляшова Н.Н., Михайлова Н.Н. Система мониторинга ядерных испытаний НЯЦ РК: развитие и возможности // Вестник НЯЦ РК. 2007. Т. 2(30). С. 5–8.

  5. Богатиков О.А., Залиханов М.Ч., Карамурзов Б.С. Природные процессы на территории Кабардино-Балкарии. М.: ИГЕМ РАН, 2004. 438 с.

  6. Володин И.А. Нелинейность и многомасштабность в сейсмоакустике // Проблемы геофизики ХХI века. М.: Наука, 2003. С. 5–36.

  7. Глинский Б.М., Ковалевский В.В., Хайретдинов М.С. Взаимосвязь волновых полей мощных вибраторов с атмосферными и геодинамическими процессами // Геология и геофизика. 1999. Т. 40. № 3. С. 431–441.

  8. Глинский Б.М., Собисевич А.Л., Хайретдинов М.С. Опыт вибросейсмического зондирования сложно построенных геологических структур (на примере грязевого вулкана Шуго) // Докл. РАН. 2007. Т. 413. № 3. С. 398–402.

  9. Глинский Б.М., Собисевич А.Л., Фатьянов А.Г. и др. Математическое моделирование и экспериментальные исследования грязевого вулкана Шуго // Вулканология и сейсмология. 2008. № 5. С. 69–77.

  10. Глинский Б.М., Караваев Д.А., Ковалевский В.В. и др. Численное моделирование и экспериментальные исследования грязевого вулкана “Гора Карабетова” вибросейсмическими методами // Вычислительные методы и программирование. 2010а. Т. 11. С. 95–104.

  11. Глинский Б.М., Ковалевский В.В., Хайретдинов М.С. Экспериментальные исследования грязевого вулкана “Гора Карабетова” вибросейсмическими методами // Труды VI Международного конгресса “ГЕО-Сибирь-2010”. Новосибирск: Изд-во СГГА, 2010б. Т. 4. С. 139–143.

  12. Глинский Б.М., Иванова И.Н. О связи спектров волнового вибросейсмического поля с геологическим строением грязевого вулкана “Гора Карабетова” // Гео-Сибирь. 2011. Т. 4. С. 100–106.

  13. Глинский Б.М., Мартынов В.Н., Сапетина А.Ф. Технология суперкомпьютерного 3D-моделирования сейсмических волновых полей в сложно построенных средах // Вестник Южно-Уральского гос. ун-та. Сер. Вычислительная математика и информатика. 2015. Т. 4. № 4. С. 101–116.

  14. Горбатиков А.В., Собисевич А.Л., Собисевич Л.Е. и др. Технология глубинного зондирования земной коры с использованием естественного низкочастотного микросейсмического поля. Изменения окружающей среды и климата природные и связанные с ними техногенные катастрофы. Т. 1. Сейсмические процессы и катастрофы. М.: ИФЗ РАН, 2008. С. 223–236.

  15. Еманов А.Ф., Кузьменко А.П., Селезнев B.C. Результаты изучения волнового поля от мощного центробежного виброисточника // Излучение и регистрация вибросейсмических сигналов. Новосибирск: ИГиГ СО АН СССР, 1986. С. 105–120.

  16. Канарейкин Б.А., Мальцев А.И., Харламов А.С. Изучение зоны проявления грязевого вулканизма на Керченском полуострове инженерно-сейсмическими методами // Геология и минерально-сырьевые ресурсы Сибири. 2019. № 1. С. 35–46.

  17. Ковалевский В.В., Глинский Б.М., Хайретдинов М.С. и др. О возможности вибросейсмических исследований грязевых вулканов // Вестник НЯЦ РК. 2012. Т. 2. № 50. С. 55–66.

  18. Ковалевский В.В. О характеристиках подземной сейсмической группы в Приэльбрусье // Вестник НЯЦ РК. 2013. Т. 2. № 54. С. 18–23.

  19. Ковалевский В.В., Белоносов А.С., Авроров С.А. и др. Локализация сейсмических событий в Приэльбрусье подземной сейсмической группой // Вестник НЯЦ РК. 2014. Т. 2. № 58. С. 123–128.

  20. Лаверов Н.П., Богатиков О.А., Гурбанов А.Г. и др. Геодинамика, сейсмотектоника и вулканизм Центрального Кавказа // Глобальные изменения природной среды и климата. М.: Наука, 1997. С. 109–130.

  21. Лаверов Н.П., Добрецов Н.Л., Богатиков О.А. и др. Новейший и современный вулканизм на территории России. М.: Наука, 2005. 604 с.

  22. Мартынов В.Н., Глинский Б.М., Караваев Д.А. и др. Моделирование вибросейсмического мониторинга вулканических структур // Интерэкспо ГЕО-Сибирь. Новосибирск: СГУГиТ, 2018. Т.4. № 2. С. 122–132

  23. Нечаев Ю.В. Космические технологии в задачах изучения локальных неоднородностей земной коры // Сборник научных трудов “Геофизика на рубеже веков”. М.: ОИФЗ РАН, 1999. С. 276–290.

  24. Овсюченко А.Н., Собисевич А.Л., Сысолин А.И. О взаимосвязи современных тектонических процессов и грязевого вулканизма на примере горы Карабетова (Таманский п-ов) // Физика Земли. 2017. № 4. С. 118–129.

  25. Проблемы нелинейной сейсмики // Сборник статей / Ред. А.В. Николаев. М.: Наука, 1987. 288 с.

  26. Преснов Д.А., Жостков Р.А., Лиходеев Д.В. и др. Новые данные о глубинном строении грязевого вулкана Джау-Тепе // Вулканология и сейсмология. 2020. № 3. С. 34–45.

  27. Собисевич А.Л., Горбатиков А.В., Овсюченко А.Н. Глубинное строение грязевого вулкана горы Карабетова // Докл. РАН. 2008. Т. 422. № 4. С. 542–546.

  28. Собисевич А.Л. Избранные задачи математической геофизики, вулканологии и геоэкологии. Т. 1. М.: ИФЗ РАН, 2012. 510 с.

  29. Собисевич А.Л., Тверитинова Т.Ю., Лиходеев Д.В. и др. Глубинное строение грязевого вулкана Джарджава в пределах Южно-Керченской антиклинальной структуры // Вопросы инженерной сейсмологии. 2015. Т. 42. № 2. С. 73–80.

  30. Соловьев В.М., Кашун В.Н., Елагин С.А. и др. Активный вибросейсмический мониторинг Алтае-Саянского региона // Интерэкспо ГЕО-Сибирь-2016. 2016. Т. 2. С. 229–233.

  31. Тверитинова Т.Ю., Собисевич А.Л., Собисевич Л.Е. и др. Структурная позиция и особенности строения и формирования грязевого вулкана горы Карабетова // Геология и полезные ископаемые Мирового океана. 2015. Т. 2. № 40. С. 106–122.

  32. Фатьянов А.Г. Полуаналитический метод решения прямых динамических задач в слоистых средах // Докл. АН СССР. 1990. Т. 310. № 2. С. 323–327.

  33. Фатьянов А.Г. Прямые и обратные задачи для тензора сейсмического момента в слоистых средах // Докл. АН СССР. 1991. Т. 317. № 6. С. 1357–1361.

  34. Фатьянов А.Г., Мирошников В.В. Энергетический метод расчета функции Грина в многомерно-неоднородных средах // Докл. РАН. 1996. Т. 351. № 2. С. 264–266.

  35. Шнюков Е.Ф., Шереметьев В.М., Маслаков Н.А. и др. Грязевые вулканы Керченско-Таманского региона. Краснодар: ГлавМедиа, 2006. 176 с.

  36. Юшин В.И., Геза Н.И., Юн Ен Дин. Об экспериментальной оценке возможности корреляционного накопления вибросейсмических сигналов для целей глубинного сейсмического зондирования// Геология и геофизика. 1981. № 8. С. 123‒126.

  37. Aki K., Ferrazzini V. Seismic monitoring and modeling of an active volcano for prediction // J. Geophysical Research: Solid Earth. 2000. V. 105(B7). P. 16617–16640.

  38. Aki K., Fehler M., Das S. Source mechanism of volcanic tremor: Fluid-driven crack models and their application to the 1963 Kilauea eruption // J. Volcanology and Geothermal Research. 1977. V. 2. P. 259–287.

  39. oki Y., Takeo M., Ayoama H. et al. P-wave velocity structure beneath Asama volcano, Japan, inferred from active source seismic experiment // J. Volcanology and Geothermal Research. 2009. V. 187. P. 272–277

  40. Alekseev A.S., Glinsky B.M., Kovalevsky V.V. et al. A Multidisciplinary Mathematical Model for Earthquake Predication Studies and Vibroseismic Monitoring of Seismic Prone Zones // Proceeding of 2-nd International Conference on Seismology and Earthquake Engineering. Tehran: JJEES, 1995. V. 1. P. 97–104.

  41. Alekseev A.S., Glinsky B.M., Kovalevsky V.V. et al. Active vibromonitoring: experimental systems and fieldwork results // Handbook of Geophysical Exploration: Seismic Exploration. 2010. V. 40. P. 105–120.

  42. Aster R.C., Slad G., Henton J. et al. Differential analysis of coda Q using similar microearthquakes in seismic gaps. Part 1. Techniques and application to seismograms recorded in the Anza Seismic gap // Bulletin of the Seismological Society of America. 1996. V. 86. P. 868–889.

  43. Benz H.M., Chouet B.A., Dawson P.B. et al. Three-dimensional P and S wave velocity structure of Redoubt Volcano, Alaska // J. Geophysical Research. 1996. V. 101. P. 8111–8128.

  44. Brenguier F., Clarke D., Aoki Y. et al. Monitoring volcanoes using seismic noise correlations // Comptes Rendus Geoscience. 2011. V. 343. Iss. 8–9. P. 633–638.

  45. Cardaci C., Coviello M., Lombardo G. et al. Seismic tomography of Etna volcano // Journal of Volcanology and Geothermal Research. 1993. V. 56. P. 357–368.

  46. Chouet B. Excitation of a buried magmatic pipe: A seismic source model for volcanic tremor // J. Geophysical Research: Solid Earth. 1985. V. 90. P. 1881–1893.

  47. Dawson P.B., Chouet B.A., Okubo P.G. et al. Three-dimensional velocity structure of the Kilauea caldera, Hawaii // Geophysical Research Letters. 1999. V. 26. P. 2805–2808.

  48. Fehler M., Aki K. Numerical study of diffraction of plane elastic waves by a finite crack with application to location of a magma lens // Bulletin of the Seismological Society of America. 1978. V. 68(3). P. 573–598.

  49. Fehler M., Roberts P., Fairbanks T. A temporal change in coda wave attenuation observed during an eruption of Mount St. Helens // J. Geophysical Research. B.: Solid Earth and Planets. 1988. V. 93(5). P. 4367–4373.

  50. Garcia-Yeguas A., Koulakov I., Ibanez J.M. et al. High resolution 3D P wave velocity structure beneath Tenerife Island (Canary Islands, Spain) based on tomographic inversion of active-source data // J. Geophysical Research: Solid Earth. 2012. V. 117. B09309.

  51. Glinskii B.M., Kovalevskii V.V., Khairetdinov M.S. Vibroseismic monitoring of earthquake-prone areas // Volcanology and Seismology. 2000. V. 21(6). P. 723–730.

  52. Glinskii B.M., Martynov V.N., Sapetina A.F. 3D modeling of seismic wave fields in a medium specific to volcanic structures // Yakutian Mathematical Journal. 2015. V. 22(3). P. 84–98.

  53. Glinskiy B., Sapetina A., Martynov V. et al. The hybrid-cluster multilevel approach to solving the elastic wave propagation problem // Parallel Computational Technologies. PCT 2017. Communications in Computer and Information Science / Eds L. Sokolinsky, M. Zymbler. Cham: Springer, 2017. V. 753. P. 261–274.

  54. Gurbanov A.G., Gazeev V.M., Bogatikov O.A. et al. Elbrus active Volcano and its geological history // Russian J. Earth Sciences. 2004. V. 6(4). P. 257–277.

  55. Karavaev D.A., Glinsky B.M., Kovalevsky V.V. A technology of 3D elastic wave propagation simulation using hybrid supercomputers // CEUR Workshop Proceedings of the 1st Russian Conference on Supercomputing – Supercomputing Days 2015. P. 26–33.

  56. Kasahara J., Korneev V., Zhdanov M. Active Geophysical Monitoring // Handbook of Geophysical Exploration, Seismic Exploration. Elsevier Science, 2010a. V. 40. P. 572.

  57. Kasahara J., Hasada Y., Tsuruga K. Seismic imaging of time lapse for CCS and oil and gas reservoirs using ultra-stable seismic source (ACROSS) // SEGJ Fall Meeting Abstract. 2010b. P. 74–77.

  58. Kasahara J., Hasada Y., Tsuruga K. Imaging of ultra-long-term temporal change of reservoir(s) by accurate seismic sources(s) and multi-receivers // Extended abstract of EAGE workshop on Permanent Reservoir Monitoring (PRM). Trondheim, Norway, 2011a. P. 40–44.

  59. Kasahara J., Ito S., Hasada Y. et al. Generation of vertical and horizontal vibrations by a synthetic method using an ultra-stable and continuous seismic source for the time lapse measurements // American Geophysical Union, Fall Meeting Abstract. 2011b. P. T43B–2310.

  60. Kasahara J., Ito S., Fujiwara T. et al. Real time imaging of CO2 storage zone by very accurate stable-long term seismic source // Energy Procedia. 2013. V. 37. P. 4085–4092.

  61. Kasahara J., Tsuruga K., Hasada Y. et al. Active monitoring of Earthquake Focal zone using ultra-stable seismic source // International Conference on Seismology Earthquake Engineering (IIEES). 2015. P. 1–7.

  62. Komatitsch D., Martin R. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation // Geophysics. 2007. V. 72(5). 1SO–Z83.

  63. Koulakov I., Jaxybulatov K., Shapiro N.M. et al. Asymmetric caldera-related structures in the area of the Avacha group of volcanoes in Kamchatka as revealed by ambient noise tomography and deep seismic sounding // J. Volcanology and Geothermal Research. 2014. V. 285. P. 36–46.

  64. Koulakov I., Gordeev E., Dobretsov N. et al. Rapid changes in magma storage beneath the Klyuchevskoy group of volcanoes inferred from time-dependent seismic tomography // J. Volcanology and Geothermal Research. 2013. V. 263. P. 75–91.

  65. Kovalevsky V.V. Estimation of sensitivity of the method of active monitoring by mono-frequency signals // Bulletin of the Novosibirsk Computing Center, Series: Mathematical Modeling in Geophysics. 2006. V. 11. P. 43–52.

  66. Kovalevsky V.V., Glinskiy B.M., Khairetdinov M.S. et al. Active vibromonitoring: experimental systems and fieldwork results // Active Geophysical Monitoring (Second Edition) / Eds J. Kasahara et al. Elsevier Science, 2020. P. 43–65.

  67. Kristek J., Moczo P., Archuleta R.J. Efficient Methods to Simulate Planar Free Surface in the 3D 4th-Order Staggered-Grid Finite-Difference Schemes // Studia Geophysica et Geodaetica. 2002. V. 46. P. 355–381.

  68. Kumazawa M., Kunitomo T., Yokoyama Y. et al. ACROSS: Theoretical and technical developments and prospect to future applications // JNC Tech. Rev. 2000. V. 9. P. 115–129.

  69. Kunitomo T., Kumazawa M. Active monitoring of the Earth’s structure by the seismic ACROSS – Transmitting and receiving technologies of the seismic ACROSS // Proc. 1st International Workshop Active Monitoring in the Solid Earth Geophysics. Mizunami, Japan, 2004. P. S4–04.

  70. Landru M., Solheim O., Hilde E. et al. The Gullfaks 4D seismic study // Petroleum Geoscience. 1999. V. 5(3). P. 213–226.

  71. Levander A.R. Fourth-order finite-difference P-SV seismograms // Geophysics. 1988. V. 53(11). P. 1425–1436.

  72. Londono B., Sanchez A., Toro E. et al. Coda Q before and after the eruptions of 13 November 1985, and 1 September 1989, at Nevado del Ruiz Volcano // Colombia. Bull Volcanol. 1988. V. 59. P. 556–561.

  73. Mazzini A., Etiope G. Mud volcanism: An updated review // Earth-Science Rev. 2017. V. 168. P. 81–112.

  74. McNutt S.R. Seismic Monitoring // Encyclopedia of Volcanoes / Eds H. Sigurdsson, B. Houghton, S. McNutt, H. Rymer, J Stix. San Diego, CA: Academic Press, 2000. P. 1095–1119.

  75. Milyukov V., Myasnikov A., Mironov A. Monitoring the State of the Magmatic Structures of Elbrus Volcano Based on Observation of Lithosphere Strains // AIP Conference Proceedings. 2008. V. 1022. P. 405–408.

  76. Paulatto M. Minshull T.A., Baptie B. et al. Upper crustal structure of an active volcano from refraction/reflection tomography, Montserrat, Lesser Antilles // Geophys. J. Int. 2010. V. 180. P. 685–696.

  77. Sapetina A.F. Supercomputer-aided comparison of the efficiency of using different mathematical statements of the 3D geophysical problem // Bulletin of the Novosibirsk Computing Center, Series: Numerical Analysis. 2016. V. 18. P. 57–66.

  78. Sapetina A.F., Glinskiy B.M., Martynov V.N. Numerical modeling results for vibroseismic monitoring of volcanic structures with different shape of the magma chamber // J. Physics: Conference Series. 2021. V. 1715. 012057.

  79. Seleznev V.S., Alekseev A.S., Goldin S.V. et al. Vibration geotechnologies in III millennium // 1st International Workshop on Active Monitoring in the Solid Earth Geophysics (IWAMO4). Proceedings, Japan. 2004. P. 39–42.

  80. Virieux J. P-SV wave propagation in heterogeneous media; velocity-stress finite-difference method // Geophysics. 1986. V. 51(4). P. 889–901.

  81. Zandomeneghi D., Barclay A.H., Almendros J. et al. Crustal structure of Deception Island volcano from P wave seismic tomography: Tectonic and volcanic implications // J. Geophys. Res. 2009. V. 114. B06310.

  82. Zollo A., D’Auria L., De Matters R. et al. Bayesian estimation of 2-D P-velocity models from active seismic arrival time data: Imaging of the shallow structure of Mt. Vesuvius (Southern Italy) // Geophys. J. Int. 2002. V. 151. P. 566–582.

Дополнительные материалы отсутствуют.