Причиной написания этой статьи стала распространенная рядом СМИ информация о прогнозе «мощнейшего» землетрясения, которое может произойти в ближайшие 30 лет в Японии и на Курилах с вероятностью до 40%. Ссылались журналисты на японских ученых. Нам удалось найти оригинал статьи, откуда, по всей видимости, и взяты сведения. Она была опубликована 20 декабря в Japan News и сейчас находится в платном архиве, но у нас есть замечательный ресурс wayback machine.

Ниже мы постараемся разобрать, о чем же действительно первоисточник, погружаясь в дебри математики и основы вероятностного анализа сейсмической опасности. Если в двух словах, то японские сейсмологи дали не прогноз мега-землетрясения на Курилах, а описали модельные характеристики сейсмических источников, которые учитываются при составлении карт сейсмического районирования со сроком службы зданий и сооружений на ближайшие 30 лет. Попытаемся произвести вычисления. Без суперкомпьютеров…

Итак, главная задача сейсмического районирования – определение наиболее вероятных сейсмических воздействий. Термин «наиболее» мы ниже расшифруем. Для начала запишем простую математику и введем некоторые термины, понятные для большинства читателей даже более-менее знакомых с теорией вероятностей и математической статистикой. Для определенности сейсмические воздействия будем измерять в физических параметрах ускорения. Скажем, пиковое ускорение грунта 0.7g (где g=9.81 м/с^2 – ускорение свободного падения) соответствует интенсивности сотрясений около IX баллов, а IX баллов – это катастрофа.

У каждого землетрясения есть своя магнитуда M. Это величина, которая характеризует силу (размер разрыва, энергию и т.д.) в очаге. Не путать с баллом! Балл – это уже проявление сотрясаемости на поверхности, то есть ощущаемость. Чем ближе к эпицентру, тем интенсивнее колеблется поверхность Земли, тем выше уровень сейсмических воздействий, а удаляясь от эпицентра, сейсмические сигналы затухают. Это фундаментальное свойство, присущее всем сигналам в сплошной среде. Назовем его законом затухания сейсмических воздействий. Очевидно, что разрушает здание не магнитуда, а пиковые ускорения грунта, которые превышают некоторый предельный уровень, заложенный при проектировании и строительстве (рис. 1).

Рисунок 1. Последствия Невельского землетрясения 2007 года

Модель затухания пикового ускорения грунта (peak ground acceleration) задается посредством функции g(m,r), которая определяет зависимость среднего значения (натурального) логарифма пикового ускорения грунта ln⁡ PGA от события с магнитудой m и на расстоянии r. Данная функция представляется регрессионным соотношением (рис. 2), построенным на основе региональной базы акселерограмм, полученных при помощи сети акселерографов. Обычно она имеет вид:

$g(m,r) = ln\overline{PGA}(m,r)=c_1+c_2m+c_3ln⁡(r+c_4 )-c_5 r+c_6 F+c_7 S+σ, (1)$

где c_1, c_2, c_3, c_4, c_5, c_6, c_7 – регрессионные коэффициенты, а F и S описывают соответственно зависимость от типа разлома (Fault) и характеристик грунта (Soil).

Рисунок 2. Кривая затухания пиковых ускорений грунта с расстоянием для землетрясения с магнитудой M=5.8: 1 – измеренные значения, 2 – региональная регрессия, 3 – доверительный интервал ±σ, 4 – плотностная функция стандартного нормального распределения натурального логарифма пикового ускорения на расстоянии 10 км от очага, 5 – уровень пиковых ускорений 0.7g.

Изменчивость (от события к событию и от пункта к пункту), наблюдаемая в данных по сильным движениям грунта, описывается нормальным распределением величины ln ⁡PGA(m,r) в каждой точке (m,r) посредством нулевого среднего и стандартной ошибки σ (рис. 2). К этому все же стоит относиться как к допущению, но вполне рабочему [Cornell, 1968].

Тогда из сделанных выше предположений относительно логнормального распределения пикового ускорения грунта следует, что из всех событий магнитуды m, происходящих на расстоянии r от рассматриваемого пункта, доля тех событий, что не вызовет сейсмические воздействия, превышающие ускорение a, есть

$Ф(\frac{lna-g(m,r)}\sigma), (2)$

где Ф – функция распределения стандартной нормальной величины. В приложении дана таблица с численными значениями нормального распределения. Мы ее позаимствовали из [Baker, 2008].

Таким образом, вероятность непревышения заданного уровня сейсмических воздействий a от сейсмического события с магнитудой m на расстоянии r будет полностью задаваться функцией распределения стандартной нормальной величины. Обозначим эту вероятность как Pr(A≤a│m,r). Тогда в соответствии с определением (2) можно записать:

$Pr(A≤a│m,r)= \int\limits_{-\infty}^a\frac{1}{\sigma\sqrt{2\pi}}\exp(-\frac{1}{2}(\frac{lnx-g(m,r)}{\sigma})^2)dx. (3)$

Графическая интерпретация формул (1)-(3) показана на рис. 2. Видно, что регрессия описывает измеренные значения пиковых ускорений грунта. Все точки относятся к одному землетрясению. Точнее говоря, сама регрессия построена по большому набору эмпирических данных разных землетрясений, но для простоты изложения приведены графики для землетрясения с магнитудой M=5.8. Скажем, вблизи Невельска (рис. 1) произошла серия землетрясений по магнитуде близких к M~6, вызвавшая разрушения в городе.

На рис. 2 виден разброс значений пикового ускорения, который определяет стандартную ошибку σ=0.783 из (1). Натуральный логарифм пикового ускорения грунта для точек m=5.8 и r=10 км равен g(m,r)=4.868 (PGA=125 см/сек^2). Сейсмические воздействия с пиковым ускорением 0.7g (687 см/сек^2) соответствую ln⁡ a=6.532. Тогда вероятность того, что сейсмические воздействия A, вызванные землетрясением с магнитудой m=5.8 на расстоянии r=10 км, не превысят 0.7g (IX баллов) равна (см. Приложение):

$Pr(A≤0.7g│m=5.8,r=10)=Ф(\frac{6.532-4.868}{0.783})=0.983$

Здания не строятся на бесконечное время эксплуатации. У всего есть свой срок службы, и он определяется временем T.

Рисунок 3. Схематическая иллюстрация сейсмических источников.

В каком же случае воздействия от землетрясения с магнитудой M_i, очаг которого располагается в точке k (рис. 3), не превысят в заданном пункте величину ускорения a в течение последующего промежутка времени T? Первый ответ на такой длинный вопрос очевиден: если вообще не произойдет ни одного землетрясения за рассматриваемый временной интервал. Второй случай, менее очевиден: если произойдет одно землетрясение, но по определенному в (1) закону затухания оно не вызовет ожидаемых сейсмических воздействий. Вообще за время T может произойти хоть два землетрясения или три или вообще Ns. Перефразируем сказанное в терминах вероятности.

Примем упрощенную модель, в которой некоторая точка генерирует землетрясения только одной магнитуды. Посмотрим на табл. 1. В ней показаны различные исходы, при которых сейсмические воздействия не будут превышены: ни одного сейсмического события (0), одно (1), два (1 + 1), три (1 + 1 + 1) и т.д. Вероятность возникновения одного события в течение времени T равна p1, двух – p2, трех – p3 и т.д., ни одного – соответственно p0. Естественно, рассматриваемое множество исходов является полным, т.е. других вариантов нет:

$p0+p1+p2+p3+\dots=1. (4)$

Вероятность того, что колебания от землетрясения заданной магнитуды не превысят заданное сейсмической воздействие в заданном пункте, обозначим как Pr, которая полностью задается через (3). Тогда вероятность того, что колебания от двух землетрясений не превысят заданное сейсмической воздействие, равна Pr*Pr, трех – Pr*Pr*Pr и т.д. При отсутствии землетрясений, очевидно, что вероятность непревышения равна 1. Таким образом, полная вероятность P того, что в течение времени T сейсмические воздействия в заданном пункте не будут превышены:

$P=p0\ast(1)+p1 \ast (Pr)+ p2 \ast (Pr \ast Pr)+p3 \ast (Pr \ast Pr \ast Pr)+\dots. (5)$

Таблица 1. Вероятностные модели источника землетрясений заданной магнитуды в заданной точке.

Исходы(количество событий)	Вероятность возникновения	Вероятность непревышения
0	p0	1
1	p1	Pr
1+1	p2	Pr*Pr
1+1+1	p3	PrPrPr
...	...	...

После того как выражение (5) переварено окончательно, можно перейти к рассмотрению общего случая. Для этого нам необходимо задать вероятность того, что в течение последующего времени T произойдет s событий, где s=0,1,2,…,Ns. Обозначим вероятность возникновения s событий с магнитудой M_i в некоторой точке k в последующее время T как P_k(s,M_i,T).

Тогда вероятность непревышения заданного уровня сейсмических воздействий a от сейсмического события с магнитудой M_i в точке k (рис. 3) в последующий промежуток времени T будет задаваться по аналогии с (5):

$P(A \leq a \mid M_i,T,k)=\sum_{s=0}^{N_s}P_k(s,M_i,T)[Pr(A\leq a \mid M_i,R_k)]^S, (6)$

где R_k – расстояние от очага землетрясения в точке k до пункта, в котором мы ожидаем сейсмические воздействия.

Дальше – проще. Необходимо рассмотреть независимые реализации магнитуд M_i. Конечно, в большой сейсмологии есть исключения. Например, когда сильное землетрясение вызывает серию афтершоков (или другими словами повторные события с магнитудой чуть меньше, чем у главного события). Но для этого производят так называемую декластеризацию каталога землетрясений, т.е. удаление афтершоков и других зависимых событий.

С учетом условия независимости магнитуд из (6) получим вероятность непревышения заданного уровня сейсмических воздействий a от серии землетрясений в точке k в последующий промежуток времени T:

$P(A≤a│T,k)=\prod_{i=1}^{N_m}Pr(A≤a│M_i,T,k). (7)$

Разумно предположить, что сейсмические источники независимы, и каждый из них живет своей жизнью. Опять-таки, в большой сейсмологии встречаются случаи, когда тектоническая активность одного разлома вызывает активизацию сейсмичности на другом. Тем не менее, из предположения независимости набора N сейсмических источников получим вероятность непревышения заданного уровня сейсмических воздействий a в нашем пункте «ожидания» в последующий промежуток времени T:

$P(A≤a│T)=\prod_{i=1}^{N}\prod_{i=1}^{N_m}Pr(A≤a│T,k). (8)$

Перейдем от вероятности непревышения в (8) к вероятности превышения:

$display$

Выражение (9) является базовым для производства карт сейсмического районирования в вероятностном анализе сейсмической опасности.

Разберем вкратце вероятность возникновения землетрясения P_k (s,M_i ). В сейсмическом районировании часто используется экспоненциальная Пуассоновская модель, то есть вероятность того, что в последующие T лет произойдет s землетрясений равна:

$P_k (s,M_i,T)=\frac{[λ_k (M_i )T]^s exp[-λ_k (M_i )T]}{s!}, (10)$

где λ_k (M_i ) – частота возникновения землетрясений с магнитудой M_i в точке k. Как и положено сумма (10) по всем s от нуля до бесконечности равна единице!

Японские ученые по данным о палеоцунами сделали реконструкция исторических землетрясений. Как известно, цунами вызываются довольно сильными сейсмическими событиями, а их заплески сохраняются в геологической и морфологической истории побережья. По оценкам японцев в районе южных Курильских островов магнитуда максимально возможного землетрясения может быть не менее M=8.8. Как они считают, рассматриваемые события происходят не случайно, а периодично. Межсейсмический интервал лежит в диапазоне 340-380 лет, т.е. средний период повторяемости составляет 360 лет. Тогда частота возникновения землетрясений λ будет оцениваться как 1/360 год^-1.

Давайте снова рассмотрим упрощенную модель сейсмического источника – в заданной точке заданной магнитуды (рис. 3). Нормативный интервал времени примем T=30 лет. Тогда очевидно, что λT << 1. Это означает, что в формуле (10) можно ограничиться «первыми» вероятностями p0 и p1:

$display$

Из (11) следует, что в каждый момент времени возникновение землетрясения равновероятно. Например, по этой модели нашумевшее «мощнейшее» землетрясение возникнет в последующие 30 лет с вероятностью 0.083.

Вероятность непревышения заданных сейсмических воздействий по формуле (5) будет оцениваться как

$display$

Искомая вероятность превышения согласно (9) и (12) задается как

$display$

Раскладывая (13) по степеням λT и оставляя только первые члены разложения, получим:

$display$

Таким образом, при некоторых допущениях мы получили простую формулу для оценки нормативных сейсмических воздействий вероятностным методом. Осталась последняя деталь. Это закон затухания сильных движений грунта для Курильских островов.

Для Курильских островов близким по геологическим и тектоническим условиям является соседний регион Японии. Здесь, как и в Калифорнии, на современном уровне разработаны уравнения затухания пиковых ускорений грунта. Наиболее унифицированной моделью, учитывающей разделение сейсмичности на субдукционный и коровый тип, является базовая модель Ши и Мидорикава [Si, Midorikawa, 1999]. Более поздние модификации этой модели во многом связаны с уточнением коэффициентов после детального анализа сейсмических воздействий мега-землетрясения Тохоку 2011 г. (Mw 9.0).

$\begin{cases} log⁡A=b-log⁡(X+c)-kX,D≤30 км \\ log⁡A=b+0.6 log⁡(1.7D+c)-1.6 log⁡(X+c)-kX,D>30 км, (15) \end{cases}$

где

$b= αMw + hD + \sum d_i S_i + e, (16)$

A – пиковое ускорение в см/с^2; D – «средняя» глубина плоскости разрыва (глубина центроида) в км; X – наикротчайшее расстояние до разрыва в км; Mw – моментная магнитуда;

+ другие коэффициенты

$display$

$c=0.006*10^{0.5Mw};$

$display$

$\sigma_{log⁡PGA} =0.27.$

Для коровых землетрясений S_1=1, S_2=0, S_3=0; для межплитовых землетрясений S_1=0, S_2=1, S_3=0; для внутриплитовых S_1=0, S_2=0, S_3=1.

Рисунок 4. Модель сейсмического источника на южных Курилах (источник)

«Наше» мега-землетрясение ожидается на стыке двух плит, т.е. тип землетрясения – межплитовый. Японское мега-землетрясение Тохоку 2011 г., известное по трагедии на атомной станции Фукусима, произошло на глубине около 11 км. Поэтому и в нашем случае примем глубину очага D=11 км, а магнитуду Mw=8.8 (как предполагают японцы). На рис. 4 изображена пространственная модель ожидаемого на Курилах мега-землетрясения. Так, например для о. Шикотан наикротчайшее расстояние от погружающегося в сторону Курильских островов разлома до центра острова будет около 40 км, поэтому примем X=40 км. Теперь, когда все входные параметры определены, определим вероятность Pr(A≤a│m=8.8,r=40) из (14) согласно закону затухания (15). Для этого надо немного повозиться, выражение (15) привести к натуральному логарифму и т.д. В результате получим: g(m,r)=7.22, σ_ln⁡ PGA =0.62. Ниже, в табл. 2, представлены посчитанные значения вероятности.

Таблица 2. Вероятности превышения заданных сейсмических воздействий

Ускорение `a, g`	Ускорение `a, см/сек^2`	Ln a	Вероятность `Pr(A≤a│m,r)`	Вероятность `Q(A>a│T)`
0.30	294.30	5.68	0.99	0.083
0.40	392.40	5.97	0.98	0.081
0.50	490.50	6.20	0.95	0.079
0.60	588.60	6.38	0.91	0.076
0.70	686.70	6.53	0.87	0.072
0.80	784.80	6.67	0.82	0.068
0.90	882.90	6.78	0.76	0.063
1.00	981.00	6.89	0.71	0.059
1.20	1177.20	7.07	0.60	0.050
1.40	1373.40	7.23	0.50	0.042

Целью такого анализа является оценка вероятности превышения уровня движения грунта, а главный результат состоит в определении зависимости вероятности превышения от уровня движения, которая называется кривой опасности. Кривая опасности наглядно приведена на рис. 5 согласно данным из табл. 2. Из рис. 5 мы видим, что вероятность превышения 0.06 в течение ближайших 30 лет соответствует уровню сейсмических воздействий около 0.9g. Это больше 9 баллов. Такая оценка для Курильских островов является волне естественной и ожидаемой, если вспомнить, что здесь в 1994 году произошло крупное землетрясение, именуемое Шикотанским. Оно сопровождалось разрушениями сооружений, гибелью людей, волнами цунами и многочисленными оползнями. Интенсивность сотрясений на о. Шикотан составила от VIII до IX баллов по шкале MSK-64.

Шикотанское землетрясение

Источник

Рисунок 5. Кривые опасности для о. Шикотан по модели сейсмического источника на южных Курилах по данным этого источника.

Кроме Пуассоновской модели существуют модели возникновения землетрясений с «памятью», которые учитывают историю предшествующей сейсмичности. К ним относится Броуновская модель (Brownian Passage Time). Ее применяют японские коллеги при построении карт сейсмического районирования для хорошо изученных территорий. По их оценкам, вероятность возникновения мега-землетрясения в ближайшие 30 лет на южных Курилах составила от 0.07 до 0.4 (рис. 4).

Представьте себе, что максимальное землетрясение происходит строго периодично – один раз в 1000 лет. И вот оно произошло, скажем, 20 лет назад, а мы собираемся в этом месте построить жилые здания, которые прослужат следующие 50 лет. За это время сильные события уже не произойдут, так как были только что, а следующее повторится не раньше, чем через 980 лет. Если такой факт периодичности установлен, то используется Броуновская модель. Тогда мы получим карты сейсмического районирования с более низкими значениями расчетной интенсивности сейсмических воздействий, и это на выходе удешевит строительство. Какую модель выбирать, решает не один специалист, а группа экспертов.

Для расчета кривой опасности по Броуновской модели требуется больше времени и материалов. Однако мы примем для простоты, что вероятность возникновения одного землетрясения в последующие 30 лет нам известна и она равна p1, согласно донесению японских коллег. Тогда вероятность того, что ни одного события в ближайшие 30 лет не произойдет примем p0=1-p1. Окончательно вероятность превышения будет задаваться как

$display$

Для максимальной вероятности возникновения мега-землетрясения 0.4 мы построили кривую опасности (рис. 5). Она показывает запредельные значения проектных сейсмических воздействий (рис. 5), более X баллов. В реальности мы не учитывали много моментов, касающихся грунтовых условий, характера затухания сейсмических воздействий от мега-землетрясений и т.д. Известно, что сейсмические воздействия от землетрясения с магнитудой Mw=9.0 не намного больше тех, что вызваны событиями с Mw=8.3. А это значит, что мы явно переоценили сейсмические воздействия по формуле (15). Тем не менее, представленный расчет позволяет почувствовать масштаб проблемы.

Последнее. Почему японские сейсмологи дали вероятностную модель сейсмических источников на ближайшие 30 лет? В японских нормах приняты эксплуатационные сроки «жизни» зданий и сооружений 30 и 50 лет, в российских строительных правилах – 50 лет. Технология сейсмического районирования во всем мире (в том числе и у нас в стране) сводится к оценке сейсмических воздействий, которые будут превышены в течение 30 или 50 лет с заданной вероятностью. Для объектов нормальной ответственности (а это наши с вами жилые дома и офисы) строительными нормами задается вероятность превышения 0.1 в течение 50 лет. Для 30 лет задается вероятность 0.06 (это Япония и некоторые другие страны).

Таким образом, японские сейсмологи приняли консервативную оценку сейсмической модели, которая будет использована при актуализации карт сейсмического районирования. Это подчеркивает высокую ответственность и культуру сейсмологических изысканий японских ученых. Это также говорит о том, что Земля не знает границ, и исследование природных явлений требует объединения усилий ученых из разных стран.

Анализ выполнил Алексей Коновалов (a.konovalov@geophystech.ru), заместитель директора по научной работе в ООО «ГЕОФИЗТЕХ».

Список литературы

1. Baker, J.W. An introduction to probabilistic seismic hazard analysis (PSHA) / J.W. Baker // Report for the US Nuclear Regulatory Commission, Version 1.3. 2008. Section 1. P. 5-27.
2. Cornell, C.A. Engineering seismic risk analysis / Cornell C.A. // Bull. Seism. Soc. Am. 1968. Vol. 58. P.1583-1606.
3. Si H., Midorikawa S. New Attenuation Relationships for Peak Ground Acceleration and Velocity Considering Effects of Fault Type and Site Condition // Journal of Structural and Construction Engineering, A.I.J. 1999. V. 523. P. 63-70, (in Japanese with English abstract).
4. japan-forward.com/magnitude-9-earthquake-in-kuril-trench-highly-imminent-government-scientists
5. www.skr.su/news/279459
6. ria.ru/science/20171220/1511300109.html
7. www.the-japan-news.com/news/article/0004136947

Приложение №1

Baker, J.W. An introduction to probabilistic seismic hazard analysis (PSHA) / J.W. Baker // Report for the US Nuclear Regulatory Commission, Version 1.3. 2008. Section 1. P. 5-27.

Вероятностный анализ сейсмической опасности