Application of the nearest neighbor method to the analysis of volcanic swarms using data from Iceland's Bárðarbunga and Fagradalsfjall volcanoes

Full Text

ВВЕДЕНИЕ

Метод ближайшего соседа Бен-Зиона‒Заляпина [Zaliapin et al., 2008; Zaliapin, Ben-Zion, 2013] позволяет устанавливать уровни статистических связей между сейсмическими событиями по схеме родитель-потомок, родителем события является его ближайший сосед, произошедшей ранее по времени. Такой метод широко применяется для задач декластеризации в тектонической сейсмичности [Баранов, Шебалин, 2019]. Известно, что распределение обобщенных расстояний до ближайшего соседа имеет бимодальную форму в большинстве случаев тектонической сейсмичности. Принято считать, что левая мода, соответствующая группе связанных событий с более близкими обобщенными расстояниями, представляет афтершоковую активность, а правая мода отвечает несвязанным событиям (иногда называют фоновыми). Стоит отметить, что это, по сути, является лишь интерпретацией наличия двух мод в распределении, хотя и вполне подтвержденной для тектонической сейсмичности. В более общем смысле различные моды такого распределения просто отражают скорее неоднородность рассматриваемой сейсмичности, которая может, например, свидетельствовать о различном механизме инициации той или иной группы, или об их различной природе. То есть в общем случае такое распределение может иметь и больше двух мод. Например, в работе [Малютин, 2023] наблюдалось распределение с тремя модами при исследовании сейсмичности Калифорнии, было показано, что одна из мод соответствовала сейсмичности района гейзеров. Также стоит отметить, что такое разделение является относительным, и если бы мы наблюдали одномодальное распределение, то однозначно интерпретировать его было бы сложно, поскольку его было бы не с чем сравнивать.

Вопрос же — почему в тектонической сейсмичности зачастую присутствует именно две группы событий — на самом деле неоднозначен. Одной из гипотез является разделение событий по природе на экзогенные/эндогенные, в этом случае часть событий считается индуцированными внешними по отношению к системе источниками, а другая часть является внутренними, самоиндуцированными [Sornette, Helmstetter, 2003]. Второй процесс в таком случае зависит от некоторого внутреннего состояния среды, ее “подготовленности”, иногда предполагается, что степень этой “подготовленности” можно охарактеризовать числом потомков, которые инициируются другими землетрясениями. Например, в работе [Nandan et al., 2021] исследуется параметр “ветвистости” процесса n, который в контексте сейсмичности имеет смысл среднего числа прямых потомков на одного родителя. Схожим по смыслу параметром является дельта продуктивность землетрясений [Баранов, Шебалин, 2019].

С этой точки зрения интересно рассмотреть вулканическую сейсмичность, так как она имеет другую природу инициации, интересно продемонстрирует ли такой сейсмический режим также две группы событий, схожие с теми, что наблюдаются для тектонической сейсмичности, или же можно будет пронаблюдать несколько мод, соответствующих различным вулканическим процессам. Науке уже известен ряд сейсмических сигналов, связанных с вулканической активностью [Minakami, 1960; Гордеев, 2007]. Стоит отметить, что далеко не все из таких сигналов являются разрывами среды аналогичными тектоническим событиям. Особое внимание стоит обратить на вулкано-тектонические сигналы — высокочастотные тектонические сигналы с гипоцентрами, локализованными на глубине в несколько десятков километров под вулканом, связанные с разрушением среды под давлением магмы [Гордеев, 2007].

Уже были предприняты попытки проанализировать вулканическую сейсмичность в обозначенном контексте, например, в работе [Traversa, Grasso, 2010]. В ней оценивалась доля фоновой сейсмичности в периоды затиший и во время интрузий магмы по дайкам для вулканов Этна и Везувий в предположении наличия двух групп, аналогичных тектонической сейсмичности. Для этого использовались распределения времен между событиями, которое для тектонической сейсмичности, как считается, аппроксимируется гамма распределением [Traversa, Grasso, 2010], причем один из параметров этого распределения характеризует долю несвязанных событий [Molchan, 2005]. Для периодов затиший на вулканах также, как и для “эталонной” тектонической сейсмичности Калифорнии было получено одинаковое гамма распределение, причем его форма не зависела ни от длительности рассматриваемого периода, ни от изменений в уровне активности, а доля фоновой сейсмичности колеблется от 20 до 40%. Однако, для участков интрузии магмы по дайкам оказалось, что распределение времен между событиями вообще не соответствует гамма распределению [Traversa, Grasso, 2010]. Этот эффект удалось затем воспроизвести на модельных данных. Эффект наблюдался, как при повышении уровня фоновой активности, так и при искусственном понижении разрешающей способности по времени.

В этой работе для той же задачи мы попытались применить другой метод — метод ближайшего соседа Бен-Зиона–Заляпина [Zaliapin et al., 2008; Zaliapin, Ben-Zion, 2013].

ДАННЫЕ

В работе анализируются данные локальной сейсмической сети Исландии (каталог доступен на сайте метеорологической службы Исландии hraun.vedur.is/ja/viku/). В частности, исследуются участки каталога, приуроченные к вулканическим процессам двух извержений вулканов Баурдарбунга (2014 г.) и Фаградальсфьядль (2021 г.).

В первом случае выбран набор событий, который представляет собой облако землетрясений, мигрирующее от центральной кальдеры к вулканическому плато Холюхрейн (Holuhraun). Это было крупное трещинное извержение, сопровождавшееся обрушением центральной кальдеры. Однако тут мы будем рассматривать события связанные именно с миграцией облака сейсмичности, предположительно вызванное интрузией магмы в породы (рис. 1а). К моменту, когда облако достигло плато, там началось истечении лавы на поверхность, след же облака совпадает с расположением системы трещин вулкана. Это свидетельствуется отчетами об извержении с сайта исландской метеорологической службы (IMO, https://en.vedur.is/earthquakes-and-volcanism/articles/nr/2947). Также подтверждения этому можно увидеть в работах [Einarsson, Brandsdóttir, 2021; Sigmundsson et al., 2015]. В работе [Sigmundsson et al., 2015] подробно проанализировано формирования дайки, выделены сегменты пути, каждый из которых является результатом пробития барьера под давлением магмы. Формирование дайки завершилось примерно 4 сентября 2014 г., а извержение началось 29 августа 2014 г. [Sigmundsson et al., 2015].

 

Рис. 1. Распределение изучаемых землетрясений в пространстве.

а — события выбранного участка каталога с 2014 по 2021 г. для извержения вулкана Баурдарбунга на плато Холюхрейн 2014‒2015 гг.; б — события каталога с 2020 по 2022 г. для извержения вулкана Фаградальсфьядль (2021 г.). Бирюзовые точки — землетрясения, красный треугольник — местоположение вулкана. Красными прямоугольниками выделены исследуемые участки.

 

Во втором случае отчетливая миграция не наблюдается, однако, сейсмичность также сконцентрирована в области системы трещин вулкана (см. рис. 1б), а на распределении событий во времени наблюдаются схожий паттерн в обоих случаях (рис. 2). Также следует отметить, что на пространственном распределении сейсмичности выражены две “ветви” сейсмичности (см. рис. 1б). Согласно исследованию [Fischer et al., 2022] одна из них соответствует сейсмичности, связанной с рифтовой зоной, другая соответствует сейсмичности в вулканической системе каналов и трещин. После извержения 2021 г. также произошло короткое извержение в 2022 г. Эти два события рассматриваются как два эпизода одного и того же вулканического процесса. Всплеск сейсмичности, начавшийся 24 февраля 2021 г., видимо, также связан с интрузией магмы в породы [Fischer et al, 2022]. Авторы связывают события примерно с 24 февраля по 19 марта 2021 г. с интрузией магмы по дайке, предшествующей извержению.

 

Рис. 2. Распределение изучаемых землетрясений во времени.

а — события выбранного участка каталога для извержения вулкана Баурдарбунга на плато Холюхрейн 2014–2015 гг.; б — события выбранного участка каталога для извержения вулкана Фаградальсфьядль (2021 г.).

Синие точки — землетрясения, салатовая линия — сейсмическая активность (число событий в сутки), вертикальными серыми линиями отмечены времена начала и конца нескольких из рассматриваемых периодов в табл. 1 и табл. 2. Красными областями выделены периоды извержений.

 

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

Границы первого участка: с 64.56 по 64.92 по широте, с –17.25 по –16.6 по долготе. Диапазон магнитуд в выбранном участке: от –0.72 до 4.63.

Границы второго участка: с 63.8 по 63.95 по широте, с –22.34 по –22.16 по долготе. Диапазон магнитуд в выбранном участке: от –0.8 до 5.72.

Далее будем ссылаться на эти два рассматриваемых случая, как BAR для вулкана Баурдарбунга и FAG для вулкана Фаградальсфьядль.

МЕТОДЫ

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

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

Оценка представительной магнитуды для вулканической сейсмичности представляет особые сложности из-за наличия нелинейных графиков повторяемости, иногда “округлых подковообразных”, иногда с двойным наклоном, иногда бимодальных (например, похожее в работе [Jacobs, Mcnutt, 2010]). Предположительно, это вызвано высокой неоднородностью данных. При этом самые популярные методы оценки представительной магнитуды (например, MAXC — метод максимальной кривизны) рассчитаны на графики “классического” вида — линейный участок с одним “острым” изломом. Для описанных же выше аномалий такие методы как правило дают заниженную оценку.

В работе реализовано три метода оценки представительной магнитуды: MAXC (Maximum curvature), GFT (Goodness-of-Fit Test), MBS (M_C by b-value stability). Эти методы подробно описаны в работе [Mignan, Woessner, 2012]. Однозначно выбрать лучший метод нельзя, поэтому надежнее всего, вероятно, использовать все три оценки, а затем выбрать из них максимальную.

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

Оценка вариаций M_C во времени проводится следующим образом.

1.  Задается некоторая начальная ширина временного окна (например, W = 10 сут).
2.  Далее это окно скользит во времени по событиям, оценивается M_C и число событий, выше порога представительности, в текущем окне.
3.  Если число событий в окне оказалось меньше некоторого заданного порога (например, N_ther = 100), тогда окно расширяется вправо, пока условие не выполнится.
4.  Следующее окно начинает скользить с расширенной правой границы предыдущего.

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

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

Далее для выбранных интервалов с полученными оценками M_C оценивается значение параметра наклона магнитудно-частотного распределения (b-value). Для этого применяется метод Бендер [Bender, 1983]. Значение фрактальной размерности, предположительно, не сильно варьируется во времени, к тому же, для ее корректной оценки требуется большое количество событий. Поэтому для всех интервалов используется одно и то же значение, рассчитанное по всем периодам в совокупности. Для случая BAR получено значение $d_f=1.2$, для случая FAG $d_f=1.3$.

Далее с полученными оценками M_C, b-value и $d_f$ строятся распределения расстояний до ближайших соседей в выделенных интервалах.

Метод ближайшего соседа

Суть метода заключается во введении метрики — обобщенного расстояния между событиями в пространстве расстояние–время–магнитуда (функция близости). Для каждого события рассчитывается это обобщенное расстояние и находится ближайший сосед, то есть наиболее вероятно связанное с ним событие [Zaliapin et al., 2008; Zaliapin, Ben-Zion, 2013].

В качестве функции близости выбрано следующее выражение [Shebalin et al., 2020]:

${\eta }_{ij}=\left\{\begin{array}{c}{t}_{ij}{r}_{ij}^{d_f}{10}^{-b{m}_i}, t_{ij}>0\\ +\infty , t_{ij}\le 0 \end{array}\right.,$

где $t_{ij}$ — время между событиями i и j, $r_{ij}$ — пространственное расстояние между ними, $m_i$ — магнитуда события i. b — параметр закона Гутенберга-Рихтера; $d_f$ — фрактальная размерность распределения эпицентров землетрясений. Более раннее из пары событий, называется “родителем”, а более позднее — его “потомком”.

Теоретическое обоснование, почему можно выбрать такое выражение в качестве функции близости было приведено в работе [Baiesi, Paczuski, 2004].

В большинстве случаев для тектонической сейсмичности события разбиваются на две группы по обобщенным расстояниям; первую группу с меньшими расстояниями, ассоциируют с кластеризованными (связанными) событиями, второю — с фоновыми (несвязанными) [Баранов, Шебалин, 2019]. Это также было проиллюстрировано модельным экспериментом в работе [Zaliapin et al., 2008].

В реальных каталогах тектонической сейсмичности присутствуют, как фоновые события, так и кластеризованные, поэтому, как правило, распределение расстояний до ближайшего соседа бимодальное и является совокупностью двух распределений. Чтобы разделить эти распределения используется немного модифицированный метод из работы [Баранов, Шебалин, 2019]. Сперва строится распределение расстояний до ближайшего соседа для полного каталога ${\rho }_{real}\left(\eta \right)$, затем левый (кластеризованный) пик грубо обрезается с варьируемым порогом (например, с шагом 0.5), события оставшейся правой части распределения перемешиваются: для каждого времени землетрясения случайным образом выбираются координаты гипоцентра и магнитуды другого события из каталога. Далее строится распределение ${\rho }_{random}\left(\eta \right)$ для ближайших соседей в полученном каталоге. Предполагается, что в случае тектонической сейсмичности такое распределение аппроксимирует фоновый пик реального распределения. Для каждого предварительного порога, процедура перемешивания применяется несколько раз и выбирается лучший вариант, поскольку перемешивание является случайной процедурой. Далее из всех предварительных порогов также выбирается лучший. Критерием, по которому выбирается лучший вариант, служит минимизация суммы квадратов разностей между столбцами гистограммы реального распределения и перемешанного распределения, умноженного на коэффициент k, который находится с помощью линейной регрессии по правому склону двух распределений. Минимизация суммы квадратов проводится в диапазоне от правого максимума реального распределения вправо. Далее можно провести декомпозицию бимодального распределения на две части:

${\rho }_{real}\left(\eta \right)=(1-k){\rho }_{clustered}\left(\eta \right)+k{\rho }_{random}\left(\eta \right)$.

Чтобы оптимизировать вес k, ищется наилучшее совпадение $k{\rho }_{random}\left(\eta \right)$ с правой ветвью ${\rho }_{real}\left(\eta \right)$ по правому склону распределений.

Пороговое значение ${\eta }_0$ находится из условия равенства интенсивностей (число событий в единицу времени) потоков кластеризованных землетрясений с ближайшими соседями $\eta \le {\eta }_0$ и некластеризованных событий с ближайшими соседями $\eta >{\eta }_0$:

$\begin{array}{c}1 - F_{clustered}\left({\eta }_0\right)= 1 -\frac{F_{real}\left({\eta }_0\right)- k{F}_{random}\left({\eta }_0\right)}{1 - k} =\\ = F_{random}\left({\eta }_0\right).\end{array}$

Входными данными для метода являются: сейсмический каталог, содержащий информацию о времени, координатах и магнитуде событий; представительная магнитуда M_C; значение наклона графика повторяемости b; фрактальная размерность каталога $d_f$.

Более подробно шаги описаны в книге [Баранов, Шебалин, 2019].

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

Критерий асимметрии выделения сейсмической энергии

Дополнительно для отличия роевой активности от крупных афтершоковых серий применяется критерий асимметрии выделения сейсмической энергии [Roland, McGuire, 2009; Passarelli et al., 2018], это особенно важно для коротких временных периодов, которые могут включать такую серию целиком, давая всплеск событий похожий на роевую активность, особенно если не наблюдается миграции роя в пространстве. Метод заключается в расчете сейсмического момента для каждого события с помощью моментной магнитуды по формуле из работы [Kanamori, 1977]:

$M_0~2*{10}^4*{10}^{1.5M_W+11.8}$.

Затем вычисляется время центроида [Jordan, 1991]:

$t_{centroid}=\overline{t}=\frac{{\displaystyle {\sum }_i^n{t}_i*{M_0}_i}}{{M_0}_{tot}}$,

где $t_i$ — время события, $M_{0i}$ — сейсмический момент события, n — число событий, $M_{0tot}$ — суммарный сейсмический момент всех событий.

Далее рассчитывается коэффициент асимметрии [Roland, McGuire, 2009]:

$skewness= \frac{{\displaystyle {\sum }_i^n{(t_i-\overline{t})}^3*{M_0}_i}}{{M_0}_{tot}*{\sigma }^3} .$

Показано [Roland, McGuire, 2009], что этот коэффициент имеет значения много более 8 для афтершоковых серий, то есть имеет левую асимметрию, так как большая часть сейсмического момента выделяется мейншоком в начале последовательности. Для сейсмических роев же высвобождение энергии более распределено во времени, давая небольшие положительные либо отрицательные значения параметра, в работе [Passarelli et al., 2018] для шести роев значения лежат в диапазоне от –4 до 4.

РЕЗУЛЬТАТЫ

Случай BAR

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

 

Рис. 3. Распределения обобщенных расстояний до ближайшего соседа в выбранных временных периодах для случая BAR.

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

 

Рис. 3. Продолжение.

 

Таблица 1. Параметры сейсмичности выделенных временных периодов для случая BAR

	Период	MC	NC	b	Активность, сутки‒1	k	Energy skewness
1	01.01.1995 ‒ 17.08.2014	2.39	91	1.28	–	0.73	–
2	17.08.2014 ‒ 21.08.2014	2.29	142	1.63	35.500	0.75	2.869
3	21.08.2014 ‒ 29.08.2014	2.29	403	1.12	50.375	0.74	0.032
4	29.08.2014 ‒ 10.09.2014	1.23	662	0.81	55.167	0.95	–0.627
5	10.09.2014 ‒ 01.12.2014	0.9	863	1.37	10.524	0.86	1.885
6	01.12.2014 ‒ 27.02.2015	1.19	138	2.19	1.568	0.16	5.671
7	27.02.2015 ‒ 01.08.2015	0.95	497	2.20	3.207	0.99	0.548
8	01.08.2015 ‒ 01.01.2016	1.39	16	–	–	–	–
9	01.01.2016 ‒ 01.01.2018	1.22	69	1.48	–	0.80	–
10	01.01.2018 ‒ 01.01.2020	1.37	44	–	–	–	–
11	01.01.2020 ‒ 01.01.2022	0.64	357	1.36	–	0.69	–
12	01.01.2022 ‒ 01.01.2023	0.7	90	1.07	–	0.79	–

Примечание. M_C — представительная магнитуда; N_C — число событий выше порога представительности; b — наклон графика повторяемости; k — доля фоновой сейсмичности; energy skewness — коэффициент асимметрии высвобожденной энергии.

 

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

В следующем периоде, начиная с 17.08.2014 начинается процесс интрузии, и форма распределения приобретает скорее одномодальную форму. Как уже упоминалось выше, при наблюдении одномодального распределения возникают трудности, поскольку масштаб расстояний оказывается не с чем сравнивать. Однако в таком случае для сравнения можно использовать как раз перемешанный каталог, его распределение показывает, какими были бы обобщенные расстояния, если бы события происходили независимо во времени (в таком случае перемешиваются все события без предварительного порога). Сравнивая наблюденное распределение с перемешанным, можно сделать выводы о наличии или отсутствии кластеризации. В данном периоде наблюденное распределение отличается от формы распределения перемешанного каталога, позволяя отделить небольшой “хвост” кластеризованной сейсмичности в левой части. Еще более выражена такая картина в следующем периоде (см. рис. 3(3)), в продолжении процесса интрузии вплоть до начала истечения лавы 29.08.2014.

С началом же извержения (см. рис. 3(4), 3(5)) сейсмичность концентрируется в области истечения лавы, а форма распределения меняется и начинает больше совпадать с распределением после перемешивания, коэффициент k при этом лежит в диапазоне от 85 до 95%.

Резкое изменение формы графика на бимодальную и падение доли фоновой сейсмичности во время извержения BAR в период № 6 (см. рис. 3(6)), видимо, объясняется одиночным сильным событием (рис. 4). Поскольку оно имеет магнитуду значительно большую, чем все остальные события в этом периоде, многие события оказываются связанными с ним, образуя крупную афтершоковую серию и обеспечивая бимодальную форму графика.

 

Рис. 4. Афтершоковая серия, являющаяся предположительной причиной выброса в доли фоновой сейсмичности во время извержения в случае BAR.

а — афтершоковая серия в пространстве; б — афтершоковая серия во времени; серые кружки — фоновые несвязанные события, синие кружки — афтершоки, красный круг — мейншок афтершоковой серии, голубые линии — связи, установленные методом ближайшего соседа.

 

Далее в период сразу после извержения (см. рис. 3(7)) сейсмичность все еще концентрируется в области, где происходила интрузия, а распределение показывает все еще низкую степень кластеризации. После чего с течением времени (см. рис. 3(9), 3(11), 3(12)) форма распределения постепенно возвращается к “классической” бимодальности.

Для периода № 3, вероятно, соответствующего интрузии магмы по дайке, также рассчитаны значения дельта продуктивности [Баранов, Шебалин, 2019]: ${\Lambda }_{0.5}=0.2$ — дельта продуктивность с $\Delta M=0.5$ и ${\Lambda }_{1.0}=0.87$ — дельта продуктивность с $\Delta M=1.0$.

Случай FAG

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

 

Рис. 5. Распределения обобщенных расстояний до ближайшего соседа в выбранных временных периодах для случая FAG.

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

 

Рис. 5. Продолжение.

 

Таблица 2. Параметры сейсмичности выделенных временных периодов для случая FAG

	Период	MC	NC	b	Активность, сутки‒1	k	Energy skewness
1	01.01.1995 ‒ 01.01.2019	1.3	1971	0.851	–	0.179	–	0.316	0.938
2	01.01.2019 ‒ 01.01.2020	0.91	714	0.776	–	0.107	–	–	–
3	01.01.2020 ‒ 01.07.2020	1.17	223	1.114	–	0.296	–	–	–
4	01.07.2020 ‒ 22.02.2021	1.61	535	0.678	–	0.198	–	–	–
5	22.02.2021 ‒ 19.03.2021	2.84	473	0.901	18.92	0.533	0.826	0.176	0.415
6	19.03.2021 ‒ 18.09.2021	0.8	1207	1.055	6.596	0.822	2.649	0.129	0.511
7	18.09.2021 ‒ 20.12.2021	0.89	1179	0.767	12.677	0.044	5.697	0.331	0.990
8	20.12.2021 ‒ 27.12.2021	2.6	153	0.898	21.860	–	–	–	–
9	27.12.2021 ‒ 29.07.2022	0.82	626	1.304	2.930	0.678	–	–	–
10	29.07.2022 ‒ 03.08.2022	3.17	43	–	–	–	–	–	–
11	03.08.2022 ‒ 21.08.2022	0.75	231	0.818	12.830	0.728	22.841	–	–
12	21.08.2022 ‒ 01.01.2023	0.54	339	1.288	2.550	0.578	–	–	–

Примечание. M_C — представительная магнитуда; N_C — число событий выше порога представительности; b — наклон графика повторяемости; k — доля фоновой сейсмичности; energy skewness — коэффициент асимметрии высвобожденной энергии; ${\Lambda }_{0.5}$ — дельта продуктивность с $\Delta M=0.5$; ${\Lambda }_{1.0}$ — дельта продуктивность с $\Delta M=1.0$.

 

Видно, что в периоды задолго до извержения (см. рис. 5(1)–5(4)) сейсмичность в основном концентрируется внутри области, выделенной зеленым эллипсом на вставках рис. 5, наблюдается хорошо выраженная бимодальность распределения обобщенных расстояний до ближайшего соседа с высокой степенью кластеризации (от 70 до 90%). Эта сейсмичность, вероятно, имеет тектоническую природу [Fischer et al., 2022].

К моменту же извержения 2021 г. события начинают концентрироваться в области вулканической системы трещин (см. рис. 5, выделено красными эллипсами на вставках). Переходным периодом является № 5 (см. рис. 5(5)), поэтому часть событий, которая может относится к тектонической сейсмичности была удалена из рассмотрения (см. рис. 5(5), удаленные события выделены зеленым цветом на вставке). В этом периоде, соответствующем процессу интрузии перед извержением, наблюдается одномодальное распределение, схожее с тем, что наблюдалось для случая BAR (см. рис. 3(2), 3(3)). Далее же во время извержения (см. рис. 5(6)) также наблюдается одномодальное распределение, но с выраженным “хвостом” в левой части, в отличие от случая BAR, однако степень кластеризации в этом периоде также падает по сравнению с другими периодами.

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

Далее в промежутке между двумя эпизодами извержения (см. рис. 5(9)) распределение снова приобретает бимодальность. Затем во время второго эпизода извержения (см. рис. 5(11)) имеет скорее одномодальную форму, однако этот эпизод слишком короткий, чтобы делать какие-то однозначные выводы.

В период же после извержения распределение снова стремится к бимодальной форме.

В целом в обоих случаях наблюдается различие в режиме сейсмичности непосредственно перед извержением, во время извержения и за долго до или после извержения (см. рис. 3(2)–3(5), 3(7)), 5(5), 5(6), 5(11)). При этом для тектонической сейсмичности в тех же регионах наблюдается “классический” бимодальный график.

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

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

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

Два типа одномодальных распределений

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

Стоит отметить, что первый тип одномодального распределения, судя по всему, также соответствует афтершоковым сериям (например, рис. 6). Параметр асимметрии энергии для этих событий: 11.10.

 

Рис. 6. Пример афтершоковой серии из региона для случая FAG.

a — распределение событий серии во времени; б — распределение обобщенных расстояний до ближайшего соседа афтершоковой серии.

 

Неоднозначность декомпозиций

Следует упомянуть, что для первого типа одномодальных распределений сложно однозначно провести декомпозицию на две части, и оно носит частично субъективный характер. Например, можно представить альтернативную декомпозицию для периода № 3 случая BAR (рис. 7).

 

Рис. 7. Пример альтернативного варианта декомпозиции на фоновую и кластеризованную части для периода № 3 (см. табл. 1).

Голубая гистограмма — фактическое реальное распределение, красная гистограмма — распределение рандомизированного каталога (аппроксимация фонового пика), желтая гистограмма — полученное распределение кластеризованной части.

 

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

Значения дельта продуктивностей для альтернативной декомпозиции получаются такие: ${\Lambda }_{0.5}=0.3$ и ${\Lambda }_{1.0}=1.13$, а доля фоновой сейсмичности — около 40%.

Чтобы попробовать решить проблему с одномодальными распределениями первого типа, можно разбить выбранный период на еще меньшие (табл. 3) и проследить изменения распределения во времени (рис. 8).

 

Таблица 3. Параметры для набора последовательных пересекающихся окон

	Период	Активность, сутки‒1	Energy skewness
1	18.08.2014 14:37:25 ‒ 22.08.2014 08:06:39	75.631	0.262
2	19.08.2014 13:04:54 ‒ 24.08.2014 10:52:52	50.120	–1.385
3	21.08.2014 00:56:32 ‒ 24.08.2014 11:36:40	70.547	–2.725
4	22.08.2014 06:39:52 ‒ 25.08.2014 19:12:23	56.493	–0.031
5	23.08.2014 13:25:44 ‒ 26.08.2014 19:25:45	103.384	0.222
6	24.08.2014 12:15:30 ‒ 28.08.2014 14:30:14	54.964	–0.399
7	25.08.2014 22:41:39 ‒ 30.08.2014 21:41:18	42.557	1.409
8	27.08.2014 23:13:11 ‒ 01.09.2014 10:00:54	57.531	–0.618
9	29.08.2014 18:53:29 ‒ 03.09.2014 01:12:21	92.913	2.693

Примечание. Energy skewness — коэффициент асимметрии высвобожденной энергии.

 

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

Голубая гистограмма — фактическое реальное распределение, красная гистограмма — распределение рандомизированного каталога (аппроксимация фонового пика), желтая гистограмма — полученное распределение кластеризованной части. Номера графиков соответствуют номерам строк из табл. 3.

 

Рис. 8. Продолжение.

 

На рис. 8 прослеживается переход формы распределения от первого типа ко второму. Графики, охватывающие период с 21 по 30 августа 2014 г. в совокупности, представляют период № 3 (см. табл. 1), на всех из них доля фоновой сейсмичности k колеблется примерно от 60 до 70%. Исходя из этого, можно отдать предпочтение первому варианту декомпозиции для этого участка (см. рис. 3(3)), который также дает значение k = 74%, хотя это и не является строгим доказательством.

В целом с точки зрения формы распределения расстояний до ближайшего соседа наблюдаются три типа ситуаций.

1. “Классическое” бимодальное распределение, легко разделяется на две составляющие.
2. Одномодальное распределение, которое после рандомизации практически полностью совпадает само с собой.
3. Одномодальное распределение с небольшим “горбом” или хвостом на левом склоне распределения. При тщательном подборе параметров мода рандомизированного распределения приближается к моде исходного распределения, а также хорошо аппроксимирует его правый склон, позволяя произвести декомпозицию на две составляющие.

Практика показывает, что первый случай наблюдается для “классической” тектонической сейсмичности в большинстве регионов Земли. Второй случай также наблюдается во некоторых ситуациях, например, для гейзерной сейсмичности [Малютин, 2023], также в лабораторных экспериментах [Маточкина, 2023]. Третья ситуация, видимо, наблюдается в случаях высокой активности, из-за чего обе компоненты оказываются слишком близко друг к другу и становятся плохо различимы, примеры такого также получены в лабораторных экспериментах [Маточкина, 2023] и для техногенной сейсмичности [Баранов и др., 2020]. Предположительно, с этим же эффектом столкнулись авторы статьи [Traversa, Grasso, 2010].

Таким образом распределение, наблюдаемое в период интрузии, вероятно, проявляет некоторую внутреннюю неоднородность сейсмичности.

ОСНОВНЫЕ ВЫВОДЫ

1) Для сейсмичности, приуроченной к извержениям, наблюдаются одномодальные формы распределения расстояний до ближайшего соседа, в отличие от “классического” бимодального распределения, известного для тектонической сейсмичности.

2) Наблюдаются два типа одномодальных распределений:

а) одномодальное распределение, которое после рандомизации практически полностью совпадает само с собой;

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

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

4) Одномодальность распределения в случае вулканических роев можно гипотетически объяснить тем, что из-за высокой сконцентрированности событий в пространстве и времени фоновый и кластеризованный пики находятся слишком близко друг другу и становятся плохо различимы. В таком случае попытки разделить распределение на два дают оценки доли фоновой сейсмичности от 60 до 90%. Такая форма распределения, вероятно, говорит о внутренней неоднородности таких сейсмических режимов. Также стоит отметить, что очень похожее распределение наблюдается в лабораторных исследованиях для акустической эмиссии в образцах горных пород [Маточкина, 2023], а также было наблюдено для техногенной сейсмичности [Баранов и др., 2020].

5) Первый же тип одномодальных распределений проявляется уже во время извержения. Стоит отметить, что очень низкая степень кластеризации наблюдается, например, для сейсмичности в геотермальных зонах, в долинах гейзеров [Малютин, 2023]. И такой тип одномодальных распределений (полное совпадение рандомизированного распределения с реальным) наблюдается для событий, инициированных закачкой воды в горные породы [Малютин, 2023].

6) Были предприняты попытки провести декомпозицию событий на фоновую и кластеризованную части для одномодальных распределений вулканических роев, предшествующих извержениям. И далее на основе этого разделения рассчитать дельта продуктивность, в виду малого числа событий и их низкой магнитуды, не удается выбрать значения $\Delta M$ больше 1. Значения получается похожими на известные для тектонической сейсмичности (около 0.5), например, [Баранов, Шебалин, 2019]. Однако, для случая BAR перед извержением (см. табл. 1, период № 3) значение аномально большое 0.87, либо 1.13 в зависимости от выбранного варианта декомпозиции.

ЗАКЛЮЧЕНИЕ

Была проанализирована вулканическая сейсмичность двух извержений вулканов региона Исландии на предмет группирования и однородности. В данной работе для этого был применен метод ближайшего соседа Бен-Зиона–Заляпина. Показано, что распределение обобщенных расстояний до ближайшего соседа, которое обычно в тектонической сейсмичности имеет бимодальную форму, в периоды, приуроченные к вулканической активности, становится одномодальным. Причем продемонстрировано, что эта особенность связана именно с периодами извержений, а не с выбором рассматриваемых регионов. Однако и среди одномодальных распределений путем сравнения с распределением рандомизированного каталога было выявлено два типа. Есть основания полагать, что один из типов, наблюдающийся преимущественно в периоды интрузий перед извержением, является скорее композицией двух близких распределений, свидетельствуя о неоднородности сейсмичности, при этом стоит иметь в виду, что природа наличия такой неоднородности может отличаться от той, что наблюдается в тектонической сейсмичности. Во время же извержений это неоднородность уменьшается или почти исчезает. В целом оба типа распределений уже наблюдались в других исследованиях: в лабораторных экспериментах [Маточкина, 2023], для техногенной сейсмичности [Баранов и др., 2020], при изучении событий, инициированных закачкой воды в горные породы [Малютин, 2023]. Установление природы неоднородности в периоды интрузий может быть важным направлением дальнейшего исследования.

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

Автор выражает благодарность исландской метеорологической службе IMO (https://en.vedur.is) за открытый доступ к каталогу землетрясений локальной сейсмической сети.

ФИНАНСИРОВАНИЕ РАБОТЫ

Работа выполнена при поддержке проекта Российского Научного Фонда № 20-17-00180П.

КОНФЛИКТ ИНТЕРЕСОВ

Автор данной работы заявляет, что у него нет конфликта интересов.

About the authors

E. M. Grekov

Institute of Earthquake Prediction Theory and Mathematical Geophysics Russian Academy of Sciences (IEPT RAS); Faculty of Physics, Department of the Physics of the Earth, Lomonosov Moscow State University

Author for correspondence.
Email: grekov.em16@physics.msu.ru
Russian Federation, Profsoyuznaya str., 84/32, Moscow, 117997; Leninskie Gory, 1, bld. 2, Moscow, 119991

References

Supplementary files