Полный текст
ВВЕДЕНИЕ
Эта статья является продолжением серии по проблеме прогнозирования опасности повторных толчков после сильных землетрясений. В данной работе ставится задача оценить распределение вероятности времени, в течение которого после сильного землетрясения могут продолжаться опасные афтершоки. Под опасными понимаются афтершоки с магнитудой Mm 2 или выше (Mm магнитуда основного толчка). Однако в конце статьи приводится формула для пересчета оценок длительности периода афтершоков с магнитудой выше произвольного порога.
Эта задача на первый взгляд кажется эквивалентной задаче оценивания длительности периода афтершоков. В такой задаче порог по магнитуде обычно устанавливается равным представительной магнитуде каталога для данного региона. Но распределение длительности афтершоков, как будет показано в данной работе, сильно зависит от рассматриваемой минимальной магнитуды, поэтому просто оценка длительности афтершоков не может использоваться для оценки периода опасных афтершоков. Кроме того, хорошо известно, что число афтершоков в единицу времени убывает со временем по степенному закону (закон ОмориУтсу [Utsu, 1961]):
|
| (1) |
Очевидно, что при последовательность афтершоков формально является бесконечной, что невозможно физически, поэтому такая зависимость возможна лишь на конечном отрезке времени. Возможны и весьма растянутые во времени реальные серии афтершоков длительностью десятки и более лет, например после серии Нью-Мадридских землетрясений 18111812 гг., M = 7.07.4 [Stein, Liu, 2009]. Столь длительные периоды вряд ли представляют существенный практический интерес. Соответственно мы вводим ограничение по времени T и рассматриваем задачу в этом условном варианте, оценивая вероятность того, что на отрезке (t, T) не произойдет афтершоков магнитуды выше заданной.
Афтершоки часто сами сопровождаются сериями вторичных афтершоков. Включение этих афтершоков в общую последовательность приводит к отклонению убывания общего числа афтершоков в единицу времени от монотонного спадания, описываемого законом ОмориУтсу (1). Для описания такого ветвящегося процесса известный японский математик И. Огата, ученик Т. Утсу, работающий в области статистической сейсмологии, предложил модель, получившую название ETAS (Epidemic-type aftershock sequence Последовательность афтершоков эпидемического типа) [Ogata, 1989; 1999]. Эта модель получила широчайшее распространение, особенно в западной научной литературе, в том числе благодаря популяризации в большой серии публикаций с участием Д. Сорнета и А. Хельмстеттер (см., например [Sornette, Helmstetter, 2002]. ETAS-модель, а также ее вариант, учитывающий пространственное распределение [Zhuang et al., 2004], широко используются для оценки сейсмической опасности афтершоков, в том числе для оценивания продолжительности серии афтершоков [Hainzl et al., 2016]. В оригинальном варианте модель ETAS представляет собой суперпозицию независимых последовательностей, подчиняющихся закону ОмориУтсу, при этом новое событие «порождает» свою последовательность с интенсивностью, зависящей от магнитуды:
|
| (2) |
где: и время и магнитуда i-ого события; μ частота событий так называемой фоновой сейсмичности; K1, α константы; c, p параметры закона ОмориУтсу.
Для прогноза активности афтершоков после сильных землетрясений по модели ETAS обычно используется «сценарный» подход. Параметры μ, K1, α, c и p оцениваются по данным за какой-либо период времени (это может быть начальная часть афтершоковой серии одного сильного землетрясения) из некоторой области пространства. Затем производится многократная симуляция модели с полученными значениями параметров, и распределение искомой величины (например, времени последнего афтершока заданной магнитуды) оценивается по наблюденной в реализациях модели частоте. При этом, однако, возникает ряд проблем, требующих введения в модель дополнительных условий и ограничений, часто мало обоснованных. К наиболее существенным для рассматриваемой в данной работе задачи следует отнести следующие проблемы.
Во-первых, при оценке параметров не учитывается локальная неполнота каталога сразу после землетрясения. Чем больше разность магнитуды инициирующего толчка и минимальной магнитуды рассматриваемых афтершоков, тем дольше интервал времени, в течение которого каталог неполон [Helmstetter et al., 2006; Shebalin, Baranov, 2017; Баранов и др., 2019].
Для оценки параметров модели (2) в принципе можно исключать начальные интервалы времени после каждого события, однако оценка параметра c при этом, очевидно, окажется весьма неточной. При численной симуляции модели также придется исключать периоды времени после событий. Именно так оценивается длительность серии афтершоков, например, в работе [Hainzl et al., 2016]. С учетом того, что именно для малых времен после событий, как следует из модели (2), интенсивность потока генерируемых событий велика, а каждое событие в исключаемом интервале также должно порождать свои афтершоки, такой подход представляется нам лишь техническим решением, позволяющим избежать при численной симуляции появление расходящихся последовательностей.
Во-вторых, как следует из модели (2), количество непосредственных афтершоков (то есть не являющихся афтершоками афтершоков) любого события зависит только от разности магнитуды этого события и рассматриваемой минимальной магнитуды, так как параметры K1 и α считаются константами. Вместе с тем, хорошо известно, что количество афтершоков может варьировать в пределах нескольких порядков [Marsan, Helmstetter, 2017] и, более того, эта величина, при фиксированной разности магнитуды основного толчка и минимальной рассматриваемой магнитуды имеет экспоненциальный вид распределения с максимумом в нуле [Шебалин и др., 2018]. В силу значительной нелинейности временных зависимостей в (2) использование этой модели неизбежно приведет к завышенным оценкам длительности серии афтершоков.
В-третьих, в модели ETAS предполагается, что интенсивность порождаемого потока афтершоков после главных толчков и после афтершоков зависит лишь от магнитуды «порождающего» события. Это удобная для стохастической модели схема, представляющая афтершоковый процесс как прямой каскад разрушений [Смирнов и др., 2010]. Вместе с тем, как главный толчок, так и афтершоки происходят в результате высвобождения накопленной в течение длительного времени упругой энергии и, таким образом, имеют общую причину. Поэтому представление о том, что каждое сейсмическое событие «порождает» свои афтершоки вряд ли правомерно. Каждое событие лишь локально ускоряет процесс релаксации напряжений и поэтому с течением времени после основного сброса напряжений количество событий в таких всплесках при равных магнитудах уменьшается [Шебалин, 2018]. Применимость представлений об афтершоковых процессах тектонических землетрясений как о прямом каскаде разрушения ставится под сомнение также недавним наблюдением независимости магнитуд таких афтершоков от времени [Баранов, Шебалин, 2019].
Таким образом, модель ETAS, несмотря на общепринятый в западной литературе статус парадигмы в решении задач прогнозирования активности афтершоков, для поставленной в данной работе задачи не применима, во всяком случае, в существующих в настоящее время модификациях.
Другой парадигмой в западных публикациях по исследованию афтершоков является модель, основанная на нелинейном трении (rate-and-state friction) [Dietrich, 1994]. Интенсивность потока афтершоков в этой модели определяется скачком напряжений
|
| (3) |
где: μ по-прежнему интенсивность потока фоновых событий; сопротивление трения; период затухания афтершоков, который обратно пропорционален скорости тектонической деформации Выражение (3) близко совпадает с законом ОмориУтсу (1) при p = 1, и и экспоненциальным затуханием при [Cocco et al., 2010]. Из модели следует, что длительность афтершоков тем больше, чем больше величина с. Недавно было установлено, что значение c зависит от глубины очага [Shebalin, Narteau, 2017]. Таким образом, в рассматриваемой в данной работе задаче необходимо учитывать глубину очага. Еще один параметр, определяющий масштаб времени афтершокового процесса, в соответствии с моделью (3) скорость тектонических деформаций. В работе [Toda, Stein, 2018] утверждается, что в районах медленных тектонических движений длительность процесса афтершоков может достигать сотен и даже тысяч лет, поэтому в таких районах афтершоковая сейсмичность преобладает над фоновой. В этой связи, при оценке длительности афтершоков важным элементом является алгоритм выделения афтершоков.
В настоящее время существует значительное число алгоритмов выделения кластеров сейсмических событий (фор- и афтершоков). До недавнего времени, в основном, использовались так называемые «оконные» методы, среди которых наиболее известным был метод работы [Gardner, Knopoff, 1974]. В США, начиная с середины 80-х годов, популярным стал алгоритм П. Ризенберга [Reasenberg, 1985], связывающий землетрясения с кластерами в соответствии с зонами пространственного и временного взаимодействия. В этом алгоритме кластеры обычно растут в размерах при обработке все большего количества землетрясений. К существенным недостаткам метода следует отнести частую идентификацию ложных афтершоков. Метод Ризенберга часто применяется без оценивания результатов и с повсеместным использованием значений параметров, оптимизированных для сейсмичности Центральной̆ Калифорнии.
Важным шагом в развитии методов идентификации афтершоков стал метод, описанный в работе [Molchan, Dmitrieva, 1991; 1992]. В этом алгоритме отличие афтершоковой последовательности от потока фоновых землетрясений определяется различием функций распределения во времени и пространстве (двумерная Гауссова функция в пространстве и степенная во времени закон Омори для афтершоков и Пуассоновский поток с равномерной плотностью для фоновой сейсмичности). При идентификации афтершоков возможны ошибки двух типов: отнесение афтершока к группе фоновых событий и, наоборот, идентификация фонового события как афтершока. Алгоритм строится таким образом, чтобы уравнять эти ошибки, что обеспечивает равенство математического ожидания количества идентифицированных афтершоков их истинному значению. Недостаток этого метода в редких случаях проявляется в виде неоднозначности решения. Это происходит при специфической конфигурации афтершоковой области, когда она состоит из двух разделенных между собою областей c сильными афтершоками, приуроченными к противоположным концам очага. Другой недостаток при локально небольшом числе событий участвующих в анализе, очевидные афтершоки могут идентифицироваться как фоновые события. Алгоритм стал популярен у российских исследователей благодаря свободно распространяемому компьютерному коду, подготовленному В. Б. Смирновым [2009].
Широкое применение модели ETAS привело к появлению стохастических методов декластеризации [Zhuang et al., 2002; Marsan, Lengline, 2010]. Ядром метода декластеризации, основанном на пространственно-временном варианте модели ETAS [Zhuang et al., 2002], является расчетная интенсивность фона. Эта интенсивность является функцией пространства, но не времени, и параметров, связанных с кластерными структурами. Алгоритм строится на оценке параметров модели и генерировании с помощью датчика случайных чисел цепочек событий, в каждой из которых каждый «предок» может иметь несколько «потомков», но каждый «потомок» имеет лишь одного «предка». Для выбора конкретного «предка» из всех вероятных, используется датчик случайных чисел и, таким образом, нет единственного решения. Оценка вероятности того, что конкретное событие является фоновым или инициированным проводится по соответствующей частоте проявления в большом количестве реализаций. Метод работы [Marsan, Lengline, 2010] предусматривает использование более широкого класса моделей. В этом алгоритме декомпозиция каталога на кластеры также осуществляется с помощью датчика случайных чисел, поэтому для практических целей необходимо использовать множество реализаций. В обоих алгоритмах оценки параметров модели весьма неустойчивы из-за большого числа параметров (выше уже говорилось о проблемах оценки параметров более простой модификации модели ETAS). Поэтому, несмотря на кажущуюся возможность более точного выделения непосредственных афтершоков какого-либо землетрясения, в рассматриваемой в данной работе задаче эти алгоритмы вряд ли применимы.
В кластерном подходе для идентификации афтершоков многими исследователями предпринимались попытки ввести определение обобщенного расстояния в пространствевремени. Удачной оказалась метрика η, предложенная в работе [Baiesi, Paczuski. 2004], учитывающая фрактальные свойства пространственного распределения сейсмичности и повторяемость событий разной магнитуды в соответствии с законом ГутенбергаРихтера [Gutenberg, Richter, 1954]:
|
| (4) |
где: время событий с индексами i и j, годы; расстояние между эпицентрами (гипоцентрами), км; фрактальная размерность пространственного распределения эпицентров (гипоцентров) в рассматриваемой области; магнитуда события i; b параметр закона ГутенбергаРихтера. Мы будем называть эту величину функцией близости.
Авторами, однако, не было предложено формальное правило для нахождения порогового значения. В работе [Zaliapin et al., 2008] были введены определения перемасштабированных времени и расстояния, произведение которых равно величине η. На графиках в этих координатах для пар событий в каталоге землетрясений Калифорнии и синтетическом каталоге на основе модели ETAS четко разделяются популяции независимых и взаимозависимых событий. В дальнейшем в работах [Zaliapin, Ben-Zion, 2013; 2016] сформулировали простой и четкий алгоритм построения кластеров связанных событий. В этом алгоритме каждое событие может иметь несколько «потомков» («offspring» в терминологии авторов), но каждый «потомок» может иметь только одного «предка» («parent» в терминологии авторов). Отметим, что такая иерархическая цепочка позволяет определить «ранг» афтершока, что может быть полезно, например, для выделения только непосредственных афтершоков, т. е. афтершоков первого ранга, для которых «предки» не являются «потомками» других событий. Главный недостаток этого метода состоит в том, что он не учитывает протяженность очага землетрясения. Кроме того, серия связанных между собой событий, выделяемых по этому методу, может оказаться значительно растянутой в пространстве и тогда эту серию уже неправомерно характеризовать как серию афтершоков какого-либо одного из сильнейших событий в этой серии.
Время возникновения последнего афтершока заданной магнитуды может существенно зависеть от определения афтершоков, то есть от алгоритма их идентификации. Мы рассмотрим несколько алгоритмов, включая алгоритм МолчанаДмитриевой и модификации метода ЗаляпинаБен-Зиона.
Задача оценки опасного периода афтершоков может рассматриваться и в несколько другой формулировке. Многими авторами рассматривалась задача оценки времени сильнейшего афтершока [Saichev, Sornette, 2005; Tahir et al., 2012; Shcherbakov et al., 2018]. Эта задача может решаться одновременно с задачей оценки магнитуды сильнейшего афтершока [Shcherbakov et al., 2018]. С практической точки зрения, на наш взгляд, более актуальна задача в нашей формулировке, поскольку не только наиболее сильный афтершок после землетрясения может представлять серьезную опасность.
Теоретическое решение поставленной задачи строится на модели афтершокового процесса, в соответствии с законом ОмориУтсу и с распределением Пуассона для числа афтершоков. Магнитудные пороги для прогноза и для оценки параметров могут различаться, в этом случае для экстраполяции мы используем закон ГутенбергаРихтера. Сравнение теоретических оценок с фактическими усредненными распределениями в разных регионах и для разных интервалов глубины очага показало хорошее совпадение теоретических и фактических распределений.
С учетом возможного практического использования мы предлагаем проводить оценку длительности опасного периода афтершоков в два этапа. На первом этапе, непосредственно после сильного землетрясения, предлагается использовать усредненную модель, зависящую только от разности магнитуды основного толчка и рассматриваемой минимальной магнитуды и от глубины очага. На втором этапе, через несколько часов после землетрясения (оптимальным представляется интервал в 12 часов), предлагается использовать оценки параметров закона ОмориУтсу по уже произошедшим афтершокам, если их число достаточно для оценки параметров.
Целью данной работы является, главным образом, разработка методики для получения чисто практических результатов, которые, в том числе, могут быть использованы в разрабатываемой автоматизированной системе прогнозирования опасности афтершоков AFCAST (www.afcast.org). Вместе с тем, впервые построенная модель распределения длительности опасного периода афтершоков, дающая хорошее совпадение с наблюдениями, служит и уточнению физических представлений о сейсмогенезе.
ИСХОДНЫЕ ДАННЫЕ
Для сопоставления теоретических и фактически наблюденных выборочных распределений вероятности длительности опасного периода афтершоков мы используем следующие данные: глобальный каталог землетрясений ANSS ComCat геологической службы США (USGS) [ANSS…]; каталог землетрясений Камчатки и Командорских островов Камчатского филиала ФИЦ ЕГС РАН [Каталог землетрясений Камчатки…], каталог землетрясений Кавказа ФИЦ ЕГС РАН, в который были добавлены данные сейсмологического бюллетеня Кавказа за 19711986 [Сейсмологический…] и исключены взрывы [Каталог землетрясений Кавказа…]; каталог землетрясений Байкала и Забайкалья Иркутского филиала ФИЦ ЕГС РАН [Каталог землетрясений Байкальского…]. Для построения эмпирических закономерностей и их сопоставления с теоретическими расчетами в работе рассматриваются серии афтершоков от землетрясений мира с магнитудой 6.5 и выше за период с 1980 по 2017 гг.; Камчатки и Курильских островов с 1962 по октябрь 2018 гг. с магнитудой 6 и выше; Кавказ и Крым с 1960 по май 2016 гг. с магнитудой 5 и выше; Байкал и Забайкалье с 1960 по 2014 гг. с магнитудой 5.5 и выше.
МЕТОД ОЦЕНКИ
Мы рассматриваем афтершоки, имеющие магнитуду выше некоторого заданного порога. Порог может задаваться в зависимости от того, насколько опасны в данном регионе события такой магнитуды. Мы предполагаем, что изменение интенсивности потока афтершоков подчиняется закону ОмориУтсу [Utsu, 1961] (1). Плотность и функция распределения времени t произвольного афтершока для условного случая, когда рассматриваются времена в пределах некоторого интервала T (везде далее мы используем значение T = 365 сут.) после основного толчка, имеют вид [Holshneider et al., 2012]:
|
| (5) |
где
|
| (6) |
Мы предполагаем также, что для каждой серии число афтершоков в интервале (0, T) имеет распределение Пуассона со средним Λ:
|
| (7) |
Вероятность того, что все k афтершоков произойдут ранее момента τ, при условии τ < T, равна Учитывая распределение (7) для величины k, по формуле полной вероятности получаем распределение вероятности времени τ последнего афтершока заданной магнитуды (при условии τ < T):
|
| (8) |
Дифференцирование (8) дает выражение для плотности вероятности:
|
| (9) |
Отметим, что выражение (8) не обращается в 0 при τ = 0. Это связано с тем, что всегда существует вероятность того, что у землетрясения не будет ни одного афтершока рассматриваемой магнитуды, которая тем больше, чем меньше величина Λ. Эта величина при прочих равных условиях может варьировать от землетрясения к землетрясению на два порядка и более [Marsan, Helmstetter, 2017]. При фиксированном пороге рассматриваемых магнитуд относительно магнитуды соответствующего основного толчка глобальное и региональные распределения величины Λ имеют экспоненциальный вид [Шебалин и др., 2018]:
|
| (10) |
Оценкой Λ0 является среднее число афтершоков в серии. Согласно работе С.Л. Соловьева и О.Н. Соловьевой [1962] Λ0 убывает с глубиной основного толчка. Это свойство является дополнительным аргументом необходимости учета глубины в настоящем исследовании.
Таким образом, величина Λ в конкретных сериях часто имеет значения, близкие 0. Для таких серий уже через короткое время после основного толчка вероятность сильных повторных толчков невелика. Для усредненных оценок, очевидно, необходимо учитывать распределение (10). Исходя из (8) и (9) получаем усредненные функцию распределения и плотность величины τ (при условии τ < T):
|
| (11) |
|
| (12) |
Ключевым параметром в распределении (11) является величина Λ0, которая определяет значение функции вероятности при τ = 0. Чем меньше значение Λ0, тем выше вероятность вообще не иметь афтершоков соответствующей магнитуды. Величина Λ0 зависит от выбранного порога магнитуды относительно магнитуды основного толчка. Удобной величиной этого порога является 2.0. Это существенно ниже средней разности магнитуды сильнейшего афтершока и основного толчка (закон Бота [Bath, 1965; Баранов, Шебалин, 2018]), но достаточно высоко, чтобы в большинстве случаев обеспечить превышение порога представительной магнитуды. Везде в дальнейшем мы будем использовать именно такой относительный порог по магнитуде. Соответственно, для величин Λ0 и τ в (11) и (12) будут использоваться обозначения Λ2 и τ2. Для параметра Λ, определяемого для каждой серии афтершоков, мы сохраним при данном пороге то же обозначение везде далее кроме специально оговоренных случаев. Мы предполагаем [Баранов, Шебалин, 2019], что в афтершоковой последовательности магнитуды не зависят от времени и распределены по закону ГутенбергаРихтера, поэтому величины Λ0 и Λ могут быть пересчитаны для любого относительного порога магнитуды, а распределение длительности опасного периода τ затем пересчитано по формулам (8), (11).
ВЫБОР АЛГОРИТМА ИДЕНТИФИКАЦИИ АФТЕРШОКОВ
Фактические значения величины τ зависят от того, какие события считаются афтершоками определенного основного толчка. Здесь будет рассмотрено четыре разных алгоритма и построены эмпирические распределения величины τ для серий афтершоков от землетрясений мира с магнитудой выделенных этими способами.
Метод 1. Алгоритм МолчанаДмитриевой [1991], кратко описанный во Введении. Этот метод, по нашему мнению, является оптимальным при достаточно большом количестве афтершоков. Однако при малом числе афтершоков алгоритм может идентифицировать очевидные афтершоки как независимые (фоновые) события. В результате число землетрясений без афтершоков оказывается завышенным. Другой недостаток этого подхода идентификация статуса каждого события (афтершок или фоновое событие) зависит не только от предыстории, но и от будущих событий.
Метод 2. «Двойной оконный» метод. Большинство «оконных» методов являются изотропными, в них не учитываются ни протяженность очага основного толчка, ни его ориентация в пространстве. Здесь мы предлагаем очень простую модификацию «оконного» метода, в которой хотя бы частично эти недостатки устраняются. Мы полагаем, что афтершоки в течение первых 12 часов после основного толчка происходят в пределах его очага. Поэтому идентификация афтершоков производится в два этапа. На первом шаге стандартным «оконным» методом выделяются «базовые» афтершоки, которые произошли в течение 12 часов. На втором шаге тем же «оконным» методом уже в полном временном «окне» отыскиваются остальные афтершоки, при этом расстоянием до гипоцентра основного толчка считается расстояние до гипоцентра ближайшего к рассматриваемому кандидату базового афтершока.
В качестве пространственного окна принимается шар радиуса В данной задаче в качестве временного окна принимается интервал В дальнейшем мы предполагаем найти способ задания временного окна для более универсального применения этого метода.
На обоих шагах исключаются события, являющиеся афтершоками другого землетрясения с гипоцентры которых расположены на меньшем относительном расстоянии до его гипоцентра. Относительное расстояние определяется путем деления фактического расстояния на соответствующую величину
Метод 3. Афтершоки «первого ранга» модификация метода ЗаляпинаБен-Зиона [Zaliapin, Ben-Zion, 2013; 2016]. Афтершоками землетрясения i магнитуды считаются события с величиной (4), не превышающей значение Для исключения случайной идентификации в качестве афтершоков далеких событий вводится дополнительное ограничение Исключаются также события с индексом j и величиной если есть такое событие k, не являющееся афтершоком события i, что и Этот подход уже использовался в работах [Shebalin, Narteau, 2017; Шебалин и др., 2018].
Метод 4. Афтершоки с учетом протяженного очага основного толчка модификация метода ЗаляпинаБен-Зиона [Zaliapin, Ben-Zion, 2013; 2016]. В этом методе мы используем идею «двойного оконного» метода. Отличие состоит в том, что вместо пространственного и временного окон используется порог для величины (4). На первом шаге, тем не менее, вводится ограничение по времени 12 часов. На обоих шагах, как и в методе 3, вводится дополнительное ограничение Для исключения афтершоков других землетрясений вместо относительного расстояния до ближайшего базового афтершока используется функция близости (4). При этом в качестве потенциальных основных толчков исключаемого события, так же как и в методе 3, могут быть любые предшествующие события с магнитудой большей или равной его магнитуде и не являющиеся рассматриваемым основным толчком, ни его афтершоками.
На рис. 1 для каждого из четырех методов показаны эмпирические функции распределения времени последнего на интервале [0, 365 сут.] афтершока магнитуды а также суммарные функции распределения всех времен t и гипоцентральных расстояний r афтершоков относительно соответствующих основных толчков и кумулятивные графики повторяемости числа афтершоков K2 с магнитудой Параметры α и для методов 2, 3 и 4 выбраны таким образом, чтобы распределения времен t и расстояний r минимально отличались от соответствующих распределений для метода 1.
Как видно из рисунка, распределение величины для афтершоков, выделенных «двойным оконным» методом, значительно отличается. При сут. начинает увеличиваться наклон графика, что, скорее всего, означает преобладание фоновых событий для этих времен. Таким образом, метод 2 требует доработки и пока не может использоваться для целей данной работы. Методы 3 и 4 дают очень близкие распределения величины При этом обращает на себя внимание огромное различие значений параметра Значения параметров для метода 3 совпадают со значениями по работе [Zaliapin, Ben-Zion, 2016], а в методе 4 значение почти на два порядка меньше. Поскольку в методе 4 неявным образом учитываются размеры и форма очага основного толчка, в данной работе мы отдаем предпочтение именно этому методу, и все дальнейшие результаты получены для афтершоков, выделенных этим методом.
Рис. 1. Эмпирические функции распределения длительности опасного периода афтершоков τ2 (а), времен (б) и гипоцентральных расстояний (в) для всех рассмотренных афтершоков относительно соответствующих основных толчков и кумулятивный график повторяемости числа афтершоков магнитуды M ≥ Mm -2 для всех рассмотренных основных толчков (г). Сплошная линия – афтершоки выделены по методу 1, штриховая линия – по методу 2 (α = 0.015 сут.), пунктирная – по методу 3 (η0 = 10–4, α = 0.015 сут., df = 1.3) штрих-пунктирная – по методу 4 (η0 = 0.000002, α = 0.02 сут., b = 1, df = 1.3).
УСРЕДНЕННАЯ МОДЕЛЬ ДЛИТЕЛЬНОСТИ ОПАСНОГО ПЕРИОДА
Как было обосновано во введении, в усредненной модели необходимо учитывать глубину очага. Усредненная модель для афтершоков с определяется формулой (11) и, с учетом распределения времен афтершоков по закону ОмориУтсу (5), зависит от трех параметров: Ключевым является параметр так как именно этим параметром определяется вероятность полного отсутствия афтершоков с (значение функции (11) при Параметры c и p влияют на форму распределения: чем больше p и/или меньше с, тем быстрее значение функции распределения (11) приближается к 1 и тем меньше наиболее вероятное значение времени последнего на интервале (0, Т) афтершока. Стоит отметить, что такая зависимость τ от параметра c соответствует модели Дитриха (3) и представлениям о том, что параметр c определяет масштаб времени афтершокового процесса.
В работе [Шебалин и др., 2018] было показано, что параметр может иметь значительные региональные отличия. Хорошо известно также, что количество афтершоков у более глубоких землетрясений обычно меньше по сравнению с менее глубокими очагами [Соловьев, Соловьева, 1962]. Мы построили (рис. 2) зависимость от глубины очага параметра для афтершоковых последовательностей от сильных землетрясений мира и Камчатки и Курильских островов В качестве оценки величины принималось среднее по всем сильным землетрясениям количество афтершоков магнитуды Действительно, оказалось, что значения параметра могут меняться в зависимости от глубины в десятки раз. При этом в обоих случаях оказалось, что эта зависимость имеет экспоненциальный характер в интервале глубин от 10 до 100 км.
Рис. 2. Зависимость параметра Λ2 от глубины очага для афтершоковых последовательностей: (а) от землетрясений мира с Mm ≥ 6.5 1980–2018 гг.; (б) от землетрясений с в Курило-Камчатском регионе. Сплошная линия – оценки параметра Λ2 по основным толчкам, упорядоченным по возрастанию глубины очага в скользящем окне 50 событий с шагом 5 событий (в качестве значения глубины принимается среднее по 50 событиям); штриховая линия – кусочно-линейная аппроксимация в логарифмическом масштабе (см. текст).
Эти зависимости в логарифмическом масштабе глубин могут быть аппроксимированы кусочно-линейными соотношениями (13) для глобального каталога и (13а) для Камчатки и Курильских островов (рис. 2).
|
| (13) |
|
| (13 a) |
Для проверки, насколько модель (11) совпадает с эмпирическими данными по глобальному каталогу мы рассмотрели 4 интервала глубины и оценили для них значение (табл. 1). Для Байкальского и Кавказского регионов данных для оценки зависимости параметра от глубины недостаточно, поэтому в этих регионах мы используем полученные общие по всем глубинам оценки (табл. 1).
Оценить реальные значения параметра с для афтершоковых серий после сильных землетрясений не представляется возможным, но возможно получить косвенные оценки параметра, используя афтершоки небольших землетрясений, рассматривая совместно большое число серий [Narteau et al., 2009]. В работе [Shebalin, Narteau, 2017] была установлена зависимость параметра с от глубины основного толчка. Более того, оказалось, что параметр с мало зависит от магнитуды как основного толчка, так и афтершоков. Для разломов сдвигового типа в Калифорнии значения c варьируют в пределах 104101 сут., демонстрируя устойчивую тенденцию уменьшения с увеличением глубины в интервале 215 км, в котором в этом регионе происходит основная масса землетрясений. В процессе подготовки публикации была получена аналогичная тенденция для зоны субдукции в Японии (рис. 3а) в интервале до 40 км, в котором параметр c уменьшается от 102 до 103 сут. Здесь мы повторили процедуру анализа этой работы по данным каталога землетрясений Курило-Камчатского региона [Каталог землетрясений Камчатки …] (рис. 3б). Оказалось, что в интервале глубин до 40 км параметр с аналогичным образом убывает от значения 102 сут. до значения 103 сут., а в интервале от 40 до 100 км снова несколько увеличивается. Можно предположить, что в других зонах субдукции параметр с изменяется с глубиной аналогичным образом примерно в таких же пределах. Землетрясения в зонах субдукции составляют основную массу землетрясений мира, поэтому в среднем по всем сильным землетрясениям величина с должна иметь тот же порядок значений. С учетом этого мы приняли значения параметра c для четырех интервалов глубины, используемых для проверки модели (табл. 1).
Таблица 1. Усредненные значения параметров для оценки длительности опасного периода афтершоков
Интервал глубин для мирового каталога или регион | Λ2 | с, сут. | p | Число основных толчков |
0–10 км | 5.0 | 0.01 | 1.25 | 233 |
10–30 км | 8.5 | 0.005 | 1.14 | 674 |
30–50 км | 6.6 | 0.001 | 1.06 | 410 |
50–100 км | 1.75 | 0.01 | 1.00 | 401 |
Курилы-Камчатка, Mm≥6.0 | 12.0 | 0.002 | 1.20 | 218 |
Кавказ, Mm≥5.0 | 2.8 | 0.001 | 1.05 | 152 |
Байкал, Mm≥5.5 | 3.2 | 0.001 | 1.11 | 443 |
Рис. 3. Зависимость оценок параметра c закона Омори–Утсу от глубины очага основного толчка, полученных по методике работы [Shebalin, Narteau, 2017]: (а) – землетрясения в зонах субдукции вблизи Японии; (б) – землетрясения Курило-Камчатского региона. На графиках отмечены оценки максимального правдоподобия (кружки) и доверительные интервалы на уровне 95%.
Для Кавказа и Байкальского региона имеющиеся данные не позволяют оценить зависимость параметра с от глубины, поэтому по методике работы [Shebalin, Narteau, 2017] мы получили лишь усредненные оценки (табл. 1).
Для оценки параметра p закона ОмориУтсу (1) использование афтершоковых серий от сильных землетрясений вполне допустимо при условии исключения из рассмотрения начальной части серии. Для глобального каталога и Курило-Камчатских данных мы получили зависимости параметра p от глубины очага (рис. 4). Для глобального каталога, кроме того, были получены аналогичные оценки для четырех рассматриваемых интервалов глубины (табл. 1). Рассматривались афтершоки с относительным порогом магнитуды При таком пороге каталог афтершоков можно считать полным, начиная с момента сут. (см. ниже). Оценки параметра p были получены с помощью Байесовского подхода [Holschneider et al., 2012] в интервале [0.01, 365] сут. Оценки проводились для стека всех афтершоков отобранной группы основных толчков [Narteau et al., 2009].
Аппроксимация зависимостей параметра p от глубины очага кусочно-линейной функцией в логарифмическом масштабе определяется соотношениями (14) для глобального каталога и (14а) для Камчатки и Курильских островов (рис. 4).
|
| (14) |
|
| (14a) |
Рис. 4. Зависимость параметра p от глубины очага для афтершоковых последовательностей (а) от землетрясений мира с Mm ≥ 6.5 1980–2018 гг.; (б) от землетрясений с Mm ≥ 6.0 в Курило-Камчатском регионе. Кружками отмечены оценки параметра p по основным толчкам, упорядоченным по возрастанию глубины очага в скользящем окне 50 событий с шагом 5 событий (в качестве значения глубины принимается среднее по 50 событиям); штриховая линия – кусочно-линейная аппроксимация в логарифмическом масштабе (см. текст).
Для Кавказа и Байкальского региона были получены усредненные для всего диапазона глубин оценки параметра p. Аналогичная усредненная оценка была получена и для Курило-Камчатского региона (табл. 1).
Проверку модели распределения величины мы провели путем сравнения теоретических и эмпирических функций распределения этой величины для четырех интервалов глубины по глобальному каталогу и по трем региональным каталогам (рис. 5). Теоретические распределения построены по формуле (11) со значениями параметров из табл. 1. Эмпирические распределения в целом значительно различаются между собой, но при этом теоретические распределения отличаются от соответствующих эмпирических в значительно меньшей степени. Наилучшее совпадение достигается для h ≥ 50 км по глобальному каталогу (доверительный уровень 0.99 по критерию Колмогорова), наихудшее для h = 3050 км (доверительный уровень 0.52). С учетом того, что в каждой серии афтершоков возможны локальные по времени отклонения от закона ОмориУтсу (1) и что зависимости модели (11) от параметров с и p этого закона существенно не линейны, полученные результаты можно считать хорошим подтверждением модели. Хорошее совпадение теоретических и эмпирических распределений также подтверждает правомерность использования косвенных оценок параметра c, поскольку значения этого параметра существенно влияют на форму распределения.
Рис. 5. Эмпирические и теоретические функции распределения длительности опасного периода афтершоков τ2 для глобального каталога, Mm ≥ 6.5 (а) и (б): интервалы глубин основного толчка км (кружки), км (треугольники), км (знаки «+»), 50 км ≤ h (крестики) и региональных каталогов (в): Курило-Камчатского региона, Mm ≥ 6.0 (жирная линия), Байкальского региона, Mm ≥ 5.5 (тонкая линия), Кавказского региона, Mm ≥ 5.0 (линия средней жирности). Теоретические функции распределения (11) с параметрами из табл. 1 показаны штриховыми линиями.
ОЦЕНКИ ДЛИТЕЛЬНОСТИ ОПАСНОГО ПЕРИОДА С УЧЕТОМ ИНФОРМАЦИИ О ПЕРВЫХ АФТЕРШОКАХ
Можно предположить, что информация о первых афтершоках может существенно улучшить оценки длительности опасного периода. Время tstart, начиная с которого каталог можно считать полным, зависит от разности магнитуды основного толчка Mm и порога магнитуды афтершоков Mc [Helmstetter et al., 2006; Баранов и др., 2019]. В работе [Баранов и др., 2019] для глобального каталога ANSS получено соотношение:
|
| (15) |
Оценкой параметра Λ в (7) при любом пороге магнитуды служит число афтершоков в исследуемом интервале [Holschneider et al., 2012]. Эта оценка может быть пересчитана для произвольного интервала. Величина Λ в соотношении (8) относится к интервалу [0, T]. Для оценки параметра Λ могут использоваться данные об афтершоках представительной магнитуды. В этом случае необходимо оценить (или задать) также параметр b закона ГутенбергаРихтера и пересчитать в соответствии с этим законом величину Λ для интересующей магнитуды. Таким образом, если известно число афтершоков магнитуды в интервале то оценка величины Λ в соотношении (8) для относительного порога определяется по формуле:
|
| (16) |
Параметры c и p при наличии достаточного количества данных на интервале также могут быть оценены [Баранов и др., 2019]. Однако, как указывалось выше, оценка параметра с в данной задаче не применима, и значения этого параметра следует задавать по косвенным данным. Значение параметра b оценивалось по методике работы [Bender, 1983] (краткое описание метода приведено в работе [Vorobieva et al., 2013]) c использованием разных вариантов априорного распределения аналогично работе [Баранов и др., 2019]. Мы провели серию экспериментов с различными способами оценки или задания параметров c, p, b. Для глобального и региональных каталогов оптимальными оказались варианты с оценкой параметра Λ по формуле (16), но с использованием заданных значений с, p и b. Это связано с тем, что оценки параметров для отдельных серий недостаточно устойчивы.
Оценивание результатов с использованием данных о первых афтершоках, аналогично работе [Баранов и др., 2019], мы проводили по критерию информационного выигрыша LG (N) [Schorlemmer et al., 2007] оцениваемого метода по отношению к референц-модели, определяемого как отношение правдоподобия полученных реализаций для двух моделей в пересчете на один прогноз:
|
| (17) |
В формуле (17) в качестве референц-модели (обозначена используется усредненная модель (12), зависящая от глубины очага сильного землетрясения. Мы рассмотрели ретроспективные прогнозы величины для одного фиксированного времени t0 = 0.5 сут. Это время представляется нам оптимальным с практической точки зрения. Для глобального каталога и для Курило-Камчатского региона параметр в референц-модели и параметр p для обеих моделей определялись из соотношений (13) и (14) или (13а) и (14а) соответственно, а значения параметра с для обеих моделей в обоих случаях брались из табл. 1. Для регионов Байкал и Кавказ для обеих моделей значения параметров с и р и для референц-модели параметра определялись из табл. 1. В тестируемой модели, аналогично работе [Баранов и др., 2019], оценка параметра проводилась при условии Информационный выигрыш (17) оценивался для совокупности только таких серий афтершоков. Результаты приведены в табл. 2.
Как следует из табл. 2, выигрыш при использовании данных об афтершоках за первые 0.5 сут. значителен по сравнению с усредненной моделью, он составляет в среднем около 50 %.
Предложенная здесь методика позволяет оценить функцию распределения величины На практике при прогнозе в реальном времени может быть удобнее оперировать конкретными значениями. Следуя подходу работы [Баранов и др., 2019], на основе функций распределения можно определить «мягкий», «жесткий» и «нейтральный» прогнозы как 10%, 90% и 50% квантили.
Таблица 2. Результаты ретроспективного теста оценок длительности опасного периода афтершоков с использованием информации об афтершоках в интервале [0, 0.5 сут.]
Каталог | Число серий N с | LG (N) |
Глобальный каталог ANSS, Mm≥6.5 | 181 | 1.48 |
Курило-Камчатский регион, Mm≥6.0 | 32 | 1.59 |
Кавказ, Mm≥5.0 | 10 | 1.39 |
Байкал, Mm≥5.5 | 4 | 1.64 |
ЗАКЛЮЧЕНИЕ
В работе разработана методика оценки длительности опасного периода возможного возникновения афтершоков магнитуды При необходимости оценка может быть пересчитана для произвольного порога магнитуды. Оценка проводится в 2 этапа. На первом этапе, сразу после сильного землетрясения, используется усредненная модель распределения величины Входной информацией для этой модели может являться только глубина очага. Параметры модели заранее определены. На втором этапе прогноз величины может быть уточнен за счет информации о первых афтершоках, например, в течение первых 12 часов после землетрясения. Параметры модели для серий с большим числом афтершоков в принципе могут быть оценены, но оказалось, что результаты лучше совпадают с реальными значениями времени последнего афтершока, если использовать заранее заданные значения параметров, по возможности зависящие от глубины очага. Поэтому фактически для оценок на втором этапе используется только число афтершоков представительной магнитуды.
Для практического применения методики оценены параметры модели по глобальному каталогу землетрясений и отдельно для Курило-Камчатского, Байкальского и Кавказского регионов. Для глобального каталога и Курил, и Камчатки параметры зависят от глубины очага. Ретроспективная проверка первого этапа методики показала хорошее совпадение теоретической модели с эмпирическими распределениями величины Уточнение прогноза на втором этапе, как показал второй ретроспективный тест, действительно целесообразно, поскольку позволяет получить информационный выигрыш в среднем около 50%.
Целью данной работы являлось, главным образом, создание методики для получения чисто практических результатов, которые, в том числе, будут использованы в разрабатываемой автоматизированной системе прогнозирования опасности афтершоков AFCAST (www.afcast.org). Вместе с тем, в работе получены новые результаты, которые могут послужить уточнению физических представлений о сейсмогенезе. В частности, установлена сильная зависимость от глубины очага среднего количества афтершоков с относительным порогом магнитуды. Эта величина экспоненциально уменьшается с глубиной, изменяясь более чем на порядок.
Важным отрицательным результатом работы является вывод о неправомерности использования широко распространенной, особенно в западной научной литературе, модели ETAS для прогноза длительности опасного периода афтершоков. Главная причина в том, что количество афтершоков у каждого события определенной магнитуды считается параметром модели, однако в природе эта величина варьирует в очень широких пределах, подчиняясь экспоненциальному распределению, имеющему максимум в нуле.
Значительное внимание в работе уделено проблеме идентификации афтершоков. Была проведена модификация метода ЗаляпинаБен-Зиона [Zaliapin, Ben-Zion, 2013; 2016], в которой неявным образом учитываются размеры и направление простирания очага сильного землетрясения. Этот алгоритм позволяет выделять последовательности афтершоков, близко совпадающие при большом числе афтершоков c последовательностями, выделенными с методом МолчанаДмитриевой [1991], при этом алгоритм лишен его недостатков при малом числе афтершоков. Оказалось, однако, что пороговое значение функции близости примерно на два порядка меньше порога в стандартном методе. Это фактически означает, что «иерархическая» схема выделения афтершоков в методе ЗаляпинаБен-Зиона в виде ветвящейся цепочки «родительпотомки» фактически имеет лишь чисто технический смысл, без фактичексих причинно-следственных связей. Таким образом, методика ЗаляпинаБен-Зиона и иерархические модели афтершокового процесса вообще требуют переосмысления.
Финансирование работы
Исследование выполнено при поддержке грантов Министерства образования и науки 14. W03.31.0033 и РФФИ 17-05-00749 в рамках Госзаданий по темам НИР АААА-А19-119011490 127-6 и 0152-2019-0010.
БЛАГОДАРНОСТИ
Авторы выражают благодарность двум анонимным рецензентам за ценные замечания.