О моделировании сейсмического режима в задачах оценки сейсмической опасности

Полный текст

ВВЕДЕНИЕ

Проблема достоверной оценки сейсмической опасности и ее применения на практике обострилась после катастрофических землетрясений 6 и 7 февраля 2023 г. на границе Турции и Сирии. В Турции, как и в России и в большинстве государств, расположенных в сейсмоопасных регионах, оценки сейсмической опасности применяются при строительстве с тем, чтобы минимизировать число жертв и материальный ущерб в случае возникновения разрушительного землетрясения. К сожалению, как это бывало в большинстве случаев катастрофических землетрясений [1], сейсмическая опасность в эпицентральных зонах землетрясений 6 и 7 февраля 2023 г. была недооценена [2]. С другой стороны, как отмечалось в прессе, большие разрушения и огромное число жертв от этих землетрясений частично объясняются несоблюдением строительных норм, учитывающих сейсмическую опасность.

Оценки сейсмической опасности условно можно разделить на детерминистские и вероятностные. В российской литературе используются сокращения ДАСО (детерминистский анализ сейсмической опасности) и ВАСО (вероятностный анализ сейсмической опасности). И в том, и в другом случае оценивается возможность воздействий от землетрясений определенной степени в заданной точке на поверхности Земли. Степень воздействия традиционно определяется в баллах шкалы интенсивности (наиболее распространенными являются шкалы MCK-64 и Mercalli) или в значениях пикового ускорения грунта (PGA, peak ground acceleration). Детерминистские оценки отвечают на вопрос “каково максимально возможное воздействие от землетрясения в данной точке?”, а вероятностные оценки отвечают на вопрос “с какой вероятностью в данной точке возможны воздействия заданной степени?”. Обычная форма представления оценок сейсмической опасности – это карты сейсмического районирования. Первые применявшиеся на практике карты сейсмического районирования опирались на детерминистский подход, поэтому традиционно такие карты строятся в виде зон максимальной ожидаемой степени воздействия. В случае вероятностного подхода строится несколько вариантов карт с расчетами ожидаемого максимального воздействия на разные интервалы времени (например, 500, 5000, 10 000 лет).

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

Несмотря на кажущуюся сложность проблемы, наибольший прогресс в мировой науке достигнут именно в решении проблемы точного моделирования воздействия от землетрясений. Наоборот, известные ошибки в оценках сейсмической опасности связаны, как правило, с неточной оценкой местоположения очагов и их магнитуды. Известно много случаев землетрясений, которые произошли там, где землетрясения такой магнитуды никто не ожидал. Наиболее яркий недавний пример – землетрясение Тохоку магнитудой 9.2 вблизи берегов Японии 11 марта 2011 г. В СССР и затем на территории СНГ неожиданностью с точки зрения карт ОСР, действовавших в тот период, стали Газлийские землетрясения 1976 г. и в Узбекистане, Спитакское 1988 г. в Армении, Зайсанское 1990 г. в Казахстане, Рачинское 1991 г. в Грузии, Сусамырское 1992 г. в Киргизии, Хаилинское 1991 г. и Нефтегорское 1995 г. в России (в Корякии и на Сахалине). После создания комплекта карт ОСР-97 для территории России произошли Олюторское землетрясение 2006 г. в Корякии, Тувинское землетрясение 2012 г., Илин-Тасское (Абыйское) землетрясение в Якутии, землетрясение 2018 г. на Урале и др. Все указанные землетрясения произошли в неожиданных для соответствующих магнитуд местах.

Другой вид ошибок в оценках сейсмической опасности – это переоценка опасности на больших территориях. Оценка суммарных фактических воздействий от землетрясений, которые произошли в период после создания карт ОСР-97, основанных на вероятностном подходе, показала, что в большинстве регионов России оценка сейсмической опасности в среднем завышена как минимум в 10 раз [3]. Основной вклад в такие ошибки – завышенная на больших территориях оценка повторяемости землетрясений заданной магнитуды. Таким образом, главный источник ошибок в картах ОСР как пропусков цели, так и ложных тревог, – это неточная модель сейсмического режима – распределения очагов землетрясений определенной магнитуды в пространстве и во времени. Традиционно, со времен первых карт ОСР, основанных на детерминистском подходе, в российской литературе модель сейсмического режима называется моделью зон ВОЗ (возможных очагов землетрясений).

В современной мировой практике для целей общего сейсмического районирования, как правило, используется модель сейсмического режима, включающая комбинацию сейсмичности вдоль линейных структур – “линеаментов” (активных тектонических разломов), сейсмичности на площадных структурах – “доменах” и отдельных потенциальных очагов землетрясений. Предполагается, что внутри структур сейсмические события имеют равновероятное положение эпицентров и одинаковую повторяемость, одинаково зависящую от магнитуды. Для перехода от модели зон ВОЗ к картам ОСР с помощью моделей распространения воздействия землетрясения от очага к точке на поверхности Земли используется либо “классический” аналитический подход, предложенный Корнеллом [4], либо формируется синтетический каталог землетрясений, и в каждой рассматриваемой точке подсчитывается повторяемость воздействий разной величины по расчетам воздействий от каждого события из синтетического каталога.

Именно такой подход был применен при создании карт общего сейсмического районирования России ОСР-97, ОСР-2015, ОСР-2016, применявшихся для строительства в сейсмоопасных регионах в качестве нормативных. Модель была создана большим коллективом авторов под руководством В. И. Уломова и получила название линеаментно-доменно-фокальной (ЛДФ) модели [5]. Эта модель, несомненно, стала существенным шагом вперед в практике районирования сейсмической опасности. Впервые стали применяться вероятностные карты, предназначенные для расчетов на сейсмостойкость строительных объектов разных категорий ответственности и сроков службы.

Вместе с тем подход построения карт ОСР на основе ЛДФ-модели лишь условно можно назвать вероятностным. Построение линеаментов и границ доменов, а также дополнительное выделение потенциальных очагов во многом носит субъективный характер, вследствие чего положение и повторяемость эпицентров в синтетическом каталоге в значительной степени предопределены этим субъективным фактором. Разделение сейсмоопасных зон на линейные и площадные структуры в ЛДФ-модели впервые было предложено основоположником ВАСО С. А. Корнеллом [4] и исторически появилось как дополнение к детерминистской модели, включавшей только линейные структуры – активные разломы. Частые разногласия исследователей по поводу границ линейных и площадных структур привели во многих странах к использованию “сглаженной модели сейсмичности”, где границы зон не нужны (например, [6]). Подход сглаженной модели сейсмичности позволяет учесть изменчивость между пространственно-временным распределением сейсмичности, согласно историческим и современным каталогам, и обеспечивает гораздо более точные в пространственном отношении оценки. Но сглаженные модели сейсмичности требуют значительно большего количества зарегистрированных землетрясений для детальной и точной оценки параметров сейсмического режима, в том числе в широком диапазоне магнитуд.

С момента создания карт ОСР-97 прошло четверть века, и за этот период накоплен огромный объем данных о землетрясениях, позволяющий перейти от зональных моделей к сглаженным моделям сейсмичности. Другой аргумент в пользу перехода к моделям сглаженной сейсмичности – радикальное улучшение вычислительных возможностей со времени разработки ВАСО. Теперь нет необходимости проводить вычисления для отдельных зон как в классическом варианте ВАСО с использованием формул Корнелла, так как расчеты выполняются посредством генерации синтетических каталогов землетрясений и прямого моделирования воздействия от каждого отдельного землетрясения.

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

В принципе в качестве модели сейсмического режима могут использоваться зональные модели типа МДФ. В этом случае параметры сейсмического режима оцениваются на основе фактического каталога для отдельных зон модели. При использовании сглаженной модели сейсмичности параметры модели обычно оцениваются по фактическому каталогу на регулярной сетке с определенным шагом. К распространенным методам сглаживания относятся те, в которых используется постоянный радиус сглаживания (например, [7]), и те, в которых радиус сглаживания увеличивается по мере уменьшения плотности наблюдаемых землетрясений (например, [6]). В обоих подходах значение оцениваемого параметра приписывается центру круга. Это всегда приводит к существенному снижению контрастности результатов оценки: на наиболее сейсмически активных участках оценка активности неизбежно оказывается заниженной, а по краям этих участков – завышенной. Дополнительное снижение контрастности вызвано тем, что при нормировании оценок активности на единицу площади не учитывается фрактальная структура распределения сейсмичности [8], и нормирование производится на всю площадь круга. В данной работе мы предлагаем высококонтрастный метод сглаживания, в котором оценки параметров приписываются не центру круга, а среднему положению эпицентров землетрясений, использованных в расчетах, и в котором площадная нормировка производится с учетом фрактальности.

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

Практически все модели сейсмического режима основаны на применении закона Гутенберга-Рихтера. Согласно этому закону, график зависимости от магнитуды числа землетрясений в логарифмическом масштабе имеет в широком диапазоне магнитуд прямолинейный вид с наклоном, тангенс которого получил название “b-value”. Это позволяет оценивать повторяемость более сильных землетрясений по повторяемости более слабых. В области малых магнитуд прямолинейность графика сохраняется до представительной магнитуды, обозначаемой Mc. В области больших магнитуд распределение ограничено максимальной возможной магнитудой Mmax; обычно используется двойное усеченное распределение. В зональных моделях величина Mmax обычно определяется на основе баз данных разломов (например, [9]). По геологическим, палеосейсмическим или геодезическим данным может оцениваться и повторяемость сильнейших землетрясений на конкретных разломах (например, [10]).

Еще одна часто используемая особенность применения закона Гутенберга-Рихтера – это так называемая характеристическая модель [11], по которой узкий диапазон магнитуды (например, очень большие магнитуды) имеет значительно большую повторяемость, чем определяемая по закону Гутенберга-Рихтера. В данной работе в отдельном разделе мы обсудим пределы применимости закона Гутенберга-Рихтера и проблему совместного использования сейсмологических и геологических данных для построения модели сейсмического режима.

Практически все используемые варианты ВАСО опираются на два важных предположения, которые связаны друг с другом. Первое – предположение о стационарности сейсмичности; иными словами, предполагается, что повторяемость землетрясений является постоянной в разные периоды времени и что период времени доступных наблюдений репрезентативен для того, какой будет сейсмическая активность в долгосрочной перспективе. Второе – возникновение землетрясений считается независимым от времени, и каждое землетрясение следует пуассоновскому процессу во времени [4]. Ключевым следствием этого является то, что 50-летний прогноз ВАСО справедлив для любого 50-летнего временного окна и не является конкретным только для следующих 50 лет. Однако хорошо известно, что из-за группирования землетрясений в пространстве и во времени распределение Пуассона допустимо лишь относительно каталогов землетрясений, прошедших процедуру декластеризации [12]. Декластеризация – это метод, при котором краткосрочные группы землетрясений (например, афтершоки и форшоки) удаляются из каталога, и оставляется самое сильное землетрясение каждой последовательности. Обычно предполагается, что если конструкции не повреждены в результате основного толчка, то они не будут затронуты афтершоками, которые, как предполагается, вызывают колебания грунта меньшей интенсивности. Однако недавние повторные сильные землетрясения в Турции, Кентербери, Новая Зеландия, центральной Италии и Кумамото, Япония, продемонстрировали, что землетрясения, которые можно было бы исключить путем декластеризации, могут вызывать жертвы и значительный ущерб. Кроме того, даже не очень сильные афтершоки могут увеличивать суммарный ущерб, так как в результате частичного разрушения конструкций при основном толчке резко снижается их прочность. Использование в ВАСО декластеризованных каталогов, таким образом, может привести к значительной недооценке риска.

Декластеризация основана на бинарной характеристике, согласно которой землетрясение либо является основным толчком, либо нет. Разделение каталога на основные и связанные землетрясения в значительной степени условно, поэтому разные методы декластеризации дают существенно разные результаты [13]. Дополнительным и часто упускаемым из виду последствием простых техник декластеризации (например, [12]) является то, что они могут вносить значительные искажения в оценки параметров закона Гутенберга-Рихтера, поскольку в результате их применения удаляется неслучайная выборка магнитуд землетрясений [13].

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

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

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

Модель сейсмического режима должна хорошо соответствовать реальной сейсмичности. Для этого необходимо с достаточной детальностью по пространству (для построения карт ОСР-97, ОСР-2015, ОСР-2016 использовалась сетка с шагом 0.25 градуса) параметризовать реальную сейсмичность. Главная проблема, которая здесь возникает, состоит в том, что необходимы локальные, точечные оценки, а подсчет числа событий возможен только в некотором пространственном объеме. Проблема решается путем использования различных способов сглаживания. Из-за значительной неравномерности пространственного распределения сейсмичности результаты сильно зависят от выбранного способа сглаживания.

Для целей оценки сейсмической опасности в различных странах приняты разные процедуры сглаживания [13], в которых используются разные варианты пространственных областей, в которых ведется подсчет числа событий. Их можно разделить на три основных типа: 1) модели типа ЛДФ [5], 2) оценки в кругах фиксированного радиуса, 3) так называемые адаптивные оценки. В моделях типа ЛДФ исследуемая территория делится на непересекающиеся области (“линеаменты” – протяженные структуры небольшой ширины вдоль активных разломов и “домены” – площадные структуры между линеаментами). Внутри каждой из структур положение эпицентров считается равновероятным. В моделях второго типа подсчет числа событий ведется в кругах фиксированного радиуса, а значения приписываются центрам кругов. В моделях третьего типа радиус круга определяется заданным числом событий, поэтому такой метод часто называется также методом к-ближайших соседей [6], а значения также приписываются центру круга.

О недостатках моделей типа ЛДФ говорилось во Введении. Адаптивные оценки дают более контрастную картину распределения сейсмичности по сравнению с другими способами сглаживания, однако излишняя детальность может быть не вполне оправданной из-за значительного пространственно-временного группирования землетрясений, которое сохраняется даже после декластеризации каталога землетрясений. Еще один недостаток такой модели состоит в том, что оценка интенсивности потока событий в единицу времени на единицу площади опирается на предположение об однородности, хотя метод нацелен именно на выявление деталей пространственной неоднородности распределения сейсмичности. Использование сглаживания с помощью кругов фиксированного радиуса, с одной стороны, позволяет рассматривать области достаточно большого размера (это необходимо для применения закона Гутенберга-Рихтера, см. следующий раздел), но, с другой стороны, может приводить к неоправданному “размазыванию” сейсмической активности.

Здесь мы предлагаем модифицированный вариант метода сглаживания кругами фиксированного радиуса, дающий значительно более контрастные результаты по сравнению со стандартным методом и не содержащий описанных выше недостатков метода k-ближайших соседей. Регион исследования сканируется кругами с постоянным радиусом R. Центры кругов расположены в узлах регулярной сетки с заданным шагом D по широте и по долготе. Хорошо известно, что сейсмичность имеет фрактальную структуру по пространству [8], которую необходимо учитывать при нормировании интенсивности потока событий. Это учитывается в нашем методе. В каждом рассматриваемом круге производится оценка $\nu \left(M\right)$ числа землетрясений с магнитудой m ≥ M землетрясений в пространственной ячейке , нормированного на длительность каталога

$\nu \left(M\right)=N\left(M\right)\frac{S_{cell}^{d_f}}{S_{circle}^{d_f}}$. (1)

Здесь N(M) – число землетрясений с магнитудой m ≥ M в круге, $S_{cell}^{d_f^{}}$ и $S_{circle}^{d_f^{}}$ – площади ячейки и круга в $d_f^{}$-мерном пространстве:

$S_{circle}^{d_f^{}}=R^p\frac{{\pi }^{\raisebox{1ex}{$d_f^{}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{\Gamma \left(1+\raisebox{1ex}{$d_f^{}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)};$ $S_{cell}^{d_f^{}}=D^{d_f^{}}{\cos }^{\raisebox{1ex}{$p{d}_f^{}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\left(\phi \right),$

где – это гамма-функция, R – радиус круга, $\phi$ – географическая широта центра круга, $d_f^{}$ – фрактальная размерность пространственного распределения эпицентров.

На рис. 1 приведена карта оценки вариаций сейсмической активности в регионах юга Сибири. Карта построена по данным каталога землетрясений за 1982–2021 гг., полученного объединением каталогов из ежегодников “Землетрясения в СССР” за 1982–1991 гг. (оцифрованная версия доступна по ссылке http://www.wdcb.ru/sep/seismology/cat_USSR.ru.html), ежегодников “Землетрясения в Северной Евразии” (1992–2021) (http://zeus.wdcb.ru/wdcb/sep/hp/seismology.ru/ cat_N_Eurasia.ru.html) и каталога Международного сейсмологического центра ISC (http://www.isc.ac.uk). Дубли в каталоге выявлены и удалены по методике работы [14]. Проведена также декластеризация каталога по методике работы [15]. Результирующий каталог доступен по ссылке https://www.itpz-ran.ru/wp-content/rnf/rnf-20-17- 00180/S_Siberia_ForUpload_MainShocks.txt. Локальные оценки представительной магнитуды для каталога (карта доступна по ссылке https://www.itpz-ran.ru/wp-content/rnf/rnf-20-17-00180/Представительная_магнитуда.pdf) позволяют принять значение M_c = 3.5 в качестве региональной оценки. Регион исследования сканируется кругами с постоянным радиусом R = 100км. Шаг сканирования D=0.1 град. Величина $a={\log }_{10}\nu \left(3.5\right)$, где $\nu \left(3.5\right)$ – оценка числа землетрясений в ячейке с магнитудой M ≥ 3.5, вычисленная по формуле (1). Значения активности привязаны к центрам кругов сканирования.

 

Рис. 1. Вариации сейсмической активности $a={\log }_{10}\nu$, v – оценка числа землетрясений с магнитудой M ≥ 3.5, вычисленная по формуле (1). Землетрясения показаны черными кружками. Значения привязаны к центрам кругов сканирования.

 

Как видно из рис. 1, зоны высокой сейсмической активности оказались значительно расширенными по сравнению с распределением эпицентров. Это связано с тем, что центры кругов сканирования могут оказаться в стороне от землетрясений, попавших в круг. Пример для узла с координатами (106^° в. д., 53^° с. ш.) показан на рис. 2.

 

Рис. 2. Карта эпицентров землетрясений в круге с координатами центра (106° в. д., 53° с. ш.). Синим крестом показан центр круга, красным крестом – среднее положение землетрясений (106.8° в. д., 52.6° с. ш.), которое сдвинуто относительно центра круга на 70 км.

 

В предлагаемой здесь модифицированной методике мы будем привязывать значение не к центру круга, а к среднему положению землетрясений выборки. В результате этой операции в некоторые ячейки ${0.1}^{°}\times {0.1}^{°}$ попадет несколько значений активности. Для каждой такой ячейки мы выберем единственное значение v, соответствующее максимальной оценке активности. В тех ячейках, куда не попало ни одного значения, мы восстановим активность с помощью интерполяции.

Интерполяция криволинейными сплайнами выполнена с помощью встроенной процедуры “Surface” пакета Generic Mapping Tool [16], являющейся модификацией известного метода минимальной кривизны. Основной проблемой методов интерполяции при картировании данных является получение необоснованных выбросов в результатах, которые плохо вписываются в реальную картину (например, резкие перепады высот на равнине в картах рельефа и т. п.). Это происходит из-за того, что значения в интервалах между точками с известными данными, ничем не ограничены. Проблема обычно решается с помощью современных методов числовых поверхностей, которые можно разделить на две группы: методы с весами и методы конечных разностей. В первом случае интерполируемое значение определяется взвешенным средним значений окружающих точек. Чем дальше находится точка, тем меньше ее вес. Однако, этот метод дает плохие результаты, когда плотность данных сильно варьируется. Области с высокой плотностью будут сглаживаться, а области с низкой плотностью интерполироваться довольно неравномерно. Для сглаживания интенсивности потока сейсмических событий применим скорее метод конечных разностей, в котором предполагается, что искомая поверхность подчиняется некоторому дифференциальному уравнению. Это уравнение затем аппроксимируется конечными разностями и решается итеративно. В самом простом случае в качестве уравнения можно использовать уравнение Лапласа:

$\frac{{\delta }^{2}z}{{\delta }^{2}x}+\frac{{\delta }^{2}z}{{\delta }^{2}y}=0.$ (2)

В таком случае значение в точке $z_{ij}$ будет определяться четырьмя значениями в соседних точках. В пакете GMT используется более сложное уравнение смещений тонкого металлического листа, изгибаемого силами, действующими в точках, так что смещение в этих точках равно наблюдаемому значению [16]. В таком случае значение в точке $z_{ij}$ будет определяться двенадцатью значениями в соседних точках. Показано, что такое уравнение обеспечивает условие минимизации квадрата кривизны поверхности. Суммарная запасенная энергия упругой деформации в согнутой пластине примерно пропорциональна кривизне. Из всех дважды дифференцируемых поверхностей, интерполирующих данные, поверхность минимальной кривизны сохраняет наименьшую энергию деформации [16]. Таким образом, накладываются ограничения на интерполируемые значения между точками с известным данными и этим ограничениям дается физическое обоснование.

Процедура “Surface” пакета GMT восстанавливает значения z в узлах регулярной сетки по заданным значениям (x, y, z) на произвольном множестве точек внутри региона с помощью решения уравнения:

$\left(1 - t\right)\Delta \left(\Delta \left(z\right)\right) + t \Delta \left(z\right) =0.$ (3)

Здесь $\Delta$ – оператор Лапласа, t – параметр натяжения (tension factor). Мы использовали рекомендованное значение t = 0.25 [16]. В нашем случае x, y – это географические координаты (долгота и широта), z – сейсмическая активность $a={\log }_{10}\nu$.

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

 

Рис. 3. Вариации сейсмической активности $a={\log }_{10}\nu$, v – оценка числа землетрясений с магнитудой M ≥ 3.5, вычисленная по формуле (1). Значения привязаны к среднему положению землетрясений выборки. Землетрясения показаны черными кружками.

 

В качестве стандарта предлагается использовать значение R = 100 км, что составляет порядка 10 размеров очага землетрясений с магнитудой М = 6.0, что позволяет получать корректные локальные оценки параметров закона повторяемости. С другой стороны, эта величина, как правило, значительно меньше размера региона исследования, что позволяет выявить локальные вариации сейсмической активности.

Аналогичный подход мы предлагаем использовать для локальных оценок параметра b-value. Оценки предлагается проводить в кругах фиксированного радиуса (мы предлагаем использовать значение R = 200 км, чтобы обеспечить достаточное количество событий для устойчивых оценок), а значения приписывать среднему положению землетрясений выборки. Значения в узлах сетки формируются с помощью метода интерполяции, описанного выше.

ПРЕДЕЛЫ ПРИМЕНИМОСТИ ЗАКОНА ГУТЕНБЕРГА-РИХТЕРА, ГИПОТЕЗА ХАРАКТЕРИСТИЧЕСКИХ ЗЕМЛЕТРЯСЕНИЙ

В соответствии с законом Гутенберга-Рихтера повторяемость более сильных землетрясений меньше повторяемости более слабых, и для оценки повторяемости более сильных событий, из-за ограниченности каталогов землетрясений по времени, ни один из методов сглаживания не может дать адекватные вероятностные оценки. Поэтому практически во всех моделях сейсмического режима используется закон Гутенберга-Рихтера, позволяющий экстраполировать оценки повторяемости слабых землетрясений на оценки для более сильных землетрясений. При этом параметр b-value закона Гутенберга-Рихтера также определяется в областях, определяемых одним из трех описанных в предыдущем разделе способов. При использовании закона Гутенберга-Рихтера и оценке его параметров необходимо соблюдение двух важных условий. Во-первых, размер рассматриваемой области должен быть много больше линейного размера очагов максимальной рассматриваемой магнитуды. Например, если мы должны оценить повторяемость землетрясений магнитуды 6, имеющих линейный размер очага порядка 10 км, то необходимо рассматривать области с линейным размером не менее 100 км; для землетрясений магнитуды 7 этот размер должен составлять не менее 500 км. Во-вторых, рассматриваемые области не должны быть слишком узкими. Закон Гутенберга-Рихтера отражает фрактальную структуру системы разломов, в которой формируются очаги землетрясения. Эта структура охватывает сложную систему разломов, ориентированных по разным направлениям. Если ограничить рассматриваемую область узкой полосой, то неизбежно будут учитываться только землетрясения, очаги которых ориентированы вдоль основной структуры. Такое искусственное ограничение может сказаться на форме закона повторяемости. Например, если рассматривается узкая область, в которой уже произошло сильное землетрясение, это распределение может иметь видимый выброс в области больших магнитуд просто потому, что для данной области рассматривается временной интервал, в котором сильное землетрясение уже точно произошло. Для корректного статистического анализа необходимо, чтобы выбор области и интервала времени исследования не зависел от данных о событиях, которые рассматриваются. Именно в нарушении этих двух важных условий может состоять реальная причина наблюдаемого явления “характеристических землетрясений”, состоящего в отклонении закона повторяемости от линейной формы в виде появления преимущественных значений магнитуды.

Гипотеза характеристических землетрясений [11] часто используется в оценках сейсмической опасности. В частности, эта гипотеза применялась для карт ОСР-97, ОСР-2015, ОСР-2016. Применимость этой модели уже давно обсуждается в литературе, но статистические тесты по-прежнему не позволяют сделать вывод для отдельных разломов.

Поскольку мы считаем нарушение двух описанных выше условий применимости закона Гутенберга-Рихтера наиболее вероятной причиной наблюдаемого эффекта, называемого гипотезой характеристических землетрясений, здесь мы предлагаем простое решение. Мы считаем, что закон Гутенберга-Рихтера выполняется в больших областях. При моделировании сейсмического режима в виде синтетического каталога землетрясений мы предлагаем генерировать времена и магнитуды землетрясений в соответствии с законом повторяемости для всей такой области. А вот координаты эпицентра могут генерироваться в соответствии с пространственным распределением локальных оценок повторяемости событий данной магнитуды, вычисляемым по результатам определения сглаженных значений $a={\log }_{10}\nu$ и b-value.

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

МЕСТА СИЛЬНЕЙШИХ ЗЕМЛЕТРЯСЕНИЙ

Слабые землетрясения в сейсмически активных районах происходят практически повсеместно. Места возможного возникновения сильных землетрясений значительно более ограничены. Обычно предполагается, что сильные землетрясения происходят на активных тектонических разломах. Эта гипотеза восходит к детерминистскому подходу к оценке сейсмической опасности и в настоящее время лежит в основе практически всех методов сейсмического районирования, применяемых на официальном уровне в мире. Такой подход имеет ряд очень существенных недостатков. Во-первых, землетрясения довольно часто происходят вне известных активных разломов. Обычно это интерпретируется тем, что разлом, на котором произошло землетрясение, ранее просто не был выявлен. Это значит, что каталоги разломов заведомо не полны. Во-вторых, картирование активных разломов часто является весьма субъективным. Этот факт хорошо демонстрирует база данных активных разломов, разработанная в Геологическом институте РАН [9]. Эта база данных представляет собой несомненную очень высокую ценность, поскольку составлена на основе 612 опубликованных источников и, таким образом, содержит наиболее полную имеющуюся информацию об активных разломах. Вместе с тем в разных источниках разные разломы картируются по-разному, а объективная возможность выбора наилучшего варианта фактически отсутствует. В-третьих, активными фактически считаются разломы, на которых проявляется или проявлялась в прошлом сейсмическая активность. Таким образом, информация об активных разломах по факту не является или, во всяком случае, не всегда является независимой от информации о землетрясениях. В результате при использовании такого подхода очень редко удается верно спрогнозировать места будущих сильных землетрясений, в которых землетрясения такой силы еще не наблюдались. Это привело к значительной недооценке ожидаемого разрушительного воздействия, рассчитанного в рамках программы GSHAP (Global Seismic Hazard Assessment Program) для большинства сильных землетрясений мира, которые произошли после построения карт в рамках этой программы [1, 3]. Пришедшая на смену GSHAP программа GEM (Global Earthquake Model), построенная на схожих принципах, не устранила указанные недостатки. Яркими примерами значительной недооценки эффекта от землетрясений являются, в частности, парное землетрясение 6 и 7 февраля 2023 г. в Турции и Сирии c (M_w 7.8 и 7.7) и землетрясение 8 сентября 2023 г. в Марокко c M_w 6.9.

Альтернативный подход к определению мест возможного возникновения сильных землетрясений давно разрабатывается с применением элементов искусственного интеллекта на основе морфоструктурного районирования (МСР) территорий [17]. МСР основано на совместном анализе современного рельефа, отображенного на топографических картах и спутниковых снимках, и геологической структуры региона, представленной на геологических картах. МСР направлено на определение иерархической линеаментно-блоковой структуры изучаемого региона и выявление наиболее тектонически активных структур – морфоструктурных узлов. Предполагается, что сильные землетрясения приурочены к морфоструктурным узлам. Основная задача этого подхода – определение мест будущих сильных землетрясений, в которых они еще не были зарегистрированы. Эта задача решается выявлением опасных морфоструктурных узлов методами распознавания образов по набору геолого-геофизических признаков.

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

МАКСИМАЛЬНАЯ ВОЗМОЖНАЯ МАГНИТУДА

В стандартных процедурах ВАСО с использованием моделей типа ЛДФ для каждого элемента модели используется вариант закона Гутенберга-Рихтера с ограничением максимальной возможной магнитуды. Обычно используется простейший вариант такого ограничения: до значения M_max принимается линейная форма закона Гутенберга-Рихтера (логарифм числа землетрясений пропорционален магнитуде с отрицательным коэффициентом b-value), а для событий большей магнитуды вероятность событий принимается равной 0. Для линейных структур значение M_max определяется исходя из линейного размера структуры. Для площадных структур обычно принимаются экспертные оценки.

Оценки M_max по размерам структуры основаны на предположении, что разломы разрываются повторяющимся и сегментированным образом; другими словами, длина разрыва ограничена границами сегментов разлома, и маловероятно, что землетрясение приведет к разрыву нескольких сегментов. Это предположение в основном было основано на геологических наблюдениях [10] и сделано по практическим соображениям упрощения расчетов. Геологические исследования некоторых исторических землетрясений показали необходимость учета многосегментных разрывов. Землетрясение 2016 года в Каикоуре, Новая Зеландия, M_w 7.8, является убедительным примером многосегментного разрыва с участием более 20 отдельных разломов [18]. Таким образом, оценки максимальной возможной магнитуды на основе информации об активных разломах оказываются весьма неоднозначными. Причина многих ошибок в оценке сейсмической опасности в рамках программы GSHAP лежит именно в недооценке максимальной возможной магнитуды по размерам структур.

Другой подход к оценке максимальной возможной магнитуды основан на статистике повторяемости землетрясений. Такой подход применим, аналогично применимости закона Гутенберга-Рихтера, для больших территорий, имеющих линейный размер много больший по сравнению с линейными размерами очага землетрясения магнитудой M_max. Обычно оценка M_max составляет наблюденное значение максимальной магнитуды с некоторой добавкой, величина которой зависит от параметров закона Гутенберга-Рихтера [19].

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

УЧЕТ ГРУППИРОВАНИЯ ЗЕМЛЕТРЯСЕНИЙ В ПРОСТРАНСТВЕ И ВРЕМЕНИ

Значительная часть землетрясений происходит как следствие предшествующих землетрясений, что проявляется в виде возникновения групп близких по времени и расстоянию событий. В каждой группе событий можно выделить одно, например, сильнейшее или первое по времени. Такие события обычно называются фоновыми. Как правило, принимается, что повторяемость фоновых землетрясений в окрестности каждой точки в рассматриваемом регионе не меняется во времени [12]. Каждое землетрясение может инициировать связанные с ним сейсмические события – афтершоки. Афтершоки составляют обычно более половины всех сейсмических событий, после сильных землетрясений, и они могут продолжаться в течение многих лет.

Как было отмечено во введении, распределение Пуассона допустимо лишь относительно каталогов землетрясений, прошедших процедуру декластеризации. Это важно, поскольку афтершоки сконцентрированы в пространстве и времени вблизи отдельных землетрясений за относительно короткий период, охватываемый фактически каталогом землетрясений (порядка 100 лет). Поэтому использование не декластеризованного каталога приведет к локальному завышению опасности вблизи эпицентров сильных землетрясений. Однако использование декластеризованных каталогов занижает общую сейсмическую опасность. Поэтому, с одной стороны, важно максимально точно разделить сейсмические события на фоновые и зависимые события. С другой стороны, при моделировании необходимо учесть и вклад афтершоков в общую сейсмичность. Обычно для стохастического моделирования афтершоков используется модель ETAS [20], в которой сейсмический режим моделируется суперпозицией последовательностей, убывающих во времени по известному закону Омори-Утсу. В пространственно-временном варианте модели ETAS учитывается степенной закон убывания интенсивности афтершоков с расстоянием от эпицентра основного толчка. В модели ETAS постулируется, что число афтершоков у землетрясений одинаковой магнитуды одинаково. Однако установленный недавно закон продуктивности землетрясений [15] опровергает этот постулат, одновременно давая объяснение факту завышения прогнозируемой вперед по модели ETAS интенсивности потока землетрясений. Закон продуктивности устанавливает, что количество афтершоков является случайной величиной, имеющей экспоненциальный вид распределения. Это распределение определяется единственным параметром – средней продуктивностью. В рамках нашего подхода мы предлагаем использовать модель ETAS-e [15], в которой дополнительно учтен закон продуктивности землетрясений.

Для процедуры декластеризации мы предлагаем использовать метод ближайшего соседа [21]. В методе постулируется однозначное определение событий-триггеров и зависимых событий по значению функции близости. Для каждого события триггером считается более раннее событие, значение функции близости с которым минимально и при этом меньше некоторого порогового значения. Это значение определяется эмпирически по распределению значений функции близости для всех пар триггер – зависимое событие (“ближайших соседей”) [15]. Продуктивность определяется как число зависимых событий в интервале магнитуды фиксированной ширины относительно магнитуды триггера.

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

СИНТЕТИЧЕСКИЙ КАТАЛОГ ЗЕМЛЕТРЯСЕНИЙ

В качестве примера применения предложенных принципов сгенерирован синтетический каталог землетрясений региона Алтай-Саяны-Прибайкалье. На первом этапе сгенерирован каталог фоновых событий, а затем для каждого фонового события по модели ETAS-e сгенерированы серии афтершоков, включающие (при наличии) серии афтершоков от афтершоков, афтершоков этих вторичных афтершоков и т. д. В большинстве случаев, в соответствии с законом продуктивности [15], число афтершоков равно нулю. Благодаря этому такая процедура быстро сходится и включает, как правило, не более 3 итераций, но возможны случаи (для сильных землетрясений), когда число итераций достигает 10. Сгенерированный на условный период 40 000 лет синтетический каталог доступен по ссылке https://www.itpz-ran.ru/wp-content/rnf/rnf-20-17-00180/baikalSynthCat.txt.

ВЕРИФИКАЦИЯ МОДЕЛЕЙ СЕЙСМИЧЕСКРОГО РЕЖИМА

Для того чтобы можно было проверить, действительно ли модель сейсмического режима выполняет прогнозную функцию, необходимо установить критерии соответствия модели реальности. Возникает также проблема сравнения с этой точки зрения разных моделей. Здесь мы предлагаем использование одного из возможных подходов, основанного на функции правдоподобия, L-тест [22]. Модель сейсмического режима определяет вероятность сейсмических событий в определенных ячейках пространства-времени-магнитуды. Функция правдоподобия определяет совместную вероятность рассматриваемых реальных событий в соответствующих ячейках. Правдоподобие определяется как произведение по всем элементам пространства-магнитуды вероятностей реализации в них фактически наблюденного числа событий. Обычно используется логарифмическая функция правдоподобия:

$\begin{array}{c}L=\underset{}{\sum _{i,j}^{}}\ln [f_{ij}\left(\omega \left(i,j\right)\right],\end{array}$ (4)

где $\omega \left(i,j\right)$ обозначает число реализаций землетрясений в интервале магнитуды i и в элементе пространства j за период времени T, $f_{ij}\left(k\right)$ указывает вероятность k реализаций землетрясений за период T.

В некоторых случаях можно предположить, что количество землетрясений в каждой ячейке имеет распределение Пуассона. В этом случае уравнение (4) имеет вид:

$\begin{array}{c}L=\sum \left[-\lambda \left(i,j\right)+\omega \left(i,j\right)\ln \lambda \left(i,j\right)-\ln (\omega \left(i,j\right)!)\right].\end{array}$ (5)

Если каталог землетрясений был использован для построения модели, то результаты L-теста могут говорить только о соответствии модели тем данным, по которым она построена: необходимое, но недостаточное условие для проверки модели. Но L-тест можно использовать также при тестировании метода по вновь происходящим землетрясениям или по части фактического каталога, не использованной для определения параметров модели. В этом случае результаты теста покажут, насколько точно модель может предсказывать будущую сейсмичность территорий.

При использовании синтетического каталога оценить, достаточно ли велико значение L, чтобы сделать вывод о соответствии модели реальным данным, можно путем многократного применения формулы (4) для непересекающихся частей синтетического каталога землетрясений длительностью T каждый. Для каждой k-й части рассчитывается логарифмическая функция правдоподобия L_k по формуле (4) аналогично расчетам с использованием реального каталога. Затем рассчитывается доля γ случаев, в которых значение правдоподобия по синтетическому каталогу меньше, чем по реальному. Очень маленькое значение γ указывает на то, что наблюдение не согласуется с моделью (на доверительном уровне 100(1-γ)%). На рис. 4 представлены эмпирическое распределение и гистограмма реализаций L_k в сравнении со значением, полученным по фактическому каталогу. Значение L находится вблизи медианы эмпирического распределения, что означает очень хорошее совпадение модели с фактическими данными.

 

Рис. 4. Результаты верификации синтетического каталога землетрясений региона Алтай-Саяны-Прибайкалье по фактическому каталогу с 1982 по 2021 гг., M – L-тест. Эмпирическая функция распределения (синяя кривая слева) и гистограмма (справа) – величины логарифмического правдоподобия. Вертикальной линией показано значение L.

 

Для сравнения разных моделей достаточно сопоставить значения L каждой из них. Модель с большим значением L лучше соответствует данным. Возможны вариации L-теста. Например, вместо набора интервалов магнитуды можно рассматривать один интервал магнитуд, соответствующий сильным землетрясениям. Тогда L-тест редуцируется к суммированию в формуле (4) только по пространственным ячейкам – S-тест. S-тест оценивает качество моделирования пространственного распределения сильных землетрясений. В М-тесте все рассматриваемое пространство редуцируется к одному элементу; этот тест оценивает качество моделирования закона повторяемости.

ОБСУЖДЕНИЕ

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

Практически все используемые в настоящее время в мире модели ВАСО в значительной степени субъективны и, кроме того, все модели всегда в той или иной степени приводят к ошибкам [13]. Эти два фактора, в случае существенной недооценки опасности для какого-то землетрясения, неизбежно приводят к снижению доверия общества к системе оценки сейсмической опасности, хотя в большинстве случаев дело не в ошибочности решений, а в недостатке данных, по которым они были построены. Единственный путь избежать такого снижения доверия – это сделать методику построения моделей ВАСО максимально объективной. В идеале желательна полная воспроизводимость результатов. В любом случае, по нашему мнению, это должна быть открытая, как минимум для научного и инженерного сообщества, система, включающая базы всех исходных данных (включая конкретные субъективные решения) и алгоритмы расчетов.

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

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

Вероятностные модели оценки сейсмической опасности изначально разрабатывались, главным образом, для внедрения системы страхования в решение проблемы снижения и компенсации общего ущерба от землетрясений. В России созданные в рамках ВАСО карты ОСР пока не привели к повсеместному внедрению страхования ущерба от землетрясений, хотя соответствующий закон уже действует с 2019 г. При совершенствовании концепции по снижению опасности землетрясений и пересмотре соответствующей нормативной базы целесообразно внедрять риск-ориентированные подходы. Система ВАСО и карты ОСР могут быть нацелены на оценки сейсмического риска, то есть возможного ущерба от землетрясений. Это позволит, например, более эффективно распределять федеральный бюджет между субъектами РФ на проведение антисейсмических мероприятий. Также это позволит выявить регионы, где землетрясения могут вызывать значительный материальный ущерб и жертвы среди населения. В таких регионах, а также в местах строительства важных объектов, целесообразно проводить более детальные исследования. В таких местах необходимо сохранить обязательность сейсмоизоляции при строительстве и соответствующую нормативную базу. При этом нормативы должны опираться на результаты детального сейсмического районирования (ДСР), которое должно не просто обеспечивать более детальный по сравнению с ОСР масштаб, но и использовать максимум доступной информации, включая карты активных разломов, палео- и архео-сейсмологические данные. Особую роль в последнее время приобретают высокоточные геодезические наблюдения с помощью спутниковых технологий. По данным Глобальной навигационной спутниковой системы (GNSS) можно оценивать скорости деформационных процессов в литосфере и выявлять наиболее опасные места, например, по балансу долговременных оценок скорости деформации по геологическим данным и краткосрочных по данным GNSS [23].

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

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

ИСТОЧНИК ФИНАНСИРОВАНИЯ

Работа выполнена в рамках гранта Российского научного фонда, проект № 20-17-00180-П.

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

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

Об авторах

П. Н. Шебалин

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

Автор, ответственный за переписку.
Email: shebalin@mitp.ru

член-корреспондент РАН

Россия, Москва; Москва

С. В. Баранов

Институт теории прогноза землетрясений и математической геофизики Российской Академии наук; Кольский филиал ФИЦ Единая геофизическая служба Российской Академии наук

Email: shebalin@mitp.ru
Россия, Москва; Москва

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

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

Email: shebalin@mitp.ru
Россия, Москва; Москва

Е. М. Греков

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

Email: shebalin@mitp.ru
Россия, Москва

К. В. Крушельницкий

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

Email: shebalin@mitp.ru
Россия, Москва

А. А. Скоркина

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

Email: shebalin@mitp.ru
Россия, Москва

О. В. Селюцкая

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

Email: shebalin@mitp.ru
Россия, Москва

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

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