Сейсмические волновые поля в сферически-симметричной Земле с высокой детальностью. Аналитическое решение

Полный текст

ВВЕДЕНИЕ

В работе рассмотрена задача построения аналитического решения для сейсмических волновых полей в слоистом шаре. Эта задача рассматривалась во многих работах (например, [1]). В данной работе показано, что эту задачу нельзя считать окончательно решенной для высокочастотных волновых полей.

Решение этой задачи строится следующим образом. На первом этапе применяется преобразование Фурье-Лежандра [1]. В итоге исходная постановка для распространения сейсмических волн сводится к двухпараметрическому семейству краевых задач для системы обыкновенных дифференциальных уравнений. Следуя [2], на втором этапе в каждом сферическом слое искомая краевая задача второго порядка сводится к двум задачам Коши первого порядка в спектральной области. Это позволило получить в явном виде аналитическое решение. В него входят функции Бесселя. Поскольку они быстро стремятся к нулю и бесконечности, в решении возникает неустойчивость. Для устранения неустойчивости нужно использовать асимптотику. В работе используется оригинальная асимптотика, разработанная ранее [3, 4].

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

В настоящее время активно развиваются конечно-разностные методы для расчета волновых полей в Земле. Несмотря на достижения в области высокопроизводительных вычислений, численное моделирование распространения высокочастотных (например, 1 Гц или выше) сейсмических волн в глобальном масштабе по-прежнему невозможно [5]. Чтобы преодолеть эту трудность, развивается гибридный (комбинированный) метод для расчета телесейсмических волновых полей. Однако при гибридном моделировании возникают погрешности [6].

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

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

$V_p=\sqrt{(k_S+\frac{4}{3}\mu )/\rho }$,

то в данном случае отличие от нуля модуля сдвига существенно влияет на характер изменения скорости продольных сейсмических волн, особенно в низах внешнего ядра Земли [7, 8]. Здесь k_S – адиабатический модуль всестороннего сжатия, ρ – плотность.

В настоящей работе приведены результаты расчета для дискретной (слоистой) модели Земли [7, 8] с несущей частотой в 1 герц. Модель Земли характеризуется тем, что в ней внешнее ядро обладает эффективным модулем сдвига, отличным от нуля. В результате аналитического расчета для дискретной модели Земли [7, 8] обнаружено, что впереди PKP-волн возникают высокочастотные колебания небольшой амплитуды, так называемые “предвестники”.

Аналитический расчет показал, что теоретические сейсмограммы для модели Земли [7, 8] во многом похожи на экспериментальные сейсмограммы, полученные мировой сетью сейсмических станций. Это сравнение, с одной стороны, подтверждает справедливость выводов о природе предвестников сейсмических волн, полученных в [4]. С другой – подтверждает правильность идей, положенных в основу построения модели Земли в [7, 8].

ПОСТАНОВКА ЗАДАЧИ

Математическая постановка задачи моделирования сейсмических волн формулируется в сферической системе координат $(0\le r\le R\mathrm{,0}\le \theta \le \pi \mathrm{,0}<\phi \le 2\pi )$ следующим образом: определить вектор смещения $\widehat{u}=(u_r,u_{\theta },u_{\phi })$ из хорошо известных динамических уравнений упругости (например, [1]).

В качестве краевого условия приложим в точке $r=R_1, \theta =0$, сосредоточенное воздействие типа радиальной осесимметричной силы:

${\sigma }_{rr}=\frac{\delta \left(\theta \right)}{R_1^2\sin \theta }f\left(t\right),$ ${\sigma }_{r\theta }=\mathrm{0,}$

$f\left(t\right)=e^{-{\left(\frac{\pi f_{0}t}{2}\right)}^2}\sin \left(2\pi f_{0}t\right).$ (1)

В (1) ${\sigma }_{rr},{\sigma }_{r\theta }$ – напряжения, $f\left(t\right)$ – функция источника по времени, $f_0$– его несущая частота, $R_1$ ‒ радиус Земли. В этом случае возбуждается поле смещений $\widehat{u}=u_r(r,\theta ,t){\widehat{e}}_r+u_{\theta }(r,\theta ,t){\widehat{e}}_{\theta }$, не зависящее от координаты φ.

В начальный момент времени ставятся нулевые начальные данные:

${\left.u_r=\frac{\partial u_r}{\partial t}=u_{\theta }=\frac{\partial u_{\theta }}{\partial t}\right|}_{t=0}=0.$ (2)

Кроме того, добавляется условие ограниченности решения в центре шара

${\left.\widehat{u}\right|}_{r=0}<\infty$. (3)

На границах $r=R_j$, где скорости продольных $V_p\left(r\right)$, поперечных $V_s\left(r\right)$ волн и плотности $\rho \left(r\right)$ терпят разрыв, вводятся известные условия сопряжения (непрерывности) [1]:

${\left.\left[u_r\right]=\left[u_{\theta }\right]=\left[{\sigma }_{rr}\right]=\left[{\sigma }_{r\theta }\right]\right|}_{r=R_j}=0.$ (4)

В (4) ${\sigma }_{rr},{\sigma }_{r\theta }$ – напряжения, $u_r,u_{\theta }$ – смещения.

ПОСТРОЕНИЕ И УСТОЙЧИВОЕ ВЫЧИСЛЕНИЕ АНАЛИТИЧЕСКОГО РЕШЕНИЯ БЕЗ ПОГРЕШНОСТИ

На первом этапе решение краевой задачи (1)–(4) ищется в виде разложения Фурье-Лежандра по переменным $(\theta ,t)$ [1, 2].

$\begin{array}{c}{u}_r\left(r,\theta ,t\right)=\\ =\frac{1}{2T}\sum _{m=-\infty }^{\infty }\sum _{k=0}^{\infty }{u}_r\left(r,k,{\omega }_m\right)\exp \left(i{\omega }_{m}t\right)P_k\left(\cos \theta \right),\end{array}$

$$$\begin{array}{c}{u}_r\left(r,\theta ,t\right)=\\ =\frac{1}{2T}\sum _{m=-\infty }^{\infty }\sum _{k=0}^{\infty }{u}_{\theta }\left(r,k,{\omega }_m\right)\exp \left(i{\omega }_{m}t\right)P_k^1\left(\cos \theta \right).\end{array}$ (5)

Здесь $P_k\left(x\right)$, $P_k^1\left(x\right) (-1\le x\le 1)$ – многочлены Лежандра [1], ${\omega }_m=m\cdot \pi /T$. Далее для сокращения записи несущественные переменные опускаются.

После применения (5) постановка (1)–(4) сведена к двухпараметрическому семейству (k, ${\omega }_m$) уравнений в каждом сферическом слое $R_{j+1}<r<R_j$.

$\begin{array}{c}\frac{d^2{u}_{\theta }}{d{r}^2}+\frac{2}{r}\frac{d{u}_{\theta }}{dr}+\frac{\lambda +\mu }{\mu }\frac{1}{r}\frac{d{u}_r}{dr}-\\ -\frac{n}{r^2}\frac{\lambda +2\mu }{\mu }{u}_{\theta }+\frac{{\omega }^2}{V_s^2}{u}_{\theta }+\frac{2}{r^2}\frac{\lambda +2\mu }{\mu }{u}_r=\mathrm{0,}\end{array}$

$\begin{array}{c}\frac{d^2{u}_r}{d{r}^2}-\frac{n}{r}\frac{\lambda +\mu }{\lambda +2\mu }\frac{d{u}_{\theta }}{dr}+\frac{2}{r}\frac{d{u}_r}{dr}+\frac{n}{r^2}\frac{\lambda +3\mu }{\lambda +2\mu }{u}_{\theta }-\\ -\frac{u_r}{r^2}\left(2+n\mu /\left(\lambda +2\mu \right)\right)+\frac{{\omega }^2}{V_p^2}{u}_r=0.\end{array}$ (6)

На поверхности Земли ($r=R_1$) и в центре Земли ($r=0$) из (1) и (3) получим краевые условия в спектральной области:

${\left.{\sigma }_{rr}\right|}_{r=R_1}=2F\left(\omega \right)\cdot \nu /R_1^2, {\left.{\sigma }_{r\theta }\right|}_{r=R_1}=\mathrm{0,}$

$\nu =k+\mathrm{0.5,}$ ${\left.u_r\right|}_{r=0}<\infty ,$ ${\left.u_{\theta }\right|}_{r=0}<\infty .$(7)

В (6) и (7) $n=k(k+1)$, λ, µ ‒ упругие параметры Ламе, ‒ плотность.

На втором этапе осуществляется переход к потенциалам продольных u_p и поперечных u_s волн. Известно классическое представление полного поля через потенциалы [9]. В данной работе оно берется в виде:

$u_r=\frac{d}{dr}\left(\frac{u_p}{r^{0.5}}\right)+\frac{k\left(k+1\right)}{r^{1.5}}{u}_s,$

$u_{\theta }=\frac{d}{dr}\left(\frac{u_s}{r^{0.5}}\right)+\frac{1}{r^{1.5}}{u}_s+\frac{1}{r^{1.5}}{u}_p.$ (8)

Представление (8) позволяет свести уравнения (6) сразу к уравнениям Бесселя для потенциалов продольных u_p и поперечных u_s волн:

$\frac{d^2{u}_p}{d{r}^2}+\frac{1}{r}\frac{d{u}_p}{dr}+\frac{{\omega }_n^2}{V_p^2}{u}_p-\frac{{\nu }^2}{r^2}{u}_p=\mathrm{0,}$

$\frac{d^2{u}_s}{d{r}^2}+\frac{1}{r}\frac{d{u}_s}{dr}+\frac{{\omega }_n^2}{V_s^2}{u}_s-\frac{{\nu }^2}{r^2}{u}_s=\mathrm{0,} \nu =k+\mathrm{0.5.}$ (9)

Краевая задача второго порядка (6), (7) с учетом (9) и (4) сводится к двум задачам Коши первого порядка. Для этого по аналогии с [2] вводятся вспомогательные функции $a_{11}\left(r\right), a_{22}\left(r\right), a_{12}\left(r\right), a_{21}\left(r\right)$ следующим образом:

$\frac{d{u}_p}{dr}=\frac{a_{11}}{r}{u}_p+\frac{a_{12}}{r}{u}_s,$ (10)

В итоге для слоистого шара получен рекуррентный процесс для нахождения функций $a_{pq}$ (p, q = 1,2) при $r=R_1$. Далее из (7) и (8) получим искомые смещения в спектральной области. Окончательно суммированием в (5) получаем решение в физической области $u_r(r,\theta ,t)$ и $u_{\theta }(r,\theta ,t)$.

Но поскольку функции Бесселя $J_{\nu }\left(z\right)$ и $Y_{\nu }\left(z\right)$ из (9) при возрастании пространственной частоты ν быстро стремятся к нулю и бесконечности [1], в аналитическом решении возникают особенности. В этой ситуации вычисление на компьютере становится неустойчивым. Происходит выход за границы числового диапазона для любой вычислительной платформы.Поэтому для расчета телесейсмических волн для несущих частот в 1 герц и выше используется асимптотика. Пусть $J_{\nu }\left(x\right)~0$, $Y_{\nu }\left(x\right)~\infty$ и $x\ne 0$. Тогда и их производные $J_{\nu }^{/}\left(x\right)~0$, $Y_{\nu }^{/}\left(x\right)~\infty$. В этом случае в [3, 4] получена новая асимптотика:

$\frac{J_{\nu }^{/}\left(x\right)}{J_{\nu }\left(x\right)}~\frac{\sqrt{{\nu }^2-x^2}}{x},$ $\frac{Y_{\nu }^{/}\left(x\right)}{Y_{\nu }\left(x\right)}~-\frac{\sqrt{{\nu }^2-x^2}}{x}.$ (11)

Использование (11) позволяет устойчиво вычислять аналитическое решение для телесейсмических волновых полей для Земли в случае несущей частоты в 1 герц и выше.Может возникнуть вопрос. Чем плохи классические асимптотики? Почему бы их не использовать? Рассмотрим классические асимптотики [10]. Они верны при $\nu \gg \left|x\right|$.

$J_{\nu }\left(x\right)=\frac{1}{\sqrt{2\pi \nu }}{\left(\frac{ex}{2\nu }\right)}^{\nu }, Y_{\nu }\left(x\right)=-\sqrt{\frac{2}{\pi \nu }}{\left(\frac{ex}{2\nu }\right)}^{-\nu }$ . (12)

Из (12) следует, что

$\frac{{J^{\text{'}}}_v\left(x\right)}{J_{\nu }\left(x\right)}=\frac{\nu }{x}$ и $\frac{{Y^{\text{'}}}_{\nu }\left(x\right)}{Y_{\nu }\left(x\right)}=-\frac{\nu }{x}$. (13)

Если $\nu \gg \left|x\right|$ то можно положить $\sqrt{{\nu }^2-x^2}~\nu$. Тогда из (11) и (13) следует, что новая асимптотика будет совпадать с классической.

На рис. 1 приведены волновые поля радиальной (вертикальной) компоненты u_r сейсмического поля с несущей частотой 1 герц для однородного шара с радиусом Земли R₁ = 6371 километров. Однородный шар имеет параметры: Vp = 5.8, Vs = 3.45, = 2.72. По горизонтали приведены расстояния в градусах. По вертикали – время в секундах (возрастает вниз). При этом на рис. 1 а используется новая асимптотика (11). А на рис. 1 б – классическая асимптотика (13). Символы P и PP означают прямую и однократную продольные волны. В этом случае кинематика для угла φ для прямой волны P будет даваться известным выражением $t_p=2R_1\sin (\phi /2)/Vp$. Соответственно кинематика однократной продольной волны PP – $t_{pp}=4R_1\sin (\phi /4)/Vp$. Для угла в 30 градусов будет $t_p\approx \mathrm{568,6}$ секунд, а $t_{pp}\approx \mathrm{573,5}$ секунд. Приход этих волн обозначен символами P и PP на рис. 1. Далее все кратные волны распространяются по граням вписанных в круг многоугольников. В пределе многоугольники стремятся к окружности. Длина дуги окружности для угла φ в 30 градусов будет $R_1\cdot pi\cdot \phi /180\approx \mathrm{3335,85}$. Поэтому для времен больших, чем $\mathrm{3335,85}/Vp\approx \mathrm{575,1}$, кратных продольных волн P быть не может. Таким образом, волновая фаза, приходящая на рис. 1 б за P- и PP-волнами, заведомо является помехой, так как она регистрируется на времени (возрастает вниз) большем 580 секунд. Она обозначена символом Noise на рис. 1 б. Из литературы известны и другие асимптотики (например, [11]). Но асимптотики из [11] приводят к формулам (13), то есть результаты ее применения к волновым задачам будут совпадать с классической асимптотикой.

 

Рис. 1. Компонента U_r сейсмического поля с несущей частотой 1 герц для однородного шара земных размеров. В алгоритме расчета используется новая (а) и классическая (б) асимптотики. По вертикали – время в секундах (возрастает вниз), по горизонтали – расстояние в градусах. Буквы P и PP обозначают прямую и однократную продольные волны. Noise – помеха при использовании классической асимптотики.

 

Таким образом, в случае высоких частот использовать классическую асимптотику нельзя. Ее использование приводит к погрешности в решении даже для однородного шара. Причина этого в том, что классическая асимптотика верна при $\nu \gg \left|x\right|$. А при расчетах аналитического решения для шара земных размеров, например, функция Бесселя $J_{\nu }\left(x\right)$становится близка к нулю, когда ν сравнима с $\left|x\right|$. Это и приводит к невозможности использования классической асимптотики. А новая асимптотика как раз применима, когда $J_{\nu }\left(x\right)~0$.

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

РЕЗУЛЬТАТЫ АНАЛИТИЧЕСКОГО РЕШЕНИЯ

На рис. 2 приведена радиальная (вертикальная) компонента для дискретной (слоистой) модели Земли из [7, 8] в диапазоне эпицентральных расстояний 130–160 градусов. Импульс f(t) в источнике взят в виде функции Гаусса в (1) с несущей частотой в 1 герц. По горизонтали приведено расстояние в градусах. По вертикали – время (возрастает вверх). Из рис. 2 видно, что впереди PKP-волн возникают высокочастотные колебания небольшой амплитуды (предвестники). Спектральный анализ показал, что их частота примерно равна 1.25 герца.

 

Рис. 2. Компонента U_r сейсмического поля. По вертикали приведено время в секундах (возрастает вверх), по горизонтали расстояние в градусах. Возникновение прекурсоров для модели Земли из [7, 8].

 

Аналогичные явления происходят и для известной модели AK135 [12]. В работе [4] показано, что для модели AK135 впереди PKP-волн также возникают высокочастотные колебания небольшой амплитуды. Отметим, что для моделей Земли с непрерывным изменением скорости распространения сейсмических волн во внешнем ядре интерференционной картины не возникает [4]. Таким образом, явление возникновения небольших высокочастотных колебаний (предвестников) впереди PKP-волн полностью объясняется сферической геометрией слоистой Земли.

На рис. 3 представлен монтаж сейсмограмм, полученный по данным мировой сети сейсмических станций [4]. Сравнение рисунков 2 и 3 показывает, что теоретические сейсмограммы во многом похожи на экспериментальные сейсмограммы. Это подтверждает правильность идей, положенных в основу построения модели Земли в [7, 8].

 

Рис. 3. Фрагмент монтажа сейсмограмм, полученного на сейсмических станциях мировой сети.

 

ВЫВОДЫ

В работе получено аналитическое решение для волновых полей сейсмических волн в сферически-симметричной Земле. В аналитическом решении за счет быстрого убывания/возрастания бесселевых функций возникают особенности. Это приводит к неустойчивости решения на любой вычислительной платформе. Для устранения не- устойчивости нужно использовать асимптотику. Исследовано применение классической и новой асимптотики для случая первого вступления продольных P-волн. Показано, что в случае шара земных размеров классическая асимптотика для высоких частот дает погрешность в решении, то есть ее использовать нельзя.

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

По созданной программе проведен расчет с высокой детальностью для дискретной (слоистой) модели Земли из [7, 8] с несущей частотой в 1 герц. Она характеризуется тем, что в ней внешнее ядро обладает модулем сдвига, отличным от нуля. В результате аналитического расчета обнаружено, что впереди PKP-волн возникают высокочастотные колебания небольшой амплитуды, так называемые “предвестники”.

Аналитический расчет показал, что теоретические сейсмограммы для модели Земли из [7, 8] во многом похожи на экспериментальные сейсмограммы, полученные мировой сетью сейсмических станций. Это сравнение говорит о правильности идей, положенных в основу построения модели Земли в работе [7, 8].

Источники финансирования

Работа выполнена в соответствии с Государственным заданием ИФЗ РАН № 0144-2019-0011 и ИВМиМГ СО РАН № 0251-2021-0004.

Авторы заявляют, что у них нет конфликта интересов.

Об авторах

А. Г. Фатьянов

Институт вычислительной математики и математической геофизики, Сибирское отделение Российской Академии наук

Автор, ответственный за переписку.
Email: fat@nmsf.sscc.ru
Россия, Новосибирск

В. Ю. Бурмин

Институт физики Земли им. О. Ю. Шмидта Российской Академии наук

Email: burmin@ifz.ru
Россия, Москва

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

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