Вулканология и сейсмология, 2023, № 2, стр. 3-11

ВВЕДЕНИЕ

Рассмотрим некоторый региональный или глобальный каталог землетрясений, магнитуды которого подчиняются распределению F(m). Относительно вида распределения F(m) в диапазоне самых больших магнитуд среди сейсмологов нет единого мнения, но в области малых и средних магнитуд общепринятым считается закон Гутенберга-Рихтера [Gutenberg, Richter, 1954], который после нормализации приводит к экспоненциальному распределению

Это распределение дает прямолинейный график повторяемости с наклоном b для функции lg(1 ‒ F(m)) и для логарифма плотности вероятности lg(f(m)). Отметим, что закон записан, используя десятичные логарифмы, но мы будем использовать натуральные из-за упрощения формул при дифференцировании и других математических операциях. Однако, поскольку приведенные выше рассуждения включают в себя неопределенности типа “в диапазоне малых и средних значений”, “в диапазоне самых больших магнитуд” и т.п., их трудно использовать на практике. В настоящей заметке мы будем исследовать вопросы, связанные с выделением прямолинейного участка графика повторяемости землетрясений и с оценкой параметра b на выделенном участке. Параметр b важен как физическая характеристика самоподобия сейсмического процесса (см. например, [Ризниченко, 1966; Смирнов и др., 1995; Соболев, 1993; Scholz, 1968]): на одно землетрясение магнитуды m в среднем приходится exp(b) землетрясений магнитуды (m ‒ 1), независимо от величины магнитуды m. В задачах сейсмического риска параметр b играет важную роль при оценивании максимальной региональной магнитуды M_max.

Будем считать, что закон Гутенберга-Рихтера справедлив на некотором конечном интервале [m₀; m₁], границы которого являются неизвестными, а вне этого интервала он может нарушаться. В диапазоне самых малых магнитуд он может нарушаться из-за неполной регистрации слабых толчков (особенно в ранние времена инструментальной сейсмологии и для исторической регистрации землетрясений), а в диапазоне самых больших толчков из-за региональных особенностей сейсмического режима и геологического строения данного региона. Таким образом, мы предполагаем, что в интервале магнитуд [m₀; m₁] справедлив усеченный закон Гутенберга-Рихтера:

Функция распределения:

Плотность распределения:

Усеченный закон Гутенберга-Рихтера (УГР) широко используется в современной сейсмологической практике из-за своей простой формы и робастности. Часто его используют для оценки максимально возможной региональной магнитуды M_max (см. например, [Cosentino et al., 1977; Holschneider et al., 2011; Kagan, Schoenberg, 2001; Kijko 2004; Kijko, Graham, 1998; Kijko, Sellevoll, 1989, 1992; Kijko, Singh, 2011]).

В настоящей работе рассмотрена задача выделения наибольшего возможного интервала [m₀; m₁], на котором УГР соблюдается с достаточной точностью, а вне этого интервала он нарушается. Параметр m₁ не следует путать с параметром M_max ‒ максимальной региональной магнитудой. Интервал [m₀; m₁] – это максимальный интервал, на котором соблюдается линейность закона повторяемости землетрясений. Как правило, m₁ существенно меньше, чем M_max. Нижняя граница интервала прямолинейности m₀ может меняться со временем и, как правило, уменьшается по мере совершенствования сейсмических сетей и методов регистрации землетрясений. Ниже мы уделим специальное внимание этому вопросу. Вопрос об учете цензурирования распределения землетрясений по магнитудам сверху при оценке наклона графика повторяемости имеет давнюю историю. По-видимому, этот эффект “открывался” разными исследователями несколько раз в период 1966‒2011 гг. (см. библиографию на этот счет в работе [Смирнов, Завьялов, 2012]. “Поправка” за учет цензурирования сверху используется, в частности, при анализе переходных сейсмических режимов, а также и в лабораторных исследованиях [Смирнов, Пономарев, 2020].

ОЦЕНКА ПАРАМЕТРА b ПРИ ИЗВЕСТНОМ ИНТЕРВАЛЕ ЛИНЕЙНОСТИ

Для выделенного интервала [m₀; m₁] оценка параметра b проводится методом максимального правдоподобия. Логарифмическая функция правдоподобия L(b, m₀, m₁| $\bar {x}$) равна

где $\bar {x}$ = (x₁, …, x_n) – выборка магнитуд, подвергаемая анализу. При известных m₀, m₁ оценка максимального правдоподобия находится как значение b, максимизирующее правдоподобие (3). Это делается с помощью стандартных алгоритмов поиска экстремума функции и обычно не вызывает затруднений, так как зависимость функции правдоподобия (3) от b имеет регулярный характер.

Изложенный выше стандартный способ оценки максимального правдоподобия для параметра b, строго говоря, справедлив для непрерывных распределений. Реальные каталоги магнитуд всегда даются с некоторой дискретизацией. Как будет показано ниже, различия в оценках параметров для непрерывного и дискретного подходов существенны, когда интервал дискретизации больше 0.05‒0.10 магнитуды и практически исчезают при дискретизации с шагом 0.01 и менее. В настоящее время встречаются каталоги с шагом дискретизации 0.01 магнитуды, например, каталог моментных магнитуд GCMT. Однако, во многих региональные каталогах (в частности, в исторических региональных каталогах) шаг дискретизации больше, иногда его значения доходят до 0.5. В этих случаях необходимо учитывать те особенности, которые вносит дискретизация в процесс статистического оценивания параметров.

Для неограниченного закона Гутенберга-Рихтера (экспоненциальное распределение магнитуд) эффективной оценкой величины b является оценка максимального правдоподобия $\hat {b}$

Эта оценка была предложена впервые в работе [Aki, 1965] для непрерывных распределений. Для дискретных распределений с шагом дискретизации Δ в формулу (4) предлагается вносить поправку

см. например, [Utsu 1966; Bender, 1983; Taroni et al., 2021].

Мы предлагаем для дискретных распределений при оценивании параметра b использовать классическую оценку максимального правдоподобия для дискретного распределения [Cramer, 1946]. Такие оценки применялись в работах [Pisarenko, Rodkin, 2010; Pisarenko et al., 2010; Писаренко и др., 2022]. Как и все оценки максимального правдоподобия в регулярном случае они являются асимптотически наиболее эффективными.

Пусть магнитуды каталога с непрерывной функцией распределения F(m |b, m₁, m₀) дискретизированы с шагом Δ. Обозначим границы интервалов дискретизации магнитуд $m_{j}^{d}$:

Ради простоты изложения мы будем считать, что на интервале [m₀; m₁] укладывается целое число r интервалов Δ.

Вычисляем дискретные вероятности p(j | b)

Пусть выборка магнитуд (х₁, …, х_n) распределилась по ячейкам (6) с числами $\bar {n}$ = (n⁽¹⁾, …, n^(r)). Логарифм функции правдоподобия имеет вид:

Также как и в непрерывном случае, оценка максимального правдоподобия для параметра b находится как значение, максимизирующее правдоподобие (8). Обозначим эту оценку b_d.

Сравним оценки b_Δ и b_d на синтетических каталогах с известным ответом. В табл. 1 приведены стандартные отклонения (STD), смещения (BIAS) и среднеквадратичные отклонения (Mean Square Error, MSE = $\sqrt {~ST{{D}^{2}}~ + BIA{{S}^{2}}} $) для синтетических каталогов, имеющих распределение (1). Параметры каталогов указаны в примечании к табл. 1, для усреднения бралось 10 000 каталогов. Для сравнения в 4-м и 7-м столбцах приведены оценки параметра $\hat {b}$, полученные при максимизации правдоподобия непрерывно распределенных магнитуд, из которых получены дискретные магнитуды (эти значения, конечно, неизвестны исследователю, но мы приводим их для сравнения с реальными оценками).

Из табл. 1 можно сделать вывод, что оценка b_d, основанная на правдоподобии дискретных магнитуд, практически не уступает по MSE идеальной (недоступной) оценке $\hat {b}$, использующей непрерывное распределение магнитуд. Для малых интервалов прямолинейности (m₁ – m₀) < 1.5 наблюдается явное превосходство по MSE оценок b_d над оценками b_Δ. Оно достигается, в основном, за счет большого систематического смещения последних. В табл. 1 красным цветом выделены те случаи, в которых MSE оценки b_Δ превышает MSE оценки b_d более, чем в 1.75 раза. В целом, можно сделать вывод об эффективности дискретных оценок правдоподобия, которые, к сожалению, не имеют столь широкого распространения как оценки с поправками b_Δ. Впрочем, следует заметить, что при очень малых выборках (n < 20‒30) оценки b_Δ могут оказаться несколько эффективнее, поскольку оценки правдоподобия (в том числе и для непрерывного случая) при малых n иногда уступают по эффективности другим оценкам [Pisarenko, Rodkin, 2017].

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

Будем считать, что распределение выборки землетрясений (каталога) подчиняется закону (1), где m₀, m₁ неизвестны, и нам нужно по возможности точнее указать такой интервал [m₀; m₁], внутри которого закон Г-Р соблюдается, а вне которого начинаются отклонения от него. Методику определения концов интервала линейности мы будем иллюстрировать на примере глобального каталога GСМТ за период 1976‒2022. Магнитуда m вычислялась по известной формуле Канамори

где M₀ ‒ скалярный сейсмический момент в единицах дин см. Поскольку lg(M₀) задается в каталоге с точностью до 0.01, можно считать, что полученные значения магнитуд могут характеризоваться такой же точностью дискретизации. Каталог был предварительно декластеризован с помощью метода, описанного в работе [Писаренко, Родкин, 2019]. График повторяемости декластеризованного каталога показан на рис. 1.

Обозначим минимальную наблюденную магнитуду каталога $m_{{{\text{min}}}}^{{{\text{obs}}}}$, а максимальную наблюденную магнитуду $m_{{{\text{max}}}}^{{{\text{obs}}}}$. Для рассматриваемого каталога $m_{{{\text{min}}}}^{{{\text{obs}}}}$ = 4.30, $m_{{{\text{max}}}}^{{{\text{obs}}}}$ = 9.08. Учитывая поведение графика повторяемости, выберем с запасом предварительный интервал линейности закона повторяемости [m; $\bar {m}$]: m = 6.0; $\bar {m}$ = 7.2. Такой выбор можно рассматривать как экспертное решение. Затем мы уточним предварительный интервал и получим окончательный интервал линейности, по которому будет получена финальная оценка максимального правдоподобия параметра b. Оценку, полученную по данным на интервале [m; $\bar {m}$], обозначим b₀ = 2.26.

Рассмотрим оценку максимального правдоподобия b_left(m) на изменяющемся интервале [m; $\bar {m}$] для $m_{{{\text{min}}}}^{{{\text{obs}}}}$ < m < $\bar {m}$. Попробуем проследить, как b_left(m) приближается к значению b₀ = 2.26 при возрастании m. Предполагая, что оценка b₀= 2.26 не сильно отличается от истинного наклона графика повторяемости, можно приближенно считать, что удвоенное логарифмическое отношение правдоподобия

имеет асимптотически χ²-распределение с одной степенью свободы (см. Линейные статистические методы и их применения. Глава 6 [Рао, 1965]).

Используя формулу (9), можно построить график функции P_left(m)

позволяющей судить о том, можно ли считать данное значение m принадлежащим интервалу линейности или нет. Для интервала линейности величины P_left(m) не должны быть малы, а когда линейность начинает нарушаться, вероятность P_left(m) уменьшается до значений 0.1 и менее. Таким образом, надо выбрать порог, за которым вероятность P_left(m) становится близкой к единице, а оценка b_left(m) достигает значения 2.26, принимаемого за истинное.

На рис. 2 показаны графики оценки b_left(m) и вероятности P_left(m), а также указана оценка b₀ = 2.26.

Из рис. 2 видно, что точка окончательного выхода Р_left(m) и b_left(m) из зоны малых значений находится вблизи магнитуды 5.72. Это подтверждается графиком функции P_left(m), которая достигает в этой точке максимального значения 1, что позволяет нам выбрать в качестве оптимального для левого конца значение m₀ = 5.72.

Итак, в результате нашего анализа мы получили для левого конца интервала линейности графика повторяемости значение m₀ = 5.72.

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

Рассмотрим убывающий с ростом m интервал [m; $m_{{{\text{max}}}}^{{{\text{obs}}}}$]. Оценку b, полученную по данным на интервале [m; $m_{{{\text{max}}}}^{{{\text{obs}}}}$] обозначим b_right(m), а в качестве истинного оставляем значение b₀ = 2.26, являющееся оценкой максимального правдоподобия на интервале [6.0; 7.2]. При сделанных предположениях удвоенное логарифмическое отношение правдоподобия

имеет асимптотически χ²-распределение с одной степенью свободы.

Используя формулу (11), можно построить график функции P_right(m)

позволяющей судить о том, можно ли считать данное значение m принадлежащим интервалу линейности или нет. Когда значение m выходит за интервал линейности, P_right(m) начинает убывать, а соответствующая оценка b_right(m) начинает отличаться от значения 2.26, принимаемого за истинное. На рис. 3 показан график оценки b_right(m), а также график вероятности P_right(m). Надо определить ту точку m₁, после которой функция P_right(m) начинает убывать и уже окончательно опускается ниже 0.1.

Из рис. 3 видно, что перелом поведения функции P_right(m) начинается вблизи магнитуды 7.25. На основании вышесказанного в качестве окончательного интервала прямолинейности графика повторяемости мы принимаем [5.72; 7.25]. Оценка максимального правдоподобия для параметра наклона b на этом интервале равна 2.19. Мы видим, что предварительная оценка максимального правдоподобия, полученная по интервалу [6.0; 7.2], отличается от окончательной оценки, полученной по интервалу [5.72; 7.25], на 0.07, что свидетельствует об устойчивости предложенной процедуры определения интервала прямолинейности и оценивания параметра b. Интересно сравнить оценку 2.19, полученную на оптимальном интервале магнитуд с оценкой максимального правдоподобия на интервале, который можно было бы выбрать просто из рассмотрения графика повторяемости на рис. 1. В качестве такого интервала можно взять, например, [4.9; 7.25], поскольку на этом интервале линейность графика повторяемости кажется на вид приемлемой. Оценка максимального правдоподобия на интервале [4.9; 7.25], равная 1.79, отличается от оценки на оптимальном интервале на 0.40, так что правильный выбор интервала для оценки параметра b может иметь существенное значение.

Для оценки разброса оценки параметра b можно было бы воспользоваться асимптотическими формулами метода максимального правдоподобия [Cramer, 1946]. Согласно этому методу, величина

имеет в пределе при n $ \to \infty $ стандартное нормальное распределение. При этом константа с находится из равенства:

Из уравнения (14) видно, что константа с зависит от неизвестного параметра b. Чтобы воспользоваться уравнением (13) для практических целей обычно подставляют в уравнение (13) вместо b оценку максимального правдоподобия, что для больших выборок вносит лишь небольшие изменения. Мы предлагаем другой подход. Используя тот известный факт, что любая непрерывная случайная величина, подставленная в свою функцию распределения, дает случайную величину, равномерно распределенную на отрезке [0 1], можно смоделировать случайные каталоги с заданной функцией распределения (см. подробнее [Pisarenko et al., 2021]). В качестве такого распределения возьмем усеченный закон Гутенберга-Рихтера (1) с интервалом линейности [5.72; 7.25] и параметром b = 2.19. Произведем большую выборку (скажем 10 000) синтетических каталогов с указанными параметрами и того же размера n = 10 125, что и наш каталог в данном интервале магнитуд. По такой большой выборке имитируемых каталогов можно с достаточной точностью оценить стандартное отклонение оценок максимального правдоподобия параметра b. Оно оказалось равным 0.03:

ОПРЕДЕЛЕНИЕ ИЗМЕНЕНИЙ ВО ВРЕМЕНИ НИЖНЕГО ПРЕДЕЛА ПРЕДСТАВИТЕЛЬНОЙ РЕГИСТРАЦИИ ЗЕМЛЕТРЯСЕНИЙ

При использовании инструментальных каталогов, покрывающих большие интервалы времени, а также исторических каталогов, нижний предел представительной регистрации землетрясений изменяется, как правило, понижаясь со временем в основном за счет совершенствования системы регистрации землетрясений. В ряде работ такие изменения отслеживают и в соответствии с ними вносят поправки в оценки параметров закона повторяемости (см. [Taroni et al., 2021]). По большей части такие изменения учитываются сейсмологами субъективно, с учетом имеющегося набора наблюдений. Мы попробуем, насколько это возможно, формализовать этот процесс, используя некоторые статистические методы.

На рис. 4 показан каталог ISC-GEM за период 1904‒2014 гг. Мы будем исследовать выборочные квантили этого каталога в скользящем временном окне.

Из рис. 4 можно видеть, что для магнитуд 5.0‒5.59 наблюдаются нерегулярные пропуски регистрации. Поэтому мы возьмем для анализа магнитуды m ≥ 5.6. В качестве характеристики интервала представительной регистрации в каждом скользящем временном окне возьмем полуинтервал наибольших величин магнитуды [Q_t ∞), имеющий вероятность q. Значение Q_t можно рассматривать как нижнюю границу интервала представительности, имеющую уровень доверия q:

Нижняя граница Q_t теоретически является квантилем уровня (1 – q) случайной магнитуды m, распределение которой зависит от времени t. Практически для его оценки берется выборочный квантиль (процентиль уровня (1 ‒ q) ) в скользящем временном окне [(t – Δ_t/2), (t + Δ_t/2)] длительности Δ_t. Выбирая q достаточно большим, скажем q = 80, 90, 95%, можно описать изменения нижней границы представительности со временем. Брать слишком близкие к единице значения q не рекомендуется, поскольку их статистические оценки становятся при этом неустойчивыми (это особенно явно проявляется при малом объеме исследуемого каталога). Наиболее подходящим для практической характеристики нижней границы интервала представительности, на наш взгляд, является уровень q = 90%. Следует отметить, что при реальной оценке квантилей в плавающем окне мы использовали для уменьшения случайных флуктуаций процедуру предварительного сглаживания по времени, которая состояла в том, что перед определением квантилей множество моментов времени (τ₁, …, τ_k) событий подвергалось k раз случайным возмущениям с равномерным распределением в интервале [–δ, +δ] и полученные по возмущенным выборкам квантили усреднялись. В нашем конкретном случае мы брали δ = 5 лет, но эту величину следует подбирать, учитывая особенности анализируемого каталога и его объем. Число усредняемых значений достаточно брать не очень большим, мы брали k = 1000. На рис. 5 в скользящем окне [(t – 1.8), (t + 1.8)], Δ_t = 3.6 (время в годах) показаны 3 выборочных сглаженных квантиля Q_t(0.20), Q_t(0.10), Q_t(0.05), соответствующие уровням представительности q = 80, 90, 95%.

Из рис. 5 видно, что после 1960‒1965 гг. значения квантилей стабилизируются. В более ранние годы они существенно больше современных. Повышение нижней границы интервала представительности в 1941‒1948 гг., по-видимому, обусловлено ухудшением системы регистрации землетрясений во время второй мировой войны 1939‒1945 гг. Изменения уровня представительности со временем необходимо учитывать при оценивании параметра b.

ЗАКЛЮЧЕНИЕ

Предложена статистически обоснованная процедура выбора интервала линейности для закона Гутенберга-Рихтера и оценивания параметра b на этом интервале (непрерывный и дискретный варианты). Как показывают примеры анализа каталогов, этот интервал [m₀; m₁] заметно меньше всего интервала регистрации землетрясений. Это следует учитывать при оценке параметра b. Включение в интервал оценки магнитуд m < m₀ приводит к занижению оценки b, а включение магнитуд m > m₁, как правило, к завышению (но возможны и случаи занижения для некоторых региональных каталогов). Следует отметить, что использование магнитуд в нижнем интервале m₀ ‒ а < < m < m₀ сказывается на оценке максимального правдоподобия сильнее, чем магнитуд в таком же по длине верхнем интервале m₁ < m < m₁+ а, так как оно сопровождается добавлением в правдоподобие гораздо большего количества наблюдений из-за экспоненциального убывания плотности. Поэтому для оценки параметра b точность оценки левого конца m₀ гораздо существенней точности оценки правого конца m₁.

Для каталогов с дискретными магнитудами вместо обычно используемых поправок на дискретность при оценивании параметра b предлагается использовать классический метод максимального правдоподобия для дискретных распределений. Показано, что для не слишком малых интервалов магнитуд этот метод дает меньшую средне-квадратичную ошибку. Кроме того, этот метод применим не только к равномерной дискретизации, но и к произвольным дискретным распределениям (см. [Писаренко и др., 2022]), энергетическим классам, агрегированным интервалам магнитуд и т.п.

Предложена статистически обоснованная нижняя граница Q_t(q) интервала представительной регистрации, имеющая уровень доверия q. Она является выборочным квантилем уровня (1 – q) случайной магнитуды m в скользящем временном окне [(t – Δ_t/2), (t + Δ_t/2)]. На примере глобального каталога ISC-GEM 1904‒2014 показано, как эта граница изменялась со временем.