Физика Земли, 2021, № 5, стр. 17-29

Механическая добротность в верхней части внутреннего ядра Земли по данным коды волн РКiКP

Д. Н. Краснощеков 1, В. М. Овчинников 1*, О. А. Усольцева 1**

1 Институт динамики геосфер имени академика М.А. Садовского РАН
г. Москва, Россия

* E-mail: ovtch1@yandex.ru
** E-mail: kriukova@mail.ru

Поступила в редакцию 24.02.2021
После доработки 22.03.2021
Принята к публикации 01.04.2021

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

Аннотация

На основе анализа методом компьютерной графики альфа-шейп огибающей коды отраженных от поверхности внутреннего ядра Земли волн PKiKP, возбужденных ядерными взрывами, получены оценки величины механической добротности для девяти областей отражения под Арктикой, Центральной и Юго-Восточной Азией со средним значением Q = 447 ± 43. Механизм формирования коды может быть обусловлен неоднородностями в верхней части внутреннего ядра с характерным линейным размером 1–4 км и вариациями скорости продольных волн порядка 1–3%.

Ключевые слова: Q-фактор, внутреннее ядро, альфа-шейп, волна PKiKP.

ВВЕДЕНИЕ

Затухание сейсмических волн во внутреннем ядре обусловлено неупругими характеристиками среды, к которым относятся мелкомасштабные неоднородности типа дислокаций в кристаллах, внутреннее трение и другие [Cormier, Li, 2002; Li, Cormier, 2002; Cormier 2009]. Эти факторы принято называть внутренним (собственным) поглощением. С другой стороны, затухание связано также с эффектами рассеивания, которые являются упругими процессами перераспределения энергии за счет отражения, преломления и обмена на границах раздела в среде. Эти рассеивающие элементы, расположенные на трассе распространения, при их малом размере, приводят к уменьшению амплитуды сейсмической волны [Vidale et al., 2000; Koper et al., 2004; Poupinet, Kennet, 2004; Krasnoshchekov et al., 2005; Leyton, Koper, 2007; Овчинников и др., 2007]. В этих исследованиях анализируется форма коды, механизмы ее формирования и связь со структурными особенностями твердого ядра и вышележащих оболочек Земли. Установлено, что кода, сформированная во внутреннем ядре, имеет форму веретена [Leyton, Koper, 2007], а огибающая – форму арки [Овчинников и др., 2007], т.е. кода имеет фазы роста и спада. Если же огибающая коды монотонно спадает (без фазы роста), то такая форма может быть объяснена механизмом ревербераций в переходной зоне внутреннее–внешнее ядро и, возможно, неоднородностями вышележащих оболочек [Poupinet, Kennett, 2004].

Количественной оценкой затухания сейсмической волны служит относительная потеря энергии за один цикл колебаний $q = \frac{{\Delta E}}{{2\pi E}}~$ или обратная ей величина $q = \frac{1}{Q} = \frac{1}{{Qa}} + \frac{1}{{Qc}}$. (Здесь Qa – неупругое поглощение, Qc – фактор рассеяния). Неупругое поглощение очень чувствительно к изменению давления и температуры среды. Поэтому в условиях внутреннего ядра, где давление меняется слабо, пространственные вариации Q могут служить источником особенностей теплового режима ядра.

За 20 лет после первой публикации [Vidale, Earle, 2000] c исследованием коды докритически отраженных от границы внешнее–внутреннее ядро волн PKiKP и возможностью их использования для определения Q во внутреннем ядре было опубликовано не более десятка работ (например, [Poupinet, Kennett, 2004; Koper et al., 2004; Leyton, Koper, 2007; Peng et al., 2008; Wu, Irving, 2017]), посвященных этому аспекту.

Величины добротности внутреннего ядра по коде, полученные из аппроксимации ее формы модельной функцией, варьируют в широком диапазоне 350–600 [Vidale, Earle, 2000; Leyton, Koper, 2007; Peng et al., 2008], а разброс индивидуальных измерений превышает два порядка [Leyton, Koper, 2007]. Значительная дисперсия при относительно небольшом общем количестве опубликованных оценок не позволяет надежно интерпретировать эти данные в терминах региональных особенностей структуры и сейсмического затухания во внешней части твердого ядра, которые установлены из других данных [Wu, Irving, 2017; Krasnoshchekov et al., 2005; Attanayake et al., 2014; Oreshin, Vinnik, 2004]. Источники неопределенности опубликованных значений Q по коде PKiKP могут быть связаны с нечеткими критериями оценки параметров моделируемой коды (длительность, время начала, и т.д.), которые обусловлены применением нестрого определенных процедур обработки сейсмической записи (например, сглаживания).

В данной статье представлены девять новых оценок Q во внутреннем ядре Земли для областей отражения волн PKiKP под Арктикой, Центральной и Юго-Восточной Азией, полученные с использованием для анализа коды нового робастного метода компьютерной графики альфа-шейп [Krasnoshchekov, Polishchuk, 2008; Nikkila et al., 2014] и численного моделирования коды.

ИСПОЛЬЗОВАННЫЕ ДАННЫЕ И МЕТОДЫ ОБРАБОТКИ

Для оценки Q были использованы сейсмограммы групп взрывов, проведенных на Семипалатинском полигоне (STS) в Казахстане в период с 1977 по 1988 гг. и зарегистрированных на одиночных сейсмических станциях BRVK (Казахстан), KEV (Финляндия), COL и LON (США), а также сейсмограммы двух взрывов на полигоне Лобнор (LNTS) в Китае, зарегистрированных на сейсмических станциях группирования ASAR (Австралия), FINES (Финляндия), KURK (Казахстан), PDAR (США), YKA (Канада). Расположение станций и полигонов показано на карте (рис. 1).

Рис. 1.

Карта взаимного расположения полигонов и сейсмических станций с дугами большого круга (пунктир): STS – Семипалатинский полигон, LNTS – полигон Лобнор, BRVK, KEV, COL, LON – трехкомпонентные сейсмические станции, ASAR, FINES, KURK, PDAR, YKA – станции сейсмического группирования; черные кружки – проекции точек отражения волн PKiKP от внутреннего ядра Земли на дневную поверхность.

Преимущество использования взрывного источника состоит в значительно более простом, чем у землетрясений, механизме очага, а также в возможности снизить неопределенность результирующих оценок за счет привлечения несейсмологической информации о координатах и времени в очаге. Чувствительность использованных широкополосных и короткопериодных каналов регистрации с плоской АЧХ в полосе частот 1–5 Гц и энергетический класс взрывов (5.9 < mb <6.1) достаточны для обнаружения коды волны PKiKP c амплитудой, составляющей единицы нанометра [Vidale et al., 2000; Адушкин, Овчинников, 2004; Krasnoshchekov et al., 2005; Leyton, Koper, 2007; Овчинников и др., 2007]. Области зондирования внутреннего ядра находятся под Арктикой, Центральной Азией и Юго-Восточной Азией (рис.1).

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

Таблица 1.  

Список станций и взрывов с датами их проведения

Код сейсмической станции или группы станций Дата взрыва,
дд.мм.гг
Эпицентральное расстояние, град Количество пунктов наблюдения или событий в группе
BRVK Перечень взрывов из [Krasnoshchekov et al., 2005] 6.2 31
KURK 17.08.1995 11.5 19
KEV 18.10.1981, 25.04.1982, 12.06.1983, 14.07.1984, 27.10.1984, 6.12.1984, 28.12.1984, 15.06.1985, 30.06.1985, 20.07.1985 31.3 10
FINES 08.06.1996 41.8 16
COL 05.12.1982, 25.04.1982, 12.06.1983,
14.07.1984, 24.05.1985, 25.04.1984,
20.07.1985,16.12.1984, 27.101984,
06.10,1983, 26.10.1983, 19.02.1984,
10.02.1985
59.9 15
YKA 08.06.1996 73.8 21
ASAR 08.06.1996 77.2 19
LON Перечень взрывов из [Krasnoshchekov et al., 2005] 82.1 24
PDAR 08.06.1996 94.2 15

Обработка взрывов, зарегистрированных станциями сейсмического группирования, включала несколько последовательно выполняемых шагов, направленных на обнаружение волны PKiKP и ee коды. Они демонстрируются на рис. 2 обработкой данных станции группирования FINES. Вначале к каждому из 16 каналов системы группирования применялась полосовая фильтрация нуль-фазовым фильтром Баттеруорта с граничными частотами 1 и 5 Гц с крутизной склонов 3-го порядка. Отфильтрованные сейсмограммы Si (i = 1 : 16) показаны на рис. 2а. Пространственные особенности сейсмического волнового поля определялись путем вычисления в скользящем временном окне двумерного спектра мощности P(k) (F–K-анализ):

$P({\mathbf{k}}) = \frac{1}{m}\mathop \sum \limits_{{{f}_{1}}}^{{{f}_{m}}} \left( {\mathop \sum \limits_{j - 1}^{j = n - 1} \mathop \sum \limits_{i = 1}^{i = n} {{F}_{{ij}}}(f){{e}^{{i{\mathbf{k}}\left( {{{{\mathbf{r}}}_{i}} - {{{\mathbf{r}}}_{j}}} \right)}}}} \right),$
где: ${{F}_{{ij}}}(f) = \frac{2}{T}\sum\nolimits_{{{k}_{\alpha }}} {S_{i}^{*}{{S}_{j}}} $ – взаимный спектр в частотном интервале [f1; fm]; ${\mathbf{k}} = \left( {{{k}_{x}}{{k}_{y}}} \right)$ – волновой вектор в горизонтальной плоскости; * – знак комплексного сопряжения, и определялось значение k0, соответствующее максимуму спектра мощности (рис. 2б) в ожидаемом интервале времени вступления волны PKiKP. Далее по отфильтрованным каналам группы вычислялась суммарная сейсмограмма (рис. 2в):
$b\left( {t:{{{\mathbf{k}}}_{0}}} \right) = \frac{1}{N}\mathop \sum \limits_{i = 1}^N {{S}_{i}}\left( {t + {{\tau }_{i}}} \right),\,\,\,\,{{\tau }_{i}} = {{ - {{{\mathbf{r}}}_{i}}{{{\mathbf{k}}}_{0}}} \mathord{\left/ {\vphantom {{ - {{{\mathbf{r}}}_{i}}{{{\mathbf{k}}}_{0}}} \omega }} \right. \kern-0em} \omega },$
где τi – временной сдвиг сейсмограммы на i-ой станции относительно опорной станции.

Рис. 2.

(а) – Фрагменты сейсмограмм, отфильтрованные в полосе 1.3–5 Гц, сейсмической группы FINESS c волной PKiKP; (б) – F–K-спектр мощности; (в) – суммарная сейсмограмма (beam) с временными задержками; (г) – статистическая характеристика сейсмического процесса Fstat (максимумы соответствуют вступлениям волн P, PcP, PKiKP); (д) – огибающая Гильберта.

Чтобы оценить значимость уровня сейсмического сигнала на фоне шума в рассматриваемом временном интервале рассчитывалась статистическая характеристика (рис. 2г) сейсмического процесса [Blandford, 1974]:

${{F}_{{stat}}}\left( {{{{\mathbf{k}}}_{0}}:T,t} \right) = \,\,~\frac{{N - 1}}{N}\frac{{\sum\limits_t^{t + T} {{{{\left\{ {\left( {b\left( {t:{{{\mathbf{k}}}_{0}}} \right)} \right)} \right\}}}^{2}}} }}{{\sum\limits_t^{t + T} {\sum\limits_{i = 1}^N {{{{\left| {{{S}_{i}}\left( t \right) - ~b\left( {t:{{{\mathbf{k}}}_{0}}} \right)} \right|}}^{2}}} } }}.$

Физический смысл Fstat – насколько полная энергия колебаний (числитель) для заданного волнового числа k0 в интервале наблюдения [t; t + T] отличается от остаточного шума (знаменатель). Заключительным шагом этого предварительного этапа обработки является вычисление огибающей Гильберта hb(t) = H(b(t : k0)) (рис. 2д), H – символ преобразования Гильберта.

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

Временная синхронизация была проведена по первому вступлению волны P путем определения временных сдвигов корреляционным анализом. Далее по сейсмическому шуму, предшествующему вступлению волны P для каждой сейсмограммы, был определен нулевой уровень по участку фона, предшествовавшему вступлению Р-волны, а сами записи преобразованы из дискретных в непрерывные функции fk(t) с помощью сплайнов. По пересечению с нулевой линией для каждой k‑ой записи было найдено время вступления Р-волны tk и длина ${{L}_{k}}$ ее первого периода. Далее были найдены максимумы коэффициентов корреляции:

${{R}_{{kj}}}\left( {\delta {{t}_{{kj}}}} \right) = \frac{1}{B}\mathop \smallint \limits_{{{t}_{k}}}^{{{t}_{k}} + {{L}_{k}}} {{f}_{k}}\left( \tau \right){{f}_{j}}\left( {\tau + {{\delta }_{{kj}}}} \right)d\tau {\text{,}}$
где
$B = \sqrt {\int\limits_{{{t}_{k}}}^{{{t}_{k}} + {{L}_{k}}} {{{f}_{k}}{{{(\tau )}}^{2}}} d\tau \int\limits_{{{t}_{k}}}^{{{t}_{k}} + {{L}_{k}}} {{{f}_{i}}{{{(\tau + {{\delta }_{{kj}}})}}^{2}}} d\tau } ,$
и соответствующие им временные сдвиги ${{\delta }_{{kj}}}$ для всех пар записей (kj; k, j = 1, N), где N – общее их число. По найденным матрицам ${{R}_{{kj}}}$ и ${{\delta }_{{kj}}}$ были вычислены взвешенные относительные сдвиги:
${{{{\theta }}}_{k}} = {{\left( {\frac{{\sum\limits_j {{{\delta }_{{kj}}}R_{{kj}}^{2}} }}{{\sum\limits_j {R_{{kj}}^{2}} }} - \frac{{\sum\limits_j {{{\delta }_{{jk}}}R_{{jk}}^{2}} }}{{\sum\limits_j {R_{{jk}}^{2}} }}} \right)} \mathord{\left/ {\vphantom {{\left( {\frac{{\sum\limits_j {{{\delta }_{{kj}}}R_{{kj}}^{2}} }}{{\sum\limits_j {R_{{kj}}^{2}} }} - \frac{{\sum\limits_j {{{\delta }_{{jk}}}R_{{jk}}^{2}} }}{{\sum\limits_j {R_{{jk}}^{2}} }}} \right)} 2}} \right. \kern-0em} 2},\,\,\,\,k = {\text{ }}1,N,$
которые и были использованы для окончательной синхронизации сейсмограмм различных взрывов по Р-волне.

Для выравнивания амплитуд использовались среднеквадратичные амплитуды:

${{A}_{k}} = \sqrt {\frac{{\int\limits_{{{t}_{k}}}^{{{t}_{k}} + {{L}_{k}}} {{{f}_{k}}{{{(\tau )}}^{2}}d\tau } }}{{{{L}_{k}}}}} ,\,\,\,\,k = 1,N,$
и соответственно скорректированные записи:
${{f}_{k}}\left( t \right) = {{f}_{k}}\frac{{{{A}_{{{\text{max}}}}}}}{{{{A}_{k}}}},$
где ${{A}_{{{\text{max}}}}} = {\text{max}}\left( {{{A}_{k}}} \right)$.

Заключительный этап обработки как данных станций группирования, так и групп взрывов состоял в построении формы огибающей с использованием α-шейпа k-го порядка. Процедура состояла из двух шагов. На первом по реализации сейсмического шума, предшествующего первому вступлению (длительность 1–2 мин), определяются радиус α и порядок k. На втором шаге восстанавливалась гладкая форма огибающей Гильберта. Последовательность оценки α-шейпа k-го порядка временного ряда (или облака точек) [Krasnoshchekov, Polishchuk 2008; 2014] схематически представлена на рис. 3. В каждый момент времени инициализируются два диска радиуса α верхний диск и нижний диск. Верхний диск располагается выше всех точек набора данных и двигается вниз. Диск останавливается в тот момент, когда он охватывает ровно k точек (рис. 3, левая панель). Аналогично нижний диск располагается под данными и движется вверх до момента включения в него k точек. Выходным элементом искомой формы в данный момент времени являеются координаты середины отрезка, соединяющего центры остановившихся верхнего и нижнего дисков.

Рис. 3.

Применение α-шейп k-го порядка к временны́м последовательностям (русифицированный рисунок из работы [Nikkila et al., 2014]). Левая панель: верхний диск продвигается вниз по данным, а нижний – вверх; диск останавливается, когда в нем содержится точно k = 4 точек множества. Восстановленный элемент выходной формы в этот момент времени соответствует середине отрезка, соединяющего центры остановившихся дисков. Правая панель: несколько точек (как внутри, так и снаружи дисков) изменили свои местоположения; несмотря на это, восстановленный элемент формы не изменился.

Чтобы пояснить работу алгоритма, рассмотрим применение алгоритма в двух предельных случаях: α = 0 и α = ∞. В первом случае в качестве выходной формы будет получена входная временная последовательность. Во втором случае выходной формой будет являться прямая линия на уровне среднего между k-ым максимумом и k-ым минимумом временнóй последовательности.

Выбор α и k для любой из сейсмограмм осуществляется по фрагменту сейсмограммы (материал обучения – сейсмический шум), которая предшествует первому вступлению. Вначале радиус α устанавливается равным одному стандартному отклонению сейсмического шума на интервале 1–2 мин. Для определения k случайным образом выбирается 1000 моментов записи, и в каждый момент времени располагается два α‑диска. Один с центром в 2α над средним уровнем шума, а другой с центром 2α под уровнем шума. Значение k устанавливается равным среднему количеству точек в 2000 дисках.

С полученными значениями α и k алгоритм α-шейп k-го порядка применяется к записи сейсмического шума и с помощью среднеквадратичного отклонения (RMSD) оценивается, насколько выходная форма отличается от прямой, соответствующей среднему уровню шума. Интуитивно ясно, что α-диски удаляют шум, что должно обеспечить прямую линию в качестве выходной формы. Если вариабельность выходной формы слишком высока, то есть когда RSMD > 0.05α, шкала времени масштабируется (сжимается в целое число раз), что приводит к увеличению числа точек, попадающих в пределы дисков радиуса 2α. Для новой шкалы снова вычисляется величина k с помощью вышеописанной процедуры (усреднение по 2000 дисков в случайные моменты времени). При увеличении масштабного коэффициента возрастает k, что соответствует меньшим вариациям выходной формы или меньшему RMSD. Процедура останавливается в тот момент, когда RMSD становится менее заранее установленного уровня 0.05α. Меньшие значения RMSD обычно приводят к избыточному сглаживанию выходной формы. Полученные значения α и k используются для определения формы огибающей для всей сейсмграммы. (Демонстационный скрипт для Матлаба доступен по ссылке http://www.cs.helsinki.fi/group/compgeom/seism.zip).

РЕЗУЛЬТАТЫ ОБРАБОТКИ

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

Начало коды PKiKP ожидается вскоре после вступления родительской фазы, время пробега которой для вышеуказанных эпицентральных расстояний составляет от 992 до 1084 с. Соответствующие фрагменты огибающей действительно проявляют признаки роста амплитуд на фоне быстрых изменений и присутствия экстремальных значений. Например, кода легко визуализируется на данных группы FINES (рис. 4). Мы пока не обсуждаем показанные на рис. 4 восстановленные методом альфа-шейп “гладкие” огибающие, это будет сделано ниже. Она проявляется в форме арки после резкого вступления волны PKiKP на времени 1011 с. Обнаруженная кода на суммарной трассе FINES позволяет сделать грубую оценку длительности коды порядка 50 с, а также разделение ее частей на фазу роста и фазу спада.

Рис. 4.

Огибающие Гильберта (серый цвет) на станциях KEV и FINES. Черные кривые результат сглаживания методом альфа-шейп.

Напротив, на записях станции KEV обнаружение коды PKiKP затруднительно. Традиционным средством обнаружения и более рельефного отображения коды и ее частей является сглаживание путем усреднения в скользящем временном окне. Эта процедура действительно позволяет устранить быстрые изменения огибающей, а также аномальные значения амплитуд. Вместе с тем, использование усреднения сопряжено с внесением дополнительных неопределенностей. Одним из негативных факторов является выбор длины временного окна усреднения. Его длина обычно определяется, исходя из доминирующей частоты выделяемого сигнала и его медленности. Последняя может оцениваться на основе частотно-волнового анализа соответствующего фрагмента огибающей. Усреднение с таким окном может приводить к неустойчивым оценкам в случае применения ко всей коде, так как выбор окна осуществлялся на основе свойств родительской фазы и, возможно, применим только для ограниченного диапазона времен.

На рис. 5 представлены результаты сглаживания огибающей суммарной трассы группы FINES во временных окнах длиной 1, 2.5, 5 и 50 с.

Рис. 5.

Огибающие Гильберта (серый цвет) и сглаженные методом скользящего среднего с различной длиной окна усреднения (черный цвет): сверху вниз – 1, 2.5, 5, 50 с.

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

Рассмотрим теперь результаты применения методики α-шейп k-го порядка для выявления гладкой формы огибающей Гильберта. На рис. 4 и рис. 6 показаны примеры восстановления огибающих методом альфа-шейп.

Рис. 6.

Огибающие Гильберта и восстановленные методом α-шейп гладкие формы коды. Название станции (группы) указано в правом верхнем углу.

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

Кода PKiKP на полученных трассах отражается в виде изменения градиента спада после теоретического времени прихода PKiKP. Как видно, кода, рассеянная во внутреннем ядре Земли, автономна, так как может наблюдаться как в присутствии, так и в отсутствии родительской волновой формы PKiKP. Длительность фазы спада составляет 30–50 с (рис. 4, рис. 6). При этом меньшая длительность соответствует станциям на расстояниях около 6° (BRVK) и 11° (KURK). По-видимому, это связано с увеличением длины пути распространения волны по неоднородной среде при увеличении расстояния, и, следовательно, данные наблюдений на малых эпицентральных расстояниях несут информацию о глубине области неоднородностей в верхней части внутреннего ядра. Согласно феноменологической теории рассеяния [Aki, Сhuet, 1975], форма коды PKiKP определяется энергией, переданной во внутреннее ядро, и не предполагает вклада со стороны собственно отражения от его границы (фазы PKiKP). На практике влияние родительской фазы будет сказываться на любой форме, полученной посредством скользящего среднего, кроме случаев исключительно длинных временных окон, надежно компенсирующих вклад импульса, соответствующего волновой форме PKiKP. В отличие от усреднения, представленные формы огибающей свободны от этого недостатка, так как используемый алгоритм обрабатывает импульс PKiKP, как любое другое аномальное значение огибающей. Это преимущество важно, так как форма коды PKiKP указывает на место ее происхождения. В частности, результаты моделирования в рамках приближения однократного рассеяния (Борновское рассеяние) [Leyton, Koper, 2007] только коды PKiKP с фазами роста и спада формируются на объемных неоднородностях внутреннего ядра без добавления эффектов, связанных с рассеянием в нижней мантии. Как видно, обнаруженные коды PKiKP проявляют слабую или выраженную фазы роста, за которой следует монотонный спад. Монотонно спадающие коды, начинающиеся сразу после импульсного вступления PKiKP, в анализируемом наборе данных не обнаружены.

Несколько пар наблюдений имеют близкие точки отражения от поверхности внутреннего ядра, поэтому рассмотрим их более подробно. У первой пары (рис. 7а) STS–BRVK и LNTS–KURK наблюдается наименьшее расхождение в эпицентральном расстоянии (чуть более 5°), при этом точки отражения PKiKP на поверхности внутреннего ядра расположены на расстоянии 173 км. У этой пары трасс наблюдаются идентичные коды PKiKP как в фазе роста, так и в фазе спада. У двух оставшихся пар код PKiKP наблюдается существенное различие в фазе роста, при этом градиенты спада после максимума коды PKiKP практически идентичны. У пары (рис. 7б) STS–COL и LNTS–YKA разность расстояний составляет 14°, точки отражения разнесены на 47 км, а у пары (рис. 7в) STS–LON и LNTS–PDAR расстояния различаются на 12°, а точки отражения разделяют 17 км. Эти наблюдения подтверждают модельные представления о форме коды PKiKP, в соответствии с которыми форма фазы роста определяется преимущественно геометрическим расхождением или эпицентральным расстоянием, в то время как фаза спада формируется во внутреннем ядре Земли и определяется его затуханием. Как видно из рис. 7, коды, сформированные в одном и том же регионе внутреннего ядра, характеризуются близкими градиентами спада коды PKiKP.

Рис. 7.

Коды, восстановленные с помощью α-шейп k‑го порядка для пар сейсмических трасс с близкими точками отражения PKiKP на поверхности внутреннего ядра Земли. Нулевой уровень по оси ординат соответствует среднему предшествующего сейсмического шума. Расчетные времена вступления PKiKP даны по модели PREM.

Необходимо отметить, что приведенные на рисунке коды PKiKP были получены без учета возможного влияния предыдущих, связанных со структурными особенностями в коре и мантии фрагментов коды. Это обстоятельство указывает на отсутствие необходимости в учете таких эффектов. Кроме того, представленная методика восстановления формы коды дает эффективный инструмент разделения фаз роста и спада коды PKiKP, что важно при аппроксимации обнаруженной формы коды модельной функцией. В частности, это позволяет избавить интерпретатора от волюнтаризма при выборе начала фазы спада, которое, например, было установлено “на глаз” равным десяти секундам после вступления PKiKP в работе [Leyton, Koper 2007].

Таким образом, использование α-шейп k-го порядка для восстановления формы коды PKiKP позволяет устранить часть неопределенностей при обработке данных, которые могут приводить к искажению искомой формы и, следовательно, повышать ошибку соответствующих оценок сейсмического затухания. Так, на стадии обнаружения коды PKiKP искажения ее формы могут происходить вследствие неоднозначного выбора длины окна усреднения и влияния родительской фазы. При аппроксимации обнаруженной формы модельной функцией увеличение ошибки может быть вызвано неправильным выбором времени начала фазы спада или длительности. Полагаем, что ключевым преимуществом использования α-шейп k-го порядка является отсутствие необходимости в выборе параметров обработки “на глаз”.

ОЦЕНКА Q

Оценка Q проведена в соответствии с формализмом, развитым в работе [Aki, Сhuet, 1975]. Предполагая, что кода формируется за счет однократно отраженных от случайно распределенных в среде неоднородностей, было получено уравнение, связывающее амплитуду коды A на частоте f в момент времени t с величиной затухания Q, обусловленной рассеянием:

$A\left( {f,t} \right) = S\left( f \right){{t}^{{ - n}}}{{e}^{{ - \pi ft/Q}}},$
где n – геометрическое расхождение. Для коды волн PKiKP n ≈ 1. Логарифмируя приведенное выше уравнение, получаем:
${\text{ln}}(A\left( {f,t} \right) + \ln t = \ln \left( {S\left( f \right)} \right) - {{\pi ft} \mathord{\left/ {\vphantom {{\pi ft} Q}} \right. \kern-0em} Q},$
откуда следует, что величина Q может быть получена из коэффициента наклона $b = \frac{{\pi f}}{Q}~$ линейной регрессии ${\text{ln}}(A\left( {f,t} \right) + \ln t$ на t.

Полученные значения Q по спаду восстановленных гладких форм огибающих (рис. 4, рис. 6, рис. 7) на основе регрессионного анализа со средним значением Q = 447 ± 43 приведены в табл. 2.

Таблица 2.  

Оценки Q для зондируемых областей внутреннего ядра

Код сейсмической станции или
группы станций
Широта точки
отражения, град
Долгота точки отражения, град Эпицентральное
расстояние, град
Q
BRVK 51.6 74.71 6.2 514
KURK 50.5 76.0 11.5 516
KEV 55.8 65 31.3 449
FINES 55.4 65.6 41.8 407
COL 57.4 145.5 59.9 448
YKA 75.9 116.4 73.8 455
ASAR 9.7 113.8 77.2 413
LON 74.3 118 82.1 423
PDAR 79.9 160.9 94.2 398

ВОЗМОЖНЫЕ МЕХАНИЗМЫ ФОРМИРОВАНИЯ И ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОДЫ

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

Чтобы оценить, какими свойствами среды можно объяснить полученные значения Q и длительность коды, было рассчитано распределение энергии в коде (кодограмма) на частоте 2 Гц в среде со случайными неоднородностями фононным методом [Sanborn, 2015] и синтетические сейсмограммы осесимметричным спектрально-элементным методом AxisemV1.1 [Nissen-Meyer et al., 2014].

На рис. 8 показано распределение энергии в коде волны PKiKP, рассчитанной фононным методом, для модели верхней части внутреннего ядра со случайными неоднородностями с характерным размером 1 км и вариацией скорости продольных волн 1.5%, из которого следует, что длительность коды составляет 40–60 с и имеет тенденцию к уменьшению на малых эпицентральных расстояниях. Если не принимать во внимание затухание за счет неупругих свойств среды, то длительность коды увеличивается до 150 с.

Рис. 8.

Распределение энергии в коде волны PKiKP.

Для полученных в данной работе величины Q ~ 450 и длительности фазы спада огибающей 30–50 с мы провели также расчеты синтетических сейсмограмм на эпицентральных расстояниях 5°–90° для глубинного источника взрывного типа (Н = 500 км) в стандартной модели Земли ak135 и в модифицированных с неоднородностями в кровле внутреннего ядра. Для моделирования использовался осесимметричный спектрально-элементный метод построения синтетически сейсмограмм AxisemV1.1 [Nissen-Meyer et al., 2014] с максимальной расчетной частотой 2 Гц. Вычисления проводились с использованием сверх- высокопроизводительных ресурсов Межведомственного суперкомпьютерного центра РАН на 1088 процессорных ядрах. Длительность счета для каждой тестовой модели около 18 ч с использованием процессоров с частотой 3.0 ГГц. Глубинный источник взрывного типа был выбран, чтобы исключить дополнительные колебания, связанные с влиянием свободной поверхности, земной коры и механизмом очага. Полученные сейсмограммы были отфильтрованы фильтром высоких частот и примерно соответствуют условиям расчета фононным методом для частоты 2 Гц.

Рассмотрены 4 модели с измененной по отношению к модели аk135 скоростью продольных волн путем введения неоднородностей в кровле внутреннего ядра Земли. В модели 1 (M1) случайные неоднородности, автокорреляционная функция которых не зависит от расстояния, расположены в верхних 100 км внутреннего ядра, а максимальная вариация скорости $\frac{{\delta V}}{V} < 4\% $. В модели 2 (M2) регулярные неоднородности размером 6 × 7 км расположены в шахматном порядке в верхних 100 км с максимальной вариацией скорости $\frac{{\delta V}}{V} < 1.5\% $. В модели 3 (M3) регулярные неоднородности размером 3 × 3.5 км расположены в шахматном порядке – в верхних 100 км с максимальной вариацией скорости $\frac{{\delta V}}{V} < 1.5\% $. Модель 4 (M4) аналогична М3, только толщина слоя с неоднородностями составляет 50 км. Наиболее характерные особенности коды на сейсмограммах для моделей M1 и М3 по сравнению со стандартной моделью ak135 показаны на рис. 9.

Рис. 9.

Синтетические сейсмограммы в зоне волны PKiKP (50-секундное временное окно) для стандартной модели ak135 и моделей с неоднородностями (М1 – модель 1, М3 – модель 3) для различных эпицентральных расстояний. Сейсмограммы выровнены относительно вступления волны PKiKP и масштабированы на максимальную амплитуду волны PKiKP.

В среде с крупными неоднородностями (модель М2) частота коды волны PKiKP ниже, чем в среде с более мелкими неоднородностями (модели М3, М4). Длительность коды PKiKP преимущественно составляет 20 с на всех эпицентральных расстояниях при толщине зоны с неоднородностями 100 км. При уменьшении толщины длительность коды PKiKP уменьшается. На эпицентральных расстояниях 70°–85° (рис. 9, М1) кода имеет типичную форму веретена с фазами роста и спада длительностью 20–30 с после вступления волны PKiKP. В модели М3 (рис. 9, М3) кода в течение примерно 15 с имеет слабо меняющуюся амплитуду колебаний. На эпицентральных расстояниях 6°–40° кода более выражена в моделях с регулярными неоднородностями, чем со случайными.

Наши результаты не противоречат исследованиям коды волн PKiKP от взрывов, зарегистрированных на сейсмической группе с большой апертурой (LASA) [Vidale, Earle 2000; Peng et al., 2008]. Наилучшее согласие теоретических сейсмограмм (кодограмм) [Vidale, Earle, 2000], рассчитанных фононным методом в Борновском приближении, было достигнуто для модели среды во внутреннем ядре с 1.2% вариацией плотности и параметров Ламе и характерным размером неоднородности 2 км с экспоненциальной автокорреляционной функцией спектра мощности неоднородностей и диссипативной составляющей Qа = 360. В другой работе [Peng et al., 2008] показано, что Q ~ 500 может соответствовать неоднородной среде с вариациями скорости порядка 1–3% и характерным размером рассеивающей неоднородности 1–4 км. Таким образом, затухание коды волн PKiKP происходит как за счет рассеяния на неоднородностях, так и за счет диссипации. В работе [Calvet, Margerin, 2008] в качестве механизма формирования коды рассмотрена среда со случайно ориентированными анизотропными блоками кристаллического железа, размеры которых для Q = 450 составляют 0.3–1 км.

Другая точка зрения на механизм затухания волн во внутреннем ядре обосновывается в работе [Cormier, Li, 2002], в которой уширение импульса волны PKPDF связывают только с фактором рассеивания в случайной среде с 8% вариацией скорости продольных волн и размером неоднородности 10 км. Полученные нами данные не поддерживают этот вывод и указывают примерно на равный вклад в затухание волн PKiKP как за счет рассеяния, так и за счет неупругих свойств.

ВЫВОДЫ

1. На основе анализа огибающей коды отраженных от поверхности внутреннего ядра Земли волн PKiKP, возбужденных ядерными взрывами, и с использованием технологии альфа-шейп восстановления гладких кривых получены оценки фактора затухания Q = 447 ± 43 для девяти областей под Арктикой, Центральной Азией и Юго-Восточной Азией.

2. Механизм формирования коды с полученным значением Q может быть обусловлен случайными неоднородностями в верхней части внутреннего ядра с характерным линейным размером 1–4 км и вариациями скорости продольных волн порядка 1–3%

3. Показана эффективность применения метода альфа-шейп для построения гладкой огибающей коды волны PKiKP по сравнению со скользящим средним.

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

  1. Адушкин В.В., Овчинников В.М. О мозаичности отражающих свойств поверхности внутреннего ядра Земли // Докл. РАН. 2004. Т. 397. № 6. С. 815–817.

  2. Овчинников В.М., Краснощеков Д.Н., Каазик П.Б. Новое представление о границе между внешним и внутренним ядром Земли // Докл. РАН. 2007. Т. 417. № 3. С. 389–392.

  3. Aki K., Chouet B. Origin of coda waves: source attenuation and scattering effects // J. Geophys. Res. 1975. V. 80. P. 3322–3342.

  4. Attanayake J., Cormier V., de Silva S. Uppermost inner core seismic structure – new insights from body waveform inversion // Earth Planet. Sci. Lett. 2014. V. 311. P. 49–58.

  5. Calvet M., Margerin L. Constraints on grain size and stable iron phases in the uppermost inner core from multiple scattering modeling of seismic velocity and attenuation // Earth Planet. Sci. Lett. 2008. V. 267. № 1. P. 200–212.

  6. Koper K., Franks J., Dombrovskaya M. Evidence for small-scale heterogeneity in Earth’s inner core from a global study of PKiKP coda waves // Earth Planet. Sci. Lett. 2004. V. 228. № 3. 227–241.

  7. Krasnoshchekov D., Polishchuk V. k-Order alpha-hulls and alpha-shapes // Information Processing Letters. 2014. V. 114. № 1–2. P. 76–83.

  8. Krasnoshchekov D., Kaazik P., Ovtchinnikov V. Seismological evidence for mosaic structure of the surface of the Earth’s inner core // Nature. 2005. V. 435. P. 483–487.

  9. Krasnoshchekov D., Polishchuk V. Robust curve reconstruction with k-order α-shapes // IEEE Shape modelling. 2008. V. 4. P. 279–280. https://doi.org/10.1109/SMI.2008.4548006

  10. Leyton F., Koper K. Using PKiKP coda to determine inner core structure: 2. Determination of Qc // J. Geophys. Res. 2007. V. 112(B5). B05317. https://doi.org/10.1029/2006JB004370

  11. Li X., Cormier V. Frequency-dependent seismic attenuation in the inner core. 1. A viscoelastic interpretation // J. Geophys. Res. 2002. V. 107(B12). P. 2361. https://doi.org/10.1029/2002JB001795

  12. Nikkilä M., Polishchuk V., Krasnoshchekov D. Robust estimation of seismic coda shape // Geophys. J. Int. 2014. V. 197. P. 557–565.

  13. Nissen-Meyer T., van Driel M., Stähler, S.C., Hosseini K., Hempel S., Auer L., Colombi A., Fournier A. AxiSEM: Broadband 3-D seismic wavefields in axisymmetric media // Solid Earth. 2014. V. 5. № 1. P. 425–445.https://doi.org/10.5194/se-5-425-2014

  14. Oreshin S., Vinnik L. Heterogeneity and anisotropy of seismic attenuation in the inner core // Geophys. Res. Lett. 2004. V. 31. L0261.

  15. Peng Z., Koper K., Vidale J., Leyton F. Shearer P. Inner-core fine-scale structure from scattered waves recorded by LASA // J. Geophys. Res. 2008. V. 113(B9). B09312. https://doi.org/10.1029/2007JB005412

  16. Poupinet G., Kennett B. On the observation of high frequency PKiKP and its coda in Australia // Phys. Earth Planet. Inter. 2004. V. 146. № 3. P. 497–511.

  17. Sanborn C.J. Radiative3D. 2015. http://rainbow.phys.uconn.edu/geowiki/Radiative3D

  18. Vidale J., Earle P. Fine-scale heterogeneity in the Earth’s inner core // Nature 2000. V. 404. P. 273–275.

  19. Wu W., Irving J. Using PKiKP coda to study heterogeneity in the top layer of the inner core’s western hemisphere // Geophys. J. Int. 2017. V. 209. P. 672–687.

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