Моделирование пространственного распределения областей повышенных предсейсмических деформаций

Полный текст

Введение

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

В настоящее время разработано множество моделей подготовки землетрясений: дилатантно-диффузная [1], лавинно-неустойчивого трещинообразования [2], консолидации и деструкции [3], неустойчивого скольжения [4], фазовых превращений [5], каскадного подъема и прогрессивной локализации [6], энергетическая модель [7] и другие. Хотя эти модели в различной степени объясняют природу геофизических аномалий, большинство из них носит феноменологический характер. Отсутствие развитого математического аппарата ограничивает их применение для моделирования напряженного состояния земной коры вокруг очагов реальных землетрясений. Наиболее разработанными в математическом отношении являются энергетическая модель, а также модели консолидации и деструкции.

Значительный вклад в развитие моделей консолидации и деструкции внесли работы И. П. Добровольского [3, 8]. В его модели деформации земной коры рассматриваются в приближении упругого однородного изотропного полупространства, где очаг формирующегося землетрясения представлен неоднородным сферическим включением с отличающимся модулем сдвига, на границах которого задано давление. Применение теории малых возмущений позволило свести решение к объемному интегралу от функций Грина для системы фиктивных сил внутри включения. Напряженно-деформированное состояние получено как разность состояний среды с включением и без него.

Дальнейшее развитие этот подход получил в работах Ю. Л. Ребецкого и А. С. Лермонтовой [9, 10], где учтена неоднородность строения земной коры с областями, проявляющими пластические и квазипластические свойства. Авторы представили очаг в виде неоднородного цилиндрического включения и решили задачу для плоского напряженного состояния в упругом и упруго-пластическом случаях. Показано, что в областях с аномальными включениями затухание деформационных возмущений происходит медленнее по сравнению с упруго-однородной средой.

Различные модификации энергетической модели представлены в работах А. С. Пережогина, Б. М. Шевцова, В. А. Салтыкова, Ю. А. Кугаенко. Авторы использовали различные силовые эквиваленты дислокационных процессов: точечные источники (единичная и двойная силы) для моделирования выступов в зоне субдукции [11, 12], комбинацию двойных пар сил с моментами [13], распределенные силовые источники, соответствующие классической дислокационной модели [14, 15]. Результаты моделирования сопоставлялись с данными наблюдений высокочастотного сейсмического шума, акустической эмиссии, поверхностных смещений и деформаций осадочных пород.

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

В данной работе предложен новый подход к идентификации областей с повышенной частотой возникновения предсейсмических деформаций. Метод апробирован на землетрясениях в районе Камчатского полуострова за 1976–2020 гг. Выбор региона обусловлен его высокой сейсмической активностью и расположением в Тихоокеанском огненном кольце. Камчатка давно служит естественным полигоном для изучения сейсмотектонических процессов, связанных с накоплением и релаксацией напряжений при взаимодействии литосферных плит.

1. Физическая постановка задачи

Очаг тектонического землетрясения представляет собой разрыв сплошности материала Земли, возникающий под действием упругих сдвиговых напряжений, накопленных в процессе длительной тектонической деформации. В момент землетрясения происходит полное или частичное снятие накопленных напряжений в его очаговой области. Согласно классическому определению Б. В. Кострова, основанному на теории упругой отдачи Г. Ф. Рейда, тектоническое землетрясение представляет собой разрыв скольжения, при котором взаимное перемещение его берегов в направлении, перпендикулярном поверхности разрыва, равно нулю. При этом процесс землетрясения сопровождается преобразованием части накопленной упругой потенциальной энергии в кинетическую энергию сейсмических волн.

Деформационные процессы, сопровождающие подготовку землетрясения, обусловлены приращением потенциальной энергии упругих деформаций $\Delta W$ в очаговой области. Важно отметить, что высвободившаяся при землетрясении сейсмическая энергия $E$ составляет лишь часть этой накопленной энергии. Отношение этих величин, называемое сейсмическим КПД $\eta$, является важной характеристикой эффективности процесса снятия напряжений:
\[ \begin{equation*}
\eta = {E}/{\Delta W} .
\end{equation*} \]

В рамках механики сплошных сред при изучении процессов потери устойчивости традиционно рассматривают случай неизменных свойств материала и сохранения его сплошности. В этом приближении очаг землетрясения может быть описан эквивалентной системой сил, распределенных по поверхности разрыва, что непосредственно следует из теоремы взаимности Бетти [16]. Для наиболее общего случая произвольно ориентированного разрыва смещений в изотропной упругой среде используется система, состоящая из девяти пар двойных сил с моментами [17]. Данная система, схематически изображенная на рис. 1, представляет собой минимально необходимую систему для построения адекватного силового эквивалента очага в общем случае.

Рис. 1. Схематическое изображение девяти пар сил, необходимых для построения силового эквивалента произвольно ориентированного разрыва смещений в сплошной среде
[Figure 1. Schematic representation of nine force couples required to construct a force equivalent for an arbitrarily oriented displacement discontinuity in a continuum medium]

При анализе деформационных процессов в земной коре важно учитывать, что фоновый уровень деформаций, обусловленный преимущественно приливными воздействиями Луны и Солнца, характеризуется величинами относительных деформаций порядка $10^{-8}$. В связи с этим в дальнейшем будем рассматривать только те деформации, которые превышают указанный фоновый уровень. Такие деформации будем называть повышенными. Большинство известных предвестников выявляются в сигналах, источники которых расположены на глубине много меньшей, чем соответствующий размер земной коры. Поэтому особое внимание в работе уделяется анализу деформационных процессов, происходящих в непосредственной близости от дневной поверхности Земли.

2. Математическая постановка задачи

Рассмотрим упрощенную геометрическую модель земной коры в виде упругого изотропного полупространства. Деформационные процессы в такой системе описываются системой дифференциальных уравнений Ламе:
\[ \begin{equation}
\mu u_{i jj} + (\lambda+\mu)u_{j ji} + X_i = 0,\quad i=1,2,3,
\end{equation} \tag{1} \]
где $u_i$ — компоненты вектора перемещений; $\lambda$, $\mu$ — коэффициенты Ламе; $X_i$ — компоненты вектора массовых сил; индексами после запятой обозначено дифференцирование по соответствующим пространственным координатам.

Примем, что упругое полупространство занимает область $x_3\leqslant 0$, а поверхность Земли соответствует плоскости $x_3=0$. На этой поверхности выполняются условия отсутствия напряжений в направлении оси $x_3$:
\[ \begin{equation}
\sigma_{31}\big|_{x_3=0}=\sigma_{32}\big|_{x_3=0}=\sigma_{33}\big|_{x_3=0}=0.
\end{equation} \tag{2} \]

Дополнительно потребуем, чтобы напряжения, создаваемые очагом формирующегося землетрясения, затухали на бесконечности:
\[ \begin{equation}
\lim_{x_1\to\pm\infty} \sigma_{ij} = \lim_{x_2\to\pm\infty} \sigma_{ij} =\lim_{x_3\to-\infty} \sigma_{ij} = 0.
\end{equation} \tag{3} \]

Для описания силового воздействия, соответствующего системе сил, показанной на рис. 1, введем плотность распределения массовых сил $\boldsymbol{X}$. Будем считать, что все силы приложены в точке с координатами $(\xi_1, \xi_2, \xi_3)$ внутри упругого полупространства. В случае сосредоточенной силы, направленной вдоль оси $j$ с интенсивностью $p$, вектор массовых сил принимает вид
\[ \begin{equation}
\boldsymbol{X}^j = p\delta(\boldsymbol{x} - \boldsymbol{\xi})\boldsymbol{e}_j,
\end{equation} \tag{4} \]
где $\delta(\boldsymbol{x} - \boldsymbol{\xi})=\delta(x_1-\xi_1)\delta(x_2-\xi_2)\delta(x_3-\xi_3)$ — трехмерная дельта-функция Дирака, $\boldsymbol{e}_j$ — единичный базисный вектор.

На основании соотношения (4) построим выражения для двойных сил. Рассмотрим сначала двойную силу, направленную вдоль оси $x_1$ с интенсивностью $p_{11}$. Как показано на рис. 2, a, такая система может быть представлена как комбинация двух сосредоточенных сил $\boldsymbol{X}^1_1$ и $\boldsymbol{X}^1_2$, приложенных в точках $(\xi_1+\Delta\xi_1, \xi_2, \xi_3)$ и $(\xi_1-\Delta\xi_1, \xi_2, \xi_3)$ соответственно, при $\Delta\xi_1\to 0$:
\[ \begin{equation*}
\begin{aligned}
\boldsymbol{X^1_1} &= p_{11}\delta(x_1-(\xi_1+\Delta\xi_1),x_2- \xi_2, x_3-\xi_3)\boldsymbol{e}_1,\\
\boldsymbol{X^1_2} &= p_{11}\delta(x_1-(\xi_1-\Delta\xi_1),x_2- \xi_2, x_3-\xi_3)\boldsymbol{e}_1,\\
\boldsymbol{X} &= \lim\limits_{\Delta\xi_1\to 0}\frac{\boldsymbol{X}^1_1-\boldsymbol{X}^1_2}{2\Delta\xi_1} = p_{11}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_1}\boldsymbol{e}_1.
\end{aligned}
\end{equation*} \]

Рис. 2. Схематическое изображение комбинации сосредоточенных сил, необходимой для получения двойной силы вдоль оси $x_1$ (a) и двойной силы вдоль той же оси, но с моментом относительно оси $x_3$ (b)
[Figure 2. Schematic representation of a combination of concentrated forces to obtain a double force along the $x_1$ axis (a) and a double force along the same axis, but with a moment relative to the $x_3$ axis (b)]

Аналогичным образом, для двойной силы вдоль оси $x_1$ с моментом относительно оси $x_3$ интенсивности $p_{12}$ (рис. 2, b) используем комбинацию сил $\boldsymbol{X}^1_1$ и $\boldsymbol{X}^2_1$, приложенных в точках $(\xi_1, \xi_2+\Delta\xi_2, \xi_3)$ и $(\xi_1, \xi_2-\Delta\xi_2, \xi_3)$, при $\Delta\xi_2\to 0$:
\[ \begin{equation*}
\begin{aligned}
\boldsymbol{X^1_1} &= p_{12}\delta(x_1-\xi_1,x_2- (\xi_2+\Delta\xi_2), x_3-\xi_3)\boldsymbol{e}_1,\\
\boldsymbol{X^1_2} &= p_{12}\delta(x_1-\xi_1,x_2- (\xi_2-\Delta\xi_2), x_3-\xi_3)\boldsymbol{e}_1,\\
\boldsymbol{X} &= \lim\limits_{\Delta\xi_2\to 0}\frac{\boldsymbol{X}^1_1-\boldsymbol{X}^1_2}{2\Delta\xi_2} = p_{12}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_2}\boldsymbol{e}_1.
\end{aligned}
\end{equation*} \]

По аналогичной схеме могут быть получены все остальные компоненты системы сил, необходимые для описания силового эквивалента произвольно ориентированного разрыва в сплошной среде. В общем случае компоненты вектора массовых сил $\boldsymbol{X}$ выражаются следующим образом:
\[ \begin{equation*}
\begin{aligned}
&X_1 = p_{11}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_1} + p_{12}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_2} + p_{13}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_3} ,\\
&X_2 = p_{21}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_1} + p_{22}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_2} + p_{23}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_3} ,\\
&X_3 = p_{31}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_1} + p_{32}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_2} + p_{33}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_3} .
\end{aligned}
\end{equation*} \]

Полученные выражения могут быть записаны в более компактной форме:
\[ \begin{equation*}
X_i=p_{ij}\frac{\partial\delta(\boldsymbol{x}-\boldsymbol{\xi})}{\partial\xi_j},\quad i=1,2,3,
\end{equation*} \]
где $p_{ij}$ характеризует интенсивность соответствующей пары сил.

3. Аналитические решения

Для задачи (1) с граничными условиями (2) и (3) существуют точные решения в виде функций Грина, впервые полученные Р. Миндлином. Рассмотрим подробно вид этих функций для различных ориентаций точечных сил.

3.1. Функции Грина для полупространства

Для единичной силы, приложенной в точке $(\xi_1, \xi_2, \xi_3)$ упругого полупространства и направленной вдоль оси $x_3$, компоненты функции Грина $\boldsymbol{g}^3(\boldsymbol{x})$ имеют следующий вид:
\[ \begin{equation*}
\begin{aligned}
g_1^3 = \dfrac{(x_1-\xi_1)}{16\pi\mu(1-\nu)}&\Bigl[\dfrac{(x_3-\xi_3)}{r_1^3} + \dfrac{(3-4\nu)(x_3-\xi_3)}{r_2^3} + \dfrac{4(1-\nu)(1-2\nu)}{r_2(r_2-x_3-\xi_3)} + \dfrac{6x_3\xi_3}{r_2^5}\Bigr],\\
g_2^3 = \dfrac{(x_2-\xi_2)}{16\pi\mu(1-\nu)}&\Bigl[\dfrac{(x_3-\xi_3)}{r_1^3} + \dfrac{(3-4\nu)(x_3-\xi_3)}{r_2^3} + {}\\
&\qquad\qquad+\dfrac{4(1-\nu)(1-2\nu)}{r_2(r_2-x_3-\xi_3)}+\dfrac{6x_3\xi_3(x_3+\xi_3)}{r_2^5}\Bigr],\\
g_3^3 = \dfrac{1}{16\pi\mu(1-\nu)}&\Bigl[ \dfrac{(3-4\nu)}{r_1} + \dfrac{5-12\nu+8\nu^2}{r_2} + \dfrac{(x_3-\xi_3)^2}{r_1^3} + {}\\
&\qquad\qquad+ \dfrac{(3-4\nu)(x_3+\xi_3)^2-2x_3\xi_3}{r_2^3}\Bigr],
\end{aligned}
\end{equation*} \]
где $\nu$ — коэффициент Пуассона среды, а величины $r_1$ и $r_2$ определяются так:
\[ \begin{equation*}
\begin{aligned}
&r_1=\sqrt{(x_1-\xi_1)^2+(x_2-\xi_2)^2+(x_3-\xi_3)^2},\\
&r_2=\sqrt{(x_1-\xi_1)^2+(x_2-\xi_2)^2+(x_3+\xi_3)^2}.
\end{aligned}
\end{equation*} \]

Для случая единичной силы, направленной вдоль оси $x_1$, компоненты функции Грина $\boldsymbol{g}^1(\boldsymbol{x})$ выражаются более сложными соотношениями:
\[ \begin{equation*}
\begin{aligned}
g_1^1=&\dfrac{1}{16\pi\mu(1-\nu)}\Bigl\{ \dfrac{3-4\nu}{r_1}+\dfrac{1}{r_2} + \dfrac{(x_1-\xi_1)^2}{r_2^3} + \dfrac{(3-4\nu)(x_1-\xi_1)}{r_2^3} + {}\\
&\qquad \qquad +\dfrac{4(1-\nu)(1-2\nu)[r_2^2-(x1-\xi1)^2-r_2(x_3+\xi_3)]}{r_2(r_2-x_3-\xi_3)^2} \Bigr\},\\
g_2^1=&\dfrac{(x_1-\xi_1)(x_2-\xi_2)}{16\pi\mu(1-\nu)}\Bigl[ \dfrac{1}{r_1^3} + \dfrac{(3-4\nu)}{r_2^3} - \dfrac{6x_3\xi_3}{r_2^5} - \dfrac{4(1-\nu)(1-2\nu)}{r_2(r_2-x_3-\xi_3)^2}\Bigr],\\
g_3^1=&\dfrac{(x_1-\xi_1)}{16\pi\mu(1-\nu)}\Bigl[ \dfrac{(x_3-\xi_3)}{r_1^3} + \dfrac{(3-4\nu)(x_3-\xi_3)}{r_2^3} - \dfrac{4(1-\nu)(1-2\nu)}{r_2(r_2-x_3-\xi_3)} - {}\\
&\qquad\qquad - \dfrac{6x_3\xi_3(x_3+\xi_3)}{r_2^5}\Bigr].
\end{aligned}
\end{equation*} \]

В силу симметрии задачи функция Грина $\boldsymbol{g}^2(\boldsymbol{x})$ для силы вдоль оси $x_2$ может быть получена из $\boldsymbol{g}^1(\boldsymbol{x})$ простой заменой индексов 1 и 2.

3.2. Решение для системы двойных сил

Функции Грина, соответствующие действию двойных сил, могут быть получены дифференцированием основных решений $\boldsymbol{g}^i(\boldsymbol{x})$ по пространственным координатам, то есть в виде $\partial\boldsymbol{g}^i(\boldsymbol{x})/\partial x_j$ [18]. При этом:

• при $i=j$ получаем решение для пары двойных сил, направленных вдоль соответствующей оси;
• при $i\neq j$ — решение для пары двойных сил с моментом относительно оси, отличной от $i$ и $j$.

3.3. Общее решение через формулу Вольтерра

В общем случае поле смещений в упругом полупространстве может быть выражено через интегральное представление Вольтерра [16]:
\[ \begin{equation*}
u_k(\boldsymbol{x}) = \int_{\Sigma} s_i(\boldsymbol{\xi})\sigma_{ij}^k(\boldsymbol{\xi}, \boldsymbol{x})n_j \, d\Sigma,
\end{equation*} \]
где
$s_i(\xi)$ — вектор смещения на поверхности разрыва $\Sigma$, $n_j$ — компоненты единичной нормали к поверхности $\Sigma$.

С учетом закона Гука для однородной изотропной среды формула Вольтерра принимает вид
\[ \begin{equation*}
u_k(\boldsymbol{x}) = \int_{\Sigma}\bigl[\mu(s_pn_q+s_qn_p) + \lambda s_k n_k\delta_{pq}\bigr]\frac{\partial g_k^p(\boldsymbol{x},\boldsymbol{\xi})}{\partial\xi_q}\, d\Sigma =
\int _{\Sigma}m_{pq}\frac{\partial g_k^p(\boldsymbol{x},\boldsymbol{\xi})}{\partial\xi_q}\, d\Sigma,
\end{equation*} \]
где $m_{pq} = \mu(s_pn_q+s_qn_p) + \lambda s_k n_k\delta_{pq}$ — тензор плотности сейсмического момента [17], полностью характеризующий механику очага землетрясения.

Для случая точечного источника решение задачи принимает особенно простой вид:
\[ \begin{equation*}
u_k(\boldsymbol{x}) = M_{pq}\frac{\partial g_k^p(\boldsymbol{x,\boldsymbol{\xi}})}{\partial\xi_q},
\end{equation*} \]
где $M_{pq}$ — тензор сейсмического момента. Данное выражение будет использовано в дальнейшем в качестве основы для численного моделирования.

3.4. Оценка коэффициента усиления деформаций

С учетом квадратичной зависимости между компонентами тензора деформаций $\varepsilon_{ij}$ и плотностью потенциальной энергии упругих деформаций коэффициент усиления напряженно-деформированного состояния при подготовке землетрясения принимается равным $\eta^{-0.5}$. Практическая формула для оценки этого коэффициента, предложенная И. П. Добровольским [3], имеет вид
\[ \begin{equation*}
\eta = 10^{0.26M_W-3.93},
\end{equation*} \]
где $M_W$ — моментная магнитуда землетрясения.

4. Результаты численного моделирования

Для численного моделирования использовались данные из каталога "The Global Centroid-Moment-Tensor Catalog" (https://www.globalcmt.org/), содержащего информацию о механизмах очагов землетрясений. Из каталога отобраны все сейсмические события, зарегистрированные вблизи Камчатского полуострова в период с 1976 по 2020 год. В анализ включены следующие параметры:

• дата и время события;
• координаты эпицентра (широта, долгота);
• глубина гипоцентра;
• магнитуда землетрясения;
• тензор сейсмического момента;
• скалярный сейсмический момент.

Географические координаты эпицентров находились в диапазоне от 49° до 60° с. ш. и от 149.0° до 170.0° в. д. Общий объем выборки составил $N = 877$ землетрясений.

4.1. Методика расчетов

Моделирование деформационных процессов выполнялось для поверхности земной коры ($x_3=0$). Расчетная область покрывала географические координаты [48°–60° с. ш.] и [151°–170° в. д.] с шагом сетки 0.5°. Параметры упругой среды приняты следующими:

• коэффициенты Ламе $\mu=\lambda=3.675\cdot10^{10}$ Н/м$^2$;
• коэффициент Пуассона $\nu=0.25$.

Для каждого землетрясения в узлах расчетной сетки вычислялись предсейсмические деформации по формуле
\[ \begin{equation*}
\varepsilon_{\max} = \frac{1}{2\mu}\sigma_{\max},
\end{equation*} \]
где $\sigma_{\max} = \max\{|\sigma_1-\sigma_2|, |\sigma_2-\sigma_3|, |\sigma_1-\sigma_3|\}$ — максимальное касательное напряжение; $\sigma_1$, $\sigma_2$, $\sigma_3$ — главные значения тензора напряжений.

Для каждого узла сетки $(i,j)$ определялось количество случаев $n_{ij}$, когда деформации превышали фоновый уровень $10^{-8}$, после чего вычислялись относительные частоты $\omega_{ij}=n_{ij}/N$.

4.2. Анализ результатов

Результаты моделирования представлены на рис. 3 и 4. На рис. 3 показано пространственное распределение относительных частот повышенных деформаций за весь исследуемый период (1976–2020 гг.). Наблюдаются следующие закономерности:

• неравномерное распределение частот с максимумами вдоль линии основного разлома у побережья Камчатки;
• наибольшие значения частот в центральной части расчетной области;
• соседство зон максимальных деформаций с районами плотной населенности.

Рис. 3. Распределение относительных частот повышенных деформаций (1976–2020 гг.): точки — эпицентры землетрясений; красная линия — основной разлом
[Figure 3. Distribution of relative frequencies of increased deformations (1976–2020). Dots show earthquake epicenters, red line indicates the main fault]

На рис. 4 представлены результаты моделирования для отдельных годов анализируемого периода. Из рисунка видно, что в 2013 г. (рис. 4, b) и 2016 г. (рис. 4, c) наблюдается увеличение относительных частот повышенных деформаций до 0.6–0.8 в областях с наибольшей концентрацией землетрясений. Суммарная накопленная потенциальная энергия упругих деформаций составила $1.6\times 10^{17}$ Дж (63 землетрясения) в 2013 г. и $4.6\times 10^{15}$ Дж (23 землетрясения) в 2016 г.

В 2005 г. (рис. 4, a) и 2020 г. (рис. 4, d) значения относительных частот оказались существенно ниже. По нашему мнению, это обусловлено большей пространственной разреженностью землетрясений и меньшей суммарной энергией деформаций, которая составила $1.3\times10^{14}$ Дж (22 землетрясения) и $5.8\times 10^{14}$ Дж (27 землетрясений) соответственно.

Рис. 4. Распределение относительных частот по годам: 2005 (a), 2013 (b), 2016 (c), 2020 (d): размер кружков соответствует магнитуде $M_W$; красная линия — основной разлом
[Figure 4. Distribution of relative frequencies by year: 2005 (a), 2013 (b), 2016 (c), 2020 (d): circle size corresponds to moment magnitude $M_W$; red line indicates the main fault]

Заключение

В рамках линейной теории упругости предложен новый подход к моделированию пространственного распределения областей повышенных деформаций, возникающих при подготовке землетрясений. Для каждой точки исследуемой поверхности земной коры выполняется моделирование предсейсмических деформаций с последующим подсчетом случаев превышения порогового значения $10^{-8}$. На основе полученных данных рассчитываются относительные частоты возникновения таких деформаций.

Предложенный метод апробирован на данных каталога "The Global Centroid-Moment-Tensor Catalog" для Камчатки — одного из самых сейсмически активных регионов планеты. Установлено, что наиболее часто аномальные деформации возникают вдоль линии основного разлома у побережья Камчатки, причем максимумы относительных частот граничат с густонаселенными районами.

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

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

Конкурирующие интересы. Авторы заявляют об отсутствии конфликта интересов в отношении авторства и публикации данной статьи.
Авторский вклад и ответственность. Все авторы внесли равный вклад в разработку концепции статьи и написание рукописи. Авторы несут полную ответственность за предоставление окончательной версии рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.
Финансирование. Исследование выполнено при поддержке Государственного задания Института космофизических исследований и распространения радиоволн Дальневосточного отделения Российской академии наук (ИКИР ДВО РАН), регистрационный номер темы исследования 124012300245-2.

Об авторах

Максим Игоревич Гапеев

Институт космофизических исследований и распространения радиоволн ДВО РАН

Автор, ответственный за переписку.
Email: gapeev.sci@yandex.ru
ORCID iD: 0000-0001-5798-7166
SPIN-код: 1393-2315
Scopus Author ID: 57212685382
https://www.mathnet.ru/person139714

младший научный сотрудник; лаб. акустических исследований

Россия, 684034, Камчатский край, с. Паратунка, ул. Мирная, 7

Юрий Валентинович Марапулец

Институт космофизических исследований и распространения радиоволн ДВО РАН

Email: marpl@ikir.ru
ORCID iD: 0000-0002-3030-9944
SPIN-код: 2976-5061
Scopus Author ID: 15725296200
https://www.mathnet.ru/person115685

доктор физико-математических наук, доцент; директор

Россия, 684034, Камчатский край, с. Паратунка, ул. Мирная, 7

Александра Андреевна Солодчук

Институт космофизических исследований и распространения радиоволн ДВО РАН

Email: aleksandra@ikir.ru
ORCID iD: 0000-0002-6761-8978
SPIN-код: 5162-3730
Scopus Author ID: 55533915300
https://www.mathnet.ru/person116312

кандидат физико-математических наук; старший научный сотрудник; лаб. акустических исследований

Россия, 684034, Камчатский край, с. Паратунка, ул. Мирная, 7

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

Дополнительные файлы