Estimation of Frequencies and Splitting of ₂S₁ and ₃S₁ Modes from the IGETS Superconducting Gravimeter Network Data after the 2013 Sea of Okhotsk Earthquake

Full Text

ВВЕДЕНИЕ

Моды ₂S₁ и ₃S₁ относятся к низкочастотным модам спектра собственных колебаний Земли (СКЗ), по модели PREM их частоты составляют 0.40 и 0.94 мГц соответственно. Из-за вращения и несферичности Земли обе моды должны наблюдаться в виде триплета. Уточнение параметров мод СКЗ – вырожденной частоты и параметров расщепления – необходимо для совершенствования модели внутреннего строения Земли. Оценки расщепления низкочастотных мод (с частотой ниже 1 мГц) особенно важны, т.к. они более всего чувствительны к трехмерной структуре мантии и ядра.

Дополнительный интерес к модам ₂S₁ и ₃S₁ объясняется тем, что они являются следующими обертонами моды Шлихтера ₁S₁. При этом, если мода Шлихтера связана с поступательным движением внутреннего твердого ядра относительно жидкого внешнего, то ₂S₁ вызвана колебанием ядра как целого относительно мантии. Амплитуда моды ₂S₁ при гравиметрических наблюдениях должна быть в несколько раз выше моды Шлихтера, однако из-за гораздо меньшей добротности также сложна в наблюдении. Так как мода Шлихтера до сих пор достоверно экспериментально не обнаружена, то наблюдение и оценка ₂S₁ представляются важными с точки зрения разработки и тестирования перспективных методов, пригодных для поиска ₁S₁.

Впервые об экспериментальном наблюдении ₂S₁ было объявлено в работе [Rosat et al., 2003]. Поиск моды осуществлялся по данным сверхпроводящих гравиметров после землетрясения в Перу 23.06.2001 г. (M_w 8.4), для этого использовались 164-часовые записи 5 станций и процедура “мультистанционного” анализа (стекинга). Выбор данных станций был объяснен меньшими шумами указанных гравиметров в исследуемом частотном диапазоне. Амплитуда синглетов в стекинговом спектре составила около 2 нГал, а найденные значения их частот оказались близки к модели PREM. В работе [Rosat et al., 2006] были проведены новые оценки моды по измененному алгоритму.

Большинство опубликованных результатов по моде 2S1 связано с мощнейшим землетрясением на Суматре 2004 г. В частности, она наблюдалась по спектральным данным широкополосных сейсмологических станций сети FDSN [Roult et al., 2010]. При этом было отмечено, что мода возбуждается очень слабо, особенно центральный синглет, что часто не дает возможности точно идентифицировать частоты синглетов по индивидуальным спектрам. В работе [Deuss et al., 2010] приведены спектральные оценки моды, полученные на сверхпроводящем гравиметре обсерватории Black Forest (BFO, Германия). Амплитуда моды незначительно превышала уровень шума, расщепления на индивидуальных спектрах видно не было. В работе [Ding, Shen, 2013] выполнена оценка частоты моды с использованием метода “Optimal sequence estimation”. В нашем предыдущем исследовании [Milyukov, Vinogradov, 2023] оценивались частоты и кориолисов параметр расщепления моды по результатам обработки данных сети сверхпроводящих гравиметров IGETS [Boy, 2016] после землетрясения в Чили 27.02.2010 г.

Цель данной работы – оценить вырожденные частоты и параметры расщепления мод ₂S₁ и ₃S₁, возбужденных землетрясением в Охотском море 2013 г. Как и в работе [Milyukov, Vinogradov, 2023], для анализа использованы данные сети IGETS, а для их обработки был применен оптимальный алгоритм, разработанный авторами на основе метода максимального правдоподобия (ММП) и предложенного для поиска и оценки моды Шлихтера [Виноградов и др., 2019]. Ранее алгоритм был апробирован на данных долговременных наблюдений Баксанским лазерным интерферометром-деформографом [Милюков и др., 2020а; 2020б] для обнаружения моды Шлихтера ₁S₁, а также для указанной выше оценки моды ₂S₁ [Milyukov, Vinogradov, 2023]. По сравнению с последней работой алгоритм был модифицирован для устранения неоднозначности с определением параметра расщепления a и наличием оценки времени начала возбуждения моды. Достоинствами ММП-алгоритма являются его базирование на фундаментальных принципах теории оптимального приема сигналов на фоне шумов, что обеспечивает наиболее эффективное обнаружение сигнала в соответствии с выбранным критерием оптимальности Неймана–Пирсона, возможность определения периода моды и параметров ее расщепления, а также расчет характеристик обнаружения и достоверности получаемых оценок.

Землетрясение в Охотском море 24.05.2013 г. с магнитудой 8.3 M_w и глубиной гипоцентра 611 км является самым мощным глубокофокусным землетрясением за все время сейсмических наблюдений. По величине сейсмического момента (M = 4 · 10²¹ Н · м) оно занимает 8-ю позицию в списке самых значительных землетрясений века. При этом сейсмический момент Охотоморского землетрясения примерно на порядок меньше землетрясений на Суматре 2004 г. и в Японии 2011 г., что делает наблюдение моды ₂S₁ сложной задачей. По нашему расчету (см. ниже), ее амплитуды должны составлять около 3 нГал, т.е. находятся на границе чувствительности сверхпроводящих гравиметров.

Несмотря на большое количество публикаций, связанных с Охотоморским землетрясением, число работ, непосредственно посвященных оценке мод СКЗ, достаточно мало. В частности, можно отметить статью [Молоденский, Молоденская, 2015], где был проведен анализ частот и добротностей по сейсмометрическим данным сети GSN, [Милюков и др., 2018], где были получены оценки 50 основных тонов СКЗ по деформографическим наблюдениям на длиннобазовом интерферометре, а также в работе [Кузин и др., 2019], в которой исследовалось расщепление моды ₀S₂.

Оценки мод ₂S₁ и ₃S₁ после Охотоморского землетрясения в литературе отсутствуют.

ИСПОЛЬЗУЕМЫЕ ДАННЫЕ И ОСОБЕННОСТИ МОД ₂S₁ И ₃S₁

На дату землетрясения в базе сервиса IGETS представлены 30 записей гравиметров, которых имеют данные в формате “level3” – с коррекцией инструментального дрейфа, приливов, смещения полярной оси и длительности суток [Boy, 2016]. С целью исключения влияния локальных неоднородностей Земли, в работе были использованы станции, находящиеся в центре Европы (см. рис. 1): это гравиметры bf056 (Schiltach, Германия); bh044 (Bad Homburg, Германия); co025 (Conrad, Австрия); mb021 (Membach, Бельгия); mo034 (Moxa, Германия); pe050 (Pecny, Чехия); st026 (Strasbourg, Франция); tr005 (Trappes, Франция) и два гравиметра we029 и we030 (Wettzell, Германия). Гравиметры bf056, mo034, we029 и we030 имеют по 2 разных датчика, поэтому общее количество исследуемых записей составило 14.

Рис.  1. Расположение станций IGETS, данные которых использовались для оценки мод.

Расчет амплитуд возбуждения моды от Охотоморского землетрясения выполнялся на основе теории, изложенной в работе [Dahlen, Tromp, 1998]. Амплитуда возбуждения моды на поверхности Земли для гравиметра может быть записана в виде:

$A_r=\frac{3}{4\pi }\left(1+\frac{g}{2{\pi }^{2}R{f}_d^2}\right)U\left(r\right)\cdot A\left(\Theta ,\Phi \right),$ (1)

где: R – радиус Земли; g – ускорение свободного падения; f_d – вырожденная частота моды; $U\text{}\left(r\right)$ – собственная функция соответствующей моды; $A\left(\Theta ,\Phi \right)$ – функция, зависящая от эпицентрического радиуса $\Theta$ и азимута $\Phi$ между источником и приемником, а также тензора сейсмического момента землетрясения. Второе слагаемое в скобках учитывает влияние изменения гравитационного поля Земли в месте расположения гравиметра (т.н. “free-air” эффект). В полном виде указанные формулы были приведены в нашей работе [Milyukov, Vinogradov, 2023].

На рис. 2 представлены собственные функции $U\text{}\left(r\right)$ внутри Земли в зависимости от расстояния от ее центра, рассчитанные по модели PREM (здесь и далее для расчета собственных функций использовался комплекс программ Фредерика Симонса: https://geoweb.princeton.edu/people/simons/software.html). Из графика следует, что мода ₂S₁ связана с ядром Земли, а ₃S₁ – с внутренним ядром и верхней мантией; также видна значительная разница в амплитудах мод на поверхности Земли.

Рис.  2. Зависимость собственных функций $\mathit{U}\left(r\right)$ мод ₂S₁ и ₃S₁ от расстояния от центра Земли.

На рис. 3 показаны теоретические амплитуды возбуждения мод  для 14 записей гравиметров, рассчитанные с учетом их географического расположения, координат гипоцентра и тензора сейсмического момента. Данные по тензору момента Охотоморского землетрясения были взяты из базы проекта Global Centroid-Moment-Tensor [Ekström et al., 2012].

Рис.  3. Теоретические значения амплитуды моды ₂S₁ (слева) и ₃S₁ (справа) для центральноевропейских станций сети IGETS после Охотоморского землетрясения 2013 г.

Мода ₂S₁ имеет амплитуду около 3 нГал, ₃S₁ примерно в 17 раз больше.

Интервал дискретизации данных формата “level3” составляет 1 мин. Длительность используемых данных для поиска периодических затухающих сигналов обычно выбирается из произведения периода сигнала и его добротности. Для обеих исследуемых мод это произведение составляет около 15 тыс. мин.

Особенности регистрации и постобработки данных IGETS исключают возможность использования данных сразу после землетрясений, т.к. они заполняются теоретическим приливным сигналом. Поэтому для последующей обработки использовались данные, начиная с 5000-го отсчета после землетрясения. С учетом не очень высокой добротности ₂S₁ (Q ≈ 400) это приводит к усложнению задачи обнаружения слабых мод из-за потери части энергии полезного сигнала.

Вращение и отклонение от сферической формы Земли приводит к расщеплению мод – вместо одной моды на вырожденной частоте наблюдается мультиплет из нескольких синглетов, количество которых определяется угловым номером l, и равно 2l + 1. Частоты синглетов зависят от вырожденной частоты f_d и могут быть представлены в виде:

$f_m=f_d\left(1+a+bm+c{m}^2\right),m\in \left[-l,l\right],$ (2)

где m – азимутальный номер синглета.

Параметры a, b и c зависят от модели Земли и определяются следующими выражениями:

$a=\frac{1}{3}\left(1-k^2\chi \right){\left(\frac{\Omega }{f_d}\right)}^2+\frac{1}{2}\frac{1}{4{\pi }^2{f}_d^2}\left(\upsilon -4{\pi }^2{f}_d^2\tau \right);$ (3)

$b=\chi \frac{\Omega }{f_d};$ (4)

$c=-\frac{3}{2}\frac{1}{k^2}\frac{1}{4{\pi }^2{f}_d^2}\left(\upsilon -4{\pi }^2{f}_d^2\tau \right),$ (5)

где: $\Omega$ – частота вращения Земли (в Гц); $k=\sqrt{l\left(l+1\right)};$ $\chi ,$ $\upsilon$ и $\tau$ являются сложными интегралами по радиусу Земли от комбинаций собственных функций мод СКЗ и их производных, а также от распределения плотности $\rho \left(r\right),$ гидростатической эллиптичности, ее производной, сжатия и упругости внутри Земли. В полном виде данные интегралы приведены в работе [Dahlen, Tromp, 1998, c. 608], здесь мы приведем выражение только для $\chi :$

$\chi =\frac{1}{k^2}{\int }_0^R\rho \left(r\right)\left[V^2+2kUV\right]r^{2}dr.$ (6)

Из формул (3)–(5) видно, что входящие в (2) параметры расщепления не являются независимыми друг от друга: в частности, параметр a может быть выражен через b и c, причем его первое слагаемое в (3) зависит только от b, второе – только от c.

Выражая $\chi$ через b, $\left(\upsilon -4{\pi }^2{f}_d^2\tau \right)$ через c и подставляя в (2), получим (при $k=\sqrt{2}$ для мод вида _nS₁):

$a=\frac{1}{3}\left(1-2b\frac{f_d}{\Omega }\right){\left(\frac{\Omega }{f_d}\right)}^2\text{}-\frac{2}{3}c.$ (7)

Используя (7), можно исключить из (2) один из параметров расщепления, получив альтернативное представление для расщепления на синглеты, не содержащее a:

$f_m=f_d\left[1+\frac{1}{3}{\left(\frac{\Omega }{f_d}\right)}^2+\left(m-\frac{2}{3}\frac{\Omega }{f_d}\right)b+\frac{1}{3}c\right].$ (8)

Таким образом, частоты триплета в рассматриваемом приближении латерально однородной модели Земли полностью определяются тремя параметрами: вырожденной частотой и двумя параметрами расщепления b и c. Представляет интерес относительные значения входящих в выражение (8) параметров для разных мод. Характерным показателем малости является отношение частоты вращения Земли к частоте моды $\frac{\Omega }{f_d}.$ Анализ показывает, что для моды ₂S₁ возможным отклонением параметра c от теоретического значения можно пренебречь, т.к. величина связанной с ним поправки к частоте в сотни раз меньше, чем от b. С ростом частоты влияние c увеличивается, но незначительно – для моды ₃S₁ поправка за его счет достигает лишь 8% от величины поправки, вызванной b.

Ранее при оценке мод мы предполагали, что начало их возбуждения совпадает с началом землетрясения. Однако это не совсем корректно, т.к. процесс высвобождения сейсмической энергии реального землетрясения занимает некоторое время – для землетрясений большой магнитуды оно может составлять сотни секунд [Goldberg et al., 2022]. Суммарный сейсмический момент рассчитывается как интеграл по так называемой скорости реализации момента, измеряемой в Н · м/сек. Максимум скорости момента обычно близок к середине временнóй функции источника и именно это время более правильно считать началом возбуждения моды. При наблюдении мод СКЗ это должно проявляться в форме начального временнóго (фазового) сдвига, учет которого может быть особенно важен для короткопериодных мод.

Амплитуда синглетов зависит только от широты места наблюдения, а фаза боковых синглетов – от его долготы [Cummins, 1991]. Амплитуды синглетов триплета для мод вида _nS₁ были приведены в [Milyukov, Vinogradov, 2023]. С учетом добавления $t_0,$ итоговую формулу для описания возбуждения моды можно записать в следующем виде:

$S\left(t,f_d,b,c,t_0\right)=a_{-1}\left(t\right)+a_0\left(t\right)+a_1\left(t\right),$ (9)

$a_{-1}=A_r{e}^{-\frac{\pi f_d}{Q}\left(t_0+t\right)}\frac{1}{2}\sqrt{\frac{3}{2\pi }}\sin \theta \cdot \cos \left(2\pi f_{-1}\left(t_0+t\right)-\phi \right);$ 

$a_0=A_r{e}^{-\frac{\pi f_d}{Q}\left(t_0+t\right)}\frac{1}{2}\sqrt{\frac{3}{\pi }}\cos \theta \cdot \cos \left(2\pi f_0\left(t_0+t\right)\right);$ (10)

$a_1=-A_r{e}^{-\frac{\pi f_d}{Q}\left(t_0+t\right)}\frac{1}{2}\sqrt{\frac{3}{2\pi }}\sin \theta \cdot \cos \left(2\pi f_1\left(t_0+t\right)+\phi \right),$

где: $A_r$ – амплитуда, возбуждаемая вырожденной модой без учета расщепления (см. рис. 3); $\theta ,\phi$ – географические координаты (коширота и долгота) приемника; Q – добротность моды; t – время; t₀ – максимум скорости реализации сейсмического момента для землетрясения (начало возбуждения моды).

Теоретические параметры исследуемых мод приведены в табл. 1. Частоты рассчитаны по модели PREM, параметры расщепления – взяты из работы [Dahlen, Sailor, 1979].

Таблица 1. Теоретические параметры мод

На предварительном этапе работы был выполнен спектральный анализ исходных данных с использованием метода периодограмм Уэлча. Благодаря большой амплитуде (десятки нГал) мода ₃S₁ хорошо видна в спектрах всех станций, причем с расщеплением (см. рис. 4). Однако на спектрах выделяются только боковые синглеты – по-видимому, это связано с недостаточным разрешением спектрального анализа, ограниченного достаточно быстрым затуханием моды. В диапазоне частот моды 2S1 видно большое число пиков, вызванных сейсмическим шумом, что делает достоверную оценку частот синглетов после Охотоморского землетрясения практически невозможной.

Рис.  4. Спектральная плотность мощности по 14 гравиметрам после Охотоморского землетрясения: (а) – мода ₂S₁; (б) – мода ₃S₁. Пунктиром показаны теоретические значения синглетов, рассчитанные по параметрам из табл. 1.

АЛГОРИТМ ОБРАБОТКИ

Обнаружение моды ₂S₁ по гравиметрическим данным относится к задачам обнаружения слабого сигнала на фоне помех.

Ранее авторами был предложен оптимальный алгоритм для одновременного обнаружения и оценки параметров моды Шлихтера с учетом ее свойств и свойств сейсмического шума, основанный на методе максимального правдоподобия. Так как форма сигналов всех мод вида _nS₁ одинакова, то аналогичный алгоритм может быть использован для моды ₂S₁. Подробно алгоритм был изложен в работе [Виноградов и др., 2019], основные положения применительно к использованию гравиметрических данных IGETS – в работе [Milyukov, Vinogradov, 2023], поэтому здесь мы кратко опишем основные принципы алгоритма.

Метод максимального правдоподобия сводится к построению так называемой достаточной статистики [Сосулин, 1992], равной логарифму отношения правдоподобия. Отношение правдоподобия пропорционально отношению плотности вероятности в присутствии полезного сигнала к плотности вероятности в его отсутствии. Достаточная статистика Z является функцией наблюдаемой реализации случайного процесса (исследуемого гравиметрического сигнала, содержащего как полезный сигнал – возбуждение моды, так и сейсмический шум), позволяющей найти оптимальное решающее правило для принятия решения о наличии или отсутствии сигнала посредством сравнения с пороговым значением h. При превышении статистикой порога (Z > h) делается вывод о наличии сигнала в наблюдаемой реализации, в противоположном случае – о его отсутствии.

Точное решение задачи оптимального приема возможно лишь в случае, если шум имеет гауссовское распределение. Такое решение реализуется через согласованную фильтрацию. Анализ шумов гравиметров сети IGETS в интересующем нас спектральном диапазоне низкочастотных мод СКЗ (0.4–1.0 мГц), показал, что шум может считаться гауссовским. Достаточно узкий частотный диапазон, который определяется максимально возможным разбросом боковых синглетов, также позволяет считать спектральную плотность сигнала постоянной, т.е. шум – белый.

ММП позволяет вместе с обнаружением сигнала выполнять оценку его неизвестных параметров. Для этого приемник делается многоканальным – для каждой комбинации значений неизвестных параметров (в нашем случае их три: вырожденная частота f_d, параметр расщепления b и задержка времени t₀) через соответствующий ей согласованный фильтр вычисляется величина достаточной статистики $Z\left(f_d\text{}, b\text{},t_0\right)$ и находится ее максимум. Те значения $f_d, b,t_0,$ которые максимизируют достаточную статистику, и являются наиболее вероятными (оптимальными) оценками сигнала.

Как было указано ранее, наиболее важными характеристиками моды являются вырожденная частота f_d и параметр расщепления b, который определяет удаление боковых синглетов триплета друг от друга. Параметр расщепления b и вырожденная частота моды связаны друг с другом, что делает обязательным включение обоих этих параметров в структуру оптимального обнаружителя. В то же время параметр расщепления c практически не влияет на результаты оценки, т.к. намного (в десятки раз) меньше b. Это позволяет считать его значение в задаче оптимальной оценки постоянным и известным (равным теоретическому значению). Параметр a определяется через b и c по формуле (7).

Структура обработки данных по ММП-алгоритму для оценки мод ₂S₁ и ₃S₁ приведена на рис. 5. В целом она аналогична приведенной в работе [Milyukov, Vinogradov, 2023], однако в ней не требуется сравнение статистики с порогом, т.к. сигнал заведомо присутствует во входных данных. Кроме этого, используется другой, более сложный согласованный фильтр, соответствующий формулам (10) и учитывающий взаимную связь параметров a, b и c, а также начальный сдвиг времени t₀.

Рис.  5. Схема алгоритма для оптимальной оценки параметров мод на основе ММП: $\mathit{g}\left(t\right)$ – гравиметри-ческий сигнал; СФ – согласованный фильтр; MAX – выбор максимального значения.

Диапазон возможных значений параметров задается исходя из теоретических представлений о сигнале и имеющихся экспериментальных данных. Поиск максимума в пространстве трех параметров требует значительного времени вычисления, поэтому задачу оценки можно оптимизировать, разбив на три этапа: на первом этапе с достаточно грубой детализацией по частоте и расщеплению определяется параметр t₀; затем его значение используется уже как известный параметр для оценки f_d и b со средней детализацией, позволяющей локализовать область абсолютного максимума статистики и сузить диапазон возможных параметров f_d и b; на заключительном этапе алгоритм применяется для локализованной зоны максимума с наибольшим разрешением.

Перед применением алгоритма данные проходили предварительную обработку, в частности, осуществлялась фильтрация полосовым фильтром в диапазоне 395–415 мкГц и 939–953 мкГц для мод ₂S₁ и ₃S₁ соответственно. Спектр СКЗ имеет достаточно плотную структуру, поэтому полосовая фильтрация важна с целью исключения взаимодействия исследуемой моды с другими близкими модами (₀T₂, ₀S₃, ₂S₂, ₁S₃ и ₀S₆). Кроме этого, фильтрация позволяет исключить воздействие приливных гармоник, снизить влияние техногенных и сейсмических шумов. Следует иметь в виду, что использование фильтров может приводить к появлению фазового сдвига относительно исходного сигнала. Поскольку ММП осуществляет обработку сигнала во временнóй области, фазовый сдвиг может оказать существенное влияние на результат. Поэтому для корректного преобразования фазы использовалась процедуры фильтрации с нулевой фазой – после прямой фильтрации отфильтрованная последовательность разворачивается и вновь пропускается через фильтр.

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

Описанный выше алгоритм был применен к каждой из 14 записей гравиметрических данных IGETS, в результате было получено 14 достаточных статистик. На первом этапе обработки были определены оптимальные задержки времени между началом землетрясения и началом возбуждения – они оказались равными нулю для обеих мод.

Были выбраны следующие диапазоны возможных значений: вырожденные частоты 402..406 мкГц для ₂S₁ и 942..946 мкГц для ₃S₁, параметры расщепления b 5..30 · 10^–3 для моды ₂S₁ и 0.5...3.0 · 10^–3 для ₃S₁. Основная (финальная) обработка осуществлялась с дискретизацией 0.001 мкГц по частоте и 5 · 10^–6 по b.

С целью уменьшения погрешностей, вносимых сейсмическими шумами, для оценки моды было выполнено усреднение достаточных статистик по всем 14 станциям с весовыми коэффициентами, обратно пропорциональными дисперсиям шумов на входе каждого согласованного фильтра. Эта процедура особенно важна для моды 2S1 в силу малости ее амплитуды.

Итоговые средневзвешенные достаточные статистики $Z\left(f_d,b\right)$ для обеих мод показаны на рис. 6.

Рис.  6. Средневзвешенная достаточная статистика по центральноевропейским станциям сети IGETS после Охотоморского землетрясения: (а) – для моды ₂S₁; (б) – для моды ₃S₁. Абсолютные максимумы статистики обозначены знаком “+”.

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

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

Таблица 2. Результаты оценки параметров мод и их сравнение с моделью PREM

Так как в большинстве работ по СКЗ приводятся оценки не параметров расщепления, а частоты синглетов, то для сравнения полученных нами результатов были по формуле (8) рассчитаны частоты триплета. Они приведены в табл. 3. и табл. 4 вместе с результатами, полученными другими исследователями.

Таблица 3. Значения частот моды ₂S₁ в мкГц по данным разных наблюдений

Примечания: * – из работы [Ding, Shen, 2013]; ** – нет данных.

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

Важнейшая особенность алгоритма – непосредственная оценка вырожденной частоты и параметра расщепления моды b, зависящих от внутреннего строения Земли. Это является существенным преимуществом перед другими методами, основанными на определении частот отдельных синглетов, т.к. ввиду неполноты систем уравнений, связывающих параметры расщепления, спектральные наблюдения не позволяют однозначно определить значения интегралов $\chi ,$ $\upsilon$ и $\tau$ (см. формулы (3)–(5)).

Сделанная в работе модификация алгоритма, учитывающая временной сдвиг и взаимную связь параметров расщепления a, b, c и вырожденной частоты, позволяет исключить неоднозначность оценок и возможные систематические ошибки. Как и в работе [Milyukov, Vinogradov, 2023], использование алгоритма ММП позволило уверенно определить возбуждение моды ₂S₁ при уровне сигнала 3 нГал – т.е. отношение сигнал–шум на входе составило не более 5 · 10^–3. Для моды ₃S₁ данный алгоритм позволил оценить вырожденную частоту с погрешностью, в несколько раз меньшей, чем в выполненных ранее работах.

Полученные новые значения вырожденной частоты и параметра расщепления b позволяют определить по формуле (7) параметр расщепления a. Для мод ₂S₁ и ₃S₁ он оказался равен – 0.16 и 0.269 соответственно. Найденные значения существенно отличаются от приведенных в табл. 1 [Dahlen, Sailor, 1979] по причине других значений f_d и b.

При анализе данных была проверена гипотеза о том, что возбуждение моды начинается не в момент начала землетрясения, а с некоторой временнóй задержкой t0, которая может быть связана с особенностями очага и обычно зависит от магнитуды землетрясения. Примененный при обработке данных модифицированный вариант алгоритма ММП позволил оценить эту задержку, однако она оказалась равной нулю (максимум статистики достигается, если начало возбуждения моды совпадет с началом землетрясения). С учетом времени дискретизации использованных данных можно сделать вывод о том, что задержка времени не превышает 1 мин. Этот вывод находится в согласии с результатами других исследователей – во многих публикациях отмечалось, что суммарная длительность очагового развития Охотоморского землетрясения, хотя и состоявшего из двух стадий, была относительно короткой – от 30 до 70 с [Wei S. et al., 2013; Абубакиров и др., 2015; Кузин и др., 2017; Чеброва и др., 2019]. Таким образом можно считать, что максимум функции реализации момента (временнóй функции источника) отстоял от начала землетрясения не более чем на половину минуты.

Полученные оценки вырожденной частоты и параметра расщепления моды ₂S₁ в пределах погрешностей соответствуют модели PREM, а также сделанным нами ранее оценкам по Чилийскому землетрясению 27.02.2010 г. [Milyukov, Vinogradov, 2023]. Достаточно большая погрешность при определении частот синглетов моды ₂S₁ связана с большой погрешностью определения параметра расщепления b. При этом необходимо учитывать, что практически все приведенные в табл. 3 оценки сделаны для гораздо более мощных землетрясений на Суматре 2004 г. и в Японии 2011 г., после которых мода ₂S₁ наблюдалась даже в обычных спектрах.

Что касается оценок моды ₃S₁, то они заметно отличаются от теоретических – частота на 0.5 мкГц ниже, а параметр расщепления b на 0.4 · 10^–3 больше. Достаточно малая погрешность, которую обеспечивает ММП, позволяет считать эти разницы значимыми. Для понимания возможных причин более низкой частоты моды можно воспользоваться т.н. ядрами Фреше, показывающими зависимость изменения частоты моды $\delta \omega$ при изменении параметров модели Земли. Для скоростей волн объемного сжатия  и сдвига  вариации частоты могут быть записаны в следующем виде:

$\delta \omega =\underset{0}{\overset{R}{}}\left(\delta \kappa K_{\kappa }+\delta \mu K_{\mu }\right)dr,$

где: $\omega =2\pi f$ – круговая частота моды; $\kappa$ и $\mu$ – модули объемного сжатия и сдвига; $K_{\kappa }$ и $K_{\mu }$ – соответствующие им ядра Фреше; R – радиус Земли.

Ядра Фреше могут быть рассчитаны для конкретной моды СКЗ через ее собственные функции $U\text{}\left(r\right)$и $V\text{}\left(r\right),$ их производные по радиусу $\dot{U}\text{}\left(r\right),$ $\dot{V}\text{}\left(r\right)$ и плотность $\rho \left(r\right)$ [Dahlen, Tromp, 1998]:

${\left(\frac{\partial \omega }{\partial \alpha }\right)}_{\beta ,\rho ,d}\text{}=K_{\alpha }\text{}=2\rho \alpha K_{\kappa },$

${\left(\frac{\partial \omega }{\partial \beta }\right)}_{\alpha ,\rho ,d}\text{}=K_{\beta }=2\rho \beta \left(K_{\mu }\text{}-\frac{4}{3}{K}_{\kappa }\right),$

$K_{\kappa }=\frac{1}{2\omega }{\left(r\dot{U}+2U-kV\right)}^2,$

$\begin{array}{c}{K}_{\mu }\text{}=\frac{1}{2\omega }\times \\ \times \left[\frac{1}{3}{\left(2r\dot{U}-2U-kV\right)}^2+{\left(r\dot{V}\text{}-V-kU\right)}^2\text{}+\left(k^2\text{}-2\right)V^2\right].\end{array}$

Радиальные зависимости ядер Фреше от центра Земли для исследуемых мод приведены на рис. 7.

Рис.  7. Ядра Фреше для мод ₂S₁ (а) и ₃S₁ (б).

График на рис. 7б показывает, что причиной меньшего значения частоты моды ₃S₁ могут быть меньшие значения скоростей волн сжатия во внутреннем жидком ядре и средней мантии и/или меньшие значения скоростей волн сдвига в твердом внешнем ядре и средней мантии, чем приняты в модели PREM.

Увеличенное по сравнению с PREM значение параметра расщепления может быть связано, в соответствии с формулами (4) и (6), с бóльшей величиной плотности внутри Земли либо с бóльшими значениями собственных функций моды ₃S₁. Также не следует исключать возможное влияние на оценку расщепления локальных неоднородностей в верхней мантии.

ВЫВОДЫ

В работе впервые выполнена оценка параметров мод ₂S₁ и ₃S₁, возбужденных глубокофокусным Охотоморским землетрясением 2013 г. Для анализа использовались гравиметрические данные 14 сверхпроводящих гравиметров сети IGETS, расположенных в Центральной Европе.

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

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

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

Применение алгоритма позволило значительно увеличить точность оценки вырожденной частоты и параметров расщепления моды ₃S₁, выявить ее отклонение от теоретических значений модели PREM. Представлена зависимость между изменением частоты моды относительно модели и величинами скоростей волн объемного сжатия и сдвига во внутреннем или внешнем ядре и мантии.

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

ФИНАНСИРОВАНИЕ РАБОТЫ

Работа выполнена при поддержке Российского научного фонда (проект № 23-27-00237).

About the authors

M. P. Vinogradov

Author for correspondence.
Email: vinogradovmp@my.msu.ru
Russian Federation, Moscow

V. K. Milyukov

Email: vinogradovmp@my.msu.ru
Russian Federation, Moscow

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Location of the IGETS stations whose data were used for mode estimation.

Download (330KB)

Indexing metadata

3. Fig. 2. Dependence of the eigenfunctions of the 2S1 and 3S1 modes on the distance from the Earth's center.

Download (184KB)

Indexing metadata

4. Fig. 3. Theoretical values of the 2S1 (left) and 3S1 (right) mode amplitude for the Central European stations of the IGETS network after the 2013 Okhotomorsk earthquake.

Download (396KB)

Indexing metadata

5. Fig. 4. Power spectral density from 14 gravimeters after the Okhotomorsk earthquake: (a) - mode 2S1; (b) - mode 3S1. The dotted line shows the theoretical values of singlet calculated from the parameters in Table 1.

Download (633KB)

Indexing metadata

6. Fig. 5. Scheme of the algorithm for optimal estimation of mode parameters based on MMP: - gravimetric signal; SF - matched filter; MAX - maximum value selection.

Download (141KB)

Indexing metadata

7. Fig. 6. Weighted average sufficient statistics for Central European stations of IGETS network after the Okhotomorsk earthquake: (a) - for mode 2S1; (b) - for mode 3S1. Absolute maxima of statistics are marked with "+" sign.

Download (953KB)

Indexing metadata

8. Fig. 7. Fréchet kernels for the 2S1 (a) and 3S1 (b) modes.

Download (390KB)

Indexing metadata