Физика Земли, 2021, № 5, стр. 91-103

ВВЕДЕНИЕ

Согласно современным представлениям основным механизмом землетрясений в верхней коре является скольжение по тектоническим разломам [Кочарян, 2016; Scholz, 2002]. Интенсивность и периодичность актов проскальзывания существенно варьируются и зависят от строения разломов, P–T-условий, напряженно-деформированного состояния и т.д. [Кочарян и др., 2014; Burgmann, 2018]. Долгое время считалось, что отдельный сегмент разлома может реализовывать накопленную энергию деформации либо посредством только повторяющихся медленных мод скольжения, которые сопровождаются излучением тремора или низкочастотных землетрясений, либо быстрых динамических срывов, приводящих к обычным землетрясениям, либо непрерывного крипа. Появляется все больше доказательств того, что на одном и том же сегменте разлома могут быть реализованы различные моды скольжения [Veedu et al., 2020; Ostapchuk et al., 2020b]. Незначительное изменение структуры зоны сместителя, нарушенности массива, флюидного режима и напряженно-деформированного состояния может привести к изменению режима скольжения [Kocharyan et al., 2021].

Считается, что активные разломы земной коры находятся в состоянии, близком к предельной прочности [Townend, Zoback, 2000]. При этом небольшие флуктуации напряженного состояния, связанные с самыми разными причинами, могут привести к возникновению динамической неустойчивости [Sobolev, 2011; Kocharyan et al., 2018b]. Переход участка разлома в состояние предельного равновесия проявляется, в том числе, в снижении сдвиговой жесткости [Кочарян, Остапчук, 2011], что характерно для различных динамических систем [Родкин, Рундквист, 2017]. При этом динамика скольжения определяется отношением эффективной величины жесткости разлома к жесткости вмещающего массива [Leeman et al., 2016]. Выявление области формирования очага землетрясения и контроль его развития является комплексной задачей, полного решения которой в ближайшее время ожидать не приходится. В настоящее время не существует инструментальных методов, позволяющих измерить действующие напряжения или статическую жесткость в зоне очага готовящегося землетрясения. Инструментальные наблюдения лишь иногда позволяют выявить косвенные проявления изменения режима скольжения [Frank et al., 2016; Kocharyan et al., 2018a; Ostapchuk et al., 2020b].

Лабораторные эксперименты являются одним из методов получения новой информации о закономерностях развития процессов скольжения. Начиная с работы [Brace, Byerlee, 1966], прерывистое скольжение контакта горных пород рассматривается многими исследователями как лабораторный аналог сейсмического цикла [Marone, 1998; Соболев, Пономарев, 2003, Кочарян, 2016; и др.]. Многие закономерности формирования и эволюции разломов исследуются путем контроля параметров АЭ [Соболев, Пономарев, 2003; Lherminier et al., 2019; и др.]. Эти методы можно считать определенным аналогом пассивного сейсмического мониторинга. Многие физические процессы, наблюдаемые в лабораторных экспериментах, статистически описываются соотношениями, подобными тем, что описывают сейсмический режим тектонического разлома (закон Гутенберга–Рихтера, закон Омори и др.) [Lei, 2003; Остапчук и др., 2016; Smirnov et al., 2019; Родкин, 2020]. Установлены многие закономерности проявления в регистрируемых колебаниях процессов разрушения [Соболев, Пономарев, 1999; 2003; Goebel et al., 2013; Пантелеев, 2020; и др.], инициирования подвижек внешними воздействиями [Kocharyan et al., 2018b], реакции системы на изменение напряженно-деформированного состояния [Смирнов, Пономарев, 2004; Bolton et al., 2020]. Исследованы микроструктурные процессы, сопровождающие формирование и эволюцию разломов и трещин [Веттегрень и др., 2017; Дамаскинская и др., 2018].

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

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

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

Ограничивая обсуждение относительно неглубокими разломами континентальной коры, можно полагать, что сдвиговые деформации при сейсмогенных подвижках локализуются в очень узкой центральной части разлома, которая сложена мелкодисперсным гранулированным материалом (англ. “gouge”); при этом закономерности скольжения по разлому длиной десятки/сотни километров определяются в значительной степени свойствами контакта частиц заполнителя разлома размером единицы–десятки–сотни микрон [Кочарян, 2016; Кочарян, Кишкина, 2020]. Узкая ширина слоев локализованного динамического сдвига делает эти объекты подходящими для исследования основных закономерностей их деформирования на лабораторных моделях.

Лабораторные эксперименты проводились на установке, в которой блок (1) под действием нормального и сдвигового усилий скользит вдоль поверхности раздела (2) (рис. 1). Контакт между поверхностями подвижного блока и основания заполнялся слоем гранулированного материала (3). Подвижный гранитный блок размером 8 × 8 × 3 см³ располагался в центре гранитного основания длиной 2.5 м и сечением 10 × 10 см². Контактирующие поверхности имели шероховатость глубиной 0.8–1.0 мм. Нормальное усилие составляло 500 Н и прикладывалось с помощью системы грузов. Сдвиговое усилие прикладывалось через пружинный элемент (6) с жесткостью 55 кН/м, конец которого смещали с постоянной скоростью V_s = 8 мкм/с. В ходе эксперимента регистрировались перемещения блока и действующее на блок сдвиговое усилие.

Рис. 1.

Схема установки геомеханического стенда ИДГ РАН. 1 – подвижный блок, 2 – гранитное основание, 3 – слой заполнителя, 4 – датчик перемещения, 5 – датчик силы, 6 – пружинный элемент, 7 – датчики АЭ.

В ходе экспериментов велась непрерывная регистрация АЭ с частотой дискретизации 1 МГц. Датчики АЭ VS30-V (7) (Vallen System, Germany) располагались на расстоянии 0.6 и 0.7 м от подвижного блока. Рабочий частотный диапазон датчиков АЭ 20–80 кГц. Пример записи АЭ представлен на рис. 2.

Рис. 2.

Зарегистрированная запись АЭ и ее спектрограмма в диапазоне частот 20–80 кГц. На врезке показан информационный критерий Акаике (AIC) и момент t_s вступления P-волны импульса АЭ.

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

Рис. 3.

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

Динамика скольжения разлома

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

(1)

$\delta X\left( t \right) = {{V}_{s}}t - x\left( t \right),$

где X(t) – положение блока относительно основания. По мере приближения к моменту динамического срыва скорость блока начинает постепенно увеличиваться. На стадии динамического проскальзывания скорость резко возрастает, а величина функции $\delta X\left( t \right)$ снижается до локального минимума. Наблюдается устойчивая картина вариаций параметров. Фазовый портрет явным образом демонстрирует регулярную повторяемость процесса (рис. 4). Траектория описывает предельный устойчивый цикл динамической системы, причем можно с высокой точностью определить области фазового пространства, соответствующие потере динамической устойчивости и началу динамического события. Для нерегулярного режима характерна более сложная картина изменения скорости скольжения. Наблюдаются как интервалы повторяющихся событий, так и интервалы квазистационарного скольжения, когда вариации дефицита перемещения $\delta X\left( t \right)$ относительно малы, а величина скорости близка к скорости протяжки V_s (рис. 3). Фазовый портрет такого режима представляет собой совокупность случайных круговых траекторий, охватывающих определенную область фазового пространства.

Рис. 4.

Временные вариации дефицита перемещения (δX) при регулярном (а) и нерегулярном (б) режимах скольжения. Справа представлены фазовые портреты реализованных режимов скольжения блока.

Вариация параметров АЭ

Скольжение модельного разлома сопровождается излучением большого количества импульсов АЭ. Для определения момента времени вступления Р-волны импульса использовался информационный критерий Акаике (AIC) [Carpinteri et al., 2012]. В основе метода лежит предположение о том, что сигнал АЭ можно разделить на отдельные интервалы, каждому из которых соответствует определенный авторегрессионный (АР) процесс [Sleeman, van Eck, 1999]. AIC используется для определения порядка АР процесса, описывающего временной ряд. Момент времени, где AIC принимает минимальное значение, соответствует оптимальному разделению на интервалы, характеризующие различные стационарные процессы, и интерпретируется как момент вступления P-волны. В двухинтервальной модели AIC сигнала рассчитывается для скользящего окна длиной N точек и является функцией момента k, разделяющей два интервала [Maeda, 1985]:

(2)

$\begin{gathered} AIC{\text{\;}}\left( k \right) = k\lg \left\{ {S\left( {A\left[ {1,k} \right]} \right)} \right\} + \\ + \,\,\left( {N - k - 1} \right)\lg \left\{ {S\left( {A\left[ {k + 1,N} \right]} \right)} \right\}. \\ \end{gathered} $

На врезке рис. 2 приведен пример временной вариации параметра AIC, который использовался для определения момента вступления P-волны импульса АЭ. AIC рассчитывался для окна длиной N = 500 точек. Для каждого выявленного импульса были определены амплитуда (A) и скорость нарастания амплитуды до максимума RA-value, которая определялась согласно соотношению:

(3)

$RA = \frac{{{{t}_{m}} - {{t}_{s}}}}{A},$

где t_s и t_m – моменты вступления P-волны и достижения максимальной амплитуды импульса.

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

Разломы земной коры как нелинейные динамические системы характеризуются свойством самоорганизованной критичности, и, соответственно, протекающие внутренние процессы самоорганизации характеризуются степенным законом распределения параметров [De Arcangelis et al., 2016]. На рис. 5 представлены распределения импульсов АЭ для регулярного режима скольжения.

Рис. 5.

Распределение импульсов АЭ по амплитуде (A) и скорости нарастания амплитуды до максимума (RA). Снизу представлены примеры импульсов АЭ с характерными значениями параметров (в скобках указано значение RA-value [c/B]).

На графике распределения импульсов АЭ по амплитуде выделяется две характерные области: A > 0.1 В характерна для импульсов АЭ, сопровождающих динамические срывы; A < 0.1 В – для импульсов АЭ, излучающихся на стадии подготовки динамических событий. Распределение по амплитуде в области амплитуд A < 0.1 В подчиняется степенному закону:

(4)

$\lg \left( N \right) = a - b\lg \left( {\frac{A}{{{{A}_{0}}}}} \right),$

где а, b – положительные константы и величина, $\lg \left( {\frac{A}{{{{A}_{0}}}}} \right)$ соответствует магнитуде импульса АЭ [Lei, 2003]. Параметр b является скейлинговым параметром, который характеризует процессы самоорганизации геосреды [Turcotte, 1999].

Распределение по параметру RA имеет более сложный вид. Выделяются две характерные области: при RA меньше характеристического значения ra₀, наблюдается “плато”, в то время как при RA ≥ ra₀ наблюдается степенной закон распределения. Стоит отметить, что величина ra₀ варьируется для различных условий проведения экспериментов и материального состава заполнителя. Для импульсов АЭ с RA < ra₀ характерна более высокая амплитуда и скейлинговый параметр b ≅ 1.3, для импульсов АЭ с RA ≥ ra₀ амплитуда импульсов существенно ниже и скейлинговый параметр b ≅ 2.0. Существенное различие скейлинговых параметров двух выборок может указывать на формирование двух динамических подсистем, эволюция которых сопровождается излучением импульсов АЭ с RA < ra₀ и RA ≥ ra₀ [Hadda et al., 2015; Ostapchuk, Morozova, 2020a].

На рис. 6 представлены временные вариации параметров АЭ. При регулярном прерывистом скольжении наблюдается периодичность изменения параметров. По мере приближения к моменту динамического срыва постепенно растет активность и доля высокоамплитудных импульсов.

Рис. 6.

Вариации параметров при регулярном прерывистом скольжении: δX – дефицит перемещения, определенный согласно соотношению (1); A и RA – амплитуда и скорость нарастания амплитуды; n – активность импульсов АЭ; b‑value – скелинговый параметр, определяемый согласно соотношению (4).

Для нерегулярного режима скольжения временные вариации параметров имеют более сложный вид (рис. 7). На стадии увеличения нагрузки активность АЭ практически не изменяется, резкий рост наблюдается только при переходе разлома в метастабильное состояние, когда скорость скольжения близка к V_s. После акта проскальзывания активность снижается до фонового уровня 10–15 c^–1. Особое внимание обращает на себя последовательность медленного и быстрого актов проскальзывания, произошедших во временном интервале 1030–1040 с. После медленного проскальзывания в окрестности t ~ 1030 c, несмотря на увеличение нагрузки, функция дефицита перемещения $\delta X\left( t \right)$ остается примерно постоянной – разлом скользит со скоростью, близкой к V_s – и в момент t ~ 1040 с происходит резкое проскальзывание. В отличие от остальных быстрых событий перед этим наблюдается снижение активности АЭ, что напоминает сейсмическое затишье перед землетрясениями [Wyss, Habermann, 1988]. Таким образом, при нерегулярном режиме скольжения можно наблюдать как рост, так и снижение активности АЭ перед актами проскальзывания.

Рис. 7.

Вариации параметров при нерегулярном режиме скольжения: δX – дефицит перемещения, определенный согласно соотношению (1); A и RA – амплитуда и скорость нарастания амплитуды; n – активность импульсов АЭ; b-value – скейлинговый параметр, определяемый согласно соотношению (4); δ – относительное количество импульсов АЭ кластера RA ≥ ra₀.

Выявить акустоэмиссионные признаки эволюции модельного разлома в условиях нерегулярного режима скольжения позволило разделение всего потока импульсов АЭ на два кластера. Такая обработка позволила наблюдать регулярные вариации доли импульсов АЭ с RA ≥ ra₀. По мере приближения к моменту динамического проскальзывания доля таких импульсов снижается, достигая минимума перед событием, а затем быстро возрастает.

Определение скорости скольжения и момента динамического события

В условиях, когда развитие процесса сопровождается сложной картиной вариаций параметров, применение методов анализа многомерных рядов данных, одним из которых является машинное обучение, позволяет определить параметры состояния исследуемой динамической системы. Для определения ожидаемой скорости скольжения блока и момента старта динамического события был применен алгоритм “случайный лес” [Breiman, 2001]. Программная реализация алгоритма выполнена в среде MATLAB. Алгоритм предполагает построение совокупности независимых деревьев принятий решений на основе набора параметров. Каждое дерево независимо определяет значение параметра, а результатом “случайного леса” является значение параметра, усредненное по всем деревьям. Количество деревьев в модели было принято равным 300, минимальная глубина каждого дерева составляла 13. Если одно из деревьев решений является ошибочным, то усреднение по совокупности независимых деревьев минимизирует данную ошибку. Для определения скорости скольжения была рассмотрена непрерывная запись сигнала АЭ, а для определения момента старта событий использовался сформированный каталог импульсов АЭ.

Применение алгоритма машинного обучения “случайный лес” предполагает, что на первом этапе выполняется обучение модели на известных данных. Обучение заключается в установлении связи между входными статистическими параметрами и выходными значениями относительной скорости смещения блоков и/или времени до момента динамического события (рис. 8). При обучении не учитывается информация о прошлом или будущем состоянии разлома, в модели используются данные АЭ только в пределах выбранного временного окна. Обучение алгоритма проводилось на данных двух экспериментов, в которых были реализованы регулярный и нерегулярный режимы скольжения. При обучении алгоритма определения скорости скольжения были использованы параметры сигнала АЭ, рассчитанные во временном окне 1 с с перекрытием 90%. В качестве параметров были взяты значения среднего и более высокие моменты, процентили от 1-го до 9-го и от 91-го до 99-го с шагом 1%, максимальные и минимальные значения амплитуды сигнала.

Рис. 8.

Обучение алгоритма “случайный лес” – предсказания скорости скольжения блока. Слева представлен пример одного из деревьев решений и параметры, на основании которых принимается решение в определенном узле дерева. Справа представлена экспериментальная зависимость изменения скорости (пунктирная линия) и предсказанные значения (сплошная линия). Коэффициент детерминации R² = 0.79.

На рис. 9 показан результат применения алгоритма восстановления скорости скольжения блока по данным регистрации акустического сигнала для последовательности событий скольжения, которые модель “случайного леса” никогда ранее не “видела”. По непрерывному сигналу АЭ достаточно успешно предсказывается скорость смещения блока как при регулярном, так и при нерегулярном режиме скольжения.

Рис. 9.

Определение скорости скольжения блока при регулярном (а) и нерегулярном (б) режимах. Пунктирная линия – реальные значения скорости, сплошная линия – значения, предсказанные моделью. Справа представлено сравнение истинных и предсказанных значений. Прямая линия соответствует точному определению скорости. При регулярном режиме точность (коэффициент детерминации) алгоритма составляет R² = 0.69, при нерегулярном – R² = 0.54.

При определении времени до динамического срыва был использован сформированный каталог импульсов АЭ. Статистические параметры, которые характеризуют поток АЭ, сопровождающий процесс деформирования, были рассчитаны в скользящем окне шириной 100 импульсов с шагом 20 импульсов. В качестве параметров были использованы среднее, дисперсия, асимметрия, эксцесс, максимальные и минимальные значения амплитуды, активность, b-value, а также все представленные параметры для каждого из кластеров. На рис. 10 представлены результаты определения времени до динамического события на основе алгоритма “случайный лес” для различных режимов скольжения. Для режима регулярного прерывистого скольжения алгоритм демонстрирует высокую эффективность с учетом малой выборки параметров, используемых в модели, и коэффициент детерминации алгоритма составляет R² = 0.73. На межсейсмической стадии время до срыва определяется с достаточно большой ошибкой 5–10 с, однако по мере приближения к моменту динамического срыва ошибка существенно снижается; наиболее вероятное значение ошибки определения времени динамического события составляет 0.4 с при длительности цикла порядка 30 c. В условиях нерегулярного режима скольжения алгоритм демонстрирует более низкую точность определения времени до срыва. Коэффициент детерминации R² = 0.20 и наиболее вероятная ошибка определения времени до динамического события составляет 2.9 с, однако при этом отчетливо наблюдаются сейсмические циклы лабораторного разлома.

Рис. 10.

Результаты определения времени до динамического события на основе алгоритма “случайный лес” при регулярном (а) и нерегулярном (б) режиме скольжения. Пунктирная линия – реальное время до срыва, сплошная – значение, предсказанное моделью. Справа представлено сравнение истинных и предсказанных значений. Прямая линия соответствует точному определению скорости. При регулярном режиме коэффициент детерминации R² = 0.73 при нерегулярном – R² = 0.22.

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

ОБСУЖДЕНИЕ И ВЫВОДЫ

Выполненные исследования демонстрируют возможность рассмотрения модельного разлома как двухкомпонентной динамической системы. Хотя в проведенных экспериментах отсутствовала возможность визуализировать протекающие на микроуровне процессы самоорганизации, выявленные вариации АЭ согласуются с результатами численных [Hadda et al., 2015; Gao et al., 2019] и лабораторных исследований [Hayman et al., 2011] деформирования слоев зернистого материала. Скольжение модельного разлома контролируется процессами формирования и разрушения в центральной части конгломератов нагруженных частиц так называемых силовых цепочек [Sammis et al., 1987]. Занимая относительно малый объем, данные образования являются естественными границами других, относительно разгруженных, областей.

Единственным источником информации о состоянии локального участка разлома являются упругие колебания, инициированные процессами самоорганизации на микро- и мезоуровне. При этом, возможно, основная доля информации содержится в сейсмическом шуме. Анализ непрерывных записей шума с применением новых методов обработки позволяет использовать информацию, полученную от источников в зоне разлома, физика которых не до конца известна [Ren et al., 2020]. Рассмотрение в лабораторных экспериментах статистических характеристик АЭ позволяет с неплохой точностью практически в режиме реального времени определить текущие параметры скольжения разлома, а также время и магнитуду готовящегося динамического события.

Хранение и высокопроизводительная обработка непрерывных высокочастотных записей сейсмических и/или акустических колебаний требует больших вычислительных мощностей. В этой связи подход, основанный на формировании и анализе каталогов сейсмических событий, выглядит более приемлемым. Однако классический сейсмический каталог, содержащий сведения о времени и магнитуде событий, целесообразно расширять, включая дополнительные параметры, характеризующие динамику косейсмического скольжения. В проведенных экспериментах при анализе нерегулярного режима скольжения особо информативным оказалось использование параметра RA-value. Вероятно, RA-value характеризует динамику межгранулярного скольжения аналогично тому, как это происходит при разрушении горных пород. В условиях развития процесса разрушения в образце горных пород скорость нарастания амплитуды уменьшается по мере развития трещины. Кроме того, при прогрессирующем разрушении доминирует низкочастотная компонента колебаний [Shiotani et al., 2001]. Вероятно, разрушение силовых цепочек сопровождается импульсами с резким вступлением, в то время как при изменении структурированности относительно разгруженных областей будут излучаться импульсы с большим RA-value [Michlmayr, Or, 2014; Gao et al., 2019].

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

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