Физика Земли, 2021, № 2, стр. 27-45

ВВЕДЕНИЕ

Параметр М_мах – максимальная региональная магнитуда – очень часто используется в различных задачах, связанных с оценкой сейсмического риска (см., например, [Kijko, Sellevoll, 1989; 1992; Писаренко, 1991; Ward, 1997; Kijko, Graham, 1998; Kagan, Schoenberg, 2001; Kijko, 2004; 2012; Aban et al., 2006; Kijko, Singh, 2011; Lasocki, Urban, 2011; Holschneider et al., 2011; Z¨oller, Holschneider, 2016; Bommer, van Elk, 2017; Vermeulen, Kijko, 2017; Beirland et al., 2019] и цитированную там литературу). Однако несмотря на его широкое распространение, этот параметр не имеет общепринятого формального определения. Мы будем рассматривать стационарные сейсмические режимы, хотя в реальной обстановке та или иная нестационарность всегда присутствует и ее необходимо, по возможности, учитывать. Если считать, что М_мах – это верхняя грань тех значений магнитуд, которые могут произойти в данном регионе в будущем, то сразу возникает вопрос о продолжительности этого будущего интервала времени. Часто под этим понимают неограниченный интервал времени, а величина М_мах иногда оценивается при этом, исходя также и из сейсмотектонических и геологических соображений, полагаемых постоянными на масштабе времени десятков тысяч лет и более. В большинстве случаев такая оценка практически непроверяема.

В данной статье проблема оценки М_мах трактуется как чисто статистическая; сейсмотектонические и физические факторы, определяющие величину М_мах, не рассматриваются. Мы будем трактовать М_мах как максимальную магнитуду события, которое произойдет в заданной области в будущий интервал времени Т. При такой трактовке М_мах является случайной величиной, характеристики которой могут быть статистически оценены с определенной надежностью. В целях практической оценки сейсмического риска можно выделить три разных временных масштаба. Масштаб 1–10 лет представляет интерес для страхования, масштаб 10–50 лет – для строительства и эксплуатации обычных сооружений, масштаб 50–сотни лет – для сооружений особого назначения (АЭС, иные особо ответственные предприятия, хранилища радиоактивных отходов и т.п.). Еще большие интервалы времени могут представлять интерес в исследованиях по геологии и тектонике земной коры. Естественно, значения М_мах для разных масштабов времени для одного и того же региона могут сильно отличаться (на единицу магнитуды и более). Кроме того, они могут отличаться еще и тем, насколько часто они могут нарушаться в реальной обстановке; при этом абсолютная верхняя грань, возможно, достигается на столь больших интервалах времени, что такая оценка теряет практический смысл. Все эти соображения неизбежно приводят к тому, что М_мах целесообразно рассматривать не как региональный (неизвестный) параметр, а как случайную величину М_мах(Т), зависящую от длины будущего интервала времени Т. Такой подход был предложен нами в работе [Pisarenko et al., 2010] и развит в ряде дальнейших публикаций [Pisarenko, Rodkin, 2010; 2013; 2015; Pisarenko et al., 2014]. Эта случайная величина имеет корректное вероятностное определение и функцию распределения, которую мы будем обозначать F_T(х), где x – магнитуда землетрясения (индекс Т мы часто будем опускать). Иногда удобнее пользоваться не самой F_T(х), а обратной по отношению к ней функцией, которая определяется из соотношения:

Если F_T(x) непрерывна и монотонна (что мы будем предполагать), то обратная функция Q_Т(q) определяется однозначно. Она представляет собой квантиль уровня q распределения F_T(x). Квантиль уровня q = 0.5 является медианой распределения. Для значений q близких к 1 квантилю Q_Т(q) может служить верхней границей доверительного интервала для магнитуд с уровнем доверия q. Задавая несколько значений уровня доверия, скажем, q = 0.9; 0.95; 0.975; 0.99; 0.999, и вычисляя оценки соответствующих квантилей Q_Т(q), можно получить достаточно полную характеристику поведения максимальных за время Т магнитуд. Прежде всего, укажем связь между квантилем максимальной за время Т магнитуды Q_Т(q) и квантилем одного случайного события K(q). Обозначим функцию распределения одного события Ф(x):

где m – магнитуда случайного события. Будем предполагать, что поток землетрясений является стационарным Пуассоновским потоком с интенсивностью λ. Тогда, как известно (см., например, [Писаренко, Родкин, 2007; Pisarenko et al., 2008; Pisarenko, Rodkin, 2010; 2014]):

Мы будем рассматривать ситуации, когда λТ  $ \gg $ 1 (это условие почти всегда выполняется в практических задачах). В этом случае можно с большой точностью считать, что:

Обозначим квантиль уровня q₀ распределения Ф(х) через K(q₀):

Квантиль Q_Т(q) определяется из уравнения

Перепишем (6) в виде

Из уравнения (7) видно, что уровень значимости q и среднее число событий на интервале Т, равное λТ, могут компенсировать друг друга. Увеличение λТ в w раз (w > 0) можно в точности скомпенсировать, взяв уровень значимости q^w (при этом надо следить за тем, чтобы правая часть (7) оставалась положительной, что может нарушаться из-за сделанного выше приближения (4)).

Согласно определению решение уравнения (7) является квантилем $K\left( {1 - \frac{{\lg \left( {\frac{1}{q}} \right)}}{{\lambda T}}} \right)~$ уровня 1 – $\frac{{{\text{lg}}\left( {\frac{1}{q}} \right)}}{{\lambda T}}$, т.е.

Равенство (8) устанавливает соотношение между квантилями распределения одиночной магнитуды и квантилями максимальной магнитуды за время Т для Пуассоновского потока событий. Для закона Гутенберга–Рихтера равенство (8) дает:

Для высоких уровней доверия (q близко к единице) можно приближенно считать, что lg(1/q) $ \approx $1 – q, и мы получаем вместо (8а) приближение:

Мы видим, что квантиль максимального события растет логарифмически с ростом Т и с ростом показателя 1/(1 – q). Если одиночная магнитуда имеет обобщенное распределение Парето (ОРП)

то из (8) получаем:

При малых по модулю отрицательных ξ получаем приближение, совпадающее с (8b):

Это приближение справедливо при следующем ограничении на значения λT:

Для больших значений λT нужно пользоваться точной формулой (8с).

Соотношения (8)–(8d) полезно иметь в виду при оценках сейсмического риска. Соотношение (8) представляет и самостоятельный интерес. Оно позволяет устанавливать своего рода эквивалентность между длительностью наблюдений произвольного Пуассоновского потока с заданной интенсивностью и уровнем значимости (надежностью) оцениваемой величины при анализе максимальных событий на будущем интервале времени.

Мы будем исследовать статистические оценки квантилей, их смещение, стандартное отклонение и среднеквадратичное отклонение. Обычно эти характеристики оцениваются по реальной выборке (каталогу землетрясений) тем или иным способом. Такой расчет обычно основан на некотором распределении, которое лишь приблизительно соответствует реальной выборке. В нашем подходе, по-видимому, впервые оценка смещения, стандартного отклонения и среднеквадратичного отклонения проводится на большом количестве (1000 штук) независимых выборок (искусственных каталогов) с известным законом распределения. Это дало возможность сравнить такую “истинную” точность оценок с той, которая получается по одному каталогу с неизвестным законом распределения. Сразу отметим, что полученные модельные характеристики разброса относительно известного ответа оказались несколько больше величин разброса, оцениваемых по одному (реальному) каталогу.

В качестве точного закона распределения (для генерирования на его основе 1000 искусственных каталогов) было использовано составное распределение, обозначаемое далее M2, состоящее из двух ветвей: распределение Гутенберга–Рихтера (ГР) (диапазон малых и средних магнитуд) и обобщенное распределение Парето (ОРП) (диапазон больших магнитуд), см. [Писаренко и др., 2020]. Обоснование такого выбора модели закона распределения дано ниже. В качестве прообразов для точных законов распределения использовались каталоги 6 разных регионов, выбранных из мирового каталога ISM-GEM (1904–2014 гг.). Как известно, этот каталог обеспечивает наиболее длинный числовой ряд в единых магнитудах M_w [Di Giacomo et al., 2018]. Были выбраны 6 регионов: Атлантика, Япония, Курильские острова, Новые Гебриды, Перу, Филиппины (см. рис. 1). Обоснование выбора этих регионов дано ниже. Для каждого из региональных каталогов подбирался теоретический закон распределения М2 (см. ниже), по которому генерировалась 1000 каталогов, что позволяло затем надежно оценить систематическое смещение и случайный разброс для каждого из регионов. При подборе закона учитывалось расстояние Колмогорова между функциями распределения и общие тектонофизические соображения, приводимые ниже.

Для оценки квантилей применялись 2 различных способа: метод максимального правдоподобия (ММП) и метод Байеса (МБ), введенный авторами в работах [Pisarenko et al., 1996; Pisarenko, Lyubushin, 1997; 1999] и модернизированный одним из авторов в последующих работах [Lyubushin et al., 2002; Lyubushin, Parvez, 2010]. В методе ММП в качестве параметрической модели использовалась модель М2, в Байесовском методе – усеченное распределение Гутенберга–Рихтера (УГР).

Помимо тестирования оценок квантилей на искусственных каталогах с известным ответом, мы провели также оценивание квантилей на их реальных “прототипах”.

ОЦЕНКА КВАНТИЛЕЙ НА ИСКУССТВЕННЫХ КАТАЛОГАХ

Искусственные каталоги генерировались следующим способом. Известно, что если любую случайную величину η c непрерывной функцией распределения F(x) = Pr{η ≤ x} подставить в эту функцию в качестве аргумента, то полученная случайная величина

будет иметь равномерное распределение на отрезке [0; 1]. Уравнение (9) можно разрешить относительно η, применив к обеим частям (9) функцию V, обратную по отношению к функции F:

Равенство (10) утверждает, что если в функцию V подставить в качестве аргумента случайную величину ζ с равномерным на [0; 1] распределением, то получится случайная величина с распределением F. Таким способом можно генерировать случайные величины с любым заданным законом, используя компьютерные псевдослучайные числа с равномерным на отрезке [0; 1] распределением.

В качестве модели теоретического (истинного) закона распределения нами использовалось составное распределение М2, которое задается следующим образом (cм. подробнее [Писаренко и др., 2020]). Функция распределения магнитуд F(x) состоит из двух частей:

Модель (11) содержит следующие параметры. Порог h разделяет 2 ветви модели, ГР и ОРП. Нижняя граница распределения m₀ зависит от системы регистрации землетрясений и обычно определяется границей соблюдения прямолинейности закона ГР на левом хвосте распределения. Мы полагаем m₀ известным. Параметр b – это наклон прямолинейной ветви ГР в логарифмическом масштабе; ξ – параметр формы ОРП, характеризующий остроту “клюва” на конце распределения; C₁, C₂, С₃ – константы (они зависят от указанных параметров), которые должны обеспечить нормировку функции распределения F(х) и ее непрерывность:

Выбор модели М2 обусловлен следующими соображениями. Общепризнано, что в диапазоне малых и средних магнитуд распределение хорошо описывается законом Гутенберга–Рихтера. Относительно диапазона больших магнитуд (хвоста распределения) у сейсмологов нет единого мнения. Ясно только одно: из физических соображений распределение магнитуд должно быть конечным, так что неограниченный закон Гутенберга–Рихтера не удовлетворяет этому требованию. Усеченное распределение Гутенберга–Рихтера является наиболее простой моделью, удовлетворяющей этому требованию; именно эта модель чаще всего и используется в приложениях. Имеются и другие модели, мы не будем останавливаться на их анализе. Модель (11) основана на том, что в пределе, при большом числе наблюдений, хвост распределения должен подчиняться закономерностям теории экстремальных значений. В соответствии с этой теорией распределение хвоста должно сходиться к обобщенному распределению Парето (см. [Gumbel, 1956; Embrechts et al., 1997; Pisarenko, Rodkin, 2010]). Это и реализовано в модели (11). Модель (11) задается тремя параметрами: h, b, ξ. Заметим, что в общем случае на значение параметра ξ в ОРП законе распределения ограничений не накладывается. Но при ξ ≥ 0 распределение ОРП, как и неограниченный закон Гутенберга–Рихтера, не отвечает физическому требованию конечности распределения. Поэтому в модель М2 в (11) введено дополнительное условие ξ < 0.

Мы выбрали 6 различных прообразов для реализации модели (11). Соответствующие им значения параметров h, b, ξ даны в табл. 1. Прообразами для этих вариантов послужили региональные каталоги 6 различных сейсмоактивных регионов – района Срединно-Атлантического хребта, Японии, Курильских островов, Новых Гебрид, Перу и Филиппин, далее эти регионы отмечены номерами с 1-го по 6-й. Данные по этим регионам (см. рис. 1) были выбраны из каталога ISC-GEM. Величина m₀ для всех регионов принята единообразно m₀ = 6.0. Критерии выбора этого набора регионов были следующие: сейсмотектоническая однородность региона, статистическая представительность, разнообразие в поведении графика повторяемости и, в особенности, хвоста распределения. Желательным дополнительным условием было возможное разнообразие сейсмотектонического характера регионов.

Из общих теоретических соображений наиболее ожидаемым поведением хвоста интегрального графика повторяемости можно полагать его более или менее систематический загиб вниз. Такие графики действительно наиболее типичны. В нашей выборке им отвечают область Срединно-Атлантического хребта, Новых Гебрид и Филиппин, соответствующие номера графиков 1, 4, 6. Однако эмпирические графики повторяемости таким типом поведения хвоста распределения далеко не исчерпываются. Часто встречаются, по-видимому, случайные нерегулярные вариации поведения хвоста распределения. Так, например, даже и на кумулятивных графиках повторяемости достаточно часто встречаются горбы (cм. рис. 2), которые, конечно, могут иметь случайный характер. Впрочем, они могут быть и не случайными, а соответствовать повышенной вероятности сильных землетрясений в некотором диапазоне магнитуд. В случае если такой диапазон достаточно узок и расположен близко к самому хвосту распределения, такие землетрясения называют характеристическими [Ellsworth, 1995]. Обсуждению проблемы существования и распространенности эффекта характеристических землетрясений посвящено большое количество работ [Stein et al., 2005; Jackson, Kagan, 2011; и др.]. В рамках модели нелинейного мультипликативного каскада [Родкин, 2011; Шерман и др., 2017] “горб” на графике повторяемости землетрясений возникает, если в некотором диапазоне размеров событий возникает нелинейная положительная связь. Аналогично, загиб вниз графика повторяемости на максимально возможном масштабе событий порождается возникающей здесь нелинейной отрицательной связью. Согласно анализу исторических данных по сейсмичности Азии [Шерман и др., 2017], ситуация с положительной обратной связью реализуется довольно часто.

Вышесказанное дает основание рассмотреть в качестве возможных графиков повторяемости также и случаи с систематическим завышением числа сильных землетрясений в некотором диапазоне сильных событий. Среди рассмотренных регионов таковы случаи Японии, Курил и Перу, номера регионов 2, 3, 5 (рис. 2). Напомним, что такие отклонения можно трактовать и как относительно маловероятные случайные флуктуации хвоста графика повторяемости; впрочем, со строго статистических позиций аналогично можно трактовать и загиб графика повторяемости вниз. Для нашей работы эти аспекты сейсмичности не являются центральными. Основная цель нашей работы заключается в том, чтобы оценить реальную эффективность оценки квантилей максимальных событий в будущем интервале времени Т на примерах с точно известным ответом, а также продемонстрировать устойчивость и робастность этих квантилей как информативной и важной характеристики сейсмического риска.

Графики эмпирических распределений магнитуд и их “прототипы” по модели (11) для перечисленных выше регионов приведены на рис. 2. В табл. 1 приведены значения параметров эмпирических распределений, выбранных для тестирования. Мы выбирали их близкими к оценкам параметров, полученных методом ММП, слегка изменив их в пределах стандартного отклонения. Для региона Японии мы намеренно выбрали для параметра ξ очень малое по абсолютной величине значение ξ с тем, чтобы посмотреть, как ведут себя оценки в “экстремальных” условиях (при значениях ξ, соответствующих очень большим возможным значениям М_мах). Поскольку в данной постановке задачи нас интересуют только относительно сильные землетрясения, нижний порог был выбран достаточно большим: m₀ = 6.0. Такие землетрясения надежно регистрируются в каталоге ISM-GEM на всем его протяжении.

Как видно из рис. 2, в регионах (2, 3, 5) эмпирические и теоретические функции имеют заметные расхождения в диапазоне больших магнитуд М > 7. В значительной степени эти расхождения являются следствием принятого выше условия ξ < 0 в модели М2. Мы не стали заниматься более тщательной подгонкой теоретических распределений к эмпирическим, так как в данном случае не это являлось нашей целью. Нам нужны были логически корректные распределения с известным ответом, которые можно было бы рассматривать как в принципе возможные. Заметим также, что отклонения эмпирических распределений в регионах (2, 3, 5) можно рассматривать и как порождаемые эффектом характеристических землетрясений, и как случайные отклонения, вызванные нестационарностью сейсмического режима и недостаточной статистикой в диапазоне больших магнитуд. Наш опыт показывает, что отклонения подобного рода могут нивелироваться, если есть возможность пополнить статистику сильных событий.

Заметим также, что поскольку наш подход использует формулы, выведенные в предположении, что сейсмический поток является стационарным Пуассоновским потоком, мы должны были провести декластеризацию исходных каталогов, чтобы избавиться от афтершоков. Мы воспользовались для этой цели методом, использующим обобщенное расстояние между двумя событиями, учитывающее как пространственную, так и временную компоненту (см. подробнее [Писаренко, Родкин, 2019]). В результате получился поток событий, гораздо лучше соответствующий Пуассоновскому свойству равномерного распределения событий на интервале времени при их фиксированном числе. Естественно, выборки уменьшились, что отразилось в уменьшении оценок интенсивности λ. События, оставшиеся после применения процедуры декластеризации, мы называем главными толчками. Процедура декластеризации не должна сильно искажать свойства максимальных толчков, поскольку именно максимальные события, как правило, и являются главными. Она влияет, главным образом, на снижение интенсивности λ. Декластеризация привела к следующим процентам главных толчков (по отношению ко всем толчкам):

Атлантика 98.9% (257 толчков); Япония 47.4% (245); Курилы 55.0% (236); Новые Гебриды 72.3% (413); Перу 88.1% (89); Филиппины 82.1% (377). Низкий процент главных толчков в Японии объясняется большим количеством афтершоков после мегаземлетрясения Тохоку (2011 г.).

В соответствии с теоретическими распределениями, изображенными на рис. 2, мы генерировали по 1000 каталогов для каждого региона, причем объем каталога соответствовал реальному каталогу основных событий. Для каждого искусственного каталога число событий и магнитудный диапазон соответствовали указанным выше значениям для данного региона. Оценка параметров производилась двумя методами: методом максимального правдоподобия (ММП) и Байесовским методом (БМ) с использованием усеченного распределения Гутенберга–Рихтера (УГР). При оценивании ММП параметр h полагался равным 75%-му квантилю выборки магнитуд, так что фактически искался максимум правдоподобия по двум параметрам b и ξ. Для региона Японии на область параметров, по которой находился максимум функции правдоподобия, к условию ξ < 0 добавлялось еще условие M_max < 11.5. В этом регионе за относительно краткое время регистрации произошло редкое по силе событие Тохоку, и использование лишь одного условия ξ < 0 приводило к выводу о принципиальной возможности возникновения событий неправдоподобно большой магнитуды, из-за чего мы были вынуждены использовать указанное выше дополнительное условие. Для остальных регионов действовало только первое условие

Байесовский метод оценки квантилей описан в Приложении 1.

Для каждого из 1000 каталогов рассматриваемого региона оценивались квантили максимального события в будущем интервале времени Т. Для квантилей были выбраны уровни q = 0.5, 0.9, 0.95, 0.975, 0.99, 0.999. Для интервалов Т мы рассчитывали довольно подробную сетку Т = 10, 25, 50, 100, 200, 1000 лет. В статье мы решили ограничиться представлением данных для Т = 50 лет, так как для остальных значений результаты получались качественно похожими; при этом расчеты для больших Т, согласно соотношению (8), эквивалентны расчетам для Т = 50 с высокими квантилями.

На рис. 3 для Т = 50 лет и для 6 регионов показаны 3 типа кривых: синим цветом – истинное модельное значение квантиля, черным цветом – средние значения по методу ММП (по 1000 каталогам, они фактически не отличаются от средних по ансамблю из-за большого объема выборки), и красным цветом – средние значения по методу БМ. По оси Х изменения q даны в логарифмическом масштабе для величины 1/(1 – q), которая дает квантиль q в удобном нелинейно растянутом виде, используемом в теории информации.

На рис. 4 приведены среднеквадратичные отклонения (MSE):

где: Q_k – оценка квантиля по k-му каталогу; n – число каталогов (n = 1000); Q_true – истинное значение квантиля. MSE показывает в среднеквадратичной метрике насколько близка данная оценка к истинному значению, т.е. дает показатель качества оценки. Этот показатель складывается из двух компонент: систематического смещения BIAS(Q) и стандартного отклонения от выборочного среднего STD(Q):

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

Имея 1000 каталогов, можно оценить стандартные отклонения оценок квантилей. Они показаны на рис. 6 черным (ММП) и красным цветами (БМ). Для реальной ситуации (один каталог) таких оценок получить нельзя. Для того, чтобы оценить стандартное отклонение в этих случаях, мы использовали следующий подход. В методе ММП мы брали оценки параметров (b, ξ), полученные по имеющемуся реальному каталогу, и генерировали 1000 каталогов с этими оценками, а затем оценивали стандартные отклонения квантилей по этим 1000 каталогам и принимали их за истинные. Из-за того, что оценки параметров М2 имеют погрешность, оценки стандартных отклонений квантилей также будут иметь погрешность, но она будет иметь высший порядок малости. Расчеты показали, что погрешность оценки стандартных отклонений при таком подходе не превышает 5%, так что этой погрешностью можно пренебречь. Для БМ методика оценивания стандартного отклонения оценки квантиля по одной реализации описана в Приложении 1. Эта оценка показана на рис. 6 лиловым цветом. Как показывает наш опыт, она несколько занижена и обычно составляет 0.8–0.9 по отношению к истинному стандартному отклонению, то есть это различие практически несущественно.

Мы видим, что в целом оценки квантилей довольно близки к истинным значениям. Среди 6 регионов, точнее среди 6 искусственных прототипов региональных каталогов, можно выделить 2 региона: Атлантику и Японию. Для Атлантики характерны небольшие значения максимальных магнитуд (по каталогу объема 257), как правило не превосходящие 7.7. График повторяемости круто загибается вниз (см. рис. 2), что отражается в оценках параметра формы ОРП ξ, он имеет относительно большие по абсолютной величине отрицательные значения –0.35 ~ –0.25. Значения MSE оценок квантилей Атлантики уровня q = 0.9 равны примерно 0.11, а для уровня q = 0.999 (очень высокий уровень, практически дающий максимально возможную магнитуду распределения) получаем примерно 0.16. Разница весьма мала, что свидетельствует о большой робастности квантилей при не слишком малых по абсолютной величине значениях параметра ξ. Для “экстремального” региона – Японии – значения MSE оценок квантилей уровня q = 0.9 равны примерно 0.35, а для уровня q = 0.999 получаем очень большое значение 0.95. Значения параметра ξ в искусственных каталогах малы по абсолютной величине и доходят практически до нуля.

В остальных 4 “регионах” наблюдается промежуточная ситуация. Заметим, что расхождения заметно выше в регионах 2, 3, 5 (Япония, Курилы и Перу), для которых на рис. 2 наблюдался горб на эмпирических графиках повторяемости и для которых, соответственно, значения параметра ξ наиболее близки к нулю. Для этих же регионов максимален и разброс оценок.

Обращает на себя внимание, что во всех случаях полученные оценки лежат ниже модельных “истинных” квантилей. Возможно, это обусловлено используемой моделью. В методе БМ к этому могло привести исходно задаваемое априорное ограничение магнитуд. В методе ММП к этому смещению могло приводить задаваемое при расчетах вариантов требование ξ < 0. Значения MSE оценок квантилей уровня q = 0.9 лежат (рис. 4) для всех регионов, кроме Японии, в интервале 0.1–0.5 (для Японии – 0.7), а для уровня q = 0.999 получаем 0.15–0.8 (для Японии – 2.2). Эти величины дают представление о том, насколько точны оценки разброса, полученные по одному реальному каталогу.

Рассмотрим теперь систематические смещения оценок. Для Атлантики смещения пренебрежимо малы: от –0.07 до +0.02 (см. рис. 5). Для Японии они велики: от –0.9 до –2.0, т.е. имеет место сильная недооценка истинного квантиля. Правда, модельные квантили были взяты очень большими: 9.6 для уровня q = 0.9; 11.7 для уровня q = 0.999, что вызывает сомнения в реальной возможности таких магнитуд. Но, как уже отмечалось выше, столь большие значения модельных ожидаемых магнитуд для Японии были выбраны для того, чтобы посмотреть, как поведут себя оценки в ситуации “экстремальных” условий (при остром клюве распределения магнитуд, соответствующих очень малым по абсолютным значениям параметра ξ в модели М2). В “регионах” Атлантики, Курилы, Новые Гебриды, Филиппины оценки ММП дают меньшее смещение, чем оценки БМ, но это, возможно, связано с тем, что искусственные каталоги генерировались с помощью модели М2. Для всех регионов кроме Японии смещения при уровне доверия q = 0.90 варьируют в диапазоне 0.05 ~ 0.15, что вполне приемлемо, учитывая саму точность, с которой определяются магнитуды.

Перейдем теперь к оценкам стандартных отклонений (СО) квантилей (рис. 6). Они не должны превосходить MSE и должны убывать с ростом объема выборки n как 1/n^1/2. В регионе “Атлантика” СО заметно меньше, чем в других регионах, примерно 0.07 ~ 0.15. В остальных “регионах” на уровне q = 0.9 СО несколько больше (0.2–0.3). Таким образом, MSE формируется, в большей степени, за счет случайного разброса, а не за счет систематического смещения оценок, кроме региона Японии. Большие значения СО для “региона” Перу обусловлены сравнительно малым объемом выборки (n = 89) по сравнению с другими “регионами”.

Подчеркнем еще раз связь формы хвоста распределения магнитуд со статистическими характеристиками разброса оценок квантилей. Чем острее “клюв” распределения хвоста, тем больше разброс оценок квантилей. “Острый клюв” означает, что возможны очень сильные толчки, но они могут происходить чрезвычайно редко. Это отражается в значениях параметра формы ξ предельного распределения магнитуд: “острый клюв” соответствует очень малым (отрицательным) значениям ξ $ \cong $ –0.1…–0.01 (и еще ближе к нулю), в то время как относительно большие (по модулю) значения ξ $ \cong $ –0.4…–0.2 дают “тупой клюв” и приводят к более устойчивым статистическим оценкам.

В целом получаем, что при Т = 50 лет и q ≤ 0.95 отклонения значений квантилей от модельного “истинного” результата, как правило, не превосходят 0.2 единицы магнитуды и только в случае Японии приближаются к 0.5. Такая погрешность может полагаться допустимой. Для больших значений ξ возрастают, причем в среднем отличие больше для случаев с меньшим по абсолютной величине значением параметра ξ.

ОЦЕНКА КВАНТИЛЕЙ НА РЕАЛЬНЫХ КАТАЛОГАХ

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

Сравнивая результаты расчетов для реальных каталогов и модельных случаев, прежде всего, отметим увеличение оценок квантилей для реальных каталогов. Различие наиболее существенно для регионов 2, 3 и 5, где оно в среднем достигает на максимальных квантилях 0.4–0.5 единицы магнитуды (0.8 для Японии). Максимальные различия величин квантилей при расчете двумя разными методами не превышают 0.4–0.5 единиц магнитуды; только в одном случае региона 4 различие достигает величины 0.6. Эти различия можно объяснить тем, что мы выбирали параметры искусственных каталогов, прежде всего, с целью испытать оценки в экстремальных условиях и не старались воспроизвести возможно более точно копии “прототипов”.

В целом оценки по двум методам вполне удовлетворительно согласуются. Оценка методом БМ является более “консервативной”: для уровней доверия, близких к единице (для значений магнитуд, близких к максимально возможным), они ниже оценок ММП. Оценки разброса весьма близки (за исключением одного региона – Японии, где значительное расхождение оценок разброса может быть объяснено упоминавшейся выше “остротой клюва” распределения магнитуд).

ОБСУЖДЕНИЕ

В настоящей работе, по-видимому, впервые оценка смещения, стандартного отклонения и среднеквадратичного отклонения проведена на большом количестве (1000 штук) независимых выборок (искусственных каталогов) с известным законом распределения. Это дало возможность сравнить такую “идеальную” точность оценок с той, которая получается по одному каталогу с неизвестным законом распределения. Полученные модельные характеристики разброса относительно известного ответа оказались несколько больше характеристик, получаемых по одному (реальному) каталогу. Наибольшие значения MSE получились для региона Японии, и это было ожидаемо, поскольку в этом регионе в 2011 г. произошло мегаземлетрясение Тохоку, резко нарушившее весь сейсмический режим этого региона. Появление в 2011 г. в каталоге этого мегаземлетрясения привело к значительному уменьшению по абсолютной величине параметра формы ξ, что привело к росту оценок квантилей и при больших значениях T и q и к снижению устойчивости этих оценок.

Для остальных регионов MSE квантиля максимального события за будущие Т = 50 лет при уровне доверия 0.9 принимают значения порядка 0.1–0.3. Такой точности оценок можно ожидать для регионов, в которых имеются каталоги землетрясений длительности Т_catлет с интенсивностью λ для событий М > 6.0, обеспечивающей условие n = λТ_cat > 100. Достижению существенно большей точности препятствуют статистические неопределенности в оценках поведения хвостов распределений. Пространственная разрешающая способность получаемых оценок зависит от интенсивности сейсмического потока. Чем больше интенсивность, тем выше разрешающая способность, поскольку имеется возможность уменьшить регион с сохранением необходимого объема наблюдений.

Проведенное сравнение не дает оснований для предпочтения одного из использованных методов оценки параметров. Существующие обрывочные данные по максимальным магнитудам исторических и палеоземлетрясений, например события с М > 8.5 для Курильской зоны субдукции [Yukinobu Okamura, Yuichi Namegaya, 2011], ложатся в пределы, полученные обоими методами, и также не дают основания для выбора одного из вариантов расчета.

Расчет методом БМ ожидаемо дает более консервативную оценку и меньший рост величин М_мах в области больших интервалов времени Т и доверительных пределов q. Возможно, именно это различие и будет в дальнейшем являться ключевым при выборе решения рекомендовать тот или иной подход. В пользу консервативного (БМ) варианта решения говорят соображения о возможном ограничении энергетического потенциала для реализации экстремально сильных сейсмических событий. В пользу варианта ММП говорит теоретическая общность теории экстремальных значений и предельных распределений.

В диапазоне квантилей для Т = 50 лет и q ≤ 0.90 погрешности определения величин М_мах обоими методами лежат в пределах практических требований к такого рода оценкам. Для более высоких квантилей расхождения больше, но указать предпочтительную модель затруднительно; для разных регионов меньшие погрешности дают разные модели.

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

Проведенные сравнения и использование соотношения (8) позволяют получить оценки ожидаемого разброса значений квантилей М_мах в области больших – заведомо ненаблюдаемых – значений интервалов времени Т и при весьма строгих границах доверительного интервала q. Оценки квантилей максимальных событий в будущем интервале времени получаются существенно более робастными при более сильном загибе вниз графика повторяемости. В регионах, где такой загиб не выражен или где наблюдаются отклонения от графика повторяемости, которые могут быть связаны с возникновением так называемых характеристических землетрясений, уровень и разброс оценок сейсмической опасности существенно возрастают.

Заметим в заключение, что рассмотренные вопросы касаются не только проблемы оценки сейсмической опасности, но и многих иных ситуаций, когда имеют место степенные распределения с тяжелым хвостом [Burroughs, Tebbens, 2001; Pisarenko, Rodkin, 2014; Beirlant et al., 2017].