Физика Земли, 2022, № 1, стр. 90-99

ВВЕДЕНИЕ

Закон Гутенберга–Рихтера для распределения количества землетрясений в данной области имеет вид [Gutenberg, Richter, 1954; 1956]:

где: m – магнитуда землетрясения; m₀ – нижний порог регистрации; N(m) – среднее число землетрясений, превосходящих магнитуду m на некотором интервале времени; a, b – коэффициенты, характеризующие сейсмичность в данной области (b > 0). Если рассматривать магнитуду произвольно выбранного землетрясения как случайную величину μ, то после соответствующей нормировки и введения функции распределения магнитуд F(m) мы получаем из (1) соотношение:

Функция распределения F(m) имеет вид:

Для удобства изложения мы будем пользоваться натуральными логарифмами и показательной функцией. При этом функция распределения F(m) и плотность распределения f(m) задают экспоненциальный закон:

Мы будем использовать параметр s обратный по отношению к β:

Параметр s измеряется в тех же единицах, что и магнитуда. Со статистической точки зрения параметр s несколько удобнее. Оценкой максимального правдоподобия для параметра s является среднее арифметические от магнитуд выборки. Оно имеет с точностью до множителя стандартное χ²-распределение с 2n степенями свободы (n – объем выборки), быстро сходящееся к нормальному закону. В то же время оценка обратной по отношению к s величины β имеет при небольших n скошенное распределение и ее сходимость к предельной величине может сопровождаться выбросами, т.е. является неробастной. В настоящей статье мы будем заниматься усеченным законом Гутенберга–Рихтера (УГР):

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

ОЦЕНКА ПАРАМЕТРА М

Оценка максимально возможной региональной магнитуды занимает важное место в проблеме оценки сейсмического риска (см. [Kijko, Sellevoll, 1989; 1992; Писаренко, 1991; Kijko, Graham, 1998; Pisarenko et al., 1996; Dargahi-Noubary, 2000; Kagan, Schoenberg, 2001; [Kijko, 2004; Pisarenko, Rodkin, 2010a; Zoller, Holschneider, 2016; Vermeulen, Kijko, 2017; Beirlant et al., 2019; Писаренко и др., 2021]. В нашей постановке максимальная региональная магнитуда выражается параметром М. Общая схема оценивания параметра М, предлагаемая нами, такова. Сначала рассматривается оценка параметров (М, s) стандартным методом правдоподобия. Обозначим имеющуюся выборку магнитуд (каталог) x = (x₁, … , x_n). Функция правдоподобия имеет вид:

Как известно (см. например, [Pisarenko et al., 1996]), оценкой максимального правдоподобия для параметра М при любом значении s является максимальная магнитуда наблюденной выборки μ_n:

Оценка максимального правдоподобия для параметра s находится как максимум правдоподобия для распределения УГР (9), в котором параметр М заменен на μ_n, т.е. как максимум функции от s:

Полученная оценка μ_n максимального правдоподобия для параметра М будет иметь систематическое, отрицательное смещение относительно истинного значения параметра М, особенно заметное при умеренных значениях n. В Приложении получена точная формула для этого смещения в виде конечной эффективно вычисляемой суммы. Эта формула имеет вид:

где: u = 1 – exp[–(M – m₀)/s]; W_n = u + $\frac{{{{u}^{2}}}}{2}$ + … + $\frac{{{{u}^{n}}}}{n}$; Е – символ математического ожидания. Как мы видим, смещение зависит от неизвестного параметра М. В предлагаемой оценке предлагается заменять его на μ_n. Параметр s нужно заменить на его стандартную оценку максимального правдоподобия $\bar {s}$. Таким образом, окончательно новая оценка $\bar {M}$ параметра М, в которой сделана поправка на смещение, имеет вид:

где: $\bar {U}$ = 1 – exp[–(μ_n – m₀)/$\bar {s}$]; $\overline {{{W}_{n}}} $ = $\bar {U}$ + $\frac{{{{{\bar {U}}}^{2}}}}{2}$ + … + $\frac{{{{{\bar {U}}}^{n}}}}{n}.$

Мы обобщим ниже оценку (13), адаптируя ее на случай оценки квантиля Q(q) распределения УГР, т.е. значения Q(q), определяемого уравнением

где q (уровень значимости квантиля) – число из интервала [0, 1]. Квантиль является обратной функцией по отношению к функции F(x). Если функция распределения F(x) непрерывна (что мы будем предполагать), то квантиль и функция распределения монотонно возрастают и взаимно-однозначно связаны друг с другом. Для УГР имеем: где, как и раньше, u = 1 – exp[–(M – m₀)/s]. Для квантиля (15) можно взять оценку, подставив в (15) μ_n вместо М и $\bar {s}$ вместо s. Эта оценка будет иметь смещение

Пользуясь тем же приемом, который мы использовали при выводе оценки $\bar {M}$, находим:

где $W_{n}^{q}$ = $~u{{q}^{{1/n}}}$ +$\frac{{{{{(u{{q}^{{1/n}}})}}^{2}}}}{2}$ + … + $\frac{{{{{(u{{q}^{{1/n}}})}}^{n}}}}{n}$. Вычитая из оценки Q(q|μ_n, $\bar {s}$) правую часть (17), в которой неизвестные параметры М, s заменены соответственно на μ_n, $\bar {s}$, получаем окончательную оценку квантиля Q(q|M, s):

(18)

$Q\left( q \right) = Q(\left. q \right|{{\mu }_{n}},\bar {s}) - \frac{{\bar {s}}}{{\overline {{{U}^{n}}} q\overline {\beta ~} }}[\lg (1 - \bar {U}{{q}^{{1/n}}}) + \bar {W}_{n}^{q}],$

где: $\bar {U}$ = 1 – exp[–(μ_n – m₀)/$~\bar {s}$]; $\bar {W}_{n}^{{q~~}}$ = $\bar {U}{{q}^{{1/n}}}$ + + $\frac{{{{{(\overline U {{q}^{{1/n}}}{\text{\;}})}}^{2}}}}{2}$ + … + $\frac{{{{{(\overline U {{q}^{{1/n}}}{\text{\;}})}}^{n}}}}{n}$.

При q = 1 оценка (18) совпадает с оценкой (13) для максимальной магнитуды.

СРАВНЕНИЕ РАЗЛИЧНЫХ ОЦЕНОК М

В работах [Kijko, Sellevoll, 1989; 1992; Kijko, Graham, 1998; Coles, Dixon, 1999; Kijko, 2004; Holschneider et al., 2011; Kijko, Singh, 2011; Lasocki, Urban, 2011; Zoller, Holschneider, 2016; Vermeulen, Kijko, 2017; Beirlant et al., 2019] введено большое количество оценок, полученных по принципу статистических оценок моментного типа. Для их получения выписывается выражение для среднего значения (или некоторого старшего момента) некоторой состоятельной оценки, сходящейся по вероятности к истинному значению при n $ \to \infty $. Например, с помощью интегрирования по частям получаем выражение для среднего значения максимального значения выборки μ_n:

Далее, в этом уравнении заменяют E{μ_n}, s соответственно на μ_n, $\bar {s}$, обосновывая это тем, что при n $ \to \infty ~$ обе величины сходятся к своему математическому ожиданию, и получают уравнение для неизвестного параметра М:

Решая это уравнение относительно М различными приближенными способами, можно получить статистическую оценку неизвестного параметра М. В качестве некоторых вариантов этого метода можно еще заменять в правой части (20) М на μ_n. Мы будем рассматривать вариант уравнения (20), предложенный в работах [Kijko, 2004; Kijko, Singh, 2011]. Соответствующую оценку М обозначим МK. Эта оценка определяется как решение трансцендентного уравнения (20) относительно М. Можно показать, что эта оценка МK определена не при любых возможных значениях μ_n. Для значений μ_n достаточно близких к М уравнение (20) не имеет решения, а теоретическая средне-квадратичная ошибка (MSE) оценки МK равна бесконечности. Поэтому для проведения сравнения с MSE других оценок мы вынуждены обрезать МK некоторым порогом h:

Рассмотрим еще одну оценку максимальной магнитуды М.

В работе автора [Pisarenko et al., 1996], введена несмещенная оценка максимальной магнитуды М, обладающая наименьшей дисперсией среди всех несмещенных оценок:

К сожалению, оценка (22) имеет довольно большую дисперсию. Поэтому мы рассмотрим усеченный вариант этой оценки:

Этот вариант оценки вполне конкурентоспособен, его MSE уменьшается (в зависимости от выбранного порога усечения). Но, естественно, появляется систематическое смещение (впрочем, не очень большое).

Сравним приведенные выше 4 оценки параметра М максимальной магнитуды.

На рис. 1 показаны MSE для четырех перечисленных выше оценок. Для усреднения использовалось N = 10 000 искусственных каталогов со значениями параметров m₀ = 6.0, M = 8.0, s = 0.4, порог отсечения для оценок MK_trunk и MP_trunk равен h = μ_n + 1; Оценка МK_trunk менее эффективна, чем другие оценки для всех n.

На рис. 2 показаны смещения для 4-х оценок магнитуды М. Смещения оценки $\bar {M}$ меньше по абсолютной величине, чем смещения оценки МK_trunk. Отметим, что смещения оценки МK_trunk положительны, в то время как смещения оценки $\bar {M}$ отрицательны.

На рис. 3 показаны стандартные отклонения 4-х оценок.

Сравнение оценок проведено для типичных значений параметров УГР: s = 0.4; M = 8.0; m₀ = 6.0. Мы использовали для сравнения и другие варианты параметров. Получаются несколько другие оценки, но в целом из этих разных вариантов можно сделать те же выводы, что мы получили для использованных значений параметров: наиболее эффективной (в смысле среднеквадратичного отклонения MSE) является оценка $\bar {M}$, близка к ней Байесовская оценка, наименее эффективной показала себя оценка МK_trunk.

ОЦЕНКА КВАНТИЛЕЙ МАКСИМАЛЬНОГО ЗЕМЛЕТРЯСЕНИЯ В БУДУЩЕМ ИНТЕРВАЛЕ ВРЕМЕНИ Т

Как уже отмечалось, параметр М в большинстве практических ситуаций оценивается неустойчиво из-за недостаточного количества наблюдений в диапазоне сильных землетрясений. К тому же возникают неясности, к какому интервалу времени относится параметр М: 1000 лет, 1 000 000 лет или “на все времена”. Поэтому в работах [Pisarenko et al., 2008; 2010; Pisarenko, Rodkin, 2009; 2010a; 2010b; 2013; Zoller et al., 2013], было предложено использовать для характеристики сейсмичности в диапазоне сильнейших событий квантили Q_T(q) уровня q случайной величины – максимальной магнитуды М(Т) землетрясения в будущем интервале времени Т (в данном регионе). Ниже мы будем считать, что каталог декластеризован с помощью какого-либо из известных методов декластеризации (см. например, [Писаренко, Родкин, 2019]) и его можно считать Пуассоновским случайным потоком с некоторой интенсивностью λ (число событий в единицу времени в заданном магнитудном диапазоне). Среднее число событий равно λТ, стандартное отклонение от среднего $\sqrt {\lambda Т} $. Случайная величина М(Т) корректно определена и при определенных условиях заметно устойчивее и робастнее, чем статистические оценки параметра М. К тому же она включает в себя в качестве частного случая при q = 1 и Т $ \to \infty $ оценку параметра М “на все времена”. Мы приводим ниже формулы для оценки квантилей Q_T(q) в рамках модели УГР. Они выводятся по той же методике, по которой выше получены формулы для оценки параметра М:

1. Берется оценка максимального правдоподобия для Q_T(q), обозначаемая ${{\hat {Q}}_{T}}$(q). Она зависит от максимальной магнитуды выборки μ_n: ${{\hat {Q}}_{T}}$(q) = ${{\hat {Q}}_{T}}$(q| μ_n).

2. Оценка ${{\hat {Q}}_{T}}$(q|μ_n) будет иметь заметное смещение (особенно сильное при небольших объемах выборки):

${\text{Bias }}({{\hat {Q}}_{T}}(\left. q \right|{{\mu }_{n}}))~ = E{{\{ }_{T}}({{\hat {Q}}_{T}}({q \mathord{\left/ {\vphantom {q {{{\mu }_{n}}}}} \right. \kern-0em} {{{\mu }_{n}}}}){\text{ }}\} ~\,\, - \,\,~{{Q}_{T}}\left( q \right).$

Это смещение будет функцией от неизвестного параметра М. Обозначим эту функцию G(M):

Оценкой максимального правдоподобия для смещения G(M) будет G(μ_n).

3. Результирующая оценка для квантиля Q_T(q), обозначаемая ${{\bar {Q}}_{T}}$(q|μ_n), получается с помощью вычитания из оценки ${{\hat {Q}}_{T}}$(q|μ_n) величины G(μ_n):

В Приложении выведены формулы для ${{\hat {Q}}_{T}}$(q|μ_n), G(μ_n):

где: $\bar {q}$ = $\frac{1}{{{{\lambda }}T}}$ lg[(1 – exp(–λT))q + exp(–λT)];

ПРИМЕР ОЦЕНИВАНИЯ М ПО РЕАЛЬНОМУ КАТАЛОГУ

Рассмотрим каталог СМТ (1976–2015 гг.) для региона Курильских островов и Камчатки:

42.81 °≤ широта ≤ 53.56°;

146.38 °≤ долгота ≤ 161.06°.

Для декластеризации каталога мы применили алгоритм, описанный в работе [Писаренко, Родкин, 2019]. В качестве нижнего порога представительной регистрации событий выбрана магнитуда m₀ = 5.7. В результате получилось 158 главных толчков в диапазоне магнитуд m ≥ 5.7. Интенсивность сейсмического потока составила: λ = 3.9606 $\frac{1}{{{\text{год}}}}$. Максимальное землетрясение m = 8.296 произошло 15.11.2006 г., широта 46.57°; долгота 153.29°. На рис. 4a показана пространственная конфигурация анализируемых землетрясений; на рис. 4б – график повторяемости.

Результаты для 4-х рассмотренных выше оценок параметра М представлены в табл. 1.

Мы пробовали 5 вариантов “свободных параметров” Байесовского метода. При этом получили довольно большой разброс оценок параметра М и его стандартного отклонения: от 16% до 25% относительных значений. Такой разброс говорит о том, что в данном случае Байесовские оценки параметров довольно существенно зависят от выбора априорной области для параметров. Допустимым можно было бы считать разброс, не превышающий 10%. Однако при выбранных нами значениях свободных параметров (априорные интервалы для М, s : μ_n ≤ M $ \leqslant ~$ μ_n + 1; 0.25 ≤ s ≤ 0.75) значения Байесовских оценок оказались близки к новым оценкам, введенным в этой статье, и мы оставили их для сравнения.

На рис. 5 представлен график оценки квантиля Q_T(0.95) ± std.

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

Дисперсию и стандартное отклонение приведенных выше оценок можно было бы оценить с помощью того же приема, который мы использовали для оценивания смещения. Однако при этом получаются весьма сложные формулы, и мы укажем более простой и универсальный путь оценивания. Этот метод основан на идеях бутстреп-оценок [Efron, 1979]. Тем или иным способом строится модель наблюдаемого каталога, а затем с помощью случайных чисел генерируется достаточно большое число случайных каталогов (скажем, 10000), по которым можно с нужной точностью оценить стандартное отклонение оценок. В качестве модели распределения мы предлагаем модель УГР, в которую вместо неизвестных параметров вставлены их оценки максимального правдоподобия. Конечно, эти оценки будут отличаться от неизвестных истинных значений, но это будут отличия, так сказать, “второго порядка” малости по сравнению с самим стандартным отклонением. При больших и умеренных объемах выборки они будут малы и ими вполне можно пренебречь. Однако при малых выборках применение этого метода требует тщательного контроля и дополнительных исследований. Мы рекомендуем подсчитать несколько вариантов моделей, в которых оценки параметров, полученные методом максимального правдоподобия, подвергнуты небольшим возмущениям, и посмотреть, насколько сильно будут меняться оценки стандартного отклонения.

Описанным выше способом оценивались характеристики разброса оценок, показанных на рис. 3, рис. 5 и в табл. 1.

Следует отметить, что в регионе Курильской островной дуги 04.11.1952 г. произошло сильное землетрясение (широта = 52.623°; долгота = = 159.779°), магнитуда которого оценивается по разным каталогам как 8.9–9.0. Оно не вошло в использованный нами каталог CMT, который начинается с 1976 г. Впрочем, магнитуда этого землетрясения входит в интервал разброса оценок, указанных в таблице, так что можно считать, что она не противоречит этим оценкам.

Если взять оценки максимальной возможной магнитуды Курильской островной дуги, полученные, исходя из региональных тектонических характеристик, то они в целом согласуются с оценками нашей таблицы (см. [Тараканов, 1990; Ермаков, 1997]). Так, оценка максимальной возможной магнитуды по номограмме Ю.В. Ризниченко дает значения M_max = 8.0–8.5, а по номограмме Р.З. Тараканова M_max = 8.6–9.5. Как отмечается в работе [Ермаков, 1997], оценки Тараканова являются завышенными, т.к. не учитывают разделения Среднего Курильского блока на среднюю и переходную части.

МОДЕЛИ С НЕТОЧНО ЗАДАННЫМИ МАГНИТУДАМИ

Довольно часто модель УГР используют с некоторым “дополнением” (см. например, [Kijko, Sellevol, 1989; Lyubushin, Parvez, 2010]), а именно, считают, что публикуемые в каталогах магнитуды m_к (их называют “кажущимися” – apparent) получаются из некоторой “истинной” магнитуды m_к0 путем добавления некоторой случайной ошибки ${{\Delta }}$m_к:

Случайную ошибку считают равномерно распределенной на некотором отрезке [$ - {{\Delta }};~ + {{\Delta \;}}]$. Для ${{\Delta }}$ обычно используют значения порядка 0.2–0.5. Эта модель удобна тем, что для нее плотность вероятности ${{{{\varphi }}}_{\Delta }}\left( m \right)$ очень просто выражается через исходную функцию распределения F(m):

Модель (27), (28) несколько сглаживает исходную модель УГР, особенно заметно это сглаживание на концах диапазона магнитуд. Следует принимать во внимание, что использование значений порядка $\Delta ~\,\, \approx 0.5$ предполагает, что истинный неизвестный параметр М может отличаться от наблюдаемого на величину ${{\Delta }}$. Кроме того, если наличие ошибок при измерении магнитуд несомненно, то само понятие “истинной магнитуды” весьма неопределенно. Для полноты обсуждения требовалось бы сравнить качество подгонки простой и сглаженной моделей УГР к исследуемому каталогу. Мы не будем углубляться в дискуссию по этому вопросу, она выходит за рамки настоящей статьи, а приведем лишь плотность сглаженной модели и сравним ее с несглаженной моделью (см. рис. 6).

Мы видим, что сглаженная модель весьма заметно искажает усеченный закон Гутенберга–Рихтера, особенно на левом конце диапазона. Это несомненно повлияет на оценку параметра s (или β).

ОБСУЖДЕНИЕ И ВЫВОДЫ

В настоящей работе предложен новый способ оценки максимально возможной региональной магнитуды М, полученный в рамках распространенной модели усеченного распределения Гутенберга–Рихтера (УГР). Описан также метод оценки квантиля Q_T(q) максимального землетрясения в будущем интервале времени. Получена точная формула для смещения оценок максимальной магнитуды М и квантиля Q_T(q) в виде конечной эффективно вычисляемой суммы. Эти смещения выражаются через неполную Бета-функцию, один из аргументов которой равен –1. Соответствующая полная функция равна при таком аргументе бесконечности. Для неполной функции (она конечна) удалось получить точную формулу в виде эффективно вычисляемой суммы (12), (17).

Показано, что новая оценка по эффективности не уступает известным распространенным оценкам максимальной магнитуды М, а в ряде случаев заметно превосходит их. Среди всех известных методов, пожалуй, лишь Байесовский метод можно считать сравнимым по эффективности с новым методом (хотя в приведенном примере он немного уступает новому методу по величине MSE). Однако Байесовский метод содержит два “свободных параметра” алгоритма. Это параметры априорных интервалов для β и М. Если же в Байесовском методе используется возмущение магнитуды (что делается очень часто), то добавляется еще и третий параметр – масштаб возмущения. Эти 2–3 параметра часто выбирают без достаточного обсуждения и обоснования, исходя из индивидуальных и интуитивных соображений о том, какое значение параметра лучше соответствует представлениям автора об изучаемом явлении. В то же время новая оценка не зависит от каких-либо свободных параметров и процедура ее применения определена однозначно. В Байесовском методе при оценке максимальной магнитуды М очень большое значение имеет выбор априорного интервала для М. Брать ли в качестве конца априорного интервала (μ_n + 0.2), (μ_n + 1) или (μ_n + 1.5) – от этого может сильно зависеть результат оценки М. К сожалению, подобные детали обычно не обсуждаются и выбор априорных интервалов не обосновывается.

Байесовский подход занимает большое место в задачах оценки сейсмического риска (см. [Lyubushin, 2010; Pisarenko et al., 1996; Pisarenko, Lyubushin, 1999; Lyubushin et al., 2002; Kijko, 2012; Zentner, 2020; Писаренко и др., 2021]. Сделаем одно общее замечание относительно Байесовского подхода. При этом подходе неизвестные параметры (в частности, параметр М) рассматриваются как случайные величины, в то время как в исходной задаче они предполагались неизвестными величинами. Среди статистиков нет единого мнения об обоснованности такой замены. Применение случайных величин всегда подразумевает ансамбль реализаций (генеральную совокупность), в которой эти случайные величины реализуются. Для Байесовского подхода это не всегда легко сделать. Если говорить о максимальной магнитуде М, то исследователь, вообще говоря, интересуется максимально возможной магнитудой в данном конкретном регионе, а не популяцией возможных величин М (см. дискуссию по этому вопросу в работе [Kendall, Stuart, 1961] и в статье (Pisarenko, Rodkin, 2021)¹¹.

В одной из ранних работ автора [Pisarenko et al., 1996] введена несмещенная оценка максимальной магнитуды М, обладающая наименьшей дисперсией среди всех несмещенных оценок:

Эту оценку можно вывести из общих результатов Тейта [Tate, 1959] (cм. также [Kendall, Stuart, 1961]), но в этих работах она получается сложным путем с помощью интегральных уравнений, в то время как в работе [Pisarenko et al., 1996] дан простой наглядный вывод. Ценное качество оценки (29) состоит в том, что она имеет нулевое смещение, и ее следует использовать в тех задачах оценки сейсмического риска, где систематическое смещение важнее случайного разброса оценки. К сожалению, оценка (29) имеет довольно большую дисперсию, так что ее MSE больше, чем у многих других оценок. Выше мы использовали усеченный вариант оценки (29), он вполне конкурентоспособный, его MSE уменьшается (в зависимости от выбранного порога усечения), но, естественно, появляется систематическое смещение.

По мнению автора наиболее эффективными оценками параметра максимальной возможной магнитуды М являются два метода: введенная в настоящей работе новая оценка $\bar {M}$ (13), (17) и Байесовская оценка. Можно рекомендовать к применению обе эти оценки (при обоснованном выборе “свободных параметров” Байесовского метода, причем рекомендуется предварительно пробовать несколько вариантов Байесовских свободных параметров). Дополнительную, весьма полезную информацию о максимальных возможных толчках даст оценка квантилей Q_T(q), изложенная выше.

Разумеется, методы оценивания максимальных, возможных магнитуд не ограничиваются рамками УГР, в которых мы рассматривали эту проблему в настоящей работе. Для решения этой проблемы предлагаются и более сложные модели, учитывающие различные детали поведения хвоста распределения магнитуд (см. например, работу [Писаренко и др., 2020], в которой предложена двухчастная модель закона повторяемости землетрясений с использованием теории экстремальных значений [Gumbel, 1958; Embrechts et al.,1997; De Haan, 2006]). Более сложные модели требуют для детального описания хвоста распределений большего числа сильных событий, а это не всегда выполнимо в конкретных практических задачах. Так что модель УГР, которая как раз предназначена для работы с относительно небольшим количеством сильных событий, будет востребована еще долгое время.