Forecasting aftershock activity: 4. Estimating the maximum magnitude of future aftershocks
- Authors: Baranov S.V.1,2, Pavlenko V.A.1,3, Shebalin P.N.1
-
Affiliations:
- Institute of Earthquake Prediction Theory and Mathematical Geophysics,Russian Academy of Sciences
- Kola Branch, Federal Research Center “Geophysical Survey of Russian Academy of Sciences”
- Schmidt Institute of Physics of the Earth, Russian Academy of Sciences
- Issue: No 4 (2019)
- Pages: 15-32
- Section: ARTICLES
- URL: https://journals.eco-vector.com/0002-3337/article/view/13432
- DOI: https://doi.org/10.31857/S0002-33372019415-32
- ID: 13432
Cite item
Full Text
Abstract
In this paper, we consider the problem of forecasting the magnitude of the future, starting from a certain instant of time, strongest aftershock. This problem is topical since the later strong aftershocks occur against the background of the less frequently repeating shocks, are less expected and thus pose an independent hazard. At the same time, the magnitudes of the strongest aftershocks decrease with time after the main shock. The purpose of accurate forecasting is to minimize the underestimation or overestimation of the magnitude of future risks. In this study, the aftershock process is represented by the superposition of the Gutenberg–Richter and Omori–Utsu laws whose parameters are estimated by the Bayess method using the data on the aftershocks that have already occurred to a given time point and the a priori information about the probable values of the parameters. This significantly improves the forecast compared to the estimates that are based on the magnitude of the main shock alone. The quality of forecasting is estimated relative to the Båth’s dynamic law with the use of two independent criteria. The first criterion is based on the similarity estimates, and the second, on the error diagram.
Full Text
ВВЕДЕНИЕ
В этой статье, продолжая серию по проблеме прогнозирования опасности повторных толчков после сильных землетрясений, мы разрабатываем метод оценивания для произвольного момента времени максимальной магнитуды афтершоков, которые произойдут после этого момента, с использованием данных об афтершоках, которые произошли до этого момента.
Задача оценки магнитуды сильнейшего афтершока, на первый взгляд, аналогична хорошо известной задаче оценки максимальной возможной магнитуды землетрясения [Holschneider et al., 2014; Kijko, 2004; Knopoff, Kagan, 1977; Pisarenko et al., 1996; Писаренко и др., 2017]. Эти две задачи, тем не менее, имеют принципиальное отличие. Количество афтершоков конечно, поэтому для задачи оценки магнитуды сильнейшего афтершока отклонение графика повторяемости от прямолинейной формы в области больших магнитуд не является определяющим, в отличие от задачи оценки максимально возможной магнитуды землетрясений, в которой основные усилия исследователей направлены именно на исследование и моделирование этого отклонения.
Рассмотрение проблемы, зависящей от времени оценки магнитуды сильнейшего афтершока, было начато в предыдущей статье [Баранов, Шебалин, 2018]. Мы получили теоретическую зависимость от времени для аналогичных оценок, но без использования информации об уже произошедших афтершоках. Теоретические распределения для разных моментов времени хорошо совпали с наблюденными. Эта модель необходима, прежде всего, как «нулевое приближение», относительно которого можно оценивать любую другую методику зависящих от времени оценок максимальной магнитуды афтершоков, в том числе разработанную в рамках данной работы. Кроме того, когда данных недостаточно, модель может быть непосредственно использована для прогноза. Необходимость «нулевого приближения» возникла на этапе тестирования и окончательной отладки предлагаемой здесь методики во времени, близком к реальному, в рамках разрабатываемой автоматизированной системы прогнозирования опасности афтершоков AFCAST (www.afcast.org).
В данной работе мы опираемся на важные результаты, полученные в указанной и других недавних работах. Во-первых, мы продолжаем следовать гипотезе РизенбергаДжонс [Reasenberg, Jones, 1989] о возможности использования прямой суперпозиции законов ГутенберхаРихтера [Gutenberg, Richter, 1954] и ОмориУтсу [Utsu, 1961] для моделирования афтершоковых процессов [Баранов, Шебалин, 2016]. Во-вторых, мы используем подход Д. Вере-Джонса [Vere-Jones, 1969; 2008] для расчетов распределения вероятности магнитуды сильнейшего афтершока в форме двойной экспоненты. В этом подходе распределение вероятности магнитуды сильнейшего афтершока зависит только от ожидаемого числа афтершоков магнитуды выше заданного порога и распределения магнитуд (закон ГутенбергаРихтера). В предыдущих работах нам удалось получить важные подтверждения применимости обоих подходов. В работе [Баранов, Шебалин, 2019] было показано, что магнитуды сильных афтершоков не зависят от времени, поскольку времена сильнейших, а также вторых по силе, третьих и т.д. афтершоков следуют тому же распределению, что и все афтершоки представительной магнитуды. Комбинация формулы Д. Вере-Джонса с обнаруженным нами в глобальной статистике экспоненциальным распределением числа афтершоков с магнитудой выше относительного порога (по отношению к магнитуде основного толчка) дает хорошее совпадение теоретического и эмпирического распределений относительной магнитуды сильнейшего афтершока [Шебалин и др., 2018]. Это дало возможность не только обосновать закон Бота [Bath, 1965], то есть величину среднего значения относительной магнитуды сильнейшего афтершока, но и впервые объяснить форму соответствующего распределения. В работе [Баранов, Шебалин, 2018] спадание магнитуды сильнейшего афтершока с течением времени (динамический закон Бота) хорошо моделируется комбинацией формулы Вере-Джонса и закона ОмориУтсу. Показано, что форма распределения магнитуды сильнейшего, начиная с некоторого момента времени t, афтершока не меняется со временем t, а сдвиг этого распределения в сторону от магнитуды основного толчка определяется только параметрами законов ОмориУтсу и ГутенбергаРихтера, что также подтверждает справедливость подхода Д. Вере-Джонса.
Полученные результаты позволяют в данной работе обоснованно опираться на подход Вере-Джонса, в котором распределение вероятности магнитуды сильнейшего афтершока определяется тремя величинами: ожидаемым числом афтершоков представительной магнитуды, значением представительной магнитуды и параметром b закона ГутенбергаРихтера. Ожидаемое число афтершоков представительной магнитуды оценивается на основании закона ОмориУтсу с учетом числа афтершоков, которые уже произошли. Таким образом, задача сводится к оценке по данным об уже состоявшихся афтершоков параметров законов ОмориУтсу и ГутенбергаРихтера. В данной работе для этого мы совершенствуем подход работы [Баранов, Шебалин, 2016].
Важным элементом данной работы является использование Байесовского подхода для получения оценок параметров. Работа над данной статьей была начата еще до начала подготовки предыдущей публикации серии [Баранов, Шебалин, 2018]. На первом этапе использовались точечные оценки максимального правдоподобия параметров, что позволяло напрямую использовать полученные аналитически квантильные функции для магнитуды сильнейшего афтершока. Однако тесты показали, что такой подход может приводить к значительным ошибкам из-за существенно нелинейной зависимости оценок от значений параметров. При значительном отклонении оценки какого-либо параметра от истинного значения это приводит к сильному искажению оценки искомой величины. Использование Байесовского подхода на следующем этапе исследования позволило преодолеть это препятствие благодаря использованию распределений оценок параметров, а не их наиболее вероятных значений. Однако этот подход оказался намного более трудоемким и потребовал проведения значительного количества тестов. Одновременно встал вопрос, как оценивать результаты прогнозов магнитуды сильнейшего афтершока, как сравнивать разные варианты оценок. Прежде всего, была необходима простая модель, подтверждаемая эмпирическими данными, относительно которой можно было бы проводить оценивание. Такой моделью и стал динамический закон Бота, обоснованный в работе [Баранов, Шебалин, 2019].
Для оценивания прогнозов магнитуды будущего сильнейшего афтершока здесь мы предлагаем два критерия. Один критерий основан на отношении функций правдоподобия анализируемой оценки и оценки по динамическому закону Бота при фактических значениях максимальной магнитуды [Schorlemmer et al., 2007]. Другой критерий использует диаграмму ошибок [Molchan, 1991] относительно динамического закона Бота.
Ретроспективное тестирование различных вариантов оценок с помощью разработанных критериев и динамического закона Бота в качестве референц-модели дало возможность поэтапно оптимизировать разные элементы алгоритма. Существенный скачок в качестве оценок произошел в результате нахождения оптимальной схемы учета неполноты каталога афтершоков сразу после основного толчка. Другое качественное улучшение произошло в результате регуляризации оценок параметров при использования информативных априорных распределений на основе глобальной статистики значений параметров ОмориУтсу и ГутенбергаРихтера. Дальнейшие тесты показали, что в таком варианте вариации параметра b, заложенные в Байесовские оценки максимальной магнитуды, значительно сильнее влияют на результат, чем фактические ошибки оценки этого параметра, и в окончательном варианте нашего подхода мы вернулись к использованию точечных оценок параметров, получаемых, тем не менее, с помощью Байесовского подхода.
Одновременно с нами методика оценки максимальной магнитуды будущих афтершоков, по данным об уже произошедших, разрабатывалась объединенной группой специалистов из Канады и Японии [Scherbakov et al., 2017], которые также использовали формулу Вере-Джонса, законы ГутенбергаРихтера и ОмориУтсу и Байесовский подход. Главные отличия нашего подхода состоят в следующем. Во-первых, мы оцениваем параметры с учетом неполноты каталога в начале серии, для чего оцениваем не только представительную магнитуду, но и время после основного толчка, начиная с которого эта магнитуда действительно является представительной. Во-вторых, для регуляризации оценок параметров мы используем распределения этих параметров на основе глобальной статистики. В-третьих, Байесовский подход мы используем только для регуляризации оценок параметров законов ГутенбергаРихтера и ОмориУтсу, но не для оценки распределения магнитуды сильнейшего афтершока. В работе мы проводим сравнение наших оценок и оценок по работе [Scherbakov et al., 2017].
Данная работа имеет, главным образом, практическую направленность, и ее результаты напрямую могут быть использованы для оценки сейсмической опасности после сильных землетрясений. Вместе с тем, полученные результаты расширяют физические представления о сейсмогенезисе.
ИСХОДНЫЕ ДАННЫЕ
Для ретроспективной проверки разработанной здесь методики оценки максимальной магнитуды последующих афтершоков мы используем каталог землетрясений ANSS ComCat геологической службы США (USGS) [ANSS…]. Для построения эмпирических закономерностей в работе рассматривается 777 серий афтершоков от землетрясений мира с магнитудой 6.5 и выше за период с 1980 по 2016 гг. Выделение основных толчков и их афтершоков проводилось по алгоритму из работы [Молчан, Дмитриева, 1991] с помощью программы В.Б. Смирнова [2009].
ОСНОВНОЙ МЕТОД ОЦЕНКИ
Мы предполагаем [Баранов, Шебалин, 2018], что в афтершоковой последовательности времена и магнитуды независимы и распределены, соответственно, по закону ОмориУтсу [Utsu, 1961]:
λ(t) = K/(t + c)p (1)
и закону ГутенбергаРихтера [Gutenberg, Richter, 1954]:
(2) |
где: t время от основного толчка, сут.; λ(t) число афтершоков c магнитудой в единицу времени; Mc уровень представительной магнитуды; с, p, b параметры; F(M) функция распределения вероятности магнитуды. Введем также дополнительные обозначения: Mm магнитуда основного толчка в серии афтершоков; T период рассмотрения, также от момента основного толчка, сут.; везде в этой работе мы принимаем значение T = 365 сут.; M1(t, T) магнитуда сильнейшего на интервале (t, T) афтершока.
Для прогноза значения M1(t, T) мы используем информацию о временах и магнитудах афтершоков, которые произошли на интервале (tstart, t). Задержка tstart > 0 необходима, поскольку каталог землетрясений сразу после основного толчка неполон; ввиду важности данного обстоятельства мы рассмотрим проблему выбора параметра tstart в отдельном разделе. Проинтегрировав выражение (1) на интервалах (tstart, t) и (t, T) и проведя преобразования, получим [Holschneider et al., 2012]:
| (3) |
где и ожидаемое число афтершоков с на соответствующих интервалах:
| (4) |
Оценкой величины является фактическое зарегистрированное на интервале количество афтершоков. Если параметры с и p известны, то выражения (3) и (4) позволяют рассчитать оценку ожидаемого на интервале количества афтершоков В сделанных предположениях вероятность того, что магнитуда произвольного афтершока в серии меньше значения M, определена формулой (2); вероятность того, что k афтершоков имеют магнитуду меньше M, равна Предположив, что число афтершоков на интервале имеет распределение Пуассона, и просуммировав по формуле полной вероятности все варианты значения k, получим распределение вероятности магнитуды M1(t, T) [Vere-Jones, 2003]:
(5)
Объединив (5) и (2) имеем распределение в виде двойной экспоненты:
| (6) |
Искомое оценочное распределение величины можно получить, подставив в (6) оценку по формулам (3) и (4), используя точечные оценки параметров c, p и b и фактическое число афтершоков представительной магнитуды на интервале в качестве оценки величины
Обратная функция распределения (6) имеет явное выражение, которое удобно использовать для нахождения квантилей распределения:
| (6a) |
С учетом сильно нелинейной зависимости (6) от параметров и опасаясь из-за этого значительного смещения оценок в случае больших фактических ошибок оценки параметров, первоначально мы усложнили процедуру. Вместо точечных оценок параметров c, p и b, мы стали использовать апостериорные Байесовские распределения оценок (c, p) по работе [Holschneider et al., 2012] и оценки b по работе [Bender, 1983], многократно генерируя методом Монте-Карло, в соответствии с этими распределениями, значения параметров и усредняя распределение (6). Фактически этот подход эквивалентен подходу работы [Scherbakov et al., 2017], в которой этапы получения Байесовских оценок параметров встроены в многократное интегрирование по формуле Байеса для получения апостериорного распределения оценки величины В нашем подходе затратное по времени численное многократное интегрирование заменялось без потери точности расчетами методом Монте-Карло с генерацией случайных значений параметров в соответствии с их апостериорными распределениями. Преимущество нашего подхода состоит также в возможности контроля оценок параметров c, p и b. Также как в работе [Scherbakov et al., 2017] мы на этом этапе использовали неинформативные априорные распределения параметров, то есть равномерные на заданных отрезках распределения.
Тестирование метода с помощью разработанных для этого критериев (см. следующий раздел) показало, что получаемые оценки параметров неустойчивы и часто приводят к нереалистичным оценкам величины Для регуляризации оценок мы стали использовать более подробную априорную информацию о параметрах, построив распределения оценок параметров по глобальным данным (см. ниже). Дальнейшее тестирование показало, что в слишком сложной схеме Байесовских оценок нет необходимости. Использование точечных оценок параметров (но с их регуляризацией с помощью Байесовского подхода с информативными априорными распределениями) оказалось предпочтительным.
ОЦЕНИВАНИЕ РЕЗУЛЬТАТОВ
Результаты прогнозов всегда целесообразно оценивать относительно простой модели, отражающей текущее на данный момент знание о прогнозируемом процессе, не вызывающей сомнения. Хорошо известен закон Бота [Bath, 1965], в соответствии которым средняя разность значений магнитуды основного толчка и сильнейшего афтершока составляет в среднем 1.2. В работе [Баранов, Шебалин, 2019] на основе глобальной статистики афтершоковых последовательностей исследована зависимость магнитуды сильнейшего афтершока от времени. Полученные результаты авторы охарактеризовали как динамический закон Бота. Показано, что с течением времени после основного толчка магнитуда сильнейших афтершоков в среднем уменьшается, но при этом форма распределения этой величины (сходная с распределением Гаусса) со временем не меняется и среднеквадратическое отклонение сохраняет значение примерно 0.65.
В работе [Шебалин и др., 2018] установлено, что количество k непосредственных афтершоков магнитуды выше порога определяемого относительно магнитуды основного толчка Mm, при глобальном и региональном рассмотрении подчиняется экспоненциальному закону c плотностью Это дает объяснение форме распределения магнитуд сильнейших афтершоков, похожего на распределение Гаусса и теоретически определяемого соотношением:
(7а)
где параметр экспоненциального распределения при относительном пороге магнитуды При = 2 для глобальной статистики афтершоков в течение года после землетрясений с получена оценка Как оказалось, формула (7а) может быть упрощена. Комбинация формулы Вере-Джонса (6) и экспоненциального распределения дает простое и красивое соотношение, численно мало отличающееся от (7а):
| (7б) |
Распределение (7б) справедливо и для произвольного интервала времени (t, T). В этом случае величина в (7б) должна быть заменена на зависимое от времени значение Зависимость от времени этой величины определяется законом ОмориУтсу:
(8)
где функция определена соотношениями (4). Таким образом, величина M1(t, T) в динамическом законе Бота имеет распределение:
| (9) |
Это распределение имеет моду, определяемую соотношением:
| (10) |
В работе [Баранов, Шебалин, 2019] были получены следующие оценки для усредненной афтершоковой серии: b = 1.0, с = 0.04 сут., p = 1.016. Формула (9) с этими параметрами и значением дает близкое совпадение с эмпирическими распределениями, полученными в указанной работе; это распределение также мало отличается от аппроксимации эмпирических распределений нормальным. В данной работе мы используем динамический закон Бота, определяемый соотношением (9) в качестве референц-модели для оценивания прогнозов по модели, учитывающей данные об уже состоявшихся афтершоках. В дальнейшем мы будем использовать обозначение плотности этого распределения
Аналогично работе [Баранов, Шебалин, 2016], мы используем два критерия для оценивания качества получаемых распределений величины Один из критериев основан на оценках правдоподобия [Shorlemmer et al., 2007], другой на диаграмме ошибок [Молчан, Дмитриева, 1991].
Пусть N число завершившихся прогнозов, по которым оценивается метод, M1, i, i = 1,.., N реализации прогноза, то есть фактические значения величины и плотность полученного оценочного распределения этой величины. Определим информационный выигрыш LG(N) оцениваемого метода по отношению к референц-модели как отношение правдоподобия полученных реализаций для двух моделей в пересчете на один прогноз:
| (11) |
Величина тем больше, чем относительно чаще значения попадают в область больших значений плотности вероятности в тестируемой модели по сравнению с референц-моделью и чем реже в область малых значений плотности вероятности.
Диаграмма ошибок, которую мы здесь предлагаем, несколько отличается от стандартной. Поскольку прогноз можно считать тем лучше, чем ближе к максимуму плотности (моде) прогнозного распределения оказывается реализация, то в качестве управляющего параметра мы рассматриваем здесь расстояние δM от моды. Поэтому, чтобы сравнить тестируемую модель с моделью динамического закона Бота, мы определяем величину τ (доля тревог) как вероятность попадания величины в интервал а величину ν (доля пропусков цели) число случаев, когда или где E(t) мода распределения для оцениваемой модели, а мода распределения для динамического закона Бота определена соотношением (10). Точки на диаграмме при разных значениях образуют «траекторию ошибок». Чем в среднем относительно меньше расстояние от реализаций до моды, тем дальше от диагонали будет находиться «траектория» и тем лучше тестируемый метод. Тангенс наклона прямой на диаграмме ошибок, проходящей от точки (0.1) через точку траектории ошибок, называется вероятностным выигрышем (probability gain) [Molchan, 1991; Shebalin et al., 2014]. Для удобства сравнения для разных t мы будем рассматривать эту величину при ν = 0.5, обозначив ее PG0.5. Для обоих критериев значение 1 означает, что качество прогноза величины анализируемым методом такое же, как и оценка по динамическому закону Бота, без привлечения информации об уже состоявшихся афтершоках. Значения больше 1 означают, что метод предпочтителен по отношению к референц-модели; значении меньше 1 предпочтителен, наоборот, динамический закон Бота.
Оценки, основанные на правдоподобии и на диаграмме ошибок дополняют друг друга. В работе [Shebalin et al., 2014] продемонстрировано, что в диаграмме ошибок относительно больший вес придается реализациям прогноза (успешные прогнозы и ложные тревоги) в областях пространствавремени, в которых тестируемая модель дает высокие вероятности прогнозируемых событий. Критерии, основанные на правдоподобии, наоборот, относительно больший вес придают областям низкой вероятности событий (включая пропуски цели). Для критериев, основанных на правдоподобии, кроме того часто возникает проблема нулевых или очень низких вероятностей. В случае пропуска цели ведь не так существенно насколько низко оценивалась вероятность события, важно, что она была много ниже значений, при которых события обычно реализуются. А «штраф» за один пропуск цели часто может компенсировать в критерии большое число очень успешных реализаций прогноза. Эта проблема легко компенсируется введением некоторого минимума вероятности события в модели: более низкие значения заменяются этим пороговым значением. Поскольку при этом возрастает суммарная вероятность, приходится производить перенормировку. Для критерия LG мы используем порог для плотности вероятности равный 0.001.
ОЦЕНКА ПАРАМЕТРОВ ДЛЯ СЕРИИ АФТЕРШОКОВ
В тех случаях, когда используемых данных немного, получаемые оценки параметра неустойчивы, что, с учетом нелинейности, может существенно сказываться на оценках величины даже при использовании Байесовского подхода с неинформативным (однородным) распределением параметра. Как показали ретроспективные тесты, из-за особенностей афтершоковых серий достаточно часто максимум апостериорного распределения параметра оказывается у границы задаваемых априорных пределов, что в конечном счете и приводит к ошибочным оценкам величины Для регуляризации оценок мы решили использовать более детальную априорную информацию о параметрах. Были построены эмпирические функции распределения параметров по большому числу афтершоковых серий. Из 777 серий были отобраны серии, у которых имеется хотя бы 20 афтершоков с магнитудой не ниже представительной, произошедших позднее 0.1 суток после основного толчка всего 334 серии. Представительная магнитуда Mc оценивалась методом работ [Wiemer, Wyss, 2000; Woessner, Wiemer, 2005].
Значение параметра b оценивалось по методике работы [Bender,1983] (краткое описание метода приведено в работе [Vorobieva et al., 2013]); программа, разработанная для данной работы, допускает использование как неинформативного (равномерного) априорного распределения параметра, так и произвольного информативного априорного распределения. На данном этапе нам была недоступна информация о распределении оценок параметра b, поэтому мы использовали равномерное априорное распределение на интервале [0.5, 1.5]. По 334 рассмотренным сериям было построено эмпирическое распределение полученных оценок и его аппроксимация нормальным распределением (рис. 1).
Рис. 1. Эмпирическое распределение оценок параметра b. Плотность f(b) (гистограмма) построена по 334 сериям афтершоков от землетрясений мира с Mm ≥ 6.5. Жирная линия – аппроксимация нормальным распределением со средним Eb = 1.12 и стандартным отклонением σb = 0.3.
Байесовская оценка апостериорного распределения оценок параметров с и p закона ОмориУтсу проводилась методом работы [Holschneider et al., 2012]. Следуя рекомендациям указанной работы, для параметра с поиск оптимального значения осуществлялся для величины lg(c). Для данной работы авторами разработана новая программа, допускающая использование произвольного априорного распределения параметров (в текущей версии оригинальной программы допускается использование лишь равномерного априорного распределения). Аналогично оценке параметра b, на данном этапе использовались равномерные априорные распределения параметра p в интервале [0.5, 2.5] и lg(c) в интервале [3, 1.7]. Полученные по 334 сериям эмпирическое распределение полученных оценок и его аппроксимация нормальным распределением показаны на рис. 2.
Рис. 2. Эмпирические распределения оценок параметров lg (c) (а) и p (б). Гистограммы построены по 334 сериям афтершоков от землетрясений мира с Mm ≥ 6.5. Жирная линия – аппроксимация нормальным распределением со средними Elg(c) = –1 и Ep = 1.05 и стандартными отклонениями σlg(c) = 0.74 и σp = 0.25.
На этапе оценки величины полученные аппроксимации нормальными распределениями используются в качестве информативных априорных распределений параметров b, c, p для оценки этих параметров на интервале (tstart, t) для каждой серии афтершоков.
ОПТИМИЗАЦИЯ ПАРАМЕТРА tstart
Важное значение на качество оценок величины оказывает выбор параметра tstart. Очевидно, что сразу после сильного землетрясения часть афтершоков не попадает в каталог, главным образом, из-за наложения их сейсмограмм на сейсмограмму основного толчка и взаимного наложения сейсмограмм афтершоков. Авторы работы [Helmstetter et al., 2006] установили для каталога землетрясений Калифонии приблизительную связь магнитуды основного толчка Mm, порога представительной магнитуды Mc и времени tstart, начиная с которого каталог с M ≥ Mc полон: Mc = Mm4.50.76 lg(tstart), где время выражено в сутках. Для рассмотренных случаев в ретроспективных тестах эта оценка, по нашему мнению, во многих случаях существенно занижена. В работе [Hainzl, 2016] предложен подход для оценивания пары значений (Mc, tstart), основанный на представлении о том, что каталог становится неполным, если число событий в единицу времени слишком велико. Для этого подхода необходим подбор параметров, однако универсального критерия для такого подбора нет.
Здесь мы используем две конкурирующие оценки. В обоих случаях сначала мы находим значение представительной магнитуды Mc на интервале (0.01, t) методом Maximum Curvature (MAXC) [Wiemer, Wyss, 2000; Woessner, Wiemer, 2005]. Задержка 0.01 сут. вводится для того, чтобы избежать завышения оценок Mc за счет отсутствия более слабых событий в начале серии афтершоков. Величина tstart определяется для найденной оценки Mc.
Первый подход, предложенный нами в работе [Shebalin, Baranov, 2017] и затем усовершенствованный в работе [Шебалин и др., 2018] состоит в нахождении такого значения tstart, при котором стабилизируется параметр b закона ГутенбергаРихтера. В качестве критерия стабилизации используется первое достижение краткосрочным средним значением магнитуды, получаемым усреднением магнитуды 10 последовательных афтершоков, относительно долгосрочного среднего значения, получаемого по серии афтершоков за весь рассматриваемый период (0, t). При таком подходе мы опираемся на соотношение Аки [Aki, 1965], устанавливающее однозначное соответствие средней магнитуды и параметра b, а также предполагаем, что низкие значения параметра b (и, соответственно, большие значения средней магнитуды) в начале серии обусловлены отсутствием в каталоге части более слабых событий.
Второй вариант значительно проще. Мы исходим из того, что значение параметра c для афтершоков сильного землетрясения в реальности невелико и значительно меньше получаемой оценки при типичном для региона и периода времени значениях Mc [Holschneider et al., 2012]. Фактическая оценка параметра c, таким образом, может служить определением величины tstart для соответствующего значения Mc. На интервале (0.0001, t) при M ≥ Mc оцениваются параметры с и p закона ОмориУтсу (2). Затем полагаем tstart = c и по данным из [tstart, t] с M ≥ Mc снова оцениваем параметры c, p, которые используются для расчета плотности распределения M1 (t, T). Оценки проводятся методом работы [Holschneider et al., 2012].
Мы сравнили два подхода путем интегрального сопоставления результатов ретроспективной оценки величины по данным на интервале [tstart, t] для всех 777 рассмотренных серий для 9 значений t (0.25, 0.5, 1, 2, 4, 8, 16, 32 и 64 сут.). В качестве значения T везде принималось значение 1 год. В качестве интегральной оценки принималось среднее по 9 значениям LG и 9 значениям PG0.5, соответствующих разным значениям t; обозначим эту величину LPG. Результаты по двум методам оказались малоразличимыми (табл. 1). Для сравнения мы рассмотрели также вариант, в котором значение tstart находилось по соотношению работы [Helmstetter et al., 2006]. Результат неожиданно оказался лишь незначительно хуже (табл. 1). Потребовалось провести ряд дополнительных численных экспериментов, чтобы выяснить причину малого расхождения качества оценок, полученных при использовании трех вариантов определения величины tstart. В частности, были рассчитаны варианты с несколькими фиксированными значениями tstart, не зависящими ни от Mm, ни от Mc. Все оценки интегрально оказались хуже всех трех рассмотренных выше вариантов. Но лишь для сут. различие составило более 0.01. По-видимому, влияние tstart на качество оценок величины обусловлено двумя противодействующими факторами. С одной стороны, определение tstart для каждой серии по данным об афтершоках дает не очень устойчивые оценки этого параметра, которые не демонстрируют какой-либо зависимости от разности магнитуд MmMc (рис. 3а). С другой стороны, завышенные или заниженные оценки частично компенсируются фактическими оценками параметра c, которые в данном случае отражают не физическое явление, а неполноту каталога.
Особенно существенно нестабильность определений tstart обоими способами проявляется при небольшом количестве данных. При большом числе, наоборот, оценки, сохраняя большой разброс значений, уже проявляют зависимость от разности магнитуд MmMc (рис. 3б). При этом для значительной доли серий эти оценки существенно больше значений по соотношению работы [Helmstetter et al., 2006], что представляется естественным, так как это соотношение было получено для района Калифорнии, где высока плотность сейсмических станций. Для глобальной сейсмичности, очевидно, это соотношение должно быть уточнено. С этой целью мы провели огибающую облака точек на графике рис. 3б, которая описывается соотношением:
Таблица 1. Критерии качества результатов ретроспективной оценки магнитуды сильнейшего на интервале (t, 365 сут.) афтершока, усредненные по 9 прогнозным интервалам времени, в зависимости от способа определения величины tstart
Метод определения tstart
| Среднее число серий* с минимум 5 афтершоков в (tstart, t)
| Среднее значение по 9 интервалам (t, 365 сут.), t = 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64 сут.
| ||
LG
| PG0.5
| (LG + PG0.5)/2
| ||
По работе [Shebalin et al., 2017]
| 409
| 1.288
| 1.081
| 1.184
|
По фактической оценке параметра с
| 402
| 1.282
| 1.102
| 1.192
|
По работе [Helmstetter et al., 2006]
| 401
| 1.277
| 1.121
| 1.199
|
343
| 1.273
| 1.172
| 1.223
|
* Учитывается число афтершоков представительной магнитуды, которое может меняться для разных методов определения tstart.
Рис. 3. Зависимости определений величины tstart от разности магнитуд Mm–Mc, полученные по скользящему среднему (кружки) и по фактической оценке параметра c (крестики): (а) – для всех серий; (б) – для серий с числом афтершоков магнитуды Mc и выше не менее 100. Параметр Mc определялся методом MAXC на интервале [0.01, 4]. Жирная линия – огибающая точки на графике, определяется уравнением , штриховая линия соответствует соотношению из работы [Helmstetter et al., 2006]
| (12) |
Использование соотношения (12) существенно улучшает интегральную оценку результатов (см. табл. 1). Везде в дальнейшем мы будем использовать этот способ определения
РЕТРОСПЕКТИВНЫЕ РЕЗУЛЬТАТЫ
Мы провели ретроспективные тесты разработанной здесь методики на 777 сериях афтершоков от сильных землетрясений мира с Mm ≥ 6.5 для времен t в геометрической прогрессии от 0.25 сут., до 64 сут. c двукратным изменением на каждом шаге. Оценки проводятся, если число афтершоков представительной магнитуды на интервале составляет не менее N0. В противном случае применяется оценка по динамическому закону Бота. Поскольку при оценивании качества метода в качестве референц-модели выступает динамический закон Бота, такие случаи при оценивании исключаются. Величина N0 является единственным свободным параметром метода. После тестов с вариацией параметра N0 в пределах от 4 до 20 было выбрано оптимальное значение N0 = 5. Это значение мы затем использовали во всех тестах. Результаты тестирования основного метода приведены в табл. 2. Для всех рассмотренных времен t значения критериев LG и PG0.5 оказались больше 1 и в среднем составили 1.22. Таким образом, наша методика, за счет использования информации об уже произошедших афтершоках, в среднем дает существенно лучшие, зависящие от времени, оценки магнитуды последующего сильнейшего афтершока по сравнению с оценками по динамическому закону Бота. В среднем выигрыш по вероятности составляет около 22%.
Различия наиболее вероятного значения величины по полученным оценкам и фактически наблюденного значения варьируют в довольно широких пределах, лишь немного более узких по сравнению с оценками по динамическому закону Бота (рис. 4). Преимущество используемого здесь подхода по отношению к оценкам по динамическому закону Бота состоит в том, что ширина каждого отдельного распределения (6) оценки величины (ширина, например, может быть определена как разность квантилей уровней 0.9 и 0.1) существенно меньше ширины распределения (9) динамического закона Бота. В показанном на рис. 5 примере при близких значениях максимумов распределений ширина первого распределения значительно меньше, а значение плотности вероятности в точке наблюденного значения больше. Отношение значений плотности вероятности двух оценок составляет основу критерия LG, равного среднему геометрическому значению этого отношения по всем сериям афтершоков. Для интервала (0.25 сут., 365 сут.), например, этот критерий оказался равен 1.392, что обусловлено в среднем относительно большими значениями плотности первого распределения (рис. 6а). Диаграмма ошибок для этого интервала показана на рис. 6б. Существенное отклонение диаграммы от диагонали, соответствующей случайному прогнозу, так же обусловлено более узким распределением оценок по данной работе. Точки на диаграмме соответствуют разным значениям расстояния в единицах магнитуды до максимума (моды) распределений, при этом величина ν это доля серий, для которых фактическое значение M1 отстоит от моды распределения (6) на большее расстояние, а величина τ вероятность случайной реализации значения M1 на меньшем расстоянии от моды распределения (9). И хотя, как следует из рис. 4, расстояния фактического значения M1 до моды двух распределений в среднем различаются несущественно, вероятности случайного попадания в интервалы одинаковой ширины для двух распределений различны. Это различие можно характеризовать величиной PG0.5. Для 9 вариантов на рис. 4 значения PG0.5 показаны графически на рис. 7 (незаштрихованные квадраты). Величина PG0.5, напомним, равна тангенсу наклона прямой, соединяющей точки (0.1) и (τ*, ν* = 0.5).
Таблица 2. Результаты ретроспективного тестирования методики прогноза величины M1(t, T)
t, сут.
| N
| Основной метод (по точечным оценкам параметров)
| Метод с использованием Байесовского подходадля оценки M1(t, T)
| По работе [Scherbakov et al., 2017]
| |||
LG
| PG0.5
| LG
| PG0.5
| LG
| PG0.5
| ||
0.25
| 107
| 1.392
| 1.337
| 1.222
| 1.215
| 1.081
| 0.743
|
0.5
| 180
| 1.354
| 1.292
| 1.224
| 1.006
| 1.078
| 0.756
|
1
| 257
| 1.243
| 1.218
| 1.136
| 1.007
| 1.036
| 0.812
|
2
| 318
| 1.273
| 1.146
| 1.176
| 1.034
| 1.084
| 0.839
|
4
| 362
| 1.266
| 1.155
| 1.201
| 1.103
| 1.121
| 0.929
|
8
| 407
| 1.302
| 1.093
| 1.203
| 1.056
| 1.116
| 0.923
|
16
| 447
| 1.148
| 1.096
| 1.139
| 1.045
| 1.096
| 0.961
|
32
| 486
| 1.216
| 1.055
| 1.178
| 1.033
| 1.155
| 1.014
|
64
| 523
| 1.265
| 1.159
| 1.220
| 1.147
| 1.206
| 1.062
|
Среднее
| 343
| 1.273
| 1.172
| 1.189
| 1.072
| 1.108
| 0.893
|
| 1.223
| 1.130
| 1.000
|
N – число серий, для которых на интервале (tstart, t) число афтершоков представительной магнитуды составляет 5 или более.
Рис. 4. Распределения разности фактической магнитуды сильнейшего на интервале (t, 365 сут.) афтершока и максимума распределения (моды) оценки этой величины основным методом для 9 значений времени t. Распределения показаны гистограммами, нормированными на число рассмотренных серий афтершоков. Тонкой линией для сравнения показаны аналогичные распределения для оценок по динамическому закону Бота.
Рис. 5. Пример оценок величины M1(0.25 сут., 365 сут.) для афтершоков землетрясения магнитуды 7.7 (07.07.1990 г.). Жирная линия – плотность распределения по основному методу данной работы, соответствует функции распределения (6), штриховая линия – по динамическому закону Бота, соответствует функции распределения (9). Магнитуда сильнейшего афтершока показана вертикальной линией.
Рис. 6. Оценивание качества оценок величины M1(0.25 сут., 365 сут.) относительно динамического закона Бота. Соотношение плотностей распределений оценки по данным о предыдущих афтершоках (6) и по динамическрму закону Бота (9) при фактически наблюденных значениях величины (а); диаграмма ошибок (б).
ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ
Одновременно с разработкой нами данного подхода, над аналогичной задачей работала объединенная группа коллег из Канады и Японии, которые также использовали Байесовский подход для оценки максимальной магнитуды последующих афтершоков по данным о предыдущих. Результаты их работы опубликованы в конце 2017 г. [Scherbakov et al., 2017]. Наиболее существенное отличие этого подхода от нашего состоит в том, что в работе [Scherbakov et al., 2017] не учитывается неполнота каталога в начале афтершоковой серии. Кроме того, для регуляризации оценок мы используем информативное априорное распределение параметров законов ОмориУтсу и ГутенбергаРихтера, взяв для этого сглаженные эмпирические распределения параметров, полученные по большому числу афтершоковых серий. В работе [Scherbakov et al., 2017] для этих параметров используется равномерные на отрезке априорные распределения, при этом границы задаются из общих представлений.
Первоначально мы также использовали Байесовский подход как для оценки параметров, так и для оценки величины результаты тестирования этого варианта приведены в табл. 2. Сравнение этих результатов с результатами по основному методу показывают существенное преимущество основного варианта, в котором используются точечные оценки параметров. На том же ретроспективном материале мы провели сравнение наших оценок с оценками по методу работы [Scherbakov et al., 2017]. Эти результаты также приведены в табл. 2. Как оказалось, оценки по работе [Scherbakov et al., 2017] значительно хуже и в среднем не дают выигрыша по отношению к оценкам по динамическому закону Бота.
Разные способы определения параметра tstart мы подробно исследовали в отдельном разделе данной работы. При сопоставлении результатов ретроспективного тестирования с применением разных оценок tstart оказалось, что предпочтительными являются не индивидуальные оценки по имеющимся данным для каждой серии, а значения, определяемые разностью магнитуды основного толчка Mm и представительной магнитуды Mc. Такое соотношение, предложенное в работе [Helmstetter et al., 2006] для Калифорнии, как оказалось, дает заниженные значения величины tstart. Используя оценки этой величины двумя способами при большом количестве данных, мы получили аналогичное соотношение (12) для глобальной сейсмичности. Использование значений tstart по этому соотношению привело к улучшению получаемых оценок
Другой важный элемент нашего подхода регуляризация оценок параметров с помощью Байесовского подхода с использованием аппроксимаций нормальными распределениями полученных нами глобальных эмпирических распределений параметров lg (c), p, b в афтершоковых последовательностях. Чтобы оценить роль использования более детальной априорной информации о параметрах (не только их пределов, как часто принимается в Байесовских оценках, но и наиболее вероятных значений и формы распределения), мы сравнили результаты ретроспективных оценок в основном варианте с результатами аналогичных оценок, но с заменой априорных распределений на равномерные на заданных отрезках (табл. 3). Результаты с регуляризацией оценок параметров с нормальными априорными распределениями оказались существенно лучше.
Таким образом, учет неполноты каталога в начале афтершоковой серии и регуляризация оценок параметров являются важнейшими факторами для эффективных оценок величины
Первоначально в рамках данного исследования планировалось учитывать две особенности закона ГутенбергаРихтера для афтершоковых последовательностей. Во-первых, распределение магнитуды афтершоков может быть естественным образом ограничено в области больших значений из-за эффекта конечного объема [Romanovicz, 1992]. Но наши тесты показали, что учет этого эффекта влияет на результирующие оценки незначительно по сравнению с другими факторами, и им можно пренебречь.
Вторая возможная особенность это излом графика повторяемости, также приводящий к снижению вероятности сильнейших событий. Такой излом мы обнаружили для афтершоков от землетрясений магнитуды 7.5 и выше в Новой Зеландии [Shebalin, Baranov, 2017]. Этот излом может быть связан с известным эффектом пост-сейсмического проскальзывания в очаге землетрясения [Vorobieva et al., 2016]. Для этого случая также были получены варианты формул (6) и (6а), а также разработана программа оценки параметров графика повторяемости с изломом. Однако увеличение числа оцениваемых параметров и соответствующее снижение устойчивости не дали улучшить оценки величины и от учета возможного излома графика повторяемости пришлось отказаться.
В данном исследовании формально мы учитываем и еще одну возможную особенность поведения графика повторяемости в афтершоковых последовательностях, часто упоминаемую в литературе. Хорошо известно, что наклон графика повторяемости значительно уменьшается перед сильными землетрясениями [Molchan et al., 1999; Rodkin, Tikhonov, 2016]. Соответственно естественно ожидать, что в афтершоковом процессе происходит возврат к обычным значениям. Такое поведение с высокой степенью надежности наблюдается в лабораторных экспериментах [Sobolev et al., 1996; Smirnov et al., 2010; 2018]. В реальных афтершоковых процессах также формально наблюдается возрастание параметра b в течение первых суток после основного толчка [Rodkin, Tikhonov, 2016], однако этот эффект в значительной (если не определяющей) степени связан с неполнотой каталога, которая зависит как от магнитуды, так и от времени [Helmstetter et al., 2006]. В работе [Баранов, Шебалин, 2019] было показано, что магнитуды сильных афтершоков не зависят от времени, а значит, вариации наклона графика повторяемости в афтершоковых сериях могут наблюдаться только за счет зависимости от времени магнитуд относительно слабых афтершоков, что напрямую проверить не представляется возможным. Результаты данной работы, однако, косвенно свидетельствуют о постоянстве параметра b, начиная с ранних стадий афтершоков. В пользу этого говорит относительно хорошее качество прогнозов по данным об афтершоках, зарегистрированных в течение менее суток после основных толчков (табл. 2). Если бы параметр b значительно вырос после этого периода, то оценки величины должны бы были оказаться значительно завышенными, что не наблюдается.
Таблица 3. Критерии качества результатов ретроспективной оценки магнитуды сильнейшего на интервале (t, 365 сут.) афтершока, усредненные по 9 прогнозным интервалам времени, в зависимости от вида априорных распределений для Байесовских оценок параметров с, p, b
Вид априорных распределенийдля оценки параметров
| Среднее значение по 9 интервалам (t, 365 сут.),t = 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64 сут.
| ||
LG
| PG0.5
| (LG + PG0.5)/2
| |
Нормальная аппроксимация эмпирических распределений
| 1.273
| 1.172
| 1.223
|
Равномерные априорные распределения
| 1.229
| 0.976
| 1.103
|
В качестве модели убывания интенсивности афтершоков мы использовали степенной закон ОмориУтсу [Utsu, 1961]. Вместе с тем, известно, что иногда завершение серии афтершоков происходит быстрее [Narteau et al., 2002]. Это характерно для очагов с подвижкой по «незрелому» разлому, образовавшемуся сравнительно недавно или непосредственно при землетрясении [Narteau et al., 2003], что является редким явлением. На начальных стадиях данного исследования мы в качестве вариантов рассматривали вместо закона ОмориУтсу его суперпозицию с экспоненциальным спаданием интенсивности афтершоков, а также модель работы [Narteau et al., 2002]. Оказалось, что незначительное улучшение оценок для нескольких землетрясений с лихвой компенсируется общим ухудшением за счет неустойчивости оценок параметров, число которых увеличивается на 1.
Выигрыш по вероятности полученных ретроспективных оценок по отношению к динамическому закону Бота составляет не многим более 20%. Есть ли какая-либо возможность улучшить этот показатель? Для ответа на этот вопрос мы рассмотрели вариант «идеального прогноза», в котором оценки числа событий представительной магнитуды на прогнозных интервалах заменены на фактическое их число. Таким образом, оценка параметров c и p не используется. Параметр b оценивается так же, как и в основном варианте. Результаты ретроспективных «идеальных прогнозов» приведены на рис. 7. Соответствующие значения критериев LG и PG0.5 являются тем переделом сверху, которого в принципе могут достичь оценки величины по нашей методике.
Рис. 7. Сравнение критериев качества относительно динамического закона Бота (выигрыша по вероятности) для ретроспективных оценок по основному варианту и для «идеального прогноза», в котором оценка числа событий представительной магнитуды на прогнозных интервалах (t, 365 сут.) заменена на их фактическое число. Треугольниками показаны значения критерия LG, квадратами PG0.5, кружками – (LG + PG0.5)/2. Закрашены символы, соответствующие «идеальному прогнозу».
Для «идеального прогноза» критерии качества составляют около 1.8. Это тот предел, который может быть в принципе достигнут. На этом фоне выигрыш около 20% и тем более 40% уже не представляется небольшим. Обращает на себя внимание, что оценки качества слабо зависят от времени t начала прогнозного интервала. Это связано, по-видимому, с тем, что разброс оценок как по формуле (6), так и по формуле (9) не зависит от t [Баранов, Шебалин, 2019]. Неожиданным и пока необъяснимым является лучшее качество оценок по наименьшему количеству доступных данных (для t = 0.25 и 0.5 сут.).
Рис. 7 красноречиво свидетельствует об ограничениях статистического подхода для оценки магнитуды сильнейших последующих афтершоков. Принципиальное улучшение таких оценок возможно, по-видимому, только при использовании информации о физических свойствах очага основного землетрясения. Перспективным в этом направлении представляется анализ волновых форм землетрясения и спутниковых данных о ко-сейсмических и пост-сейсмических деформациях.
Рис. 8. Гистограммы отношения оценки числа Λ (4, 365) афтершоков представительной магнитуды на интервале (4 сут., 365 сут.) по формулам (3) и (4) к наблюденному числу N (4, 365) для рассмотренных 362 серий (см. табл. 2): (а) – для каждой серии одна оценка с использованием точечных оценок параметров (максимумы Байесовских апостериорных распределений); (б) – для каждой серии 100 оценок методом Монте-Карло по апостериорным распределениям параметров.
В процессе разработки методики на одном из последних этапов мы отказались от Байесовской оценки величины Тот факт, что использование точечных оценок параметров (но, по-прежнему, полученных с помощью Байесовского подхода с использованием информации о распределении параметров в афтершоковых последовательностях глобальной сейсмичности) дает лучшие ретроспективные результаты по сравнению с полной схемой Байесовских оценок величины является контринтуитивным и требует пояснения. Действительно, оценка числа событий на прогнозном интервале в полной схеме Байесовских оценок (рис. 8б) несколько лучше, чем при использовании точечных оценок параметров (рис. 8а), среднеквадратичный разброс почти вдвое меньше (0.21 против 0.39). Однако при использовании точечных оценок ширина распределения оценки сравнительно невелика, среднеквадратичный разброс, независимо от прогнозного интервала, составляет примерно 0.45. При использование полной Байесовской схемы усредняются не оценки ожидаемого числа событий на прогнозном интервале, а распределения оценки что приводит к «размазыванию» этого распределения и приближает его форму к форме динамического закона Бота, имеющего среднеквадратичный разброс примерно 0.66.
Таким образом, использование точечных оценок параметров (при условии их регуляризации описанным выше способом) действительно предпочтительно по сравнению с полной Байесовской схемой оценок величины Дополнительное небольшое улучшение оценок может быть получено использованием промежуточного варианта с использованием Байесовских оценок величины Λ(t, T) в формуле (6). Тесты показали, что такой подход в среднем улучшает выигрыш по вероятности ретроспективных оценок еще примерно на 1%.
ЗАКЛЮЧЕНИЕ
Рассмотрена задача прогноза магнитуды сильнейшего афтершока, ожидаемого через некоторое время после основного толчка. Афтершоковый процесс моделировался суперпозицией законов ГуттенбергаРихтера ОмориУтсу (подход РизенбергаДжонс), параметры которых оценивались методом Байеса с использованием данных об уже произошедших афтершоках и априорной информации о вероятных значениях параметров.
Качество прогноза оценивалось относительно референц-модели, определяющей распределение вероятности максимальной магнитуды последующих афтершоков, по двум независимым критериям. Первым критерием являлся информационный выигрыш модели равный отношению правдоподобия полученных реализаций для оцениваемого метода и референц-модели в пересчете на один прогноз. В качестве второго критерия использовался выигрыш по вероятности оцениваемого метода, рассчитанный по диаграмме ошибок, построенной относительно референц-модели, для доли пропусков цели 0.5. В качестве референц-модели был использован динамический закон Бота [Баранов, Шебалин, 2018].
Рассмотрено несколько вариантов оценки максимальной магнитуды последующих афтершоков. Значительное внимание уделено вариантам определения представительной магнитуды, зависящей от интервала времени после основного толчка, исключаемого из рассмотрения при оценке параметров. Учет неполноты каталога в начале серии афтершоков играет ключевую роль для точных оценок. Получено оптимальное соотношение, связывающее магнитуду основного толчка, период исключения данных и представительную магнитуду, уточняющее соотношение работы [Helmstetter et al., 2006].
Ретроспективная проверка вариантов осуществлялось по данным 777 серий афтершоков от землетрясений мира с магнитудой 6.5 и выше за период с 1980 по 2016 гг., выделенных из каталога ANSS ComCat Геологической службы США по алгоритму Молчана Г.М. и Дмитриевой О.Е. [1991].
В результате тестирования была выработана оптимальная методика оценивания магнитуды сильнейшего афтершока, ожидаемого через некоторое время после основного толчка. Исходными данными для оценивания является афтершоки, произошедшие за время t = 0.25, 0.5, 1, 2, 4, 8, 16, 32 или 64 суток после основного толчка. Методика позволяет оценить магнитуду сильнейшего афтершока, ожидаемого на интервале (t, 365) суток после основного толчка. В основе методики лежат точечные оценки параметров по методу Байеса с использованием априорной информации в виде усредненных распределений параметров, аппроксимируемых нормальными распределениями.
Разработанная методика дает выигрыш по вероятности от 20 до 40% по сравнению с оценками, полученными по динамическому закону Бота, при этом идеальный прогноз (параметры модели оцениваются на интервале (0, 365) после основного толчка) дает выигрыш 80%. На этом фоне выигрыш около 20% и тем более 40% уже не представляется небольшим. По результатам тестирования по глобальным данным методика также значительно превосходит метод работы [Scherbakov et al., 2017]. Главными отличиями нашего метода является учет неполноты каталога и использование усредненных распределений параметров в качестве априорной информации для их оценок по методу Байеса.
Разработанная методика может использоваться для прогноза максимальной магнитуды последующих повторных толчков после сильных землетрясений с использованием информации об уже произошедших афтершоках. Эти оценки включены в разрабатываемую автоматизированную систему оценки опасности сильных афтершоков AFCAST (www.afcast.org/afcast).
Финансирование работы
Исследование выполнено при частичной поддержке Российского фонда фундаментальных исследований (проекты 17-05-00749 и 19-05-00812.
About the authors
S. V. Baranov
Institute of Earthquake Prediction Theory and Mathematical Geophysics,Russian Academy of Sciences; Kola Branch, Federal Research Center “Geophysical Survey of Russian Academy of Sciences”
Author for correspondence.
Email: bars.vl@gmail.com
Russian Federation, Moscow; Apatity
V. A. Pavlenko
Institute of Earthquake Prediction Theory and Mathematical Geophysics,Russian Academy of Sciences; Schmidt Institute of Physics of the Earth, Russian Academy of Sciences
Email: pavlenko.vasily@gmail.com
Russian Federation, Moscow
P. N. Shebalin
Institute of Earthquake Prediction Theory and Mathematical Geophysics,Russian Academy of Sciences
Email: shebalin@mitp.ru
Russian Federation, Moscow
References
- Баранов С.В., Шебалин П.Н. О прогнозировании афтершоковой активности. 1. Адаптивные оценки на основе законов Омори и Гутенберга–Рихтера // Физика Земли. 2016. № 3. С. 82–101. doi: 10.7868/S0002333716020034
- Баранов С.В., Шебалин П.Н. О прогнозировании активности афтершоков. 2. Оценка области распространения сильных афтершоков // Физика Земли. 2017. № 3. C. 43–61. doi: 10.7868/S0002333717020028
- Баранов С.В., Шебалин П.Н. Глобальная статистика афтершоков сильных землетрясений: независимость времен и магнитуд // Вулканология и сейсмология. 2019. № 2. С. 67–76.
- Баранов С.В., Шебалин П.Н. О прогнозировании афтершоковой активности. 3. Динамический закон Бота // Физика Земли. 2018. № 6. С. 129–136.
- Молчан Г.М., Дмитриева О.Е. Идентификация афтершоков: обзор и новые подходы // Вычислительная сейсмология. 1991. Вып. 24. С. 19–50.
- Писаренко В.Ф., Родкин М.В., Рукавишникова Т.А. Оценка вероятности редких экстремальных событий для случая малых выборок, методика и примеры анализа каталога землетрясений // Физика Земли. 2017. № 6. С. 3–17.
- Смирнов В.Б., Пономарев А.В., Бернар П., Патонин А.В. Закономерности переходных режимов сейсмического процесса по данным лабораторного и натурного моделирования // Физика Земли. 2010. № 2. С. 17–49.
- Смирнов В.Б., Пономарев А.В., Станчиц С.А., Потанина М.Г., Патонин А.В., Dresen G., Narteau C., Bernard P., Строганова С.М. Лабораторное моделирование афтершоковых последовательностей: зависимость параметров Омори и Гутенберга–Рихтера от напряжений // Физика Земли. 2019. № 1. С. 17–49.
- Шебалин П.Н., Баранов С.В., Дзебоев Б.А. Закон повторяемости количества афтершоков // Докл. РАН. 2018. Т. 481. № 3. С. 320–323.
- Aki K. Maximum likelihood estimate of b in the formula log N a – bM and its confidence level // Bull. Earthquake Res. Inst. 1965. V. 43. P. 237–239.
- ANSS Comprehensive Earthquake Catalog (ComCat) URL: https: // earthquake.usgs.gov/data/comcat/
- Bath M. Lateral inhomogeneities in the upper mantle // Tectonophysics. 1965. V. 2. P. 483–514.
- Bender B. Maximum likelihood estimation of b values for magnitude grouped data // Bulletin of the Seismological Society of America. 1983. V. 73. № 3. P. 831–851.
- Gutenberg B., Richter C.F. Seismicity of the Earth and Associated Phenomena 2nd ed./Princeton N.J. Princeton University Press. 1954. 273 р.
- Hainzl S. Rate-Dependent Incompleteness of Earthquake Catalogs // Seismological Research Letters. 2016. V. 87. № 2 A. P. 337–344. Helmstetter A., Sornette D. Båth’s law derived from the Gutenberg–Richter law and from aftershock properties // Geophys. Res. Lett. 2003. V. 30. P. 2069.
- Holschneider, M., Zöller G., Clements R., Schorlemmer D. Can we test for the maximum possible earthquake magnitude? // J. Geophys. Res. Solid Earth. 2014. V. 119. P. 2019–2028. doi: 10.1002/2013 JB010319
- Holschneider M., Narteau C., Shebalin P., Peng Z., Schorlemmer D. Bayesian analysis of the modified Omori law // Journal of Geophysical Research. 2012. V. 117. B05317. doi: 10.1029/2011 JB009054
- Kijko A. Estimation of the Maximum Earthquake Magnitude, Mmax // Pure appl. Geophys. 2004. V. 161. P. 1–27.
- Knopoff L., Kagan Y. Analysis of the Extremes as Applied to Earthquake Problems // J. Geophys. Res. 1977. V. 82. P. 5647–5657.
- Molchan G. Structure of optimal strategies in earthquake prediction // Tectonophysics. 1991. V. 193. P. 267–276.
- Molchan G., Kronrod T., Nekrasova A. Immediate foreshocks: time variation of the b-value // Phys. Earth Planet. Int. 1999. V. 111. P. 129–140.
- Narteau C., Shebalin P., Holschneider M. Temporal limits of the power law aftershock decay rate // Journal of Geophys. Res. 2002. V. 107. P. 1201–1214.
- Narteau C., Shebalin P., Hainzl S., Zoller G., Holschneider M. Emergence of a band-limited power law in the aftershock decay rate of a slider-block model // Geophys. Res. Letters. 2003. V. 30. P. 1568. doi: 10.1029/2003 GL017110
- Pisarenko V.F., Lyubushin A.A., Lysenko V.B., Golubeva T.V. Statistical Estimation of Seismic Hazard Parameters: maximum possible magnitude and related parameters // Bull. Seism. Soc. Am. 1996. V. 86. P. 691–700.
- Rodkin M.V., Tikhonov I.N. The typical seismic behavior in the vicinity of a large earthquake // Physics and Che¬mistry of the Earth. 2016. V. 95. P. 73–84.
- Romanovicz B. Strike-slip earthquakes on quasi-vertical transcurrent faults: Inferences for general scaling relations // Geophys. Res. Lett. 1992. V. 19. Is. 5. P. 481–484. doi: 10.1029/92 GL00265
- Reasenberg P.A., Jones L.M. Earthquake Hazard After a Mainshock in California // Science. 1989. V. 242. № 4895. P. 1173–1176. doi: 10.1126/science.243.4895.1173
- Schorlemmer D., Gerstenberger M., Wiemer S., Jackson D.D., Rhoades D.A. Earthquake likelihood model testing // Seismol. Res. Lett. 2007. V. 78. P. 17–29.
- Shcherbakov R., Zhuang J., Ogata Y. Constraining the magnitude of the largest event in a foreshock-mainshock-aftershock sequence // Geophys. J. Int. 2018. V. 212. P. 1–13. doi: 10.1093/gji/ggx407
- Shebalin P., Narteau C. Depth dependent stress revealed by aftershocks // Nature Communications. 2017. V. 8. № 1317. doi: 10.1038/s41467-017-01446 y
- Shebalin P., Narteau C., Holschneider M., Zechar J. Combining earthquake forecast models using differential probability gains // Earth, Planets and Space. 2014. V. 66:37. P. 1–14.
- Shebalin P., Baranov S. Long-Delayed Aftershocks in New Zealand and the 2016 M7.8 Kaikoura Earthquake // Pure and applied Geophysics. 2017. V. 174 (10). P. 3751–3764. doi: 10.1007/s00024-017-1608-9
- Sobolev G.A., Ponomarev A.V., Koltsov A.V., Smirnov V.B. Simulation of triggered earthquakes in the laboratory // Pure and Applied Geophysics. 1996. V. 147. P. 345–355. doi: 10.1007/bf00877487
- Utsu T.A. Statistical study on the occurrence of aftershocks // Geoph. Magazine. 1961. V. 30. P. 521–605.
- Vere-Jones D. A note on the statistical interpretation of Båth’s law // Bull. Seismological Soc. Amer. 1969. V. 59. P. 1535–1541.
- Vere-Jones D. A limit theorem with application to Båth’s law in seismology // Adv. Appl. Prob. 2008. V. 40. P. 882–896.
- Vorobieva I., Narteau C., Shebalin P., Beauducel F., Nercessian F., Clouard V., Bouin M.P. Multiscale Mapping of Completeness Magnitude of Earthquake Catalogs // Bulletin of the Seismological Society of America. 2013. V. 103. P. 2188–2202.
- Vorobieva I., Shebalin P., Narteau C. Break of slope in earthquake size distribution and creep rate along the San Andreas fault system // Geophysical Research Letters. 2016. V. 43 (13). P. 6869–6875.
- Wiemer S., Wyss M. Minimum magnitude of completeness in earthquake catalogs: Examples from Alaska, the wes¬tern United States, and Japan // Bulletin of the Seismo¬logical Society of America. 2000. V. 90. P. 4859–4869. doi: 10.1785/0119990114
- Woessner J., Wiemer S. Assessing the quality of earthquake catalogues: Estimating the magnitude of completeness and its uncertainty // Bulletin of the Seismological Society of Ame¬rica. 2005. V. 95. P. 2684–2698. doi: 10.1785/0120040007