Синтетический каталог землетрясений Восточного сектора Арктической зоны Российской Федерации

Полный текст

ВВЕДЕНИЕ

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

Параллельно, с развитием системы страхования от последствий землетрясений, возникла необходимость введения вероятностных оценок силы воздействия от землетрясений в заданных точках или на заданных объектах. Так возник вероятностный подход к оценке сейсмической опасности. В России такой подход получил сокращение ВАСО (вероятностный анализ сейсмической опасности), в международной литературе – PSHA (probabilistic seismic hazard assessment). Несмотря на критику этого подхода со стороны большого числа ученых, оценки сейсмической опасности, рассчитанные с помощью этого подхода, остаются нормативными в большинстве стран, в которых есть сейсмоопасные регионы, включая Россию.

Более того, этот подход повсеместно используется именно в том виде, в котором он был разработан в 1960–1970 гг. и внедрен по всему миру в 1990-е гг. в рамках программы GSHAP (Global Seismic Hazard Assessment Program) [Giardini et al., 1999] и затем GEM (Global Earthquake Model). По сути, это лишь надстройка над детерминистским подходом с добавлением элементов оценки вероятности сейсмических проявлений на заданном интервале времени. Как и в детерминистском анализе, этот подход включает решение двух независимых задач: 1) определение потенциальных очагов будущих землетрясений (в вероятностном анализе с учетом их средней повторяемости), 2) построение модели затухания сейсмических волн от очага землетрясения до конкретного объекта или точки на поверхности Земли.

По-прежнему основными сейсмогенерирующими структурами считаются активные разломы. Фактор времени добавлен путем включения главного элемента ВАСО – закона Гутенберга–Рихтера [Gutenberg, Richter, 1944], который устанавливает подобие числа землетрясений в единицу времени для землетрясений разной магнитуды. Этот закон позволяет по статистике менее сильных землетрясений оценивать повторяемость более сильных.

Фактически любые модели затухания могут применяться как в ВАСО, так и в ДАСО. Первоначально использовались простейшие модели макросейсмического поля Блейка–Шебалина [Шебалин, 2003], которые устанавливают соотношение между интенсивностью сотрясений в баллах (шкалы Mercalli, MSK-64) и магнитудой, глубиной очага и эпицентральным расстоянием. Одновременно с макросейсмической шкалой стали применяться инструментально измеряемые величины (пиковое ускорение грунта, PGA, Peak Ground Acceleration; пиковая скорость грунта, PGV, Peak Ground Velocity). В современных условиях для инженерных целей помимо пиковых значений часто используются спектральные параметры, позволяющие рассчитать отклик сооружений на воздействия сейсмических волн разной частоты.

Быстрый рост объемов и качества наблюдательной базы, совершенствование вычислительной техники привели за последние десятилетия к значительному прогрессу в решении второй задачи: воздействия на объекты от землетрясений с конкретными параметрами можно рассчитать достаточно точно. А вот решение первой задачи – определения мест и повторяемости будущих землетрясений разной силы – по-прежнему часто приводит к значительным ошибкам. Многие сильнейшие землетрясения происходят в неожиданных местах [Wyss et al., 2012; Шебалин, 2024], но при этом на больших территориях оценки вероятной силы воздействия от землетрясений за определенный период оказываются завышенными [Шебалин и др., 2022; Шебалин, 2024].

В России в настоящее время используется комплект карт Общего сейсмического районирования ОСР-2015, которые являются нормативными при строительстве в сейсмоопасных районах. На картах обозначены зоны нормативной балльности, рассчитанной для разных интервалов времени. Карта А представляет зоны ожидаемой в течение 50 лет интенсивности воздействия землетрясений по шкале МСК-64 с вероятностью 10%, карта В – с вероятностью 5%, карта С – c вероятностью 1%. Карты ОСР-2015 в значительной степени воспроизводят карты ОСР-97 [Комплект …, 1999], впервые построенные в рамках ВАСО. Главные отличия, помимо включения территории Крыма, состоят в увеличении нормативной балльности в восточной части Российской Арктики: на севере Якутии и на севере Камчатки. Эти изменения были вызваны Олюторским землетрясением 2006 г. в Корякии и Иллин-Тасским землетрясением 2013 г. на севере Якутии. Оба землетрясения можно считать пропусками на карте ОСР-97, поскольку их воздействие оказалось выше ожидаемого примерно на 2 балла [Шебалин, 2024]. Частично эти пропуски можно оправдать недостаточной изученностью малонаселенных районов, в которых произошли землетрясения. Помимо карт ОСР-97 и ОСР-2015 короткое время действовал также комплект карт ОСР-2016 [Пояснительная ..., 2016]. Для рассматриваемого региона эти карты отличаются как от карт ОСР-97, так и от карт ОСР-2015 (рис. 1). Нормативная балльность в регионе несколько снижена по сравнению с ОСР-2015. Тем не менее отличия трех поколений карт недостаточны, чтобы компенсировать общее завышение нормативной балльности на больших территориях, обнаруженное в работах [Шебалин и др., 2022; Шебалин, 2024]. Таким образом, вопрос о более точной оценке сейсмической опасности в регионе остается актуальным.

В основу трех поколений вероятностных карт ОСР легла линеаментно-доменно-фокальная (ЛДФ) модель зон ВОЗ (возможных очагов землетрясений) [Ulomov, 1999]. Идеология моделей типа ЛДФ была предложена Корнеллом в 1968 г. [Cornell, 1968]: комбинация линейных и площадных структур позволяла ему получать аналитические выражения для перехода от модели зон ВОЗ к модели сотрясаемости. Позднее появилась возможность замены аналитических выражений Корнелла расчетами с помощью метода Монте-Карло. При таком подходе формируется синтетический каталог землетрясений, реализующий повторяемость землетрясений разной магнитуды в соответствии с параметрами закона Гутенберга–Рихтера, установленными для элементов ЛДФ-модели, а величина нормативной интенсивности определяется путем многократной реализации синтетического каталога на таких интервалах и пересчетом по формулам затухания.

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

Модели типа ЛДФ давно подвергаются критике, главным образом из-за их субъективного характера [Gerstenberger et al., 2020]. Альтернативой моделям типа ЛДФ являются модели сглаженной сейсмичности [Akinci et al., 2018; Helmstetter, Werner, 2012; Stock, Smith, 2002], учитывающие закон Гутенберга–Рихтера. Этот подход обеспечивает гораздо более точные в пространственном отношении оценки параметров по сравнению с моделями типа ЛДФ, но требует значительно большего количества зарегистрированных землетрясений в широком диапазоне магнитуд для построения детальной и точной модели сейсмического режима [Ризниченко, 1980]. Главная проблема в таких моделях – “размазывание” сейсмической активности, почти неизбежное при использовании процедуры сглаживания. В результате в наиболее опасных местах активность оказывается заниженной, а на краях активных зон завышенной.

Проблема усугубляется тем, что при нормировании активности на площадь, как правило, не учитывается фрактальная структура пространственного распределения эпицентров [Kossobokov, Mazhkenov, 1994; Spada et al., 2011]. Значительный прогресс в построении моделей сглаживания был достигнут с появлением высококонтрастного метода среднего положения (МСП) [Vorobieva et al., 2024; Shebalin et al., 2024], разработанного в Институте теории прогноза землетрясений и математической геофизики РАН (ИТПЗ РАН) в рамках проекта РНФ № 20-17-00180 “Развитие сценарного подхода в задачах оценки сейсмической опасности и риска”. В отличие от практически всех моделей сглаживания сейсмичности, в МСП значения параметров приписываются среднему положению эпицентров землетрясений, использованных в расчетах (а не центру круга), а площадная нормировка производится с учетом фрактальности.

В расчетах ОСР при оценке параметров ЛДФ-модели по реальному каталогу землетрясений афтершоки исключались из рассмотрения. Оставшиеся землетрясения обычно называются фоновыми. Такой подход основан на гипотезе, что повторяемость фоновых землетрясений в окрестности каждой точки в рассматриваемом регионе не изменяется во времени [Gardner, Knopoff, 1974]. Афтершоки концентрируются в пространстве и времени вблизи отдельных землетрясений и их число обычно сопоставимо с числом фоновых землетрясений. Использование полного каталога, без удаления афтершоков (так называемой процедуры декластеризации), неизбежно приводило бы к локальному завышению опасности вблизи эпицентров сильных землетрясений.

С другой стороны, игнорирование афтершоков может приводить к недооценке сейсмической опасности [Шебалин, 2024]. В рамках проекта РНФ № 20-17-00180 в ИТПЗ РАН была предложена “симметричная” процедура декластеризации каталога и моделирования последовательности афтершоков в синтетическом каталоге землетрясений [Shebalin et al., 2024]. Декластеризация осуществляется по методу Заляпина [Zaliapin et al., 2008; Zaliapin, Ben-Zion, 2013], в котором строятся цепочки афтершоков разных поколений, реализующих модель ETAS [Zhuang et al., 2002]. При создании синтетического каталога сначала генерируются фоновые события, а затем по модели ETAS генерируются афтершоки. Используется вариант модели ETAS-e, в котором учитывается экспоненциальный закон продуктивности [Baranov et al., 2019; Shebalin et al., 2020]. Применение метода Заляпина для декластеризации каталога и модели ETAS-e в значительной степени компенсирует возможные ошибки декластеризации [Baranov et al., 2022]. В случае, если число фоновых событий при декластеризации занижено, оказывается завышенным число афтершоков. Это приводит к увеличению средней продуктивности землетрясений, благодаря чему при моделировании по модели ETAS-e генерируется большее число афтершоков, частично компенсирующее недостаток фоновых событий.

В зонах высокой сейсмичности главную опасность представляют сильные землетрясения. Они происходят очень редко, но их воздействие распространяется на большие территории. Одна из важнейших задач в моделях сглаженной сейсмичности состоит в том, чтобы максимально сузить области, в которых могут возникать сильные землетрясения [Gerstenberger et al., 2020]. Без решения этой задачи зона возможного возникновения сильных землетрясений “размазывается”, вследствие чего в зонах, где сильные землетрясения действительно могут возникнуть, их ожидаемая повторяемость оказывается заниженной, а по краям этих зон завышенной. Один из возможных способов сужения областей, в которых ожидаются сильные землетрясения, – использование части модели ЛДФ, определяющей линейные структуры. Этот подход, однако, остается весьма субъективным. Кроме того, он неизбежно приведет к ошибкам типа пропуска цели, поскольку сильные землетрясения довольно часто происходят далеко от уже известных активных разломов и тем более от линеаментов, определяющих в ЛДФ-модели места возможных сильных землетрясений. Другой возможный подход – это использование результатов распознавания мест возможного возникновения сильных землетрясений по методике EPA (Earthquake Prone Areas) [Гвишиани и др., 2020]. Именно этот подход было предложено использовать в методике создания синтетического каталога землетрясений, разработанной в ИТПЗ РАН в рамках проекта РНФ № 20-17-00180 [Shebalin et al., 2024]. Выделенные и распознанные по геолого-геофизическим признакам морфоструктурные узлы определяют в синтетическом каталоге возможные места сильных землетрясений. Эпицентры таких землетрясений сосредоточены вблизи центров узлов, но могут удаляться в пределах 30 км. Распределение направлений при этом не однородно, а определяется геометрией близлежащих активных разломов, которые находятся по базе данных Геологического института РАН (ГИН РАН) [Zelenin et al., 2022]. В данной работе мы реализуем третий возможный подход, используя результаты FCAZ-распознавания [Dzeboev et al., 2021] мест возможного возникновения сильных землетрясений в рассматриваемом регионе [Gvishiani et al., 2022b].

ПОСТАНОВКА ЗАДАЧИ И КЛЮЧЕВЫЕ ПРИНЦИПЫ ЕЕ РЕШЕНИЯ

В данной работе мы ставим задачу построения синтетического каталога землетрясений в Восточном секторе Арктической зоны Российской Федерации (АЗРФ) на условный период 20 000 лет, воспроизводящего основные свойства реального каталога. К таким свойствам относятся:

• региональная повторяемость (рассчитанная для всего региона в единицу времени) фоновых землетрясений во всех интервалах магнитуды;
• региональное соотношение числа афтершоков и фоновых событий;
• локальная повторяемость (оцененная для малой окрестности каждой точки рассматриваемой территории на единицу площади в единицу времени) фоновых землетрясений (без учета афтершоков);
• локальные соотношения числа фоновых землетрясений разной магнитуды (параметр b закона Гутенберга–Рихтера);
• локальные свойства группирования землетрясений;
• пространственное распределение эпицентров сильных землетрясений.

Методика построения синтетического каталога землетрясений с такими свойствами была недавно разработана в ИТПЗ РАН в рамках проекта РНФ № 20-17-00180 “Развитие сценарного подхода в задачах оценки сейсмической опасности и риска” [Shebalin et al., 2024]. Ключевые принципы этого подхода следующие:

• на первом этапе генерации синтетического каталога формируется каталог фоновых событий с магнитудой $M\ge M_0$ (порог генерации выбирается исходя из целей создания синтетического каталога);
• время каждого фонового события генерируется в соответствии с распределением Пуассона с параметром, равным региональной повторяемости числа событий с $M\ge M_0$;
• магнитуда каждого события генерируется в соответствии с региональным магнитудно-частотным распределением [Shebalin et al., 2024];
• координаты эпицентра каждого фонового события генерируются по модели распределения локальной повторяемости землетрясений данной магнитуды; если магнитуда события превышает заданный порог $M_{strong}$, то используется модель пространственного распределения сильных землетрясений [Shebalin et al., 2024];
• для каждого фонового землетрясения по модели ETAS-e [Shebalin et al., 2020] генерируется иерархическая последовательность афтершоков; такая последовательность может быть пустой и не содержать ни одного события.

Для проверки выполнения в синтетическом каталоге свойств каталога фактически зарегистрированных землетрясений проводится ряд тестов.

МЕТОД

Здесь мы используем методику работы [Shebalin et al., 2024]. Она включает: 1) разделение каталога землетрясений на фоновые события и афтершоки; 2) построение высококонтрастной модели параметров сейсмического режима фоновых землетрясений; 3) оценку параметров афтершоков: региональную оценку параметров c и p закона Омори–Утсу и параметра b закона Гутенберга–Рихтера, а также локальные оценки средней продуктивности; 4) генерацию синтетического каталога фоновых событий; 5) генерацию последовательностей афтершоков и формирование полного синтетического каталога; 6) верификацию синтетического каталога. Важным элементом этой методики является ограничение возможных мест сильных землетрясений. В работе [Shebalin et al., 2024] предлагается использовать для этого результаты распознавания по методике EPA [Гельфанд и др., 1972; Гвишиани и др., 2020]. В данной работе ограничения вводятся в соответствии с результатами FCAZ-распознавания [Gvishiani et al., 2022b].

Первые два и последний элемент были описаны и реализованы для рассматриваемого региона в работе [Воробьева и др., 2024], где было приведено детальное обоснование применяемых методик. Здесь мы приведем их сокращенное описание, остановившись подробнее на остальных элементах.

Разделение каталога землетрясений на фоновые события и афтершоки

Используемый здесь метод ближайшего соседа Заляпина–Бен-Зиона [Zaliapin et al., 2008; Zaliapin, Ben-Zion, 2013] для декластеризации каталога по сути является обратным относительно метода генерации афтершоков в синтетическом каталоге по модели ETAS-e. В методе устанавливается иерархия последовательности афтершоков, в которой однозначно выделяются афтершоки фоновых событий, афтершоки афтершоков, их афтершоки и т.д. В модели ETAS-e такие иерархические последовательности синтезируются. Таким образом обеспечивается максимальное соответствие “афтершоковой” части в синтетическом и реальном каталогах.

В методе Заляпина–Бен-Зиона выделяются пары связанных событий, одно из которых считается триггером (или “родителем”), а другое инициируемым событием (или “потомком”); потомок всегда возникает позже родителя. Каждое событие может иметь несколько потомков, но не более одного родителя. Степень связанности событий в пространстве‒времени‒магнитуде для событий с индексами i и j определяется значениями функции близости [Baiesi, Paczuski, 2004]:

${\eta }_{ij}=\left\{\begin{array}{c}\begin{array}{ll}{t}_{ij}{\left(r_{ij}\right)}^{d_f}{10}^{-b{m}_i},& t_{ij}>\mathrm{0,}\\ & \end{array}\\ +\infty , t_{ij}\le \mathrm{0,}\end{array}\right.$ (1)

где $t_{ij}$ – время между событиями (сутки), $r_{ij}$ – расстояние между эпицентрами, d_f – фрактальная размерность пространственного распределения землетрясений, b – параметр закона Гутенберга–Рихтера (b-value), – магнитуда события i. События считаются связанными, если ${\eta }_{ij}\le {\eta }_0.$

Пороговая величина ${\eta }_0$ оптимально разделяет независимые и связанные события. Значение ${\eta }_0$ определяется по представленному в работе [Shebalin et al., 2020] методу, в котором распределение функции близости для независимых событий моделируется путем случайного перемешивания времен событий относительно их координат гипоцентра и магнитуд. Каталог фоновых событий формируется из событий в исходном каталоге, для которых значение функции близости относительно любого события η_ij > η₀ (то есть событие не имеет “родителя”).

Построение модели параметров сейсмического режима фоновых землетрясений

В отличие от методики построения карт ОСР-97, ОСР-2015 и ОСР-2016, мы не используем данные об активных разломах, опираясь исключительно на данные каталога землетрясений. Более того, даже результаты распознавания возможных мест сильных землетрясений по методике FCAZ [Gvishiani et al., 2022b], используемые в данной работе, также основаны исключительно на данных каталога землетрясений.

Для построения локальных оценок параметров закона Гутенберга–Рихтера использован метод среднего положения (МСП) [Vorobieva et al., 2024; Shebalin et al., 2024]. Регион сканируется кругами с постоянным радиусом R с центрами в узлах регулярной сетки с шагом D по широте и долготе. В каждом круге подсчитывается число событий $N_{i,j}$ с $M\ge M_c$ ($M_c$ – представительная магнитуда) или оценивается значение параметра b закона Гутенберга–Рихтера. Для оценки параметра b диаметр кругов должен как минимум в 5 раз превышать размер очага сильнейшего события в исследуемом диапазоне магнитуды [Крушельницкий и др., 2024]. Радиус круга для подсчета числа событий с $M\ge M_c$ обычно выбирается вдвое меньшим.

Главная особенность МСП состоит в том, что значения оцениваемого параметра приписываются не центру круга, а среднему положению землетрясений из выборки. Вторая важная особенность – учет неравномерности пространственного распределения эпицентров (фрактальности). Для подсчета числа событий в пространственной ячейке в единицу времени ${\lambda }_{i,j}$ используется следующая нормировка [Vorobieva et al., 2024]:

${\lambda }_{i,j}=\frac{1}{T}{N}_{i,j}\frac{S_{cell}}{S_{circle}}$, (2)

где $S_{circle}=R^{d_f^b}\frac{{\pi }^{\raisebox{1ex}{$d_f^b$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{\Gamma \left(1+\raisebox{1ex}{$d_f^b$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)}$ и $S_{cell}=D^{d_f^b}\text{cos}\left(\phi \right)$ – площади круга и ячейки в $d_f^b$-мерном пространстве; $\phi$ – широта центра ячейки; Г – гамма-функция; $R$– радиус круга; $d_f^b$ – фрактальная размерность пространственного распределения эпицентров фоновых землетрясений. В некоторые ячейки попадет несколько значений ${\lambda }_{i,j}$ – в этом случае выбирается максимальное. В некоторые ячейки не попадает ни одно значение – в них значения ${\lambda }_{i,j}$ определяются путем интерполяции криволинейными сплайнами с помощью встроенной процедуры “Surface” пакета Generic Mapping Tool [Wessel et al., 2019].

Параметр b оценивается, если число событий в круге составляет 50 или более. В противном случае используется региональная оценка параметра b. Оценивание проводится методом максимального правдоподобия для группированных данных с ограниченным максимумом магнитуды [Bender, 1983]. Этот метод дает несмещенную оценку для выборок небольшого объема. Полученные оценки также приписываются среднему положению эпицентров, в случае нескольких оценок выбирается сделанная по максимальному числу событий.

Оценка параметров афтершоков

Данных для проведения локальных оценок параметров афтершоков, кроме продуктивности, обычно недостаточно, поэтому для параметров c и p закона Омори–Утсу и параметра $b_{aft}$ закона Гутенберга–Рихтера выполняются региональные оценки.

Для оценки параметров закона Омори–Утсу составляется общий каталог времен афтершоков относительно их “родителей” в схеме Заляпина–Бен-Зиона. Параметры определяются методом, представленным в работе [Holschneider et al., 2012]. Оценка параметра $b_{aft}$ проводится методом, представленным в работе [Bender, 1983], для всех афтершоков.

Средняя дельта-продуктивность Λ_∆M [Shebalin et al., 2020] также оценивается локально методом среднего положения. В кругах радиуса R_Λ подсчитывается среднее число “потомков” по схеме Заляпина–Бен-Зиона для каждого события с эпицентром в этом круге, а значения приписываются среднему положению эпицентров. Потомки не обязательно должны иметь эпицентры в круге. Магнитуда потомков ограничена снизу величиной, на ΔM меньшей магнитуды “родителя”. Если в ячейку (i,j) попадает несколько значений, выбирается максимальное. Пустым ячейкам приписывается региональное значение.

Генерация синтетического каталога фоновых событий

Формирование синтетического каталога производится с помощью генератора случайных чисел с заданным распределением. Для каждого фонового события сначала синтезируется время в соответствии с распределением Пуассона с параметром, равным региональному числу событий в единицу времени с $M\ge M_0$. Затем для этого события синтезируется значение магнитуды в соответствии с региональным магнитудно-частотным распределением. Стандартно принимается модель Гутенберга–Рихтера в пределах $M_0\le M\le M_{max}$. Магнитуда определяется с точностью 0.1.

Для определения координат гипоцентра используется трехмерная таблица интенсивности событий заданной магнитуды ${\lambda }_{i,j,m}$, пересчитанная из двухмерных таблиц ${\lambda }_{i,j}$ и $b_{i,j}.$ Здесь третий индекс соответствует магнитуде. Поскольку в данной работе вводятся ограничения на возможное положение эпицентров землетрясений с М ≥ 5.5, определяемые результатами FCAZ-распознавания, то значения в ячейках трехмерной таблицы, соответствующие значениям магнитуды , обнуляются. Для всех значений магнитуды интенсивность нормируется так, чтобы сумма по пространственным элементам (i,j) для каждого значения магнитуды равнялась единице. Генерируется случайная величина с равномерным распределением на интервале (0,1). По этому значению находится ячейка (i,j). Затем более точное положение определяется в предположении равномерного распределения эпицентров внутри ячейки.

Глубина очага синтезируется аналогично работе [Малютин и др., 2024] по распределению Вейбулла с параметрами, определенными для всего региона.

Генерация полного синтетического каталога с афтершоками

Генерация афтершоков в синтетическом каталоге осуществляется по модели ETAS-e [Shebalin et al., 2020] для каждого события из каталога фоновых землетрясений. Для этого последовательно генерируются непосредственные афтершоки: сначала самого события, на второй итерации – афтершоки непосредственных афтершоков и т.д. Времена между ними определяются в соответствии со свойством Пуассоновского потока событий. Генерируется случайная величина u, имеющая равномерное распределение в интервале (0,1). Время между событием i и первым непосредственным афтершоком t_i1 определяется как $t_{i1}=-\text{ln}\left(u\right)c^p/f\left(M_i\right)$, где $f\left(M_i\right)$ – интенсивность потока афтершоков в момент события i. В соответствии с законом продуктивности [Shebalin et al., 2020, 2022] величина $f\left(M_i\right)$ сама является случайной и имеет экспоненциальное распределение. Среднее значение величины $f\left(M_i\right)$ находится из соотношения

$\begin{array}{l}{\Lambda }_{\Delta M}{10}^{b_a\left(M_m-M_c-\Delta M\right)}={\int }_{t_i}^{t_i+T}\frac{f\left(M_i\right)}{{(t-t_i+c)}^p}dt=\\ ={\int }_0^T\frac{f\left(M_i\right)}{{(t+c)}^p}dt=f\left(M_i\right)\left[u\left(T+c\right)-u\left(c\right)\right]\end{array}$ (3)

где $u\left(x\right)=\left\{\begin{array}{c}\ln \left(x\right), p=1\\ \frac{x^{1-p}}{1-p}, p\ne 1\end{array}\right.$.

Для определения времени следующего афтершока генерируется новое значение u. Время второго афтершока от момента инициирующего землетрясения определяется как $t_{i2}=t_{i1}-\ln \left(u\right){\left(c+t_{i1}\right)}^p/f\left(M_i\right)$.

Итерации продолжаются (для каждого нового “родителя” генерируется новое значение функции $f\left(M_i\right)$), пока время нового афтершока не превысит значение T. Мы используем значение T = 365 сут. Часто оказывается, что у события i нет ни одного афтершока, благодаря чему вся процедура конечна.

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

S(км²) = 10^{−3.49+0.91Mi}. (4)

Затем вычисляются величины:

L1(км) = 0.5 × 10^{−2.44+0.59Mm}, (5)

L2(км) = S/L1/π.

Обычно , но возможна и обратная ситуация, поэтому далее мы оперируем величинами L_max = max(L1, L2) и L_min = min(L1, L2). Координаты эпицентра афтершока случайно генерируются как двумерная нормально распределенная случайная величина со стандартными отклонениями, соответствующими значениям большей (L_max) и меньшей (L_min) полуосям эллипса, ориентированного вдоль ближайшего разлома из базы данных активных разломов ГИН РАН [Zelenin et al., 2022], и средними, равными координатам основного толчка. Если же направления разломов равноправны, то область имеет форму круга и координаты афтершока генерируется по двумерному нормальному распределению со стандартными отклонениями, равными радиусу круга R[км] = (S/π)^0.5, и средними, равными координатам основного толчка.

Верификация моделей

Чтобы проверить, действительно ли синтетический каталог моделирует реальный каталог землетрясений, в работе [Shebalin et al., 2024] было предложено использование L-теста [Zechar et al., 2010] и его модификаций, основанных на функции правдоподобия. Функция правдоподобия определяет совместную вероятность реальных событий, произошедших в соответствующих ячейках. Правдоподобие определяется как произведение по всем элементам пространства‒магнитуды вероятностей реализации в них фактически наблюденного числа событий. Обычно используется логарифмическая функция правдоподобия:

$L=\underset{i,j,m}{}\text{ln}\left[f_{ijm}\left(\omega \left(i,j,m\right)\right)\right]$, (6)

где $\omega \left(i,j,m\right)$ – число реализаций землетрясений в интервале магнитуды m и в элементе пространства (i,j) за период времени T; $f_{ijm}\left(k\right)$ – вероятность реализаций землетрясений за период T в ячейке пространства‒магнитуды. При использовании синтетического каталога вероятность может быть рассчитана как частота реализаций значений k при разбиении каталога на интервалы длительности T. Чтобы понять по значению L, соответствует ли синтетический каталог реальным данным, можно путем многократного применения формулы (6) к непересекающимся частям синтетического каталога землетрясений длительностью каждый построить распределение величины L. По распределению оценивается доля γ случаев, в которых значение правдоподобия по синтетическому каталогу меньше, чем по реальному. Очень маленькое значение γ указывает на то, что модель не согласуется с наблюдениями (на доверительном уровне 100 (1− γ)%). Значение γ ~ 50% означает хорошее согласование модели с реальными данными, по которым она построена.

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

РЕГИОН ИССЛЕДОВАНИЯ И ИСХОДНЫЕ ДАННЫЕ

Детальное изучение сейсмичности Восточного сектора АЗРФ началось сравнительно недавно. Из-за малой плотности населения, отсутствия в досоветский период письменности у многих населяющих его народов в регионе очень мало свидетельств проявления землетрясений в прошлом. Фактически “сейсмическая история” региона начинается лишь с XX-го столетия. Наиболее детальная информация о сильных землетрясениях в регионе собрана в работах [Новый каталог …, 1977; New Catalog …, 1982] (в англоязычной версии каталога данные дополнены до 1981 г.). В более поздний период в регионе произошло около десятка землетрясений магнитудой 6.0 и более. Два из них были неожиданностью с точки зрения карт ОСР, что привело к существенной коррекции карт именно в рассматриваемом регионе (см. рис. 1). Это Олюторское землетрясение с M = 7.6 [Chebrov, 2010; Lander et al., 2010; Rogozhin et al., 2010] и Илин-Тасское (Абыйское) землетрясение с [Shibaev et al., 2020]. Совместному анализу тектоники, геодинамики и сейсмичности в регионе в последние годы посвящен целый ряд исследований [Имаев и др., 2020; Imaeva et al., 2019; Daragan-Sushchova et al., 2015; Kanao et al., 2015].

 

Рис. 1. Фрагменты карт ОСР для Востока Российской Арктики.

Показаны фрагменты карт ОСР-97А (слева), ОСР-2015А (в центре), ОСР-2016А (справа).

 

В качестве исходных данных в работе мы использовали наиболее полный и представительный каталог землетрясений с однородной магнитудной шкалой [Gvishiani et al., 2022a], созданный путем системно-аналитического алгоритмического объединения [Vorobieva et al., 2022] данных различных каталогов и исследования соотношений разных шкал магнитуды и энергетического класса. Каталог содержит информацию о 23 254 сейсмических событиях за период с 1962 по 2020 гг. Он находится в открытом доступе на сайте Мирового центра данных по физике твердой Земли (http://www.wdcb.ru/arctic_antarctic/arctic_seism.html).

РЕЗУЛЬТАТЫ

Выделение в объединенном каталоге [Gvishiani et al., 2022a] фоновых землетрясений и иерархических цепочек афтершоков по схеме Заляпина–Бен-Зиона было проведено в работе [Воробьева и др., 2024]. Для этого был выбран период каталога, в котором данные можно считать однородными. В период до 1982 г. магнитуды большинства событий были пересчитаны из целочисленных энергетических классов и имеют шаг дискретизации 0.5, поэтому для анализа был выбран период с 1982 по 2020 г. Оценка представительной магнитуды для этого периода составляет $M_{ñ}=4.0$. Для фрактальной размерности и параметра закона Гутенберга–Рихтера получены оценки d_f = 1.81 и $b=0.924$. По методике, представленной в работе [Shebalin et al., 2020], определено пороговое значение функции близости $lg{\eta }_0=-0.61$. В каталоге по схеме Заляпина–Бен-Зиона выделены пары “родитель”‒“потомок”. В отдельный каталог записаны фоновые события, не имеющие “родителя”. Каталог включает 371 событие с $M\ge M_{ñ}=4.0$.

Модель параметров сейсмического режима фоновых землетрясений также была построена в работе [Воробьева и др., 2024]. Проведена необходимая для этого оценка фрактальной размерности и параметра b закона Гутенберга–Рихтера, получены региональные значения $d_f^b=1.88$ и $b=0.9$. Для локальных оценок ${\lambda }_{i,j}$ по методу МСП использовались круги радиусом R = 200 км, для интерполяции в кругах без эпицентров принималось значение ${\lambda }_{i,j}={10}^{-5} го{д}^{-1}$. Оценки $b_{i,j}$ проводились в кругах радиусом R_b = 300 км при минимальном количестве событий 50. При меньшем числе событий принималось региональное значение $b=0.9$.

Как было показано в работе [Воробьева и др., 2024], даже простой расчет числа землетрясений с $M\ge 6$ по таблицам ${\lambda }_{i,j}$ и $b_{i,j}$ достаточно хорошо отражает возможные места таких землетрясений: большинство эпицентров землетрясений с $M\ge 6$, в том числе за период с 1900 по 1981 г., не учитывавшийся при построении модели, располагаются в локальных зонах высоких значений ${\lambda }_{i,j}\left(M\ge 6\right)$. Однако площадь таких зон слишком велика, вследствие чего из-за необходимости соблюдения региональной повторяемости сильных землетрясений локальные оценки оказываются заниженными.

В данной работе увеличение контрастности зон возможного возникновения сильных землетрясений достигается путем использования результатов FCAZ-распознавания [Gvishiani et al., 2022b]. Зоны, в которых в соответствии с этими результатами возможны эпицентры землетрясений с $M\ge 5.5$, показаны на рис. 2. При расчете трехмерной нормированной таблицы ${\lambda }_{i,j,m}$ значения вне зон FCAZ обнуляются, благодаря чему, с учетом перенормировки, оценки повторяемости сильных землетрясений внутри зон увеличиваются.

 

Рис. 2. Зоны возможного возникновения землетрясений с M ≥ 5.5, определенные методом FCAZ-распознавания [Gvishiani et al., 2022b].

Зоны FCAZ показаны красной штриховкой. Сплошной контур показывает границы рассматриваемого региона.

 

Параметры закона Омори–Утсу оценивались на интервале (0.001, 100) сут (рис. 3). Получена устойчивая оценка с = 0.0128, р = 1.09407 (см. рис. 3a); модель с такими параметрами хорошо воспроизводит реальное распределение времен непосредственных афтершоков (см. рис. 3б).

 

Рис. 3. Оценка параметров c и p закона Омори–Утсу для всех пар “родитель”–“потомок” в каталоге, размеченном по схеме Заляпина–Бен-Зиона.

а – апостериорное распределение байесовских оценок c и p: контурами с маркерами показаны линии уровня квантилей, белым кружком отмечено положение максимума правдоподобия; б – распределение времен афтершоков: кружками обозначено эмпирическое распределение по парам “родитель”–“потомок”, сплошной линией – теоретическое распределение по закону Омори–Утсу.

 

Локальные оценки продуктивности для ΔM = 1.0 проводились методом МСП в кругах радиусом R_Δ = 300 км (рис. 4). Если число событий в круге оказывалось 5 или менее, то продуктивности присваивалось среднее значение 0.6. Высокие значения продуктивности наблюдаются на севере Камчатки и в районе Чукотки.

 

Рис. 4. Локальные оценки средней продуктивности землетрясений Λ_∆M при ∆M = 1 методом СМП.

 

Синтетический каталог фоновых событий

Для оценки сейсмической опасности и, в частности, для построения карт ОСР используются модели затухания. Интенсивность воздействия от землетрясения с определенными параметрами зависит не только от расстояния до эпицентра, но и от глубины очага. Поэтому важной характеристикой очага землетрясения является его глубина. В данной работе мы моделируем распределение глубины очага распределением Вейбулла $f\left(h\right)=\frac{B}{A}{\left(\frac{h}{A}\right)}^{B-1}{e}^{-{\left(\frac{h}{A}\right)}^B}$ аналогично работе [Малютин и др., 2024]. Это распределение лучше других рассмотренных вариантов аппроксимирует эмпирическое распределение (рис. 5). Получены оценки параметров распределения Вейбулла A = 14.7, B = 1.81.

 

Рис. 5. Гистограммы функции плотности вероятности распределения глубин очага землетрясений с рассмотренными вариантами аппроксимаций. Наилучшая аппроксимация достигается распределением Вейбулла.

 

По методике, описанной в разделе “Генерация синтетического каталога фоновых событий”, сгенерирован каталог фоновых землетрясений, который содержит 189 552 землетрясения с магнитудой $M\ge 4$. Использована модель сейсмического режима в виде таблицы параметров ${\lambda }_{i,j,m}$, пересчитанной из двухмерных таблиц ${\lambda }_{i,j}$ и $b_{i,j}$, рассчитанных ранее в работе [Воробьева и др., 2024]. Для землетрясений с $M\ge 5.$ положение эпицентров ограничено областью, выделенной по результатам FCAZ-распознавания [Gvishiani et al., 2022b].

Полный синтетический каталог был сгенерирован на основе синтетического каталога фоновых событий по методике, описанной в разделе “Генерация полного синтетического каталога с афтершоками”. Всего в каталоге представлены данные о 471 416 землетрясениях, произошедших за последние 20 тыс. лет, 281 864 из которых – афтершоки. Общая доля афтершоков, таким образом, составляет 60% от общего числа событий, что близко аналогичному значению в реальном каталоге, равному 65%. Небольшое различие может быть связано с недооценкой значений продуктивности из-за временного локального увеличения уровня представительной магнитуды в области сильных землетрясений [Shebalin, Baranov, 2017].

По методике, описанной в разделе “Верификация моделей”, была проведена верификация двух вариантов синтетического каталога: только фоновых событий и полного.

В первом случае была рассчитана логарифмическая функция правдоподобия L = −3211.2 для фоновых землетрясений из реального каталога с магнитудой $M\ge 4$. Значения функции $f_{ijm}\left(k\right)$ в формуле (6) рассчитывались из распределения Пуассона с параметром, равным значению в соответствующей ячейке трехмерной таблицы λ_ijm, нормированной на общее число событий. Чтобы оценить, насколько модель соответствует реальности, синтетический каталог фоновых событий разделен на отрезки длительностью 39 лет (это длительность фактического каталога) и по каждому отрезку, как по фактическому каталогу, определено значение L. Гистограмма полученных значений L приведена на рис. 6а. Как видно из этого рисунка, значение L-теста лежит вблизи медианы распределения значений по синтетическому каталогу. Это говорит об очень хорошем соответствии модели фактическим данным.

 

Рис. 6. Результаты верификации синтетического каталога землетрясений по фактическому каталогу согласно L-тесту.

а – синтетический и реальный каталоги фоновых землетрясений, б – синтетический и реальный полные каталоги. Показаны гистограммы L-теста для отрезков синтетического каталога, вертикальной штриховой линией показано значение L-теста по реальному каталогу; гистограммы построены по реализациям L-теста одинаковой длительности (39 лет), равной длительности реального каталога, в которых реальный каталог заменяется выборкой из синтетического каталога.

 

Для полного каталога теоретические значения функции $f_{ijm}\left(k\right)$ рассчитать очень сложно, поэтому были использованы значения частоты реализации значения k в ячейках пространства‒времени–магнитуды в 512 вариантах отрезков синтетического каталога (20000/39 = 512.82). Для реального каталога значение L-теста составило −4776.1.

В качестве дополнительного теста в данной работе было проведено сравнение положения эпицентров землетрясений с $M\ge 6.0$ в реальном и синтетическом каталогах (рис. 7). Как видно из этого рисунка, лишь небольшая доля эпицентров сильных землетрясений прошлого оказалась вне зон проявления сильных землетрясений в синтетическом каталоге. Эти землетрясения сосредоточены на севере Якутии и на севере Чукотки.

 

Рис. 7. Сравнение положений эпицентров с M ≥ 6.0 в реальном и синтетическом каталогах.

В качестве подложки использованы зоны нормативной балльности по карте ОСР-2015А (обозначения цветов см. на рис. 1). Синими точками обозначены эпицентры из синтетического каталога за 20 000 лет (на карте слева – все события, на карте справа – только фоновые), красными кружками – фактические эпицентры за 1900–1981 гг. [New Catalog …, 1982], звездочками – фактические эпицентры за 1982–2020 гг. Рамкой показаны границы рассматриваемого региона.

 

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

С использованием интегрированного каталога Восточного сектора АЗРФ [Gvishiani et al., 2022a] определены параметры сейсмического режима для региона в пределах 57.5°–77° с.ш., 110° в.д.–165° з.д. Проведена оценка параметров не только фоновых землетрясений, но и афтершоков. Это позволило сгенерировать синтетический каталог землетрясений, который воспроизводит те же значения параметров. Результаты L-теста свидетельствуют о хорошем соответствии синтетического каталога реальному, как в части только фоновых событий, так и в целом.

Синтетический каталог создавался и при подготовке карт ОСР-97, ОСР-2015, ОСР-2016. Главные отличия данного синтетического каталога состоят в следующем:

1) синтетический каталог сгенерирован по данным только каталога землетрясений и их строго алгоритмической обработки без использования модели линеаментов и доменов, основанной на экспертных решениях;

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

3) синтетический каталог содержит в качестве одного из параметров время событий;

4) каталог моделирует афтершоковую активность.

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

Было бы закономерно сравнить действующие карты ОСР с аналогичными картами, построенными по сгенерированному синтетическому каталогу. Мы пока не ставили такую задачу, так как ее решение в значительной степени зависит от выбора модели затухания, а точное воспроизведение модели, использованной в ОСР-2015, не представляется возможным. Тем не менее были сделаны предварительные оценки с использованием уравнений макросейсмического поля Н.В. Шебалина с коэффициентами, приведенными в работе [Новый каталог …, 1977]. Чтобы избежать дискуссий на основе неполной информации, не касающихся сути данной работы, мы приводим здесь лишь расчеты площадей зон каждого нормативного балла, рассчитанного для вероятности превышений 10% за 50 лет (табл. 1). Для удобства сравнения в табл. 1 приводятся кумулятивные значения площадей (сумма площади зон данного и более высоких баллов). Расчеты направлены на то, чтобы выяснить, в какой степени разные элементы модели построения синтетического каталога влияют на расчеты ожидаемой балльности. Мы сравнили три варианта синтетического каталога (см. табл. 1): 1) каталог фоновых событий на основе оценки параметров сейсмического режима (этот каталог был получен в работе [Воробьева и др., 2024]); 2) каталог фоновых событий, в котором места сильных землетрясений ограничены результатами FCAZ-распознавания [Gvishiani et al., 2022b]; 3) полный каталог, описанный выше.

 

Таблица 1. Оценки доли площади зон внутри рассматриваемой территории, в которых с вероятностью 10% за 50 лет ожидается превышение заданного балла макросейсмической интенсивности, от общей рассматриваемой площади

	ОСР2015-А1	Вариант 12	Вариант 22	Вариант 32
VI	65.68	34.50	27.20	29.06
VII	41.57	5.52	7.52	9.09
VIII	13.28	1.48	1.70	2.46
IX	0.47	0.00	0.34	0.59

Примечание. 1 – доли площади нормативного балла на карте ОСР-2015-А от общей площади рассматриваемой территории; значения просуммированы от низшего балла к высшему; 2 – вероятности превышения рассчитаны по частоте расчетного воздействия (по уравнениям макросейсмического поля Н.В. Шебалина с коэффициентами из работы [Новый каталог …, 1977]) от землетрясений соответствующего варианта синтетического каталога.

 

Как и следовало ожидать, максимальное отличие наблюдается при переходе от первого ко второму варианту. Сужение области, в которой возможны эпицентры самых сильных землетрясений, перераспределяет ожидаемую повторяемость таких событий, в результате чего в “разрешенных” местах ожидаемая повторяемость возрастает. Различие двух вариантов карт значимо именно в таких местах. Переход от второго варианта к третьему дает не столь значительные различия. Однако они тоже весьма существенны, в особенности для высоких баллов. Для понимания причин такого расхождения мы приводим карты эпицентров “синтетических” землетрясений с M ≥ 6.0: всех (рис. 7, слева) и только фоновых (рис. 7, справа). Можно видеть, что афтершоки не только увеличивают суммарную повторяемость землетрясений с M ≥ 6.0, но и несколько расширяют зону распространения их эпицентров.

Получен и другой важный результат: предварительные оценки значительно уменьшают площади территорий, отвечающих каждому нормативному баллу от VI до VIII на карте ОСР-2015, что является желательным результатом с учетом обнаруженным недавно их значительным завышением [Шебалин и др., 2022].

Синтетический каталог может и должен быть использован для построения карт ОСР. При этом сгенерированный в данной работе каталог включает фактор времени. Это позволит получать на его основе более точные оценки сейсмического риска, учитывающие пространственные корреляции сейсмических воздействий [Шебалин и др., 2024].

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

Работа выполнена в рамках гранта Российского научного фонда № 21-77-30010 “Системный анализ динамики геофизических процессов в российской Арктике и их воздействие на развитие и функционирование инфраструктуры железнодорожного транспорта”.

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

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

Об авторах

П. Н. Шебалин

Геофизический центр Российской академии наук; Институт теории прогноза землетрясений и математической геофизики Российской академии наук

Автор, ответственный за переписку.
Email: p.n.shebalin@gmail.com
Россия, ул. Молодежная, 3, Москва, 119296; ул. Профсоюзная, 84/32, Москва, 117997

А. Д. Гвишиани

Геофизический центр Российской академии наук; Институт физики Земли им. О.Ю. Шмидта Российской академии наук

Email: p.n.shebalin@gmail.com
Россия, ул. Молодежная, 3, Москва, 119296; ул. Большая Грузинская, 10, стр. 1, Москва, 123242

П. А. Малютин

Институт теории прогноза землетрясений и математической геофизики Российской академии наук

Email: p.n.shebalin@gmail.com
Россия, ул. Профсоюзная, 84/32, Москва, 117997

Е. М. Греков

Институт теории прогноза землетрясений и математической геофизики Российской академии наук

Email: p.n.shebalin@gmail.com
Россия, ул. Профсоюзная, 84/32, Москва, 117997

А. О. Антипова

Геофизический центр Российской академии наук; Институт теории прогноза землетрясений и математической геофизики Российской академии наук

Email: p.n.shebalin@gmail.com
Россия, ул. Молодежная, 3, Москва, 119296; ул. Профсоюзная, 84/32, Москва, 117997

И. А. Воробьева

Геофизический центр Российской академии наук; Институт теории прогноза землетрясений и математической геофизики Российской академии наук

Email: p.n.shebalin@gmail.com
Россия, ул. Молодежная, 3, Москва, 119296; ул. Профсоюзная, 84/32, Москва, 117997

Б. А. Дзебоев

Геофизический центр Российской академии наук; Институт физики Земли им. О.Ю. Шмидта Российской академии наук

Email: p.n.shebalin@gmail.com
Россия, ул. Молодежная, 3, Москва, 119296; ул. Большая Грузинская, 10, стр. 1, Москва, 123242

Б. В. Дзеранов

Геофизический центр Российской академии наук

Email: p.n.shebalin@gmail.com
Россия, ул. Молодежная, 3, Москва, 119296

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

Дополнительные файлы