Forecasting aftershock activity: 4. Estimating the maximum magnitude of future aftershocks

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]:

$$F\left(M\right)=1-{10}^{-b\left(M-M_c\right)}, M\ge M_c,$$

(2)

где: t $–$ время от основного толчка, сут.; λ(t) $–$ число афтершоков c магнитудой $M\ge M_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]:

$$\Lambda \left(t,T\right)=\frac{D\left(t,T;c,p\right)}{D\left(t_{start},t;c,p\right)}\Lambda \left(t_{start},t\right),$$

(3)

где $\Lambda \left(t,T\right)$ и $\Lambda \left(t_{start},t\right)$ $–$ ожидаемое число афтершоков с $M\ge M_c$ на соответствующих интервалах:

$$D\left(t_1,t_2;c,p\right)=\left\{\begin{array}{c}\frac{\left(1-p\right)/c}{{\left(1+\frac{t_2}{c}\right)}^{1-p}-{\left(1+\frac{t_1}{c}\right)}^{1-p}}, p\ne 1\\ \frac{1/c}{ln\left(1+\frac{t_2}{c}\right)-ln\left(1+\frac{t_1}{c}\right)}, p=1.\end{array}\right.$$

(4)

Оценкой величины $\Lambda \left(t_{start},t\right)$ является фактическое зарегистрированное на интервале $\left(t_{start},t\right)$ количество афтершоков. Если параметры с и p известны, то выражения (3) и (4) позволяют рассчитать оценку ожидаемого на интервале $\left(t,T\right)$ количества афтершоков $\Lambda \left(t,T\right).$ В сделанных предположениях вероятность того, что магнитуда произвольного афтершока в серии меньше значения M, определена формулой (2); вероятность того, что k афтершоков имеют магнитуду меньше M, равна $F{\left(M\right)}^k.$ Предположив, что число афтершоков на интервале имеет распределение Пуассона, и просуммировав по формуле полной вероятности все варианты значения k, получим распределение вероятности магнитуды M1(t, T) [Vere-Jones, 2003]:

 

$\begin{array}{c} P\left(M_1\left(t,T\right)\le M\right)= G\left(M;t,T\right)=\\ ={\displaystyle \sum }_{k=0}^{\infty }\frac{{\left[\Lambda \left(t,T\right)\right]}^k}{k!}{e}^{-\Lambda \left(t,T\right)}{\left[F\left(M\right)\right]}^k=e^{-\Lambda \left(t,T\right)\left[1-F\left(M\right)\right]}\end{array}$ (5)

Объединив (5) и (2) имеем распределение в виде двойной экспоненты:

$$G\left(M;t,T\right)=e^{-\Lambda \left(t,T\right){10}^{-b\left(M-M_c\right)}}.$$

(6)

Искомое оценочное распределение величины $M_1\left(t,T\right)$ можно получить, подставив в (6) оценку $\Lambda \left(t,T\right)$ по формулам (3) и (4), используя точечные оценки параметров c, p и b и фактическое число афтершоков представительной магнитуды на интервале $\left(t_{start},t\right)$ в качестве оценки величины $\Lambda \left(t_{start},t\right).$

Обратная функция распределения (6) имеет явное выражение, которое удобно использовать для нахождения квантилей распределения:

$$Q\left(\alpha \right)=M_c-\frac{1}{\text{ln}10\cdot b}ln\left[-\frac{1}{\Lambda \left(t,T\right)}\text{ln}\left(\alpha \right)\right].$$

(6a)

С учетом сильно нелинейной зависимости (6) от параметров и опасаясь из-за этого значительного смещения оценок в случае больших фактических ошибок оценки параметров, первоначально мы усложнили процедуру. Вместо точечных оценок параметров c, p и b, мы стали использовать апостериорные Байесовские распределения оценок (c, p) по работе [Holschneider et al., 2012] и оценки b по работе [Bender, 1983], многократно генерируя методом Монте-Карло, в соответствии с этими распределениями, значения параметров и усредняя распределение (6). Фактически этот подход эквивалентен подходу работы [Scherbakov et al., 2017], в которой этапы получения Байесовских оценок параметров встроены в многократное интегрирование по формуле Байеса для получения апостериорного распределения оценки величины $M_1\left(t,T\right).$ В нашем подходе затратное по времени численное многократное интегрирование заменялось без потери точности расчетами методом Монте-Карло с генерацией случайных значений параметров в соответствии с их апостериорными распределениями. Преимущество нашего подхода состоит также в возможности контроля оценок параметров c, p и b. Также как в работе [Scherbakov et al., 2017] мы на этом этапе использовали неинформативные априорные распределения параметров, то есть равномерные на заданных отрезках распределения.

Тестирование метода с помощью разработанных для этого критериев (см. следующий раздел) показало, что получаемые оценки параметров неустойчивы и часто приводят к нереалистичным оценкам величины $M_1\left(t,T\right).$ Для регуляризации оценок мы стали использовать более подробную априорную информацию о параметрах, построив распределения оценок параметров по глобальным данным (см. ниже). Дальнейшее тестирование показало, что в слишком сложной схеме Байесовских оценок нет необходимости. Использование точечных оценок параметров (но с их регуляризацией с помощью Байесовского подхода с информативными априорными распределениями) оказалось предпочтительным.

ОЦЕНИВАНИЕ РЕЗУЛЬТАТОВ

Результаты прогнозов всегда целесообразно оценивать относительно простой модели, отражающей текущее на данный момент знание о прогнозируемом процессе, не вызывающей сомнения. Хорошо известен закон Бота [Bath, 1965], в соответствии которым средняя разность значений магнитуды основного толчка и сильнейшего афтершока составляет в среднем 1.2. В работе [Баранов, Шебалин, 2019] на основе глобальной статистики афтершоковых последовательностей исследована зависимость магнитуды сильнейшего афтершока от времени. Полученные результаты авторы охарактеризовали как динамический закон Бота. Показано, что с течением времени после основного толчка магнитуда сильнейших афтершоков в среднем уменьшается, но при этом форма распределения этой величины (сходная с распределением Гаусса) со временем не меняется и среднеквадратическое отклонение сохраняет значение примерно 0.65.

В работе [Шебалин и др., 2018] установлено, что количество k непосредственных афтершоков магнитуды выше порога $\Delta M,$ определяемого относительно магнитуды основного толчка Mm, при глобальном и региональном рассмотрении подчиняется экспоненциальному закону c плотностью $p\left(k\right)={\Lambda }_0{e}^{-k{/}_0}.$ Это дает объяснение форме распределения магнитуд сильнейших афтершоков, похожего на распределение Гаусса и теоретически определяемого соотношением:

$P\left(M_1-M_m<m\right)=\frac{1-e^{-\frac{1}{{\Lambda }_0}}}{1-e^{-\frac{1}{{\Lambda }_0}}\left[1-{10}^{-b\left(m-\Delta M\right)}\right]},$ (7а)

где ${\Lambda }_0$ $–$ параметр экспоненциального распределения при относительном пороге магнитуды $\Delta M.$ При $\Delta M$ = $–$2 для глобальной статистики афтершоков в течение года после землетрясений с $M_m\ge 6.5$ получена оценка ${\Lambda }_0=\mathrm{6.7.}$ Как оказалось, формула (7а) может быть упрощена. Комбинация формулы Вере-Джонса (6) и экспоненциального распределения дает простое и красивое соотношение, численно мало отличающееся от (7а):

$$P\left(M_1-M_m<m\right)=\frac{1}{1+{\Lambda }_0{10}^{-b\left(m-\Delta M\right)}}.$$

(7б)

Распределение (7б) справедливо и для произвольного интервала времени (t, T). В этом случае величина ${\Lambda }_0$ в (7б) должна быть заменена на зависимое от времени значение ${\Lambda }_0\left(t,T\right).$ Зависимость от времени этой величины определяется законом Омори$–$Утсу:

${\Lambda }_0\left(t,T\right)={\Lambda }_0\frac{{{\displaystyle \int }}_t^T{\left(t+c\right)}^{-p}dt}{{{\displaystyle \int }}_0^T{\left(t+c\right)}^{-p}dt}={\Lambda }_0\frac{D\left(t,T;c,p\right)}{D\left(\mathrm{0,}T;c,p\right)},$ (8)

где функция $D\left(t_1,t_2;c,p\right)$ определена соотношениями (4). Таким образом, величина M1(t, T) в динамическом законе Бота имеет распределение:

$$\begin{array}{c}P\left(M_1\left(t,T\right)-M_m<m\right)=\\ =\frac{1}{1+{10}^{-b\left\{m-\Delta M-\frac{1}{b}\text{lg}\left({\Lambda }_0\right)-\frac{1}{b}\text{lg}\left(\frac{D\left(\mathrm{0,}T;c,p\right)}{D\left(t,T;c,p\right)}\right)\right\}}}.\end{array}$$

(9)

Это распределение имеет моду, определяемую соотношением:

$$\begin{array}{c}Mode \left(M_1\left(t,T\right)\right)= \\ =\Delta M+\frac{1}{b}\lg \left({\Lambda }_0\right)+\frac{1}{b}\text{lg}\left(\frac{D\left(\mathrm{0,}T;c,p\right)}{D\left(t,T;c,p\right)}\right).\end{array}$$

(10)

В работе [Баранов, Шебалин, 2019] были получены следующие оценки для усредненной афтершоковой серии: b = 1.0, с = 0.04 сут., p = 1.016. Формула (9) с этими параметрами и значением ${\Lambda }_0=6.7$ дает близкое совпадение с эмпирическими распределениями, полученными в указанной работе; это распределение также мало отличается от аппроксимации эмпирических распределений нормальным. В данной работе мы используем динамический закон Бота, определяемый соотношением (9) в качестве референц-модели для оценивания прогнозов по модели, учитывающей данные об уже состоявшихся афтершоках. В дальнейшем мы будем использовать обозначение плотности этого распределения $p_{dynBath}\left(M,t\right).$

Аналогично работе [Баранов, Шебалин, 2016], мы используем два критерия для оценивания качества получаемых распределений величины $M_1\left(t,T\right).$ Один из критериев основан на оценках правдоподобия [Shorlemmer et al., 2007], другой $–$ на диаграмме ошибок [Молчан, Дмитриева, 1991].

Пусть N $–$ число завершившихся прогнозов, по которым оценивается метод, M1, i, i = 1,.., N $–$ реализации прогноза, то есть фактические значения величины $M_1\left(t,T\right),$ и $g\left(M;t,T\right)$ $–$ плотность полученного оценочного распределения этой величины. Определим информационный выигрыш LG(N) оцениваемого метода по отношению к референц-модели как отношение правдоподобия полученных реализаций для двух моделей в пересчете на один прогноз:

$$LG\left(N\right)={\left\{\frac{{{\displaystyle \prod }}_{i=1}^{N}g\left(M_{\mathrm{1,}i};t,T\right)}{{{\displaystyle \prod }}_{i=1}^N{p}_{dynBath}\left(M_{\mathrm{1,}i},t\right)}\right\}}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$N$}\right.}.$$

(11)

Величина $LG\left(N\right)$ тем больше, чем относительно чаще значения $M_{\mathrm{1,}i}\left(t,T\right)$ попадают в область больших значений плотности вероятности в тестируемой модели по сравнению с референц-моделью и чем реже в область малых значений плотности вероятности.

Диаграмма ошибок, которую мы здесь предлагаем, несколько отличается от стандартной. Поскольку прогноз можно считать тем лучше, чем ближе к максимуму плотности (моде) прогнозного распределения оказывается реализация, то в качестве управляющего параметра мы рассматриваем здесь расстояние δM от моды. Поэтому, чтобы сравнить тестируемую модель с моделью динамического закона Бота, мы определяем величину τ (доля тревог) как вероятность попадания величины $M_1\left(t,T\right)$ в интервал $\left[E_{dynBath}\text{}\left(t\right)-\delta M\text{},E_{dynBath}\text{}\left(t\right)+\delta M\right],$ а величину ν (доля пропусков цели) $–$ число случаев, когда $M_{i,1}\left(t,T\right)<E\left(t\right)-\delta M$ или $M_{i,1}\left(t,T\right)>E\left(t\right)+\delta M,$ где E(t) $–$ мода распределения для оцениваемой модели, а мода распределения для динамического закона Бота $E_{dynBath}\left(t\right)$ определена соотношением (10). Точки на диаграмме при разных значениях $\delta M$ образуют «траекторию ошибок». Чем в среднем относительно меньше расстояние от реализаций до моды, тем дальше от диагонали будет находиться «траектория» и тем лучше тестируемый метод. Тангенс наклона прямой на диаграмме ошибок, проходящей от точки (0.1) через точку траектории ошибок, называется вероятностным выигрышем (probability gain) [Molchan, 1991; Shebalin et al., 2014]. Для удобства сравнения для разных t мы будем рассматривать эту величину при ν = 0.5, обозначив ее PG0.5. Для обоих критериев значение 1 означает, что качество прогноза величины $M_1\left(t,T\right)$ анализируемым методом такое же, как и оценка по динамическому закону Бота, без привлечения информации об уже состоявшихся афтершоках. Значения больше 1 означают, что метод предпочтителен по отношению к референц-модели; значении меньше 1 $–$ предпочтителен, наоборот, динамический закон Бота.

Оценки, основанные на правдоподобии и на диаграмме ошибок дополняют друг друга. В работе [Shebalin et al., 2014] продемонстрировано, что в диаграмме ошибок относительно больший вес придается реализациям прогноза (успешные прогнозы и ложные тревоги) в областях пространства$–$времени, в которых тестируемая модель дает высокие вероятности прогнозируемых событий. Критерии, основанные на правдоподобии, наоборот, относительно больший вес придают областям низкой вероятности событий (включая пропуски цели). Для критериев, основанных на правдоподобии, кроме того часто возникает проблема нулевых или очень низких вероятностей. В случае пропуска цели ведь не так существенно насколько низко оценивалась вероятность события, важно, что она была много ниже значений, при которых события обычно реализуются. А «штраф» за один пропуск цели часто может компенсировать в критерии большое число очень успешных реализаций прогноза. Эта проблема легко компенсируется введением некоторого минимума вероятности события в модели: более низкие значения заменяются этим пороговым значением. Поскольку при этом возрастает суммарная вероятность, приходится производить перенормировку. Для критерия LG мы используем порог для плотности вероятности $g\left(M_{\mathrm{1,}i};t,T\right),$ равный 0.001.

ОЦЕНКА ПАРАМЕТРОВ ДЛЯ СЕРИИ АФТЕРШОКОВ

В тех случаях, когда используемых данных немного, получаемые оценки параметра неустойчивы, что, с учетом нелинейности, может существенно сказываться на оценках величины $M_1\left(t,T\right)$ даже при использовании Байесовского подхода с неинформативным (однородным) распределением параметра. Как показали ретроспективные тесты, из-за особенностей афтершоковых серий достаточно часто максимум апостериорного распределения параметра оказывается у границы задаваемых априорных пределов, что в конечном счете и приводит к ошибочным оценкам величины $M_1\left(t,T\right).$ Для регуляризации оценок мы решили использовать более детальную априорную информацию о параметрах. Были построены эмпирические функции распределения параметров по большому числу афтершоковых серий. Из 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.

 

На этапе оценки величины $M_1\left(t,T\right)$ полученные аппроксимации нормальными распределениями используются в качестве информативных априорных распределений параметров b, c, p для оценки этих параметров на интервале (tstart, t) для каждой серии афтершоков.

ОПТИМИЗАЦИЯ ПАРАМЕТРА tstart

Важное значение на качество оценок величины $M_1\left(t,T\right)$ оказывает выбор параметра tstart. Очевидно, что сразу после сильного землетрясения часть афтершоков не попадает в каталог, главным образом, из-за наложения их сейсмограмм на сейсмограмму основного толчка и взаимного наложения сейсмограмм афтершоков. Авторы работы [Helmstetter et al., 2006] установили для каталога землетрясений Калифонии приблизительную связь магнитуды основного толчка Mm, порога представительной магнитуды Mc и времени tstart, начиная с которого каталог с M ≥ Mc полон: Mc = Mm$–$4.5$–$0.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].

Мы сравнили два подхода путем интегрального сопоставления результатов ретроспективной оценки величины $M_1\left(t,T\right)$ по данным на интервале [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. Все оценки интегрально оказались хуже всех трех рассмотренных выше вариантов. Но лишь для $t_{start}\ge 0.2$ сут. различие составило более 0.01. По-видимому, влияние tstart на качество оценок величины $M_1\left(t,T\right)$ обусловлено двумя противодействующими факторами. С одной стороны, определение tstart для каждой серии по данным об афтершоках дает не очень устойчивые оценки этого параметра, которые не демонстрируют какой-либо зависимости от разности магнитуд Mm$–$Mc (рис. 3а). С другой стороны, завышенные или заниженные оценки частично компенсируются фактическими оценками параметра c, которые в данном случае отражают не физическое явление, а неполноту каталога.

Особенно существенно нестабильность определений tstart обоими способами проявляется при небольшом количестве данных. При большом числе, наоборот, оценки, сохраняя большой разброс значений, уже проявляют зависимость от разности магнитуд Mm$–$Mc (рис. 3б). При этом для значительной доли серий эти оценки существенно больше значений по соотношению работы [Helmstetter et al., 2006], что представляется естественным, так как это соотношение было получено для района Калифорнии, где высока плотность сейсмических станций. Для глобальной сейсмичности, очевидно, это соотношение должно быть уточнено. С этой целью мы провели огибающую облака точек на графике рис. 3б, которая описывается соотношением:

 

Таблица 1. Критерии качества результатов ретроспективной оценки магнитуды сильнейшего на интервале (t, 365 сут.) афтершока, усредненные по 9 прогнозным интервалам времени, в зависимости от способа определения величины t_start

Метод определения 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
$\lg \left(t_{start}\right)=\frac{1}{0.7}\left(M_m-M_c-3.5\right)$	343	1.273	1.172	1.223

* Учитывается число афтершоков представительной магнитуды, которое может меняться для разных методов опреде­ления t_start.

 

Рис. 3. Зависимости определений величины tstart от разности магнитуд M_m–M_c, полученные по скользящему среднему (кружки) и по фактической оценке параметра c (крестики): (а) – для всех серий; (б) – для серий с числом афтершоков магнитуды M_c и выше не менее 100. Параметр M_c определялся методом MAXC на интервале [0.01, 4]. Жирная линия – огибающая точки на графике, определяется уравнением $\mathit{l}\mathit{g}\left({\mathbf{t}}_{\mathbf{s}\mathbf{t}\mathbf{a}\mathbf{r}\mathbf{t}}\right)\mathbf{=}\frac{\mathbf{1}}{\mathbf{0}\mathbf{.7}}\left({\mathbf{M}}_{\mathbf{m}}\mathbf{-}{\mathbf{M}}_{\mathbf{c}}\mathbf{-}\mathbf{3}\mathbf{.5}\right)$, штриховая линия соответствует соотношению из работы [Helmstetter et al., 2006] $\mathit{l}\mathit{g}\left({\mathbf{t}}_{\mathbf{s}\mathbf{t}\mathbf{a}\mathbf{r}\mathbf{t}}\right)\mathbf{=}\frac{\mathbf{1}}{\mathbf{0}\mathbf{.76}}\left({\mathbf{M}}_{\mathbf{m}}\mathbf{-}{\mathbf{M}}_{\mathbf{c}}\mathbf{-}\mathbf{4}\mathbf{.5}\right)\mathbf{.}$

 

$$t_{start}={10}^{1.4\left(M_m-M_c-3.5\right)}.$$

(12)

Использование соотношения (12) существенно улучшает интегральную оценку результатов (см. табл. 1). Везде в дальнейшем мы будем использовать этот способ определения $t_{start}.$

РЕТРОСПЕКТИВНЫЕ РЕЗУЛЬТАТЫ

Мы провели ретроспективные тесты разработанной здесь методики на 777 сериях афтершоков от сильных землетрясений мира с Mm ≥ 6.5 для времен t в геометрической прогрессии от 0.25 сут., до 64 сут. c двукратным изменением на каждом шаге. Оценки проводятся, если число афтершоков представительной магнитуды на интервале $(t_{start}, t)$ составляет не менее N0. В противном случае применяется оценка по динамическому закону Бота. Поскольку при оценивании качества метода в качестве референц-модели выступает динамический закон Бота, такие случаи при оценивании исключаются. Величина N0 является единственным свободным параметром метода. После тестов с вариацией параметра N0 в пределах от 4 до 20 было выбрано оптимальное значение N0 = 5. Это значение мы затем использовали во всех тестах. Результаты тестирования основного метода приведены в табл. 2. Для всех рассмотренных времен t значения критериев LG и PG0.5 оказались больше 1 и в среднем составили 1.22. Таким образом, наша методика, за счет использования информации об уже произошедших афтершоках, в среднем дает существенно лучшие, зависящие от времени, оценки магнитуды последующего сильнейшего афтершока по сравнению с оценками по динамическому закону Бота. В среднем выигрыш по вероятности составляет около 22%.

Различия наиболее вероятного значения величины $M_1\left(t,T\right)$ по полученным оценкам и фактически наблюденного значения варьируют в довольно широких пределах, лишь немного более узких по сравнению с оценками по динамическому закону Бота (рис. 4). Преимущество используемого здесь подхода по отношению к оценкам по динамическому закону Бота состоит в том, что ширина каждого отдельного распределения (6) оценки величины $M_1\left(t,T\right)$ (ширина, например, может быть определена как разность квантилей уровней 0.9 и 0.1) существенно меньше ширины распределения (9) динамического закона Бота. В показанном на рис. 5 примере при близких значениях максимумов распределений ширина первого распределения значительно меньше, а значение плотности вероятности в точке наблюденного значения $M_1$ больше. Отношение значений плотности вероятности двух оценок составляет основу критерия LG, равного среднему геометрическому значению этого отношения по всем сериям афтершоков. Для интервала (0.25 сут., 365 сут.), например, этот критерий оказался равен 1.392, что обусловлено в среднем относительно большими значениями плотности первого распределения (рис. 6а). Диаграмма ошибок для этого интервала показана на рис. 6б. Существенное отклонение диаграммы от диагонали, соответствующей случайному прогнозу, так же обусловлено более узким распределением оценок $M_1\left(t,T\right)$ по данной работе. Точки на диаграмме соответствуют разным значениям расстояния в единицах магнитуды до максимума (моды) распределений, при этом величина ν $–$ это доля серий, для которых фактическое значение 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 – число серий, для которых на интервале (t_start, t) число афтершоков представительной магнитуды составляет 5 или более.

 

Рис. 4. Распределения разности фактической магнитуды сильнейшего на интервале (t, 365 сут.) афтершока и максимума распределения (моды) оценки этой величины основным методом для 9 значений времени t. Распределения показаны гистограммами, нормированными на число рассмотренных серий афтершоков. Тонкой линией для сравнения показаны аналогичные распределения для оценок по динамическому закону Бота.

 

Рис. 5. Пример оценок величины M₁(0.25 сут., 365 сут.) для афтершоков землетрясения магнитуды 7.7 (07.07.1990 г.). Жирная линия – плотность распределения по основному методу данной работы, соответствует функции распределения (6), штриховая линия – по динамическому закону Бота, соответствует функции распределения (9). Магнитуда сильнейшего афтершока показана вертикальной линией.

 

Рис. 6. Оценивание качества оценок величины M₁(0.25 сут., 365 сут.) относительно динамического закона Бота. Соотношение плотностей распределений оценки по данным о предыдущих афтершоках (6) и по динамическрму закону Бота (9) при фактически наблюденных значениях величины (а); диаграмма ошибок (б).

 

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

Одновременно с разработкой нами данного подхода, над аналогичной задачей работала объединенная группа коллег из Канады и Японии, которые также использовали Байесовский подход для оценки максимальной магнитуды последующих афтершоков по данным о предыдущих. Результаты их работы опубликованы в конце 2017 г. [Scherbakov et al., 2017]. Наиболее существенное отличие этого подхода от нашего состоит в том, что в работе [Scherbakov et al., 2017] не учитывается неполнота каталога в начале афтершоковой серии. Кроме того, для регуляризации оценок мы используем информативное априорное распределение параметров законов Омори$–$Утсу и Гутенберга$–$Рихтера, взяв для этого сглаженные эмпирические распределения параметров, полученные по большому числу афтершоковых серий. В работе [Scherbakov et al., 2017] для этих параметров используется равномерные на отрезке априорные распределения, при этом границы задаются из общих представлений.

Первоначально мы также использовали Байесовский подход как для оценки параметров, так и для оценки величины $M_1\left(t,T\right),$ результаты тестирования этого варианта приведены в табл. 2. Сравнение этих результатов с результатами по основному методу показывают существенное преимущество основного варианта, в котором используются точечные оценки параметров. На том же ретроспективном материале мы провели сравнение наших оценок с оценками по методу работы [Scherbakov et al., 2017]. Эти результаты также приведены в табл. 2. Как оказалось, оценки по работе [Scherbakov et al., 2017] значительно хуже и в среднем не дают выигрыша по отношению к оценкам по динамическому закону Бота.

Разные способы определения параметра tstart мы подробно исследовали в отдельном разделе данной работы. При сопоставлении результатов ретроспективного тестирования с применением разных оценок tstart оказалось, что предпочтительными являются не индивидуальные оценки по имеющимся данным для каждой серии, а значения, определяемые разностью магнитуды основного толчка Mm и представительной магнитуды Mc. Такое соотношение, предложенное в работе [Helmstetter et al., 2006] для Калифорнии, как оказалось, дает заниженные значения величины tstart. Используя оценки этой величины двумя способами при большом количестве данных, мы получили аналогичное соотношение (12) для глобальной сейсмичности. Использование значений tstart по этому соотношению привело к улучшению получаемых оценок $M_1\left(t,T\right).$

Другой важный элемент нашего подхода $–$ регуляризация оценок параметров с помощью Байесовского подхода с использованием аппроксимаций нормальными распределениями полученных нами глобальных эмпирических распределений параметров lg (c), p, b в афтершоковых последовательностях. Чтобы оценить роль использования более детальной априорной информации о параметрах (не только их пределов, как часто принимается в Байесовских оценках, но и наиболее вероятных значений и формы распределения), мы сравнили результаты ретроспективных оценок в основном варианте с результатами аналогичных оценок, но с заменой априорных распределений на равномерные на заданных отрезках (табл. 3). Результаты с регуляризацией оценок параметров с нормальными априорными распределениями оказались существенно лучше.

Таким образом, учет неполноты каталога в начале афтершоковой серии и регуляризация оценок параметров являются важнейшими факторами для эффективных оценок величины $M_1\left(t,T\right).$

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

Вторая возможная особенность $–$ это излом графика повторяемости, также приводящий к снижению вероятности сильнейших событий. Такой излом мы обнаружили для афтершоков от землетрясений магнитуды 7.5 и выше в Новой Зеландии [Shebalin, Baranov, 2017]. Этот излом может быть связан с известным эффектом пост-сейсмического проскальзывания в очаге землетрясения [Vorobieva et al., 2016]. Для этого случая также были получены варианты формул (6) и (6а), а также разработана программа оценки параметров графика повторяемости с изломом. Однако увеличение числа оцениваемых параметров и соответствующее снижение устойчивости не дали улучшить оценки величины $M_1\left(t,T\right),$ и от учета возможного излома графика повторяемости пришлось отказаться.

В данном исследовании формально мы учитываем и еще одну возможную особенность поведения графика повторяемости в афтершоковых последовательностях, часто упоминаемую в литературе. Хорошо известно, что наклон графика повторяемости значительно уменьшается перед сильными землетрясениями [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 значительно вырос после этого периода, то оценки величины $M_1\left(t,T\right)$ должны бы были оказаться значительно завышенными, что не наблюдается.

 

Таблица 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 являются тем переделом сверху, которого в принципе могут достичь оценки величины $M_1\left(t,T\right)$ по нашей методике.

 

Рис. 7. Сравнение критериев качества относительно динамического закона Бота (выигрыша по вероятности) для ретроспективных оценок по основному варианту и для «идеального прогноза», в котором оценка числа событий представительной магнитуды на прогнозных интервалах (t, 365 сут.) заменена на их фактическое число. Треугольниками показаны значения критерия LG, квадратами PG_0.5, кружками – (LG + PG_0.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 оценок методом Монте-Карло по апостериорным распределениям параметров.

 

В процессе разработки методики на одном из последних этапов мы отказались от Байесовской оценки величины $M_1\left(t,T\right).$ Тот факт, что использование точечных оценок параметров (но, по-прежнему, полученных с помощью Байесовского подхода с использованием информации о распределении параметров в афтершоковых последовательностях глобальной сейсмичности) дает лучшие ретроспективные результаты по сравнению с полной схемой Байесовских оценок величины $M_1\left(t,T\right),$ является контринтуитивным и требует пояснения. Действительно, оценка числа событий на прогнозном интервале в полной схеме Байесовских оценок (рис. 8б) несколько лучше, чем при использовании точечных оценок параметров (рис. 8а), среднеквадратичный разброс почти вдвое меньше (0.21 против 0.39). Однако при использовании точечных оценок ширина распределения оценки $M_1\left(t,T\right)$ сравнительно невелика, среднеквадратичный разброс, независимо от прогнозного интервала, составляет примерно 0.45. При использование полной Байесовской схемы усредняются не оценки ожидаемого числа событий на прогнозном интервале, а распределения оценки $M_1\left(t,T\right),$ что приводит к «размазыванию» этого распределения и приближает его форму к форме динамического закона Бота, имеющего среднеквадратичный разброс примерно 0.66.

Таким образом, использование точечных оценок параметров (при условии их регуляризации описанным выше способом) действительно предпочтительно по сравнению с полной Байесовской схемой оценок величины $M_1\left(t,T\right).$ Дополнительное небольшое улучшение оценок может быть получено использованием промежуточного варианта с использованием Байесовских оценок величины Λ(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

Supplementary files