Физика Земли, 2020, № 3, стр. 40-51

ВВЕДЕНИЕ

Продуктивность – это общее количество сейсмических событий, возникающих в результате возмущения напряженного состояния, вызванного другим, более ранним событием. Такое понятие продуктивности впервые было использовано для разработки подходящей модели возникновения афтершоков c учетом закона Омори–Утсу [Utsu, 1961; 1970]. Продуктивность землетрясения является ключевой характеристикой в исследовании сейсмического режима и имеет критическое значение для оценивания опасности повторных толчков.

С момента появления стохастических моделей сейсмичности, рассматривающих ее как ветвящийся процесс, изучение свойства продуктивности стало основной задачей, поскольку она является ключевым параметром, определяющим увеличение уровня сейсмичности после каждого землетрясения [Kagan, Knopoff, 1981; Ogata, 1989; Helmstetter, Sornette, 2002]. Во всех этих моделях число событий, вызванных землетрясением c магнитудой m, считается случайной величиной, имеющей распределение Пуассона с интенсивностью

Значения α варьируются от 0.5 до 2.3 [Felzer et al., 2004; Hainzl, Marsan, 2008; Hainzl et al., 2013; Marsan, Helmstetter, 2017; Wang et al., 2010; Werner, Sornette, 2008; Zhuang et al., 2004], но они всегда близки к наблюдаемому значению наклона графика повторяемости b. Тем не менее, эти оценки остаются неопределенными из-за сложности выделения относительного вклада последовательных событий в серии.

С точки зрения распределения магнитуд землетрясений (закон Гутенберга–Рихтера [Gutenberg, Richter, 1944]), продуктивность является постоянной величиной, которую можно рассматривать независимо от масштабирующего значения b, представляющего собой наклон графика повторяемости. Зная продуктивность и принимая во внимание закон Омори–Утсу [Utsu, 1961], можно рассчитать число повторных толчков, ожидаемых в заданном временном интервале [Shcherbakov, Turcotte, 2004; Holschneider et al., 2012; Баранов и др., 2019].

Недавние исследования показали [Шебалин и др., 2018; Баранов, Шебалин, 2019], что продуктивность тектонических землетрясений на глобальном и региональном уровнях подчиняется экспоненциальному закону распределения независимо от магнитуды землетрясения. Продуктивность систематически уменьшается с глубиной, сохраняя свою экспоненциальную форму.

В настоящей статье на примере землетрясений Хибинской природно-технической системы (ПТС) мы покажем, что продуктивность в условиях техногенной сейсмичности также распределена экспоненциально. При этом вид распределения не зависит от нижнего порога магнитуд рассматриваемых инициированных событий, магнитуд и глубин событий-триггеров. Данный результат позволяет утверждать, что закон продуктивности подобно двум другим законам сейсмологии (Гутенберга–Рихтера и Омори–Утсу) справедлив как для тектонических, так и для техногенных землетрясений. Это позволяет сделать предположение об универсальности этого закона.

КРАТКАЯ СЕЙСМОТЕКТОНИЧЕСКАЯ ХАРАКТЕРИСТИКА ХИБИНСКОГО МАССИВА И ИСХОДНЫЕ ДАННЫЕ

Хибинский массив представляет крупную щелочную интрузию палеозойского возраста (подробнее см. [Пожиленко и др., 2002] и ссылки там же) площадью 1300 км² и располагается в центральной части Кольского полуострова. Наиболее крупные и разрабатываемые в настоящее время Кировским филиалом (КФ) АО “Апатит” месторождения сосредоточены в юго-западной и юго-восточной частях рудного поля на участке протяженностью около 15 км (рис. 1).

Хибинский массив характеризуется разнообразием форм тектонических структур, которые по времени связаны с воздействием конических и радиальных разломов, относящихся к шести зонам [Онохин, 1975; Пожиленко и др., 2002]: Кукисвумчоррской, Поачвумчоррской, Суолуайвской, Ньоркпахкской, Коашкарской, Куэльпорской. Массив представляет собой сложную иерархично-блочную среду, разрушения в которой реализуются по наиболее слабым зонам, какими являются тектонические нарушения между блоками, заполненные шпреуштейнизированными водовмещающими породами. Разломы на месторождениях, близкие по ориентировкам к плоскостям действия максимальных кулоновских напряжений, могут быть наиболее сейсмоопасными [Сим и др., 2013]. Кроме того, при землетрясениях, зарегистрированных в пределах месторождений, происходили подвижки по имеющимся в массиве разрывным нарушениям, представленным окисленными породами.

Естественное поле напряжений в области разрабатываемых месторождений является результатом взаимодействия всех элементов блочной структуры. Хибинский массив является высоко тектонически напряженным, что установлено многочисленными натурными измерениями [Марков, 1977; Козырев и др., 2005, с. 123; Семенова, 2016] и подтверждено практикой ведения горных работ. Экспериментально обоснованные данные о действии в массиве, наряду с гравитационным, тектонического силового поля появились в 50–60-х годах. Величина максимальных тектонических напряжений (напряжений сжатия) горных пород в массиве достигает 40–60 МПа на глубинах в 200–600 м от дневной поверхности, что в 3–5, а в некоторых случаях в 20 раз превышает гравитационные напряжения [Марков, 1977], обусловленные весом налегающих пород. В качестве возможной причины высокого уровня напряжений горизонтального сжатия может служить механизм формирования остаточных напряжений гравитационного напряженного состояния, вызываемых денудацией [Ребецкий и др., 2017]. О продолжающемся тектоническом формировании этого района свидетельствуют современные поднятия массива со скоростью от 0.5 до 2–4 мм/год [Яковлев, 1982] и периодические землетрясения [Kremenetskaya, Triapitsin, 1995; Виноградов и др., 2016].

Хибинская ПТС включает в себя промышленные предприятия Апатитско-Кировского района, в состав которых входят объекты КФ АО “Апатит”: три рудника (Кировский, Расвумчоррский и Восточный), две фабрики по обогащению руд, хвостохранилища, энергетические и транспортные объекты; горно-обогатительный комбинат “Олений Ручей АО “Северо-Западная Фосфорная Компания”; объекты гражданского назначения.

Отработка апатитовых месторождений и выемка 4.5 млрд тонн руды и пустых пород в Хибинском массиве привела к повышению сейсмической активизации района, увеличению скорости деформирования массива на порядок в отдельных блоках [Сейсмичность …, 2002]. При отработке месторождений на рудниках КФ АО “Апатит” наблюдается весь спектр геодинамических явлений: шелушение, динамическое заколообразование, стреляние пород, микроудары и горные удары, толчки, горно-тектонические удары и техногенные землетрясения [Козырев и др., 2009; 2011; Жукова, Федотова, 2017; Zhukova et al., 2018]. В этих исследованиях выделяются следующие причины динамического проявления горного давления: природные (действие в массиве горизонтальных тектонических напряжений, высокие упругие свойства пород и наличие тектонических нарушений) и техногенные (ослабление межблоковых связей массива в сочетании с нарушением устойчивости блоков в результате ведения горных работ). Места пространственного группирования сейсмических событий в Хибинской ПТС связаны с зонами активного ведения горных работ и приурочены как к имеющимся, так и к формирующимся разрывным нарушениям в Хибинском массиве [Журавлева, 2017].

Таким образом, сейсмическая активность Хибинского массива есть результат совместного влияния целого ряда факторов: геодинамических (тектонические процессы), антропогенных (горные работы), гидрогеологических (обводненность массива пород), метеорологических (резкие перепады температуры) [Козырев и др., 2011; Жукова, Самсонов, 2014].

ИСХОДНЫЕ ДАННЫЕ

В исследовании использован каталог сейсмических событий, зарегистрированных сетью сейсмического мониторинга КФ АО “Апатит” [Корчак и др., 2014] за период с 1996 по июнь 2019 гг. (рис. 1). В настоящее время сеть состоит из 50 3-х компонентных сейсмических датчиков, расположенных на Кировском и Расвумчоррском рудниках с частотой дискретизации входных сигналов 1000 Гц. Мониторинговая сеть позволяет определять гипоцентры сейсмических событий с энергией E = 10³ Дж с точностью до 25 м в зоне повышенной точности и до 100 м в районе уверенной регистрации.

C 2010 г. мониторинг техногенной сейсмичности Хибинской ПТС ведется объединенной системой контроля сейсмичности массива с привлечением данных сети сейсмических станций Кольского филиала Единой геофизической службы РАН (КоФ ФИЦ ЕГС РАН). С этой целью в КоФ ФИЦ ЕГС РАН была разработана информационная система ЛОРС (локатор региональной сейсмичности).

При обработке данных сети сейсмического мониторинга КФ АО “Апатит” рассчитывается энергия события E, Дж. В статье пересчет энергии в магнитуду выполнялся по формуле Т.Г. Раутиан [1960] lg E (Дж) = 1.8M + 4.

Начиная с 1996 г., использованный каталог имеет представительную магнитуду не выше 0 (обозначим ее M_c), что соответствует энергии E_c ≤ 10⁴ Дж. Отметим, что такая представительность и точность расчета гипоцентров определяет уникальность использованных данных, позволяя проводить исследования для очень слабой сейсмичности, тем самым заполняя разрыв между лабораторными экспериментами и натурными наблюдениями. Это является дополнительной проверкой универсальности закономерностей, выявленных как в результате лабораторных исследований, так и в результате анализа глобальных и региональных каталогов тектонических землетрясений.

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

ВОЗМОЖНЫЕ ПОДХОДЫ К ИССЛЕДОВАНИЮ ПРОДУКТИВНОСТИ

Механизмы передачи напряжений и вероятностные модели широко исследовались для объяснения особенностей кластеризации землетрясений [Das, Scholz, 1981; Cocco et al., 2010; Hainzl et al., 2015; Смирнов и др., 2010; 2019], но эти подходы все еще далеки от того, чтобы выявить причинную связь между отдельными событиями.

Несмотря на разнообразие методов декластеризации, применявшихся ранее (см. обзорные статьи [Stiphout et al., 2012; Дещеревский и др., 2016; Шебалин, 2018; Писаренко, Родкин, 2019]), исследование причинно-следственных связей внутри каскадов инициированной сейсмичности все еще находится на начальном этапе и окончательной классификации пока не выработано. Первый подход состоит в том, чтобы отделить ветвящуюся структуру цепочек землетрясений от фоновой сейсмичности, используя итерационный алгоритм, связанный с оценками максимального правдоподобия эпидемической модели [Zhuang et al., 2002]. Другой способ состоит в том, чтобы идентифицировать прямо и косвенно инициированные события, предполагая линейный вклад каждого землетрясения в общий уровень сейсмичности без использования априорной модели [Marsan, Lengline, 2008].

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

СХЕМА РАСЧЕТОВ

В последующем анализе мы принимаем, что каждое землетрясение может инициировать несколько связанных с ним толчков, но каждый толчок может быть инициирован лишь одним определенным землетрясением. Продуктивностью события в такой схеме называется число землетрясений, инициированных этим событием-триггером. Мы будем вести подсчет числа инициированных событий магнитудой не ниже M относительно магнитуды основного толчка M_m (M ≥ M_m – ΔM). Для определения связей между событиями каталога мы воспользуемся методом ближайшего соседа [Zaliapin, Ben-Zion, 2013; 2016] с использованием функции близости [Baiesi, Paczuski, 2004]:

где: t_ij = t_j − t_i – время между событиями, которое положительно, если событие j происходит после события i и отрицательно в противном случае; r_ij ≥ 0 – пространственное расстояние между эпицентрами или гипоцентрами событий; m_i – магнитуда i-го события; b – параметр закона Гутенберга–Рихтера [Gutenberg, Richter, 1956]; d_f – фрактальная размерность распределения эпицентров или гипоцентров. Согласно работе [Писаренко, Родкин, 2019], эффективность метода ближайшего соседа для декластеризации каталога выше, чем эффективность оконных методов.

Для каждого события в каталоге землетрясение-триггер определяются по минимуму значения функции близости из всех предшествующих событий относительно рассматриваемого. Если соответствующее значение функции близости превышает заданный порог η₀, то связь разрывается, и оказывается, что данное событие не имеет “предка”.

Стандартная процедура определения значения η₀ предполагает, что распределение функции близости для ближайших соседей описывается смешанной моделью, состоящей из линейной комбинации двух гауссовых распределений [Zaliapin, Ben-Zion, 2013; 2016]. Порог η₀ находится из условия равенства плотностей с двумя оцененными гауссовыми модами. Модели смесей гамма-распределений и распределения Вейбулла можно рассматривать как альтернативу [Bayliss et al., 2019], которая может дать гораздо лучшее соответствие данным. В этом случае найденный порог η₀ зависит от модели. Отметим, что теоретические предпочтения между различными моделями отсутствуют.

Здесь для выбора порога η₀ используется модельно-независимый метод [Баранов, Шебалин, 2019]. Идея метода заключается в том, что распределение функции близости для некластеризованных землетрясений может быть смоделировано с использованием рандомизированного каталога, полученного путем случайного перемешивания времен событий относительно координат их гипоцентров и магнитуд [Писаренко, Родкин, 2019]. Рандомизированный таким образом каталог, полученный из полного, а не декластеризованного каталога, все еще может сохранять некоторую пространственно-временную кластеризацию. Чтобы преодолеть эту проблему, перед рандомизацией мы сначала грубым методом исключаем из каталога очевидные афтершоки. Полная процедура состоит из следующих шагов.

Шаг 1. Начальная декластеризация каталога. Используя только представительные данные (события c магнитудой выше представительной), для каждого события найдем его “предка” или “ближайшего соседа”, то есть предшествующее событие, имеющее минимальное значение функции близости относительно этого события. Построим распределение функции близости для ближайших соседей (серая кривая на рис. 2). Находим позицию правой моды распределения η_m и позицию половины высоты его правой ветви η_1/2 > η_m. В качестве порога для грубой декластеризации принимаем

Объединяя события, удовлетворяющие условию $\eta \leqslant \eta {\kern 1pt} '$, строим иерархические деревья попарно связанных событий. Из каждого иерархического дерева мы выбираем землетрясение с наибольшей магнитудой. Деревья могут состоять только из одного землетрясения.

Шаг 2. Перемешивание декластеризованного каталога. Для каждого времени землетрясения мы случайным образом выбираем из каталога гипоцентр и магнитуду. Строим распределение F_random(η) для ближайших соседей в полученном каталоге.

Шаг 3. Предполагая, что F_random(η) воспроизводит распределение для некластеризованных землетрясений, выполним декомпозицию распределения F_real(η) для ближайших соседей в реальном каталоге на две части [Zaliapin, Ben-Zion, 2013; 2016]:

Аналогичное соотношение справедливо для плотностей распределений p_real(η), pclustered(η), prandom(η). Толстая и тонкая черные кривые на рис. 2 показывают кластеризованную и случайную компоненты соответственно.

Чтобы оптимизировать вес κ, найдем наилучшее совпадение κp_random(η) с правой ветвью p_real(η) для значений η > η_4/5. Здесь η_4/5 > η_m такое, что p_real(η) = 4 p_real(η_m)/5.

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

Шаг 4. Находим значение η₀. Пороговое значение η₀ определяется как значение η, для которого доля кластеризованных событий с ближайшими соседями η > η₀ (ошибка I рода) равна доле некластеризованных событий с ближайшими соседями η ≤ η₀ (ошибка II рода) (рис. 2б):

Рисунок 2 иллюстрирует этот подход на примере используемого каталога и функции близости (1).

Альтернативой шагу 4 является рассмотрение вместо доли землетрясений для ошибок двух родов порогового значения η₁ такого, что:

Значение η₁ может быть более подходящим для декластеризации каталога при оценке сейсмической опасности, поскольку оно компенсирует число независимых землетрясений, которые рассматриваются как группированные эквивалентным числом группированных землетрясений, рассматриваемых как независимые [Molchan, Dmitrieva, 1992]. В подавляющем большинстве случаев полученные значения η₀ и η₁ близки друг к другу (см. следующий раздел, где в табл. 1 приведены значения η₀ и η₁, а также соответствующие им оценки средней продуктивности Λ_ΔM и $\Lambda _{{\Delta M}}^{'}$). Далее используются результаты, полученные при пороговом значении η₀.

РЕЗУЛЬТАТЫ РАСЧЕТОВ

Из используемого каталога за период с 1996 по июнь 2019 гг. 429 сейсмических событий c магнитудами М_m ≥ 1.5 были идентифицированы как события-триггеры. Они ассоциированы с 1177 инициированными событиями с относительными магнитудами M_r ≥ 1.5. Под относительной магнитудой понимается величина M_r = M_m – ΔM, где ΔM – относительный порог. Поскольку представительная магнитуда каталога M_c = 0, то использование таких относительных магнитуд является корректным. Распределение числа инициированных сейсмических событий вместе с экспоненциальным распределением и распределением Пуассона с одним и тем же параметром Λ_1.5 = = 1177/429 = 2.7 показано на рис. 3а. Сравнение эмпирических и теоретических распределений показывает, что продуктивность подчиняется экспоненциальному распределению

с плотностью а не распределению Пуассона, как обычно предполагается [Kagan, Knopoff, 1981; Ogata, 1989; Helmstetter, Sornette, 2002]. Формула (5) – закон продуктивности землетрясений.

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

Подчеркнем, что каждое иерархическое дерево кластеризации строится из основного события-триггера. На более низких уровнях иерархии инициированные события сами по себе являются вторичными событиями-триггерами и так далее до последних ветвей кластера. На рис. 3б показано, что распределения продуктивности инвариантно относительно иерархического уровня события-триггера, а среднее число инициированных событий Λ_ΔM убывает при увеличении уровня.

Распределение продуктивности остается экспоненциальным, когда относительный порог магнитуды ∆M меняется от 1.5 до 0.8 (рис. 4a). Как и предполагалось, средние значения Λ_ΔM уменьшаются в соответствии с параметром b закона Гутенберга–Рихтера (рис. 4б). Кроме того, при сохранении значения величины ∆M = 1.5, распределения числа инициированных событий и их средние значения Λ_1.5 практически одинаковы независимо от магнитуды землетрясений-триггеров M_m (рис. 5). Таким образом, экспоненциальное распределение продуктивности можно рассматривать как общее свойство всех сейсмических событий в рассматриваемой ПТС независимо от их силы.

Можно предположить, что, аналогично тектоническим землетрясениям продуктивность техногенной сейсмичности зависит от глубины очагов. На рис. 6 показано распределение продуктивности для различных глубин событий-триггеров. Приведенные расчеты показывают наличие вариаций величины Λ_1.5 с глубиной, меняющейся от –1000 до 500 м, считая от уровня Балтийского моря (отрицательные значения соответствуют глубинам выше уровня Балтики). При этом экспоненциальный вид распределения сохраняется для разных значений глубин событий-триггеров, характеризуемых различными значениями литостатического давления и горизонтальных напряжений, а также интенсивностью ведения горных работ и образованием трещин скола консоли налегающих пород [Сейсмичность…, 2002].

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

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

Экспоненциальный спад распределения продуктивности, наблюдаемый как для тектонических землетрясений по глобальным и региональным данным [Шебалин и др., 2018; Баранов, Шебалин, 2019], так и для техногенной сейсмичности по данным сети КФ АО “Апатит” для Хибинской ПТС, по-видимому, отражает общую прочность породы, независимо от величины возмущения, благодаря чему наиболее вероятное количество инициируемых событий равно нулю. Помимо индивидуального вклада каждого землетрясения в инициирование будущих событий, ветвящийся процесс в каскадах сейсмичности отвечает за свойства сейсмических кластеров, в том числе и в условиях техногенной сейсмичности. Как оказалось во всех рассмотренных случаях, распределение Пуассона неверно учитывает изменчивость числа землетрясений, вызванных каким-либо событием. Это опровергает уже сложившееся и широко используемое в моделировании сейсмического режима представление. В частности, как и в случае тектонических землетрясений [Баранов, Шебалин, 2019], закон продуктивности техногенной сейсмичности фактически опровергает эпидемическую модель ETAS [Ogata, 1989]. В основе этой модели лежит предположение о постоянстве числа инициируемых событий для данной магнитуды. Распределение Пуассона естественно моделирует отклонения от среднего для конкретных реализаций. На самом деле, как было установлено, число генерируемых событий распределено экспоненциально и, следовательно, в большинстве случаев меньше среднего значения. Именно поэтому модель ETAS, как правило, дает завышенные оценки опасности при прогнозе.

Неточность оценок с помощью модели ETAS может быть продемонстрирована на синтетическом каталоге сейсмических событий, сгенерированном по пространственной модели ETAS [Zhuang et al., 2002]. Соответствующая программа была любезно предоставлена Agnes Helmstetter, Institut des Sciences de la Terre, Гренобль, Франция. Число фоновых событий оценивалось по декластеризованному каталогу, и составило 4 события с M ≥ 0 в сутки. Остальные параметры взяты из таблицы. Продуктивность землетрясений по данным синтетического каталога (рис. 7) соответствует распределению Пуассона, а не экспоненциальному распределению, стабильно наблюдаемому по реальным данным. Такое расхождение реальных и модельных данных свидетельствует о необоснованности применения модели ETAS в существующих модификациях к оценке опасности повторных толчков.

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

Подтверждение закона продуктивности для техногенной сейсмичности является важным для оценки опасности повторных толчков. Используя распределения продуктивности и Гутенберга–Рихтера для магнитуд событий, а также, учитывая закон Омори–Утсу, можно рассчитать ожидаемое количество повторных толчков на заданном временном интервале и оценить магнитуду сильнейшего толчка и длительность опасного периода [Баранов и др., 2019; Шебалин, Баранов, 2019].

ЗАКЛЮЧЕНИЕ

Закон продуктивности (способность землетрясений вызывать последующие толчки) подтвержден для техногенной сейсмичности на примере Хибинской ПТС. Установлено, что распределение продуктивности имеет экспоненциальный вид независимо от магнитуды и глубины землетрясений. Полученный результат позволяет утверждать, что закон продуктивности, подобно двум другим законам сейсмологии (Гутенберга–Рихтера и Омори–Утсу), справедлив как для тектонических, так и техногенных землетрясений с энергией от 10⁴ Дж (M ≥ 0) Это свидетельствует об универсальном характере этого закона.

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