1. На главную
  2. Электронные версии
  3. Физика Земли
  4. № 3, 2022

Физика Земли, 2022, № 3, стр. 56-75

Оценка сейсмической ситуации в зоне восточного побережья Эгейского моря Греции перед землетрясением с М = 6.7 (30.10.2020 г.)

Г. А. Попандопуло *

Институт физики Земли им. О.Ю. Шмидта РАН
г. Москва, Россия

* E-mail: gap@ifz.ru

Поступила в редакцию 11.04.2021
После доработки 30.09.2021
Принята к публикации 01.11.2021

Полный текст (PDF)

Аннотация

Проведены детальные исследования временных вариаций сейсмичности в зоне восточного побережья Эгейского моря Греции, где в окрестности о. Самос, 30.10.2020 г. произошло катастрофическое землетрясение с М = 6.7. Для анализа сейсмичности использовались: параметр b закона Гутенберга–Рихтера, величина выделенной сейсмической энергии в форме lg E2/3 и показатель сейсмической активности в виде lg N. Результаты анализа показали, что перед исследуемым землетрясением уверенно выделяются аномалии сейсмических параметров, длительностью до 15 лет. Аналогичные аномалии продолжительностью около 4–6 лет наблюдаются перед 5-ю сильными землетрясениями с магнитудой M ≥ 5.5, произошедшими в этом же регионе в течение последних 50 лет. Показано, что амплитуды временных аномалий увеличиваются, при приближении выборки землетрясений к глубинам очагов подготовки сильных событий. Делается вывод о том, что зависимость амплитуды временных вариаций сейсмических параметров от глубины может быть использована для оценки глубины очагов будущих сильных событий. Исследования временных вариаций сейсмичности в трех независимых сейсмоактивных регионах, в окрестностях о. Лесбос, о. Самос и о. Кос, показывают, что в период с 2010 по 2020 гг. перед тремя сильными землетрясениями с магнитудами М ≥ 6.1, произошедшими в этих регионах, наблюдаются предвестники исследуемых параметров. Пространственно-временная последовательность возникновения этих событий объясняется на основе концепции образования сейсмической бреши в окрестности о. Самос. Делается предположение о детектировании деформационной волны в зоне восточного побережья Эгейского моря в период с 2005 по 2020 гг., фронт которой распространяется по направлению с ССЗ на ЮЮВ, со скоростью 40–50 км/год. Выдвигаются соображения о том, что в ближайшие годы (10–15 лет), в этом регионе на участке между о. Лесбос и о. Кос возникновение сильных землетрясений с М ≥ 6.1 маловероятно. В то же время не исключено, что в ближайшие 2–3 года южнее о. Кос, в районе о. Родос, возможно возникновение сильного землетрясения с М ≈ 6.0–6.5.

Ключевые слова: сейсмичность, временные вариации сейсмичности, параметр b закона Гутенберга–Рихтера, сейсмическая энергия, прогноз землетрясений, оценка глубины прогнозируемого землетрясения.

1. ВВЕДЕНИЕ

Разрушительное землетрясение магнитудой M = 6.7 (магнитуда Ml) произошло 30.10.2020 г. на восточном побережье Эгейского моря, приблизительно в 14 км к северо-востоку от греческого острова Самос, на глубинах в диапазоне от 2 до 8 км. Землетрясение вызвало большие разрушения как на самом острове, так и на обширной территории западного побережья Турции. Особенно пострадал г. Измир – третий по численности город Турции. Здесь было разрушено более 20 домов, сотни людей остались без крова, а вызванное этим землетрясением цунами накрыло прибрежные территории. В г. Сеферихисаре, расположенном в 50 км к югу от г. Измира, цунами затопило несколько городских кварталов. В результате землетрясения, в Турции погибли 111 человек и 1035 человек было ранено. В Греции на о. Самос погибли два подростка, еще десять человек получили ранения.

Известно, что в районе о. Самос возникновение сильных землетрясений, приводящих к значительным ущербам, не редкость, о чем свидетельствуют как исторические записи, так и инструментальные наблюдения. Эпицентры сильных землетрясений возникали как на суше, так и в прибрежной зоне Эгейского моря [Tan et al., 2014]. В качестве последнего значительного события известно землетрясение 1904 г. с М = 6.8, в результате которого был нанесен серьезный ущерб городам и деревням на северо-западном побережье Турции. Отсюда следует, что произошедшее катастрофическое землетрясение с М = 6.7 (30.10.2020 г.) является одним из наиболее значимых событий, зарегистрированных в регионе за период более чем 100 лет, и для которого имеются инструментальные наблюдения. Все это предоставляет уникальную возможность ученым, связанным с науками о Земле, исследовать проявление предвестников различных геофизических и сейсмологических параметров перед этим событием.

Настоящая работа посвящена анализу пространственно-временных вариаций сейсмичности в восточной части Эгейского моря с целью оценки сейсмической опасности в период, предваряющий это катастрофическое землетрясение. Исследования проводились в зоне, включающей окрестности о. Самос и близлежащих регионов, охватывающих территории о. Лесбос и о. Кос. Для описания сейсмичности нами использовались параметр закона Гутенберга–Рихтера b-value (или просто b) [Ishimoto, Iida, 1939; Gutenberg, Richter, 1944], величина выделенной сейсмической энергии, представленной как lg E 2/3 [Keylis-Borok, 1959], и величина сейсмической активности, определяемая как логарифм числа землетрясений lg N.

Известно, что одна из основных характеристик сейсмичности описывается уравнением закона Гутенберга–Рихтера (ГР):

(1)
$\lg N~\,\, = \,\,~a--bM,$
где: lgN — логарифм числа землетрясений с магнитудами M (или ≥ M для кумулятивного распределения); a и b – параметры. Параметр a характеризует величину сейсмической активности региона наблюдения, а b это угол наклона прямой правой части частотно-магнитудного распределения (ЧМР) землетрясений. Наибольший интерес для сейсмологов представляет параметр b.

С физической точки зрения величина b в уравнении (1) показывает соотношение числа землетрясений малых и больших магнитуд в заданном интервале времени, в фиксированном районе выборки. Высокие значения b характеризуют относительно небольшое число землетрясений больших магнитуд, а, соответственно, низкие значения этого параметра свидетельствуют об обратном, т.е. об относительно большем количестве землетрясений больших магнитуд. Отсюда вытекает одно из фундаментальных свойств этого параметра: поскольку магнитуда землетрясений пропорциональна размерам длин разрывных нарушений, то, в случае пространственного распределения, величина b характеризует степень раздробленности (неоднородности) земной коры. В случае же временных вариаций, этот параметр свидетельствует об изменениях тектонических напряжений в земной коре [Mogi, 1962; Scholz, 1968; Mori, Abercrombie, 1997; Schorlemmer et al., 2004а; 2004в; 2005; Попандопуло, Баскутас, 2011; Попандопуло, 2018; 2020]. Параметр выделенной сейсмической энергии lgE 2/3 отражает текущую сейсмотектоническую деформацию земной коры [Keilis-Borok, 1959; Садовский, Писаренко, 1983]. И, наконец, третий параметр lgN характеризует сейсмическую активность региона наблюдения и, по сути, пропорционален величине параметра a уравнения (1).

Для вычисления перечисленных параметров используется компьютерная программа FastBEE, позволяющая на цифровой карте региона наблюдения определять область исследования, задавать параметры выборки землетрясений, выводить на экран компьютера ЧМР землетрясений, строить карты эпицентров, поперечный разрез сейсмичности и проводить пространственно-временное картирование величин минимальной представительной магнитуды Мс и трех упомянутых параметров [Papadopoulos, Baskoutas, 2009; Попандопуло, Баскутас, 2011; Попандопуло, 2018; 2020].

В качестве исходных данных выбран каталог землетрясений Греции, созданный Геодинамическим Институтом Национальной Обсерваторией Афин (Geodynamic Institute of National Observatory of Athens – GINOA), который содержит более 320 тысяч событий, зарегистрированных с 01.05.1964 по 01.11.2020 гг. Данный каталог находится в свободном доступе по адресу http://www.gein.noa.gr

2. ЭЛЕМЕНТЫ ТЕКТОНИКИ РЕГИОНА ИССЛЕДОВАНИЯ

Территория Эгейского моря и Западной Анатолии является одним из наиболее сейсмоактивных регионов Восточного Средиземноморья (рис. 1). Этот регион альпийско-гималайского пояса, с его геодинамическими процессами, представляет собой современную природную лабораторию, позволяющую исследовать различные проблемы геологии и геофизики [Jackson, 1994; Papazachos, 1999; Karakostas et al., 2010].

Рис. 1.

Восточно-Средиземноморский регион с основными активными границами и относительным движением тектонических плит. Район исследования обозначен прямоугольной областью. NAF – Северо-Анатолийский разлом (North Anatolian Fault); NAT – Северо-Эгейский прогиб (North Aegean Trough); CTF – трансформный разлом Кефалонии (Cephalonia Transform Fault), Hellenic Arc (Греческая дуга) – зона погружения Африканской плиты под Эгейскую плиту. Стрелками показаны относительные движения тектонических плит по разломам. Толстые серые стрелки указывают направления движения тектонической плит. Пунктирная стрелка – направление внутренней деформации (растяжения) Эгейского региона. Звездочки – землетрясения с магнитудой М ≥ 5.5.

Большая часть Эгейского моря относится к фрагменту евразийской континентальной плиты, называемой Эгейской микроплитой (рис. 1). На севере она примыкает к главной евразийской континентальной плите, а на востоке – к Анатолийской плите. На юге Эгейская микроплита сталкивается с Африканской тектонической плитой. Доминирующей геотектонической особенностью данного региона является субдукция литосферы Восточного Средиземноморья (Eastern Mediterranean) под Эгейское море вдоль Греческой дуги (Hellenic arc) и связанной с ней зоны Беньофа [Papazachos, Comninakis, 1971].

Анатолийская плита (Anatolian Plate) движется на запад вдоль разлома Северная Анатолия (NAF – North Anatolian Fault). Это правостороннее движение вдоль NAF продолжается вдоль Северного Эгейского прогиба (NAT – North Aegean Trough) и его южных рукавов [McKenzie, 1972; 1978]. Механизмы очагов землетрясений показывают, что деформация Эгейского моря обусловлена расширением земной коры в направлении с севера на юг, с преобладанием разломов, имеющих сбросовый характер [Scordilis et al., 1985; Papazachos, Delibasis, 1969; Papazachos et al., 1998].

Непосредственно зона исследования, в окрестности о. Самос, занимает восточную часть Эгейской микроплиты, сильно деформированную область растяжения, с компонентами движения в юго-западном направлении. Тектоника этой области определяется зоной перехода Анатолийской плиты к Эгейской по границе, которая проходит вдоль центральной и южной части побережья Турции, включающей соседние греческие острова [Papazachos, 1999].

3. МЕТОДИКА И АЛГОРИТМ ВЫЧИСЛЕНИЯ ПАРАМЕТРОВ СЕЙСМИЧНОСТИ

3.1. Краткое описание алгоритма FastBEE

В программе FastBEE упомянутые параметры сейсмичности определяются следующим образом:

Параметр b вычисляется по методу максимального правдоподобия с использованием формулы, предложенной в работах [Utsu, 1965; Aki, 1965]:

(2)
$b(t) = \frac{{\lg e}}{{\bar {M}(t) - (Mc - {{\Delta M} \mathord{\left/ {\vphantom {{\Delta M} 2}} \right. \kern-0em} 2})}},$
где: e – натуральное число; Мс – минимальная представительная магнитуда; $\bar {M}(t)$ – средняя магнитуда землетрясений, полученная по выборке данных землетрясений во временном окне сглаживания w для М ≥ Мс; ΔM – величина ячейки группирования (binned) магнитуды, в нашем случае ΔM = 0.1. Среднеквадратическая погрешность величины b вычисляется согласно модифицированной формуле, предложенной [Shi, Bolt, 1982]:

(3)
${{\sigma }_{b}} = 2.30{{b}^{2}}\sqrt {\sum\limits_{i = 1}^n {{{{{{({{M}_{i}} - \bar {M})}}^{2}}} \mathord{\left/ {\vphantom {{{{{({{M}_{i}} - \bar {M})}}^{2}}} {n(n - 1)}}} \right. \kern-0em} {n(n - 1)}}} } .$

Параметр выделенной сейсмической энергии lgE2/3 представлен в виде логарифма средней сейсмической энергии, выделенной в единице времени:

(4)
$\lg {{E}^{{2/3}}} = \lg \left( {\frac{1}{n}\sum\limits_{i = 1}^n {E_{i}^{{2/3}}} } \right),$
где: n – число землетрясений во временном окне осреднения; Ei – энергия i-го землетрясения с магнитудами М ≥ Мс, определяемая по формуле E = 101.5M+ 4.8 [Gutenberg, Richter, 1956].

Сейсмическая активность lgN представляется как логарифм числа землетрясений N с магнитудами M ≥ Mс, произошедшими в единицу времени (во временном окне осреднения). Имея в виду, что число землетрясений во времени описывается распределением Пуассона и среднеквадратическая ошибка вычисления n определяется по формуле ${{\sigma }_{N}} = \sqrt N $, можно определить среднеквадратическую ошибку логарифма числа землетрясения, которая будет равна дифференциалу этой функции: ${{\sigma }_{{\lg N}}} = {{0.4343} \mathord{\left/ {\vphantom {{0.4343} {\sqrt N }}} \right. \kern-0em} {\sqrt N }}$.

3.2. Метод средней точки для определения величины минимальной представительной магнитуды Мc

Как видно из формулы (2), одним из важных и необходимых условий вычисления величины b, является знание минимальной представительной магнитуды Мс, которая, как известно, является характеристикой системы наблюдения, и отражается в данных каталога землетрясений. В литературе известны множество методов оценки Мс [Смирнов, 1997; Wiemer, Wyss, 2000; 2002; Woessner, Wiemer, 2005; Papadopoulos, Baskoutas, 2009; Mignan, Woessner, 2012]. В то же время, признано, что использование ЧМР землетрясений для оценки величины Мс, вероятно, является самым простым и надежным способом [Wiemer, Wyss, 2000].

Отметим, что практически во всех компьютерных программах процедура определения величины Мс скрыта от пользователя и процесс вычисления определяется только математическими моделями, используемыми в алгоритмах. В результате пользователь полностью зависит от выбранного алгоритма и модели вычисления параметра Мс. Более того, при изменении параметров выборки землетрясений по глубине, площади или временному интервалу, величина Мс варьирует и тем самым пользователь не уверен в сохранении этой величины, а значит и в правильности результатов исследований. В отличие от этого, в используемой программе FastBEE величина Мс определяется в диалоговом режиме, и на протяжении всего процесса исследования пользователь осуществляет контроль и корректировку этой величины.

Процедура определения Мс заключается в следующем: при начальном запуске программы на экран компьютера выводится дискретное и кумулятивное ЧМР землетрясений (рис. 2). Цифровое представление ЧМР позволяет в интерактивном режиме непосредственно на экране компьютера визуально определить величину Мс. Как известно, точка Мс соответствует точке на оси магнитуд ЧМР, при которой наблюдается максимум числа событий при загибе кривой, характеризующей начало линейного падения логарифма числа землетрясений.

Рис. 2.

Кумулятивное (круги) и дискретное (треугольники) частотно-магнитудное распределение (ЧМР) землетрясений, произошедших в районе о. Самос в период с мая 1965 г. по ноябрь 2020 г. (рис. 1). Пример визуального определения величины представительной магнитуды Мс при использовании программы FastBEE: (a) – первая итерация (по умолчанию) при величине Мс = 2.0; (б) – вторая итерация при выбранной величине Мс = 3.6; b-value- (или b) – величина угла наклона прямой графика повторяемости; r – коэффициент корреляции между экспериментальными данными и аппроксимирующей прямой, проведенной через среднюю точку $\lg \overline N $ и $\bar {M}$ с углом наклона, равным b.

Отметим, что при любой наперед заданной величине Мс величины параметров a и b закона ГР определяются согласно формулам (2) и (1). Однако для вычисления реальной величины b необходимо определить такую точку Мс, при которой правый край экспериментальных данных ЧМР землетрясений был бы наилучшим образом аппроксимирован прямой уравнения (1). Напомним, что одним из условий проведения прямой на плоскости является знание по крайней мере двух величин: угла наклона прямой и одной точки на плоскости. В нашем случае, угол наклона прямой b вычисляется по формуле (2), а в качестве второго параметра используется точка на ЧМР землетрясений с координатами, равными средним значениям $\lg \overline N $ и $\bar {M}$, которые вычисляются для магнитуд M ≥ Мс. Таким образом, вычисляя величину b и среднюю точку ($\lg \overline N ,$ $\bar {M}$) на каждом шаге перебора Мс, мы можем подобрать такое значение Мс, при которой прямая уравнения (1) будет наилучшим образом аппроксимировать правый край ЧМР. Критерием окончательного выбора точки Мс является такая точка Мск на оси магнитуды, при которой наблюдается наивысший коэффициент корреляции между экспериментальными данными и аппроксимирующей прямой уравнения (1). Конкретный пример такого построения представлен на рис. 2.

Экспериментальным способом установлено, что при правильном определении Мс, коэффициент корреляции r, как правило, больше 0.95. Обычно, для определения Мс и b для одной выборки, необходимо 2–3 итерации, что занимает около 5–10 с. Точность визуального определения величины Мс определяется величиной шага ячейки группирования (binned), используемой при построении ЧМР и равна половине ее величины. В нашем случае величина ячейки группирования ΔM = 0.1, что указывает на то, что точность вычисления параметра Мс составляет ±0.05 [Лукк, Попандопуло, 2012; Попандопуло, Лукк, 2014].

Здесь уместно провести сравнение способа оценки величина Мс методом средней точки (СТ) с широко используемым методом (GFT – goodness-of-Fittest)), предложенным в работе [Wiemer, Wyss, 2000]. Идея предложенного GFT-метода заключается в подгонке такого синтетического распределения, основанного на степенном законе, с коэффициентами (a, b и Mс), которое наилучшим образом соответствовало бы наблюденному ЧМР землетрясений. В то же время авторы упомянутой работы допускают, что в некоторых случаях ЧМР слишком изогнуты или бимодальны, чтобы получить удовлетворительную подгонку для соответствия простому степенному закону. Это может быть обусловлено сильной пространственной и временной неоднородностью данных каталога землетрясений или различными физическими процессами, происходящими в земной коре.

Сравнивая алгоритм метода СТ с методом GFT, можно заметить, что по сути это два идентичных метода, где аппроксимирующая прямая в методе СТ заменяет прямую линейного синтетического распределения степенного закона в методе GFT. Однако подход, используемый в программе FastBEE, имеет некоторое преимущество, поскольку, во-первых, более простой с точки зрения вычислительных процедур и, во-вторых, позволяет визуально контролировать форму ЧМР в каждом конкретном случае и избежать ошибок, вызванных плоскими или биноминальными распределениями наблюденных ЧМР землетрясений [Popandopoulos, Chatziioannou, 2014; Попандопуло и др., 2016]

3.3. Построение и фильтрация временных рядов

В алгоритме FastBEE процедура построения временных рядов состоит из двух этапов. На первом этапе создается временный файл каталога землетрясений, относящийся к району исследования, определяется величина Mс и вычисляются средние значения величин a и b за исследуемый период времени. Затем создаются массивы сумм месячных значений для каждого из исследуемых параметров. Так, для вычисления параметра b используется один двухмерный массив месячных значений ЧМР землетрясений:

$FMD(t,j) = \sum\limits_{i = 1}^{n(t)} {\sum\limits_{j = 1}^K {{{k}_{{i,j}}}} } ,$
где: $FMD(t,j)$ – ЧМР землетрясений в одном месяце; t – порядковый номер месяца с начала момента исследования; j – порядковый номер ячейки (binned) группирования магнитуд; $n(t)$ – число землетрясений в месяце; k = int(M/ΔM) – номер ячейки группирования, куда попало i-е землетрясение, равное целой части его магнитуды, деленной на размер ячейки ΔM( = 0.1); K – число ячеек группирования.

Для построения временных рядов выделенной сейсмической энергии lg E2/3 необходимо создание двухмесячных массивов:

$SE(t) = \sum\limits_{i = 1}^{n(t)} {E_{i}^{{2/3}}} ,\,\,n(t),$
где: $SE(t)$ сумма месячных значений выделенной энергии землетрясений; t – порядковый номер месяца; i-порядковый номер землетрясения в месяце; $n(t)$ – число землетрясений в месяце t.

Для построения массивов величины сейсмической активности lgN(t) используются массив $\lg (n(t))$ – состоящий из логарифма суммы месячных чисел событий.

Создание сумм месячных массивов позволяет при помощи рекуррентных соотношений проводить построение временных рядов с одновременным сглаживанием и статистической оценкой текущих средних значений в любом временном окне w, кратном одному месяцу.

На втором этапе проводится осреднение и фильтрация месячных значений в выбранном окне сглаживания w. Известно, что использование простой процедуры сглаживания равносильно процессу фильтрации временного ряда фильтром прямоугольной формы. С целью подавления побочных лепестков сглаживания, автоматически возникающих при такой фильтрации, производится вторичная фильтрация временных рядов фильтром треугольной формы с тем же временным окном сглаживания w. Фильтр треугольной формы, как известно, является фильтром низких частот, который без искажения пропускает сигналы с периодами, равными интервалу более половины окна сглаживания, и практически не имеет побочных эффектов [Bath, 1974].

Подчеркнем, что при первоначальном осреднении (сглаживании) значения на выходе фильтра приписываются к концу временного окна сглаживания, в то же время, при вторичной фильтрации – фильтром треугольной формы, значения на выходе приписываются уже к центру временного окна. При анализе необходимо иметь в виду, что в результате вторичной фильтрации, на графиках временных рядов возникает пустое пространство, между концом кривой и моментом окончания периода исследования.

3.4. Определение параметров района исследования и обработки данных землетрясений

Для адаптации программы FasBee к региону исследования создается векторный файл географических координат административных границ и основных геотектонических структур региона наблюдения, который выводится на экран компьютера (рис. 1). В интерактивном режиме задаются параметры выборки локального района (ЛР) исследования и параметры выборки данных землетрясений, а для построения временных рядов – временное окно сглаживания.

В процессе анализа, параметры выборки землетрясений, к которым относятся временной интервал исследования T, диапазон интересующих магнитуд Mmin и Mmax, интервал глубин очагов землетрясений Hmin и Hmax и временное окно сглаживания w варьируют.

Важно отметить, что длина временного окна сглаживания w зависит от сейсмической активности ЛР выборки, размер которой определяется путем подбора такой величины w, при котором минимальное число землетрясений N в текущем окне сглаживания на протяжении всего периода исследования, была бы больше 50. Известно, что величина N ≥ 50 является необходимым критерием для надежной оценки точности вычисления b [Marzocchi, Sandri, 2004; Sandri, Marzocchi, 2007; Woessner, Wiemer, 2005; Mignan, Woessner, 2012]. Параметр длины окна сглаживания (фильтрации) также может меняться в зависимости от цели задачи исследования, при анализе длиннопериодных или короткопериодных вариаций сейсмических параметров. Здесь возникает неопределенность, чем больший размер ЛР выбирается, тем меньшая длина окна сглаживания w может быть использована, и наоборот.

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

Таблица 1.

Список землетрясений с М ≥ 5.5, произошедших в окрестностях о. Самос (рис. 1) с мая 1965 г. по ноябрь 2020 г. в диапазоне глубин 0–45 км

№ Дата
год. мес. день 		Широта (сш.), град. Долгота (вд.), град. Глуб H, км Магн. ML
1 1969.04.06 38.400 26.800 10 5.5
2 1979.06.14 38.700 26.600 10 5.5
3 1992.11.06 38.090 27.190 39 5.7
4 2005.10.17 38.130 26.590 29 5.6
5 2005.10.20 38.150 26.250 25 5.6
6 2020.10.30 37.924 26.810 02 6.7

4. КРАТКАЯ ХАРАКТЕРИСТИКА КАТАЛОГА ЗЕМЛЕТРЯСЕНИЙ ГРЕЦИИ

Известно, что каталоги землетрясений, имеющие долговременную историю создания, являются крайне неоднородными с точки зрения минимальной представительности магнитуды Мс. Безусловно, что и каталог землетрясений Греции не является исключением. Начиная с момента его создания, с 1964 г., его характеристики постоянно улучшалось. История развития системы сейсмологических наблюдений в Греции и изменения параметров каталога землетрясений (уменьшение величины Мс, вычисление координат эпицентра и глубины и т.д.), довольно детально описаны в работах [Papanstasiou, 2011; D’Alessandro et al., 2011; Chouliaras et al., 2013; Mignan, Chouliaras, 2014; Дещеревский, Сидорин, 2011; Popandopoulos, Chatziioannou, 2014; Попандопуло и др., 2016].

Согласно работам [Papanstasiou, 2011; Mignan, Chouliaras, 2014], с момента основания наблюдательной системы в Греции с 1964 г. выделяются три основных этапа развития:

(I) Переход в 1995 г. от аналоговых приборов к цифровым сейсмическим станциям.

(II) Создание с 2007 по 2011 гг. Объединенной Сейсмологической Сети Греции (Hellenic Unified Seismological Network – HUSN), состоящей более чем из 100 сейсмических станций вместо 35 существовавших ранее.

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

Здесь необходимо вкратце рассмотреть влияние развития системы сейсмологических наблюдений Греции на качество каталога землетрясений в восточной части Эгейского моря. Локальный каталог таких землетрясений получен в результате выборки данных землетрясений из района, имеющего прямоугольную форму, размером 130 × 180 км с центром в эпицентре сильного события с М = 6.7 (30.10.2020 г.). Азимутальная ориентация района совпадает с эллипсом распределения плотности сейсмичности в этой зоне и приблизительно соответствует направлению береговой линии Турции (рис. 1).

На рис. 3 представлен временной ряд логарифма чисел землетрясений lgN с магнитудой М ≥ 1.0, произошедших в выбранном регионе с мая 1964 г. по ноябрь 2020 г. Кривая lg N получена путем осреднения данных временным окном 9 мес. с шагом 1 мес. Видно, что картина изменения чисел землетрясений lg N хорошо отражает упомянутые этапы развития системы наблюдения. На рис. 2 представлено дискретное и кумулятивное ЧМР землетрясений, произошедших в исследуемом районе за тот же период времени. Видно, что ЧМР имеет сложную структуру с явно выраженной двумодальной формой с максимумами магнитуд М = 3.2 и 3.6. Кроме того, видна пологая область в диапазоне магнитуд от 2.0 до 3.0. Причина такой формы кривой кроется в наложении трех совокупностей данных землетрясений с различными величинами представительной магнитуды Мс, возникших в результате развития системы наблюдений в Греции. Так, например, пик, наблюдаемый при магнитуде Мс = 3.6, отражает общую минимально представительную магнитуду всего периода наблюдения с мая 1965 г. по ноябрь 2021 г. Основной максимум, наблюдаемый при М $ \approx $ 3.2, скорее всего отражает сейсмичность, регистрируемую на II и III этапах (рис. 3). И, наконец, пологая часть кривой ЧМР, в диапазоне магнитуд ΔМ = 2–3 отражает сейсмичность только III этапа, регистрируемую сетью HUSN после 2007 г., где разрешающая способность сети существенно улучшилась.

Рис. 3.

График средних чисел землетрясений с магнитудой М > 1.0, построенный с временным окном сглаживания 9 месяцев с шагом 1 месяц в районе о. Самос (рис. 1); I, II и III –основные этапы развития сети сейсмических станций в Греции.

5. ВРЕМЕННЫЕ ВАРИАЦИИ СЕЙСМИЧНОСТИ В ОКРЕСТНОСТЯХ О. САМОС

Настоящая работа посвящена исследованию предвестников сейсмичности перед землетрясением с М = 6.7 (30.10.2020 г.). Анализ проводился путем построения временных рядов с различными параметрами выборки данных землетрясений путем варьирования минимальной представительной магнитуды Mc, диапазона исследуемых глубин Hmin и Hmax и временного окна сглаживания w.

5.1. Временные вариации сейсмичности при различных величинах Мс и длине окна сглаживания w

На рис. 4 представлены временные вариации b, lg E2/3 и lg N для района о. Самос (рис. 1) во временном интервале с мая 1965 г. по ноябрь 2020 г. При построении временных рядов использовались землетрясения с магнитудами М ≥ (Мс = 3.6) (рис. 2б). Временные ряды строились путем осреднения с последующей фильтрацией данных землетрясений в диапазоне глубин ΔH = 0–45 км, во временном окне сглаживания w = 25 мес. и с шагом в 1 месяц. Напомним, что в результате вторичной фильтрации возникло пустое пространство между окончанием временного ряда и моментом возникновения сильного землетрясения с М = 6.7 (см. разд. 3.3).

Рис. 4.

Временные вариации b, lg E2/3 и lg N для землетрясений с М ≥ 3.6, произошедших в диапазоне глубин 0–45 км, в окрестности о. Самос (рис. 1). Кривые построены с временным окном сглаживания данных 25 мес., с шагом 1 мес. Пунктирными линиями показаны долговременные средние величины. Вертикальными стрелками отмечены моменты возникновения землетрясений с магнитудой М ≥ 5.5. Номера (N) и магнитуды (М) землетрясений соответствуют номерам и магнитудам, указанным в табл. 1. На кривых b и lg N вертикальными отрезками нанесена величина 1σ текущих средних значений, определенных в окне сглаживания. Для кривой lg E2/3 горизонтальными линиями по обе стороны от долголетних средних значений отмечен среднестатистический разброс данных 1σ текущих величин.

На графиках вертикальными стрелками отмечены моменты возникновения сильных землетрясений с магнитудой M ≥ МПЗ, произошедших в пределах района выборки. За величину минимальной магнитуды прогнозируемого землетрясения (МПЗ), согласно ЧМР (рис. 2), принята величина M = 5.5. Номера стрелок (N) и магнитуды землетрясений (M) на графиках соответствуют данным, указанным в табл. 1, а их эпицентры показаны на рис. 1.

Для кривых b и lg N величины доверительных интервалов 1σ показаны вертикальными отрезками для каждого временного окна сглаживания (разд. 3.1). На кривых lgE 2/3 горизонтальными линиями по разные стороны от долговременных средних значений отмечен среднеквадратический разброс данных ±1σ, который соответствует 70% интервалу доверия.

Результаты исследования, полученные ранее, при использовании описанной методики показали, что перед сильными землетрясениями с магнитудой M ≥ МПЗ, наблюдаются аномалии параметра b “бухтообразной” формы выпуклые вверх (положительные аномалии). Причем за аномальные значения принимаются текущие значения b, превышающие ее средний уровень более чем на 1σ.

Важно отметить, что на первой стадии развития аномалии, т.е. при возрастании кривой b, превышающей средний уровень, возникновение сильного землетрясения крайне маловероятно, в то время как на второй стадии аномального процесса после прохождения кривой своего максимума и начала снижения, вероятность возникновения землетрясения велика [Попандопуло, 2020]. Возникновение сильных землетрясений возможно на протяжении всего периода падения кривой b, вплоть до достижения минимальных значений, которые могут выходить за рамки долговременных средних значений [Попандопуло, Баскутас, 2011; Попандопуло, 2018; 2020]. Анализ качества среднесрочного прогноза сильных землетрясений, при описанном подходе методики FastBEE, показал, что в 84% случаев появления наблюдаемых аномалий b в районах наблюдения заканчиваются сильными землетрясениями с М ≥ МПЗ, а средняя вероятность успешного прогноза, с учетом ложных тревог и пропуска целей, составляет 71% [Попандопуло, 2020].

Поведение кривых lgE2/3 перед сильными землетрясениями показывает “бухтообразные” изменения, выгнутые вниз (отрицательные аномалии), величины которых выходят за нижний предел –1σ от величины долговременного среднего. В первом приближении картину временных вариаций lg E2/3 можно представить как зеркальное отражение кривых b. Однако корреляционное уравнение связи временных вариаций b и lgE 2/3 имеет более сложный характер и хорошо аппроксимируется степенной функцией [Попандопуло, 2018]. Дальнейший анализ зависимости b и lgE 2/3 выходит за рамки настоящей работы. В то же время, в первом приближении, можно сказать, что кривые lgE 2/3 имеют противоположный характер проявления временных аномалий по сравнению с кривыми b.

Временные вариации величины lgN носят вспомогательный характер и только в особых случаях отражают процесс подготовки сильных землетрясений, выраженный в “бухтообразном” изменении кривой lgN выгнутом вниз (отрицательные аномалии). Как правило, амплитуды таких аномалий превышают величину 1σ [Попандопуло, Баскутас, 2011; Baskoutas, D’Alessandro, 2014; Попандопуло, 2018; 2020].

Данные, представленные на рис. 4, показывают, что за период исследования в районе о. Самос, в диапазоне глубин 0–45 км, произошло 6 землетрясений с M ≥ МПЗ (см. табл. 1). Видно, что в большинстве случаев перед сильными событиями наблюдаются описанные выше аномалии параметра b. Поведение кривой lgE 2/3 перед сильными землетрясениями имеет “бухтообразное” понижение кривой вплоть до нижней долговременной отметки 1σ с последующим ее восстановлением, однако выражено менее четко. Как было упомянуто выше, кривая lgN не показывает четких предвестников сильных событий.

Рассматривая кривые параметра b, можно заметить, что перед землетрясением № 2 наблюдается 4-х летняя аномалия, перед землетрясением № 3, соответственно, 5-летняя, а перед землетрясениями № 4, № 5 длительность аномалии составляет порядка 6 лет. Наблюдаемая аномалия в течение 2008–2011 гг. формально может быть охарактеризована как ложная, а перед главным событием № 6 с M = 6.7 наблюдается пятилетняя аномалия относительно небольшой амплитуды. В результате картина временных вариаций сейсмических параметров показывает, что практически перед всеми событиями с M ≥ МПЗ наблюдаются характерные аномалии b и lgE 2/3 по форме, удовлетворяющей феноменологической модели подготовки сильных землетрясений, предложенной ранее в работах [Попандопуло, Баскутас, 2011; Попандопуло, 2018; 2020]. Как было упомянуто выше, момент возникновения сильного землетрясения наблюдается во второй половине аномального периода.

Анализируя кривую lg N (рис. 4), видно, что число событий, попадающих в выбранное окно сглаживания (w = 25 мес.), в период с 1965 по 1980 гг., в среднем, равно 19, в период 1980–1996 гг. равно 35 и с 1996 г. по ноябрь 2020 г., соответственно, равно 68. Известно, что минимальное число данных, необходимых для статистически надежной оценки b, должно быть не менее 50 [Wiemer, Wyss 2002; Marzocchi, Sandri, 2003; Sandri, Marzocchi, 2007]. Из этого следует, что к кривым, отражающим временные вариации b в период с 1965–1995 гг., необходимо относиться с осторожностью. В то же время, наблюдаемые закономерности временных изменений b и lg E 2/3 перед сильными событиями дают определенную уверенность в объективности представленной картины. Более того, картина временных изменений после 1995 г., где среднее число данных в окне сглаживания превышает 50 (рис. 3), в целом, не меняет характер поведения кривых перед сильными землетрясениями, что может также косвенно свидетельствовать о надежности полученных результатов. И, наконец, видимая корреляция кривых b и lgE 2/3 также может свидетельствовать об объективности наблюдаемых временных изменений исследуемых параметров.

При анализе кривой lgN (М ≥ 3.6) возникает вопрос, требующий дополнительного исследования используемых данных. С одной стороны, для построения временных рядов использовались данные землетрясений с магнитудами М ≥ (Мс = = 3.6), для которых Мс определена по ЧМР землетрясений во всем интервале исследования 1965–2020 гг. (рис. 2б), что должно свидетельствовать о том, что в этом промежутке времени все землетрясения с М ≥ 3.6 были зарегистрированы без пропуска. В этом случае можно было ожидать некоторую стационарность кривой lgN с точки зрения долговременной составляющей временных вариаций землетрясений с М ≥ 3.6 за 55-летний период наблюдения [Box, Jenkins, 1970]. Но поведение lgN показывает некоторую корреляцию с этапами развития системы наблюдений в Греции, т.е. наблюдаемое изменение в 1995 и 2007 гг. (см. разд. 4), что, в принципе, противоречит предыдущему высказыванию о пропуске землетрясений с М ≥ 3.6.

Для разрешения этой неопределенности было проведено повторное построение временных рядов с другими параметрами выборки данных землетрясений. Во-первых, была увеличена, насколько это было возможно, минимальная представительная магнитуда до величины Мс = 3.8 (рис. 2б). Дальнейшее увеличение Мс приводило к отказу функционирования вычислительной программы FastBEE, по причине возникновения пустых окон осреднения. И, во-вторых, поскольку при увеличении Мс до 3.8 число данных землетрясений в окне осреднения (w = 25 мес.) уменьшилось, необходимо было его увеличение до w = = 53 мес. Результаты построения временных рядов с новыми параметрами осреднения представлены на рис. 5.

Рис. 5.

Временные вариации b, lgE2/3 и lg N для землетрясений с М ≥ 3.8, произошедших в диапазоне глубин 0–45 км, в окрестности о. Самос (рис. 1). Кривые построены с окном сглаживания данных 53 мес. с шагом 1 мес. Условные обозначения те же, что и на рис. 4.

Полученные кривые показывают незначительное увеличение амплитуд аномалий перед землетрясениями № 3 и № 4, № 5. Одновременно, наблюдается сглаживание амплитуды аномалии перед событием № 2 и № 6 с M = 6.7. Такая картина временных изменений является следствием увеличением длины окна сглаживания. Можно также заметить, что, в целом, картина временных вариаций чисел землетрясений lgN (М ≥ 3.8) по сравнению с lgN (М ≥ 3.6) существенно не изменилась, что может говорить о объективности полученных аномалий.

И, наконец, представляет интерес отдельное определение величины Мс для временного периода (I) с 1965 по 1995 гг. (рис. 3). График ЧМР землетрясений, представленный на рис. 6, показывает, что в упомянутом интервале времени величина Мс = 3.6, а это значит что эта величина действительна и для всего интервала исследования с мая 1965 г. по ноябрь 2020 г.

Рис. 6.

Частотно-магнитудное распределение (ЧМР) землетрясений, произошедших в районе о. Самос в период с мая 1965 г. по январь 1995 г. (рис. 1). Условные обозначения те же, что и на рис. 2.

На основании проведенного анализа можно утверждать, что при выборе землетрясений с магнитудами М ≥ 3.6 изменения, происходившие в системе сейсмологических наблюдений Греции, не влияют на результаты анализа, а наблюдаемые аномалии исследуемых параметров (рис. 4 и рис. 5) отражают процессы подготовки сильных землетрясений с М ≥ 5.5 в исследуемом регионе.

5.2. Временные вариации в различных диапазонах глубин ΔH

Ранее проведенные исследования в Гармском районе Таджикистана показали, что картина временных вариаций параметра b неоднородна с глубиной и зависит от местоположения гипоцентра очага подготовки сильного землетрясения [Попандопуло, 2020]. Рассматривая таблицу сильных землетрясений с М ≥ 5.5, произошедших в районе исследования в окрестности о. Самос, видно, что гипоцентры этих землетрясений находятся в диапазонах глубинах от 2 до 39 км. В связи с этим, интерес представляет анализ временных вариаций исследуемых параметров в различных диапазонах глубин. Для этой цели, в качестве минимальной представительной магнитуды была принята величина Мс = 3.2, которая появляется на первом максимуме в ЧМР землетрясений, построенном в период с января 1995 г. по ноябрь 2020 г. (рис. 7). Выбор этой величины Мс дает возможность увеличить число событий в выборке данных землетрясений, с одновременным делением земной коры на сейсмоактивные слои.

Рис. 7.

Определение минимальной представительной магнитуды Мс = 3.2 для землетрясений, произошедших в районе о. Самос в период с января 1995 г. по ноябрь 2020 г. (рис. 1). Условные обозначения те же, что и на рис. 2.

На рис. 8а, 8б представлены временные вариации исследуемых параметров по данным землетрясений из двух диапазонов глубин ΔH1 = 0–23 км и ΔH2 = 0–18 км. Анализ представленных кривых показывает, что временные вариации сейсмических параметров, полученные в двух диапазонах глубин, имеют устойчивую картину. Отметим, что дальнейшее уменьшение диапазона глубин не представлялось возможным по причине недостатка данных. Кривые, представленные на рис. 8а, 8б показывают, что аномалии b и lgE 2/3 в исследуемых диапазонах глубин наблюдаются только в периоды, предшествующие двум сильным событиям № 2 и № 3 (событие № 3 на рис. 8a, 8б соответствует событию № 6 на рис. 4, рис. 5). Четырехлетняя аномалия небольшой амплитуды, слегка превышающая 1σ, наблюдается также перед землетрясением, произошедшим 06.11.1992 г. с М = 5.7 на глубине 39 км (см. табл. 1), гипоцентр которого выходит за пределы выбранного диапазоне глубин землетрясений 0–23 км.

Рис. 8.

Временные вариации b, lgE2/3 и lg N для землетрясений с М ≥ 3.2, произошедших в диапазоне глубин ΔH1 = 0–23 км (a) и ΔH2 = 0–18 км (б), в период с мая 1965 г. по ноябрь 2020 г. в окрестности о. Самос (рис. 1). Временное окно сглаживания w = 25 мес., шаг 1 мес. Условные обозначения те же, что и на рис. 4.

В то же время, в период с 1995 по 2005 гг., т.е. в течение 10 лет, мы не наблюдаем существенных изменений параметров b и lgE 2/3, и одновременно с этим, в этот период времени в исследуемом диапазоне глубин не возникает сильных событий с М ≥ 5.5 (см. табл. 1). Видно, что величина максимальных амплитуд временных изменений параметра b для глубин ΔH2 (рис. 8б) больше, чем для ΔH1 (рис. 8а), так, например, максимальные значения b для ΔH2 равно 2.5, в то время как для ΔH1 равно 1.88 соответственно. Видно, что такое значительное увеличение амплитуд наблюдается в период, предшествующий катастрофическому землетрясению № 3 с М = 6.7.

Длительность аномалий параметров b перед землетрясением № 3 (рис. 8) хорошо прослеживается, начиная с 2005 по 2020 гг., т.е. на протяжении 15 лет. Кроме того, примерно в этот же период времени в этих диапазонах глубин перед землетрясением № 3 надежно выделяются аномалии на кривых lgN, которые указывают на значительное уменьшение сейсмической активности землетрясений с М ≥ 3.2 в исследуемом районе. Данный феномен хорошо известен в сейсмологии как эффект сейсмического затишья, наблюдаемый во многих регионах мира (напр. [Нерсесов и др., 1976; Тейтельбаум, Пономарев, 1983]). Формально начало аномалии lgN можно определить, начиная с 2011г. Однако видно, что в период с 1995 по 2020 гг., при улучшении наблюдательной системы Греции (рис. 3), кривые lgN находятся выше долговременных средних значений, и на самом деле процесс зарождения аномалии lgN мог начаться в 2008 г.

Большой интерес представляет сравнение кривых исследуемых параметров для диапазонов глубин ΔH = 0–45 км (рис. 4 и рис. 5), с соответствующими кривыми для диапазонов глубин ΔH1 = 0–23 км и ΔH2 = 0–18 км (рис. 8a, 8б). На рис. 4 и рис. 5 отчетливо выделяются аномалии b и lgE 2/3 перед сильными землетрясениями № 3, и № 4, № 5, однако в период, предшествующий тем же землетрясениям, на рис. 8a, 8б таких изменений не обнаружено. Аналогичную картину можно наблюдать и в период, предшествующий землетрясениям № 2 и № 6 (№ 3 на рис. 8). Здесь на рис. 4 и рис. 5 мы видим относительно небольшое изменение параметра b в диапазоне глубин 0–45 км и значительные изменение параметров b, lgE 2/3и lgN в диапазонах глубин на ΔH1 = 0–23 км и ΔH2 = = 0–18 км (рис. 8a, 8б).

Рассматривая данные, представленные в таблице, видно, что гипоцентры событий № 3 и № 4, № 5 находятся в нижних слоях земной коры на глубинах 39, 29 и 25 км соответственно. В то же время, землетрясения № 2 (Η = 10 км) и № 6 (№ 3 на рис. 8) (Н = 2–10) произошли на небольших глубинах. Все это говорит о том, что аномалии сейсмических параметров выделяются в диапазонах глубин, близких к глубинам очагов подготовки сильных землетрясений. В результате, аномалии перед землетрясениями № 3 (H = 39) и № 4, № 5 (H = 29.25) хорошо прослеживаются в диапазонах глубин 0–45 км и мало заметны в верхних слоях 0–23 км. В то же время, землетрясения № 2 и главное сильное событие № 6 с М = 6.7 находятся в приповерхностной зоне, и поэтому аномалии хорошо выделяются в верхних диапазонах глубин 0–23 и 0–18 км.

Отсюда следует, что при анализе временных вариаций выбранный диапазон глубин играет роль фильтра и при определенных условиях может или увеличивать, или уменьшать амплитуды аномалий сейсмических параметров. Так, видно, что использование общего диапазона глубин 0–45 км нивелирует эффекты, происходящие в верхних слоях земной коры, в то время как в верхних диапазонах глубин 0–18 км, где происходила подготовка сильного землетрясений с М = 6.7, наблюдается усиление амплитуд временных аномалий сейсмических параметров.

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

5.3. Временные вариации сейсмичности в трех независимых локальных районах

Особый интерес представляет анализ поведения сейсмических параметров в районах, расположенных вблизи исследуемой зоны. Такой анализ, с одной стороны, позволит проверить статистическую значимость результатов как результат эксперимента, полученный при использовании независимых выборок, а с другой стороны, при таких исследованиях можно получить представление о геодинамических процессах в обширной зоне в период, предшествующий возникновению сильного землетрясения с М = 6.7 (30.10.2020 г.).

Для проведения такого эксперимента были выбраны дополнительно два локальных района (ЛР), расположенных к северу и к югу от исследуемого района в окрестности о. Самос (рис. 9). Первый из них ЛР (а), размером 100 × 110 км, включает о. Лесбос, где 12.06.2017 г. произошло сильное землетрясение с М = 6.1; второй ЛР (в), размером 90 × 110 км расположен в окрестности о. Кос, где 20.07.2017 г. произошло также сильное землетрясение с М = 6.2. Первоначально исследуемый район (рис. 1) в окрестности о. Самос был уменьшен до размера 95 × 100 км – ЛР (б).

Рис. 9.

Районы выборки данных землетрясений в окрестностях о. Лесбос (a), о. Самос (б) и o. Кос (в) (карта справа) и, соответствующие им (слева) ЧМР землетрясений, произошедших в период с января 1999 г. по ноябрь 2020 г. Условные обозначения те же, что и на рис. 2. Звездочки на карте (справа) – сильные землетрясения с магнитудой М = 6.1 (а), М = 6.7 (б) и М = 6.2 (в). Стрелками на карте обозначены последовательность возникновения сильных землетрясений и вероятное возникновение будущего землетрясения с М ≥ 6.1 в районе о. Родос.

Для каждого ЛР определялась величина Мс по землетрясениям, произошедшим во временном интервале с января 1995 г. по ноябрь 2020 г. Уменьшение интервала исследования обусловлено стремлением к увеличению статистики данных, наблюдаемой с 1995 г. (рис. 3). Графики ЧМР землетрясений и определение Мс для трех независимых ЛР представлены на рис. 9. Видно, что величины Мс для ЛР (а), (б) и (в), близки и равны 2.9, 3.0 и 3.1 соответственно.

На рис. 10–рис. 12 представлены временные изменения b, lgE 2/3 и lgN соответственно для трех районов выборки. Кривые (а), (б), (в) на рис. 10–рис. 12 соответствуют районам (а), (б) и (в) на рис. 9. На рис. 10–рис. 12 показаны фрагменты временных изменений исследуемых параметров в интервале с января 1995 г. по ноябрь 2020 г., взятые из временных рядов общего периода наблюдения с мая 1965 г. по ноябрь 2020 г. Кривые получены по землетрясениям, произошедшим в диапазоне глубин 0–20 км, с магнитудами М ≥ Мс, определяемыми для каждого ЛР. Построение временных рядов проводилось путем осреднения данных во временном окне сглаживания 27 мес. со сдвигом в 1 мес. и с последующей их фильтрацией фильтром треугольной формы с тем же временным окном (см. разд. 3.3).

Рис. 10.

Временные вариации b, по данным землетрясений, произошедшим в трех районах в окрестностях о. Лесбос (a), о. Самос (б) и о. Кос (в), в диапазоне глубин 0–20 км, в период с января 1995 г. по ноябрь 2020 г. (рис. 9). Кривые построены с окном сглаживания данных 27 мес. с шагом 1 мес. Вертикальными стрелками отмечены моменты возникновения сильных землетрясений и их магнитуды. На кривых вертикальными отрезками нанесены величины 1σ текущих средних значений. Штрих-пунктирные линии (1, 2, 3, 4) указывают предполагаемое направление и скорость фронта деформационной волны с ССЗ на ЮЮВ.

Рис. 11.

Временные вариации выделенной сейсмической энергии lgE2/3. Горизонтальными линиями по обе стороны от долголетних средних значений, отмечен среднестатистический разброс текущих средних значений 1σ. Условные обозначения те же, что и на рис. 10.

Рис. 12.

Временные вариации сейсмичности lgN. Условные обозначения те же, что и на рис. 10.

При рассмотрении кривых временных изменений параметра b (рис. 10) видно, что перед тремя сильными землетрясениями с М ≥ 6.1, произошедшими в трех ЛР (а), (б), (в), наблюдаются значимые временные аномалии продолжительностью от 7 до 10 лет. Как и наблюдаемые ранее, кривые на рис. 4 и рис. 5, временные аномалии параметра b имеют “бухтообразную” форму, выпуклую вверх.

Для ЛР (а) наиболее ярко аномалия проявляется в течение 2010–2017 гг., для ЛР (б) аналогичная аномалия наблюдается в период с 2011 по 2020 гг., и для ЛР (в), соответственно, в период 2014–2017 гг. Подчеркнем, что признаком окончания аномалии параметра b является достижение минимальных значений кривой и их переход в плато или к возрастанию [Попандопуло, 2018]. Такая картина хорошо прослеживается на кривых (a) и (в), где в 2017 г. в этих районах произошли сильные землетрясения. В то же время видно, что наблюдаемая аномалия на кривой (б) продолжала уменьшаться до тех пор, пока не произошло сильное землетрясение с М = 6.7. Такое поведение кривых косвенно свидетельствует о том, что временные изменения изучаемых параметров отражают процессы подготовки сильных землетрясений.

При рассмотрении кривых, представленных на рис. 11 видно, что параметр выделенной сейсмической энергии lgE 2/3 имеет отрицательные “бухтообразные” аномалии, предваряющие сильные землетрясения во всех исследуемых ЛР. Аномалии lgE 2/3 представляются как их зеркальное отображение на кривых b (рис. 10). В целом, в пределах аномального периода, начиная с 2005 г., видны незначительные волнообразные колебания, синхронные с кривыми b на рис. 10. Отметим, что параметры lgE 2/3 и b с точки зрения алгоритма их вычисления являются независимыми величинами, и их видимая корреляция также свидетельствует о надежности полученных результатов.

Особый интерес представляет рассмотрение временных вариаций активности землетрясений lgN (М ≥ Мс), представленных на рис. 12. Хорошо видно, что здесь наблюдаются “бухтообразные” изменения этого параметра перед сильными землетрясениями с М ≥ 6.1, для всех выбранных районов. Однако продолжительность и начало аномалий несколько отличаются от кривых, представленных на рис. 10 и рис. 11. В целом, аномалии lgN несколько запаздывают, т.е. смещены вправо по отношению к началу возникновения аномалий параметров b и lgE2/3. Начало аномалий lgN начинается в 2012 г., и они продолжаются до момента возникновения сильных землетрясений в этих зонах. Наиболее ярко картина аномалии lgN наблюдается для кривой (а), где выделяется четкая аномалия продолжительностью 5 лет с 2012 по 2017 гг. Принимая во внимание то, что при построении временных рядов lgN для зоны ЛР (а) использовались землетрясения с М ≥ (Мс = 2.9), можно сказать, что в аномальный период времени среднее число землетрясений с М ≥ 2.9, происходящих в окне сглаживания в окрестностях о. Лесбос, составляло 15 (lgN = 1.18). В то же время, среднее число землетрясений таких магнитуд в этом районе вне аномального периода равно, примерно, 50 (lgN = 1.7). Из этого следует, что в аномальный период времени наблюдалось значительное уменьшение (затишье) количества землетрясений с М ≥ 2.9.

Подводя итог, хотелось бы отметить, что параметр lgN является одной из наиболее независимых физических величин в сейсмологии, с простым алгоритмом вычисления статистик, для которого методические ошибки вычисления маловероятны, и поэтому наблюдаемые аномалии этого параметра явно свидетельствуют о выделении реальных предвестников сильных землетрясений. Более того, корреляция трех сейсмических параметров говорит об обнаружении объективных геодинамических процессов, происходящих в земной коре перед сильными землетрясениями.

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

Анализ сейсмичности на восточном побережье Эгейского моря показал, что перед катастрофическим землетрясением с M = 6.7, произошедшим 30.10.2020 г. в окрестностях о. Самос, наблюдались значимые аномальные изменения сейсмических параметров b, lgE2/3 и lgN. На временных рядах (рис. 4, рис. 5 и рис. 8), построенных по данным о землетрясениях с магнитудой М ≥ Мс, где Мс = 3.6, 3.8 и 3.2 соответственно в диапазонах глубин ΔH = 0–45 км, ΔH1 = 0–23 км. и ΔH2 = 0–18 км, перед всеми землетрясениями с магнитудами М ≥ 5.5, включая и катастрофическое землетрясение вблизи о. Самос с М = 6.7, наблюдаются значимые аномалии параметров b и lgE 2/3 “бухтообразной” формы противоположных знаков. Обнаружено, что аномалии кривых b, lgE 2/3 и lgN, полученные в пределах глубин очагов подготовки сильных событий, имеют большие амплитуды, чем вне этих глубин. Можно утверждать, что наблюдаемое закономерное изменение амплитуд аномалий сейсмических параметров с глубиной может быть использовано для оценки глубины очага будущего сильного землетрясения и всей прогностической ситуации в регионе в целом.

Как было упомянуто выше, активная тектоника на исследуемой территории контролируется движением Эгейской плиты с северо-востока на юго-запад (рис. 1). Принимая во внимание такую геотектонику, интересно рассмотреть динамику возникновения трех сильных землетрясений с М ≥ 6.1, произошедших в исследуемом регионе за истекшее десятилетие.

Совместное рассмотрение временных вариаций параметра b для трех исследуемых районов (рис. 9), выявляют как общие, так и индивидуальные особенности, присущие всем кривым, представленным на рис. 10. Так, можно заметить наличие характерных двух горбов на всех кривых, сдвинутых по времени сверху вниз от ЛР (а) до ЛР (в). Первый “горб” можно выделить между 1 и 3 штрих-пунктирными линиями, в то время как второй “горб” хорошо прослеживается перед всеми сильными землетрясениями в ЛР-нах и отмечен штрих-пунктирной линией 4. Некоторое волнообразное искажение первого “горба” на рис. 10 (в) могло быть вызвано локальными процессами, происходящими в окрестности ЛР (в). Картина, описанная выше, наводит на мысль, что наблюдаемые временные вариации b отражают единый процесс, имеющий волнообразную форму, который начался ранее в 2005 г. и продолжался вплоть до конца 2020 г. Так, на кривой (а) выделяется характерный небольшой минимум в течение 2009 г., на кривой (б) похожий минимум можно выделить в 2011 г., и, наконец, аналогичный минимум наблюдается на кривой (в), но с относительно большей амплитудой, выделяется в 2014 г. Аналогичную картину, но менее ярко выраженную, можно заметить и на предыдущем минимуме отмеченной штрих-пунктирной линией 1. Можно предположить, что наблюдаемый фазовый сдвиг между кривыми может быть объяснен развитием деформационного процесса, происходящего в направлении с севера на юг (рис. 10).

Используя расстояние между центрами ЛР (а) и ЛР (в), равное примерно 230 км, и время распространения деформационного фронта между этими районами, равное примерно 5 годам, можно рассчитать приблизительную скорость распространения фронта деформационной волны по компоненте с ССЗ на ЮЮВ, которая в этом случае равна 40–50 км/год. Точный азимут вектора распространения фронта волны определить в данном эксперименте невозможно, по причине отсутствия второй компоненты, а детальный анализ выходит за рамки настоящей работы. Однако интересно отметить, что такая же величина скорости распространения фронта деформационной волны 40–50 км/год в направлении СВ–ЮЗ получена и для Гармского полигона в Таджикистане [Попандопуло, 2020]. Подводя итог, следует отметить, что наблюдаемая закономерность распространения аномалии b по направлению расположения трех независимых районов, которая показывает фазовый сдвиг с одновременной трансформацией ее формы в пространстве и времени, свидетельствует о реальности наблюдаемого геодинамического процесса.

Интересно также рассмотреть динамику возникновения сильных землетрясений, произошедших в исследуемых районах от острова Лесбос до о. Кос (рис. 1 и рис. 9). Кривые, представленные на рис. 11–рис. 13 , показывают, что первоначально предвестники возникают на севере в ЛР (а) и распространяются на юг в направлении к ЛР (в), где с разницей в один месяц возникли сильные землетрясения с М = 6.1 (12.06.2017 г.) и М = 6.2 (20.07.2017 г.) соответственно (рис. 9). В дальнейшем, по прошествии трех лет, произошло землетрясение с М = 6.7 (30.10.2020 г.) в ЛР (б), с эпицентром, расположенным практически посередине между упомянутыми островами. Такая последовательность возникновения сильных событий может быть интерпретирована на основе концепции сейсмической бреши (seismic gap) первого рода, которая объясняет поведения плиты, стремящейся скользить как жесткое целое вдоль границы поверхности раздела с примыкающей плитой [Федотов, 1965; Касахара, 1985; Моги, 1988]. На основе этой модели в результате произошедших последовательно первых двух событий вблизи о. Лесбос и о. Кос, образовалась сейсмическая брешь между ними, в окрестности о. Самос, которая в дальнейшем реализовалась посредством возникновения землетрясения с М = 6.7.

Интересно отметить, что временной интервал между двумя землетрясениями с М = 6.1 (12.06.2017 г.), и М = 6.2 (20.07.2017 г.), произошедшими в районах (а) и (в) и землетрясением с М = 6.7 (30.10.2020 г.), произошедшим в районе (б), составляет около 3 лет (рис. 9). Можно предположить, что в результате высвобождения тектонической энергии в зоне, где произошло сильное землетрясение с М = 6.7 (30.10.2020 г.), брешь перестала существовать. Однако, согласно концепции сейсмической бреши, это могло привести к образованию новой области с сейсмической брешью южнее о. Кос, в районе о. Родос (рис. 9), т.е. по направлению движения Эгейской плиты. Из этих рассуждений следует, что в ближайшие 2–3 года возможно возникновение сильного землетрясения с магнитудой М ≈ 6.0–6.5 в районе о. Родос, тем более, что в этой области такие землетрясения имели место быть. Последнее из них произошло 25 лет назад 20.07.1996 г. с М = 6.1.

Интересно также провести анализ длительности временного интервала аномалий, наблюдаемые перед сильными событиями. Кривые, представленные на рис. 10–рис. 12 показывают, что длительность периода аномалий сейсмических параметров перед сильными землетрясениями в исследуемой зоне с магнитудой М ≥ 6.1 составляет примерно 10–15 лет. Отсюда следует, что процесс накопления тектонических деформаций земной коры в этой зоне перед этими событиями занимал, по крайней мере, 10–15 лет. Тогда можно предположить, что на восточном побережье Эгейского моря в зоне, расположенной между о. Лесбос и о. Кос, возникновение сильных землетрясений с М ≥ 6.1, по крайней мере, на протяжении ближайших 15 лет маловероятно.

7 . ЗАКЛЮЧЕНИЕ

При помощи компьютерной программы FastBEE, проведены детальные исследования временных вариаций сейсмичности в зоне восточного побережья Эгейского моря Греции, где в окрестности о. Самос 30.10.2020 г. произошло катастрофическое землетрясение с М = 6.7. Для описания сейсмичности использовались параметр b закона Гутенберга–Рихтера, величина выделенной сейсмической энергии lgE 2/3 и показатель сейсмической активности lgN.

Результаты анализа показали, что перед исследуемым землетрясением уверенно выделяются аномалии сейсмических параметров, длительностью около 15 лет. Аналогичные аномалии продолжительностью около 4–6 лет наблюдаются перед 5-ю сильными землетрясениями с магнитудой M ≥ 5.5, произошедшими в этом же регионе в течение последних 50 лет.

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

Исследования временных вариаций сейсмичности в трех независимых сейсмоактивных регионах, в окрестностях о. Лесбос, о. Самос и о. Кос показывают, что в период с 2010 по 2020 гг. перед тремя сильными землетрясениями с магнитудой М ≥ 6.1, произошедшими в этих регионах, наблюдаются предвестники исследуемых параметров. Пространственно-временная последовательность возникновения этих событий объясняется на основе концепции образования сейсмической бреши в окрестности о. Самос [Федотов, 1965; Касахара, 1985; Моги, 1988].

Делается предположение о возможном детектировании деформационной волны в зоне восточного побережья Эгейского моря в период с 2005 по 2020 гг., фронт которой распространяется по направлению с ССЗ на ЮЮВ, со скоростью 40–50 км/год.

Выдвигаются соображения о том, что в ближайшие 10–15 лет в этом регионе на участке между о. Лесбос и о. Кос (рис. 9) возникновение сильных землетрясений с М ≥ 6.1 маловероятно. В то же время не исключено, что в ближайшие 2–3 года южнее о. Кос, в районе о. Родос, возможно возникновения сильного землетрясения с М ≈ 6.0–6.5.

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

  1. Дещеревский А.В., Сидорин А.Я. Изменения представительности каталога землетрясений Греции в пространстве и времени // Сейсмические приборы. 2011. Т. 47. № 4. С. 64–80.

  2. Касахара К. Механика землетрясений. М.: Мир. 1985. 246 с.

  3. Лукк А.А., Попандопуло Г.А. Надежность определения параметров распределения Гутенберга–Рихтера для слабых землетрясений Гармского района в Таджикистане // Физика Земли. 2012. № 9. С. 31–55.

  4. Моги К. Предсказание землетрясений. М.: Мир. 1988. 376 с.

  5. Нерсесов И.Л., Пономарев В.С. Тейтельбаум Ю.М. Эффект сейсмического затишья при больших землетрясениях. Исследования по физике землетрясений. М.: Наука. 1976. С. 149–169.

  6. Попандопуло Г.А., Баскутас Я. Закономерности временных вариаций сейсмически параметров и возможность предсказания сильных землетрясений на территории Греции // Физика Земли. 2011. № 11. С. 27–38.

  7. Попандопуло Г.А., Лукк А.А. Изменения с глубиной параметра b-value в магнитудно-частотном распределении землетрясений Гармского района (Таджикистан) // Физика Земли. 2014. № 2. С. 124–140.

  8. Попандопуло Г. А., Баскутас Я., Хатзииоанну Э. Пространственно-временной анализ представительной магнитуды и наклона графика повторяемости по данным каталога землетрясений Греции // Физика Земли. 2016. № 1. С. 1–17.

  9. Попандопуло Г.А. Детальные исследования временных вариаций параметра b-value закона Гутенберга–Рихтера по данным высокоточных сейсмических наблюдений на Гармском полигоне в Таджикистане // Физика Земли. 2018. № 4. С.79–99.

  10. Попандопуло Г.А. Пространственно временные вариации параметра b-value закона Гутенберга–Рихтера в зависимости от глубины и латерального положения в земной коре Гармском районе Таджикистана // Физика Земли. 2020. № 3. С. 52–73.

  11. Садовский М.А., Писаренко В.Ф. О зависимости времени подготовки землетрясения от его энергии // Докл. АН СССР. 1983. Т. 271. № 2. С. 330–333.

  12. Смирнов В.Б. Опыт оценки представительности каталогов землетрясений // Вулканология и сейсмология. 1997. № 4. С. 93–105.

  13. Тейтельбаум Ю.М., Пономарев В.С. Сейсмический режим и сильные землетрясения. Экспериментальная сейсмология. М.: Наука. 1983. С. 88–98.

  14. Федотов С.А. Закономерности распределения сильных землетрясений Камчатки, Курильских островов и Северо-Восоточной Японии. Сейсмическое микрорайонирования // Вопросы инженерной сейсмологии. 1965. Вып. 10. № 36(293). С. 66–93.

  15. Aki K. Maximum likelihood estimate of b in the formula lgN = a – bM and its confidence limits // Bull. Earhq. Res. Inst Univ. Tokyo. 1965. V. 43. P. 237–239.

  16. Bath M. Spectral Analysis in Geophysics. Amsterdam: Elsevier. 1974. 563 p.

  17. Baskoutas I., D’Alessandro A. Study of the seismicity temporal variation for the current seismic hazard evaluation in Val d’Agri, Italy // Nat. Hazards Earth Syst. Sci. 2014. V. 14. P. 3169–3174. https://doi.org/10.5194/nhess-14-3169-2014

  18. Box G.E.P., Jenkins G.M. Time Series Analysis: Forecasting and Control. San Francisco: Holden Day. 1970. 406 p.

  19. Chouliaras G., Melis N.S., Drakatos G., Makropoulos K. Operational network improvements and increased reporting in the NOA (Greece) seismicity catalog // Adv. Geosci. 2013. V. 36. P. 7–9. https://doi.org/10.5194/adgeo-36-7-2013

  20. D’Alessandro A., Papanastassiou D., Baskoutas I. Hellenic Unified Seismological Network: an evaluation of its performance through SNES method // Geophys J. Int. 2011. V. 185. P. 1417–1430.

  21. Gutenberg B., Richter Ch.F.Frequency of earthquakes in California // Bull. Seismol. Soc. Am. 1944. V. 34. P. 185–188.

  22. Gutenberg B., Richter C.F. Earthquake Magnitude, Intensity, Energy and Acceleration (second paper) // Bull. Seismol. Soc. Amer. 1956. V. 46. № 2. P. 105–145.

  23. Ishimoto M., Iida K. Observations of earthquakes registered with the microseismograph constructed recently // Bull. Earthquake Res. Inst. Tokyo Univ. 1939. V. 17. P. 443–478.

  24. Jackson J. Active tectonics of the Aegean region // Annu. Rev. Earth Planet. Sci. 1994. V. 22. P. 239–271.

  25. Karakostas V.G., Papadimitriou E.E., Tranos M.D., Papazachos C.B. Active seismotectonic structures in the area of chios island, north aegean sea, revealed from microseismicity and fault plane solutions. Bull. Geol. Soc. of Greece, Pros 12. Inter. Cong. 2010.

  26. Keylis-Borok V. On estimation of the Displacement in an Earthquake Source and of Source Dimentions // Ann. Geof. 1959. V. 12. № 2.

  27. Marzocchi W., Sandri L. A review and new insights on the estimation of the b-value and its uncertainty // Ann. Geophys. 2003. V. 46. P. 1271–1282.

  28. Mignan A., Woessner J. Completeness magnitude in earthquake catalogs. Community Online Resource for Statistical Seismicity Analysis. 2012. https://doi.org/10.5078/corssa-00180805

  29. Mignan A., Chouliaras G. Fifty Years of Seismic Network Performance in Greece (1964–2013): Spatiotemporal Evolution of the Completeness Magnitude // Seism. Res. Let. 2014. V. 85. № 3. P. 657–667.

  30. Mogi K. Study of the elastic shocks caused by the fracture of heterogeneous materials and its relation to earthquake phenomena // Bull. Earthq. Res. Inst. Tokyo Univ. 1962. V. 40. P. 125–173.

  31. Mori J., Abercrombie R.E. Depth dependence of earthquake frequency-magnitude distributions in California: implications for the rupture initiation // J. Geophys. Res. 1997. V. 102. P. 15081–15090.

  32. McKenzie. D. Active tectonics of the Mediterranean region // Geophys. J. Roy. Astr. Soc. 1972. V. 30. 2. P. 109–185. https://doi.org/10.1111/j.1365-246X.1972.tb02351.x

  33. McKenzie D. Active tectonics of the Alpine–Himalayan belt: the Aegean Sea and surrounding regions // Geophys. J. Roy. Astr. Soc. 1978. V. 55. № 1. P. 217–254. https://doi.org/10.1111/j.1365-246X.1978.tb04759.x

  34. Papanastassiou D. Detection-location capability of the Hellenic Unified Seismological Network (HUSN) operating by the Institute of Geodynamics, National Observatory of Athens // Hellenic J. Geosci. 2011. V. 45. P. 209–216.

  35. Papadopoulos G.A., Baskoutas I.G. New tool for the temporal variation analysis of seismic parameters // Nat. Hazards Earth Syst. Sci. 2009. V. 9. P. 859–864. https://doi.org/10.5194/nhess-9-859-2009

  36. Papazachos B.C., Comninakis P.E. Geophysical and tectonic features of the Aegean arc // J. Geophys. Res. 1971. V. 76. P. 18517–8533.

  37. Papazachos B.C., Delibasis N.D. Tectonic stress field and seismic faulting in the area of Greece // Tectonophysics. 1969. V. 7. P. 231–255. https://doi.org/10.1016/0040-1951(69)90069-9

  38. Papazachos B.C., Papadimitriou E.E., Kiratzi A.A., Papazachos C.B., Louvari E.K. Fault plane solutions in the Aegean Sea and the surrounding area and their tectonic implications // Boll. Geof. Teor. Appl. 1998. V. 39. P. 199–218.

  39. Papazachos C.B. Seismological and GPS evidence for the Aegean-Anatolia interaction // Geophys. Res. Let. 1999. V. 26. № 17. P. 2653–2656.

  40. Popandopoulos G.A., Chatziioannou E. Gutenberg–Richter Law Parameters Analysis Using the Hellenic Unified Seismic Network Data Through FastBEE Technique // Earth Science. 2014. V. 3. № 5. P. 122–131. https://doi.org/10.11648/jearth.20140305.12

  41. Sandri L., Marzocchi W. A technical note on the bias in the estimation of the b-value and its uncertainty through the Least Squares technique // Ann. Geophys. 2007. V. 50. № 3. P. 329–339.

  42. Scordilis E.M., Karakaisis G.F., Karakostas B.G., Panagiotopoulos D.G., Comninakis P.E., Papazachos B.C. Evidence for transform faulting in the Ionian Sea: The Cephalonia Island earthquake sequence // Pure Appl. Geophys. 1985. V. 123. P. 388–397.

  43. Scholz C.H. The Frequency-Magnitude Relation of Microfracturing in Rock and its Relation to Earthquakes // Bull. Seism. Soc. Am. 1968. V. 58. P. 399–415.

  44. Schorlemmer D., Wiemer S., Wyss M. Earthquake statistics at Park field: 1. Stationarity of b-values // J. Geophys. Res. 2004a V. 109. B12. https://doi.org/10.1029/2004JB003234

  45. Schorlemmer D., Wiemer S., Wyss M., Jackson D.D. Earthquake statistics at Park field: 2. Probabilistic forecasting and testing // J. Geophys. Res. 2004b. V. 109. B12. https://doi.org/10.1029/2004JB003235

  46. Schorlemmer D., Wiemer S., Wyss M. Variations in earthquake-size distribution across different stress regimes // Nature. 2005. V. 437. P. 539–542. https://doi.org/10.1038/nature04094

  47. Shi Y., Bolt B.A. The standard error of the Magnitude-frequency b-value // Bull. Seism. Soc. Am. 1982. V. 72. P. 1677–1687.

  48. Tan O., Papadimitriou E.E., Pabucçu Z., Karakostas V., Yörük A., Leptokaropoulos K. A detailed analysis of microseismicity in Samos and Kusadasi (Eastern Aegean Sea) areas // Acta Geophys. 2014. V. 62. № 6. P. 1283–1309. https://doi.org/10.2478/ s11600-013-0194-1

  49. Utsu T. A method for determining the value of b in a formula lgN = a – bM showing the magnitude frequency for earthquakes // Geophys. Bull. Hokkaido Univ. 1965. V. 13. P. 99–103.

  50. Wiemer S., Wyss M. Minimum magnitude of complete reporting in earthquake catalogs: examples from Alaska, the western United States, and Japan // Bull. Seismol. Soc. Am. 2000. V. 90. P. 859–869.

  51. Wiemer S., Wyss M. Mapping spatial variability of the frequency – magnitude distribution of earthquakes // Adv. Geophys. 2002. V. 5. P. 259–302.

  52. Woessner J., Wiemer S. Assessing the quality of earthquake catalogues: Estimating the magnitude of completeness and its uncertainty // Bull. Seismol. Soc. Am. 2005. V. 95. P. 684–698.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Физика Земли

Архивы выпусков Информация о журнале Отправить рукопись в журнал