ПРИМЕНЕНИЕ МЕТОДОВ МАШИННОГО ОБУЧЕНИЯ В ЗАДАЧЕ ОБНАРУЖЕНИЯ ПРЕДВЕСТНИКОВ ЗЕМЛЕТРЯСЕНИЙ ПО ДАННЫМ МОНИТОРИНГА СЕЙСМИЧЕСКИХ ШУМОВ В БАЙКАЛЬСКОЙ РИФТОВОЙ СИСТЕМЕ

1. ВВЕДЕНИЕ

В последние десятилетия в геофизических исследованиях наблюдается резкий рост неструктурированных массивов данных большого объема в результате существенного повышения плотности размещения и качества приборов. Вследствие этого появилась необходимость поиска эффективных способов обработки и анализа большого объема применительно к задачам геофизики. Быстрый прогресс в алгоритмах машинного обучения (МО) привлекает все большее внимание геофизического сообщества [Kong et al., 2019]. МО находит широкое распространение в сейсмических исследованиях, в частности в микросейсмическом мониторинге [Anikiev et al., 2023], в связи с возможностью обработки больших объемов данных в реальном времени. В нескольких обзорных статьях [Kong et al., 2019; Ridzwan, Yusoff, 2023; Mignan, Broccardo, 2020; Kubo et al., 2024] обобщаются приложения МО для выявления различных параметров микросейсмического режима и обнаружения скрытых закономерностей и признаков, связанных с подготовкой землетрясений.

Землетрясения являются одними из самых разрушительных стихийных бедствий, вызывающими значительные социальные и экономические потери, в связи с чем поиск решения задачи среднесрочного и краткосрочного прогноза является одной из наиболее актуальных практических целей современной сейсмологии. Несмотря на совершенствование измерительной сейсмологической аппаратуры и развитие технологий, точное прогнозирование сильных разрушительных землетрясений все еще невозможно. В ряде научных работ последних лет [Bletery, Nocquet, 2023; Seminsky et al., 2021, 2022; Bornyakov et al., 2023; и др.] приводятся примеры аномальных изменений геофизических параметров различной физической природы, которые могут выступать в качестве эффективных индикаторов краткосрочных предвестников землетрясений в пределах локального участка литосферы сейсмоактивной зоны.

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

В последние десятилетия с переходом на цифровую сейсмическую регистрацию активно развиваются методы анализа компонент волнового поля микросейсм, направленные на получение прогностических данных [Saltykov, 2017; Sobolev, 2004, 2011; Sobolev, Ponomarev, 2003; Sobolev et al., 2008; Lyubushin, 2011]. Хорошо зарекомендовала себя методика выделения предвестников землетрясений на основе анализа приливных микросейсм [Saltykov, 2017]. Примером наиболее успешного прогноза сейсмического события по вариациям энтропии сейсмического шума (с выявлением долгосрочных предвестников) является прогноз катастрофического землетрясения Тохоку (Япония) 11 марта 2011 г. с магнитудой М=9 [Lyubushin, 2011; Sobolev et al., 2008]. В частности, при оценке текущей сейсмической опасности в Японии последние измерения среднего значения энтропии показали ее максимальный уровень с 1997 г., что свидетельствует о высоком уровне сейсмической опасности согласно методу А.А. Любушина (персональное сообщение, сайт https://alexeylyubushin.narod.ru).

В Байкальской рифтовой системе (БРС) изучение вариаций микросейсмических колебаний и их возможной связи с процессами подготовки и реализации землетрясений также проводится в течение многих лет. Например, зафиксировано значимое уменьшение уровня микросейсмического шума на расстоянии 25 км за 1 ч перед Куморским землетрясением 16.09.2003 г. (энергетический класс К=14.3, моментная магнитуда Мw=5.6) с минимумом на частоте 6.6 Гц (обзор в [Korol et al., 2022]). В работе [Korol et al., 2022] для восьми событий с К=10.1–14.3 (М=3.4–5.7) было установлено понижение уровня амплитуд сейсмического шума в частотном диапазоне от 1 до 5 Гц за время от 12 мин до 6 ч 45 мин перед толчком. Для сильных землетрясений – Култукского 27.08.2008 г. (К=15.9, M=6.3) и Кударинского 09.12.2020 г. (К=14.6, Mw=5.6) – наблюдалась противоположная картина – увеличение уровня микросейсмического шума в области низких частот: в первом случае за 2 ч до главного толчка на частотах 0.5–1.4 Гц, во втором – за 10 дней на частотах 0.01–1 Гц [Dobrynina et al., 2023].

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

2. ИССЛЕДУЕМЫЙ РЕГИОН

Байкальская рифтовая система расположена в Северной Евразии и является одним из наиболее сейсмически опасных регионов Российской Федерации (рис. 1, а). Она состоит из ряда кайнозойских рифтовых впадин, разделенных горными хребтами и междувпадинными перемычками (рис. 1, б). Впадины ограничиваются крупными рифтоконтролирующими разломами преимущественно сбросового кинематического типа [Logatchev, Florensov, 1978] (рис. 1, б). В тектоническом смысле БРС приурочена к границе двух литосферных плит – стабильной Североевразийской и мобильной Амурской, относительная скорость движения между которыми по данным GPS-геодезии оценивается в 3.4±0.7 мм/год [San'kov et al., 1999]. Основная масса землетрясений БРС сосредоточена в достаточно узкой полосе, протягивающейся вдоль оси рифтовых впадин. Сибирская платформа практически асейсмична, Забайкальский блок характеризуется рассеянной умеренной сейсмичностью (рис. 1, а). В год в пределах БРС регистрируется до 10 тыс. землетрясений, большинство из которых относятся к слабым и умеренным событиям. Последние сильные землетрясения произошли в 2020–2021 гг. и локализовались на юго-западном фланге рифтовой системы (Быстринское землетрясение 21.09.2020 г., Mw=5.6 [Seminsky et al., 2021], Хубсугульское землетрясение 11.01.2021 г., M=6.7 [Emanov et al., 2023]) и Кударинское землетрясение 09.12.2020 г. в Южнобайкальской впадине Mw=5.6 [Dobrynina et al., 2022]).

3. МЕТОДИКА И ДАННЫЕ

В основе используемой в настоящей работе методики поиска краткосрочных предвестников землетрясений лежит оценка вариаций спектральной плотности мощности (СПМ) микросейсмических шумов в ближней к очагу зоне до и после землетрясения. Этот подход основан на гипотезе о том, что существует значимый скачок в напряжениях и деформациях земной коры перед землетрясением. Многочисленные исследования последних лет подтверждают возможность этого явления [Bletery, Nocquet, 2023; Bornyakov et al., 2023; Dobrynina et al., 2023; Korol et al., 2022; и мн. др.].

В самой Байкальской впадине наиболее сейсмически активной является ее южная часть (рис. 2, а), на прилегающих к ней территориях Южного Прибайкалья и Забайкалья (Иркутская область, Республика Бурятия) находятся крупные города, населенные пункты, а также опасные производства, линейные сооружения, плотина ГЭС (г. Иркутск), автомагистрали и линии Восточно-Сибирской железной дороги. В 2020 г. на западном побережье оз. Байкал в местечке Куяда был открыт пункт комплексного мониторинга опасных геологических процессов Института земной коры СО РАН «Бугульдейка», оснащенный современным геофизическим оборудованием, в том числе широкополосной сейсмической станцией (см. рис. 1, в). Пункт мониторинга расположен вблизи зоны влияния крупнейших разломов центральной части рифта – Приморского и Морского (см. рис. 1, в). Сейсмическая станция «Куяда» (региональный код KYD) оснащена трехкомпонентным широкополосным сейсмометром Trillium Compact Posthole 120 c и цифровым регистратором Centaur, фирма Nanometrics. Станция работает в непрерывном режиме, интервал дискретизации – 100 отсчетов в секунду, рабочая полоса частот от 0.0083(3) до 108 Гц. Регистрируются три компоненты вектора скорости сейсмических колебаний: север – юг, восток – запад и вертикальная.

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

Для анализа были выбраны 40 землетрясений различного энергетического диапазона (K=9.5–14.9, М=3.1–5.9), зарегистрированных широкополосной сейсмической станцией KYD на эпицентральных расстояниях от 10 до 120 км (рис. 2). Указанные землетрясения в основном локализованы в зоне влияния Обручевского, Морского и Приморского разломов (кинематический тип – сброс), ограничивающих Южнобайкальскую впадину с запада, небольшая часть землетрясений сосредоточена в пределах восточного борта впадины (см. рис. 1, а). Выбранная область характеризуется как область чистого растяжения (pure extension zone): ось напряжений сжатия близвертикальна, оси напряжений растяжения ориентированы вкрест основным рифтовым структурам [Radziminovich et al., 2019] (см. рис. 1, б). Радиус 120 км охватывает сейсмические события, происходящие в зонах влияния этих разломов, предоставляя полный набор данных для анализа закономерностей и поиска предвестников будущих землетрясений, которые могут повлиять на регион. Тестовые испытания модели МО с разными наборами данных показали, что радиус 120 км обеспечивает оптимальный баланс, гарантируя достаточность данных без ущерба для релевантности и точности модели, в то время как выбор большего радиуса может внести дополнительные помехи, включив данные из областей с различными сейсмическими характеристиками, что потенциально снизит точность прогнозирования наших моделей.

Для оценки СПМ использовался метод периодограмм Уэлча на 10-минутных участках записи микросейсмического шума. Параметры периодограмм: длительность сегментов 10 секунд с перекрытием 50 %, значения СПМ медианные. Это позволило получить оценки СПМ с шагом 0.1 Гц.

Чтобы обеспечить хорошее качество модели, необходимо эффективно сегментировать исторические данные мониторинга на периоды сейсмической активности и сейсмической неактивности. Согласно ранее полученным результатам [Korol et al., 2022; Dobrynina et al., 2023] можно сделать вывод, что для одиночных землетрясений, не предваряемых и не сопровождающихся форшоковыми и афтершоковыми событиями, время проявления значимых изменений в микросейсмических шумах не превышает нескольких часов до и после толчка. На этой основе в работе анализировались записи МСШ в указанный период времени перед землетрясениями. Количество образцов, временной интервал между ними и частотный диапазон подбирались экспериментально по наибольшей точности модели МО, для обучения которой эти данные использовались. Полученные оценки СПМ сравнивались с оценками, полученными в спокойные периоды – не менее чем через 12 часов после землетрясения. Для анализа использовался вертикальный канал сейсмометра, так как на нем менее проявлены техногенные, ветровые, озерные и другие помехи.

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

Исследование выполнялось с помощью цифровой платформы (ЦП) (https://izk.sscc.ru) для комплексных геофизических исследований в Байкальском регионе [Braginskaya et al., 2023]. На платформу поступают данные с пунктов комплексного мониторинга опасных геологических процессов ИЗК СО РАН, в том числе данные с сейсмостанции «Куяда». В составе ЦП имеются интерактивные средства для визуализации и анализа данных в режиме онлайн. В частности, можно оценивать СПМ сейсмических записей по алгоритму Уэлча и отображать результаты в виде диаграмм. Первоначально анализ и сравнение различных участков записей МСШ выполнялись с помощью ЦП в интерактивном режиме, позже были написаны программы на языке Python для автоматического отбора и анализа данных. Построение интерактивных карт с эпицентрами землетрясений обеспечивает входящая в состав платформы ГИС.

4. РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЯ И ИХ ОБСУЖДЕНИЕ

В настоящей работе проблема выявления предвестников землетрясений в данных сейсмического мониторинга формулируется как проблема бинарной классификации с категориями активных и неактивных участков сейсмического шума. Активным участком МСШ считается участок, непосредственно предшествующий землетрясению, неактивным – отстоящий от землетрясений не менее чем на 12 час. Нами использовались контролируемые модели обучения, т.е. обучение на размеченных образцах (активный или неактивный участок).

Чтобы получить представление о характеристиках и закономерностях в записях микросейсмического шума, нами был проведен разведывательный анализ данных (EDA). Согласно полученным результатам, за несколько часов перед землетрясениями наблюдается статистически значимое увеличение уровня МСШ в частотном диапазоне 0.5–2.5 Гц, в то время как на более высоких частотах – от 4 Гц до 16 Гц – картина обратная – уровень МСШ понижается. На остальных частотах уровень МСШ для активных и неактивных участков одинаковый. Это коррелирует с данными, полученными в работах [Korol et al., 2022; Dobrynina et al., 2023]. На рис. 3 представлены усредненные нормализованные спектры МСШ перед землетрясениями относительно спокойного фона.

Для каждого землетрясения строились спектры МСШ для пяти активных образцов (за 250, 200, 150, 100 и 50 мин до события) и пяти фоновых. Таким образом, результирующие спектры получены усреднением 200 спектров активных образцов и 200 неактивных. Нормализация осуществлялась путем вычитания среднего значения и нормирования на стандартное отклонение. На рис. 3 также показаны 95-процентные доверительные интервалы. Как можно видеть, разница СПМ между активными и фоновыми образцами МСШ наблюдается на низкочастотной части диапазона, приблизительно до 2.5 Гц, и на более высоких частотах, но с обратным знаком. Опытным путем было установлено, что наибольшая точность модели бинарной классификации достигается при использовании для ее обучения значений СПМ из диапазона 0.8–2.4 Гц с шагом 0.1 Гц. Аналогично, по наилучшему качеству модели подбиралось приведенное выше количество образцов МСШ и интервал между ними для каждого события. Спектры МСШ для выбранного диапазона частот отдельно приведены на рис. 4.

Наблюдаемое увеличение МСШ перед землетрясениями может быть объяснено в аспекте тектонофизической модели очага тектонического землетрясения [Bornyakov et al., 2024]. Положения модели, включая формирование предвестников, основаны на выявлении общих закономерностей структурно-динамического развития процессов формирования зон разломов в результате физического моделирования. В качестве такого предвестника выступает процесс деформационной самоорганизации активных сегментов в очаге готовящегося землетрясения. В работе [Bornyakov et al., 2024] показано, что этот процесс поддается превентивной диагностике по данным деформационного и микросейсмического мониторинга и проявляется в виде низкочастотных автоволновых колебаний непосредственно перед сейсмическим событием во временном интервале от первых дней до первых часов.

На основе тектонофизической модели очага [Bornyakov et al., 2024] и результатов исследования вариаций МСШ, полученных в работах [Korol et al., 2022; Dobrynina et al., 2023], наблюдаемые изменения спектрального состава МСШ перед близкими землетрясениями могут быть интерпретированы как отражение процесса подготовки землетрясения и связанного с ним изменения свойств среды. Консолидация (сжатие) среды перед толчком приводит к закрытию мелких трещин в области, окружающей очаг готовящегося землетрясения, и, как следствие, уменьшению уровня МСШ на высоких частотах – от 4 до 16 Гц, в то время как непосредственно в очаге перед толчком постепенно начинается самоорганизация активных сегментов, проявляющаяся в виде медленных, сходных с крипом, движений, генерирующих низкочастотные колебания, выражающиеся в виде усиления низкочастотной части спектра в диапазоне 0.5–2.5 Гц. Оба этих явления могут классифицироваться как краткосрочный предвестник близкого умеренного или сильного землетрясения.

Для автоматического распознавания, относится запись микросейсмического шума к неактивному или активному периоду, методом машинного обучения была построена модель бинарной классификации. В качестве независимых переменных (признаков) для обучения и тестирования моделей использовались значения СПМ в диапазоне 0.8–2.4 Гц с шагом 0.1 Гц, всего 17 значений. 400 образцов СПМ микросейсмического шума для обучения и тестирования делились в отношении 3:1. Обучение проводилось на моделях KNN, Gradient Busting и Random Forest из библиотеки Sklearn (https://scikit-learn.ru). Методом кросс-валидации установлено, что наилучшие результаты показала модель Random Forest («Случайный лес»). Матрица ошибок и основные метрики для этой модели приведены на рис. 5.

Матрица ошибок отображает статистику, полученную при тестировании модели. Тестирование производилось на ста образцах МСШ, не участвовавших в обучении. В двух верхних клетках приведено количество неактивных образцов, правильно определенных – 45 и неправильно определенных – 7. В нижних клетках количество активных образцов, правильно определенных – 37 и неправильно определенных – 11. Метрика accuracy равна доле правильно классифицированных образцов, precision можно интерпретировать как долю образцов, названных моделью активными и при этом действительно являющихся активными, а recall показывает, какую долю активных образцов из всех активных образцов нашла модель.

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

5. ЗАКЛЮЧЕНИЕ

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

В результате проведенных исследований получила развитие методика обнаружения признаков подготовки землетрясений на примере близких умеренных и сильных сейсмических событий по данным о вариациях спектрального состава МСШ. На примере 40 землетрясений в БРС с энергетическим классом К=9.5–14.5 на эпицентральных расстояниях от 10 до 120 км выявлены предвестники за несколько часов перед толчком. На основе полученных данных методом машинного обучения создана модель для бинарной классификации записей МСШ на активные и фоновые. Зафиксированные изменения спектрального состава МСШ могут классифицироваться как краткосрочный предвестник, связанный с консолидацией (сжатием) среды и с деформационной самоорганизацией активных сегментов в очаге готовящегося землетрясения за 3–4 ч до толчка.

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

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

7. ЗАЯВЛЕННЫЙ ВКЛАД АВТОРОВ / CONTRIBUTION OF THE AUTHORS

Все авторы внесли эквивалентный вклад в подготовку рукописи, прочли и одобрили финальную версию перед публикацией.

All authors made an equivalent contribution to this article, read and approved the final manuscript.

8. РАСКРЫТИЕ ИНФОРМАЦИИ / DISCLOSURE

Авторы заявляют об отсутствии конфликта интересов, связанного с этой рукописью.

The authors declare that they have no conflicts of interest relevant to this manuscript.