Физика Земли, 2021, № 4, стр. 133-141

ВВЕДЕНИЕ

Широкое распространение метода функции приемника (или приемных функций) для получения глубинных скоростных моделей обусловлено тем, что метод позволяет строить детальные одномерные скоростные разрезы до глубин 250–300 км по данным одиночной широкополосной сейсмической станции. Приемные функции P и S представляют собой обменные поперечную или продольную волны, образованные конверсией на контрастных границах телесейсмических P- или S-волн соответственно. Это позволяет не только определить положение глубинных сейсмических границ (зоны фазовых переходов 410–660 км, астеносферного слоя, зоны Мохо и т.д.), но и получить скоростной разрез земной коры и верхней мантии. Обзор выделения и интерпретации обменных волн приведен в статье [Винник, 2019].

Применение приемных функций существенно осложняется, если под сейсмостанцией присутствует приповерхностный слой с низкими значениями сейсмических скоростей. Такого рода слой различной мощности существует на большей части территории Земли. Как правило, он связан с осадочными породами, однако возможны и другие причины его образования, такие как трещиноватость и/или флюидонасыщенность пород [Vinnik et al., 2002; Aleshin et al., 2006]. Сложности обработки связаны с тем, что из-за отражений от верхней и нижней границ слоя образуются кратные обменные волны, сравнимые по амплитуде с первичной. Поэтому слой, а тем более несколько слоев пониженной скорости разной мощности, формируют сложную волновую картину. Это приводит к “маскировке” вступлений обменных фаз от нижележащих коровых границ, а в некоторых случаях – границ переходной зоны мантии (см. Приложение).

Мощность низкоскоростного слоя лежит в пределах от нескольких десятков метров до нескольких километров, что сравнимо с пределом разрешения метода приемных функций. В таких случаях для построения скоростной модели приходится использовать специальные приемы. Например, в работах [Санина и др., 2014; Гоев и др., 2018] при изучении центральной части Русской платформы подбор модели осуществлялся в два этапа. Сначала верхняя часть модели подбиралась по начальному участку приемной функции. Для определения скоростной структуры осадочного чехла проводился расчет “короткопериодных” функций приемника (в диапазоне частот 0.5–3.4 Гц). В результате инверсии первых 10 с восстанавливалась детальная скоростная структура приповерхностных слоев. Затем выполнялась еще одна инверсия уже по “длиннопериодным” данным (в диапазоне периодов от 3 с), чтобы построить нижнюю часть разреза. При этом рассчитанная ранее модель осадочного чехла фиксировалась. Кроме сложности реализации такой подход имеет ряд очевидных недостатков, связанных с эмпирическим характером процедур. В частности, фиксация верхней части модели, вообще говоря, искажает значения параметров глубоких слоев, причем степень этого искажения с трудом поддается оценке.

Обычно чтобы уменьшить степень неопределенности расчета скоростной модели необходимо использовать дополнительные независимые сейсмические данные. Наиболее естественно в качестве таких данных привлечь приемные функции S-волн [Винник, 2019], дисперсионные зависимости фазовой скорости поверхностных волн [Julia et al., 2000], спектральные характеристики параметра эллиптичности волн Рэлея [Мунирова, Яновская, 2001] или комбинацию дисперсии и эллиптичности [Berg et al., 2019]. Однако упомянутые данные являются длиннопериодными, поэтому они бесполезны для решения интересующей нас проблемы определения сейсмических и геометрических параметров тонких слоев верхней части коры.

Одним из хорошо известных локальных методов определения скоростных характеристик приповерхностной части разреза является подход, описанный в работе [Nakamura, 1989]. Он заключается в анализе спектрального отношения горизонтальной компоненты (H) к вертикальной компоненте (V) сейсмического шума. За последние несколько десятков лет метод спектрального отношения H/V получил значительное развитие и сегодня активно применяется, в том числе и для построения скоростного разреза под станцией (см., напр., [Arai, Tokimatsu, 2004; Sánchez-Sesma et al., 2011]). И приемная функция, и отношение H/V определяются строением среды под станцией. Причем отношение H/V зависит от свойств самой верхней части разреза, которая и вызывает сложности интерпретации начального участка приемной функции. Поэтому совместное использование этих данных позволит естественным образом учесть влияние низкоскоростного слоя при построении модели среды.

Характерный частотный спектр сейсмических шумов, используемых при наблюдении H/V, составляет 0.5–20 Гц. В то же время высокочастотная граница спектра приемных функций не превышает 10 Гц. Эти различия в частотном составе могут, вообще говоря, привести к тому, что сигналы будет невозможно описать в рамках единой модели. Предпосылкой к совместной интерпретации приемной функции и частотного спектра отношения H/V является наличие тонкого (порядка 50–100 м) слоя, в несколько раз отличающегося по упругим свойствам от подстилающей его породы. С одной стороны, это приводит к формированию соответствующей обменной фазы, а с другой стороны, положение и форма максимума спектра H/V определяются параметрами тонкого слоя. Подобная скоростная модель была получена в работе [Гоев и др., 2018] по данным приемных функций станции “Монаково” (Нижегородская область, Россия, Русская платформа). В предлагаемой работе мы попытались обосновать возможность совместной интерпретации приемной функции и частотного спектра отношения H/V сейсмического шума.

МОДЕЛЬ СРЕДЫ ПО ДАННЫМ ПРИЕМНЫХ ФУНКЦИЙ

Скоростная модель среды под группой станций “Монаково” по приемным функциям в терминах плоскослоистой модели была получена ранее в работе [Гоев и др., 2018] минимизацией функционала невязки:

здесь RF^(syn), RF^(obs) означают соответственно синтетическую и наблюденную приемные функции; $\vec {m}$ – параметры модели (совокупность сейсмических скоростей, плотностей и толщин слоев); суммирование ведется по всем моментам времени t, в которых производились измерения; вспомогательная функция $q\left( {\vec {m}} \right)$, зависящая только от параметров модели, используется для стабилизации инверсии [Kosarev et al., 1987]. Было показано, что непосредственно под станцией находится слой пониженной скорости толщиной около двух километров, сформированный осадочными породами. Этот слой имеет сложную структуру, важной частью которой является расположенный непосредственно у поверхности низкоскоростной слой толщиной несколько десятков метров (рис. 1, пунктир).

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

СПЕКТР ОТНОШЕНИЯ ГОРИЗОНТАЛЬНОЙ И ВЕРТИКАЛЬНОЙ КОМПОНЕНТ: ТЕОРИЯ И НАБЛЮДЕНИЯ

В качестве исходных данных для оценки спектрального отношения H/V были использованы записи сейсмического шума, зарегистрированного трехкомпонентной сейсмической станцией “Монаково”. Станция оборудована широкополосным датчиком RefTek 151-60 с предельным периодом 60 с. Запись трех компонент сейсмического сигнала выполнялась регистратором RefTek 130 с частотой опроса 100 Гц на канал. Для расчета искомого спектрального соотношения были отобраны 10 получасовых интервалов, записанных в ночное время и не содержащих сейсмических событий. По каждой из компонент отобранных записей был сформирован спектр Фурье амплитуд колебательных скоростей. Горизонтальный спектр $H\left( f \right)$ вычислялся как среднеквадратичное значение спектров компонент N-S $Y\left( f \right)$ и E-W $X\left( f \right)$ по формуле:

Полученные таким образом горизонтальные спектры Фурье, а также соответствующие вертикальные спектры $V\left( f \right)$ сглаживались скользящим прямоугольным окном шириной 0.1 Гц, после чего были рассчитаны спектры ${{{{\chi }}}^{{\left( {obs} \right)}}}\left( f \right) = {{H\left( f \right)} \mathord{\left/ {\vphantom {{H\left( f \right)} {V\left( f \right)}}} \right. \kern-0em} {V\left( f \right)}}$ для каждой из записей. Окончательный вариант был получен усреднением индивидуальных записей. Он представлен на рис. 2 толстой сплошной линией.

В соответствии с подходом авторов статьи [Nakamura, 1989] контрастная граница осадочного чехла и фундамента проявляется как резонанс в спектральной зависимости отношения H/V, так называемый “резонанс SH-волн”. Вообще говоря, в горизонтальных компонентах шумовых колебаний кроме фундаментальной могут присутствовать более высокие рэлеевские моды, поверхностные волны Лява и объемные волны. Поэтому для максимально полной интерпретации сейсмического шума в терминах слоистой модели следует использовать основанный на функции Грина подход (см., напр., [Sánchez-Sesma et al., 2011]), требующий данных измерительной сети, составленной из группы сейсмических станций. В настоящем исследовании мы по нескольким причинам полагаем, что основной вклад в отношение спектров H/V сейсмического шума вносит именно фундаментальная мода поверхностной волны Рэлея. Во-первых, в ряде экспериментальных и теоретических работ [Bonnefoy-Claudet et al., 2006; Malischewsky, Scherbaum, 2004; Hobiger et al., 2013] было показано, что в диапазоне частот, близком к резонансной частоте ${{f}_{0}} \approx {{{{V}_{S}}} \mathord{\left/ {\vphantom {{{{V}_{S}}} {~4h}}} \right. \kern-0em} {~4h}}$, именно амплитуды рэлеевских волн наиболее чувствительны к вариациям параметров разреза. В предыдущей формуле V_S и h означают скорость поперечных волн в осадочном слое и его толщину соответственно. Во-вторых, сейсмический шум от различных источников, расположенных вблизи поверхности [Bonnefoy-Claudet et al., 2006], сформирован в основном поверхностными волнами, причем высокочастотные моды затухают быстрее. Наконец, ниже будет показано, что использованный нами подход позволяет удовлетворительно описать наблюдаемый спектр. Это косвенно подтверждает правомерность сделанных допущений, а также позволяет оценить величину ошибки интерпретации, исходя из невязки между наблюденным и синтетическим спектрами.

Существует несколько методов описания поверхностных волн в слоистой среде. Здесь мы использовали метод пропагатора в варианте, предложенном авторами работ [Abo-Zena, 1979; Menke, 1979]. Это позволяет обойти известные проблемы оригинального подхода Томсона–Хаскелла [Haskell, 1953], возникающие в вычислениях высокочастотной части спектра. Для фиксированной фазовой скорости c из дисперсионного соотношения ${{E}_{5}}\left( {c,{{\omega }}} \right) = 0$ вычислялись собственное значение циклической частоты ${{\omega }}$, а также горизонтальная и вертикальная компоненты соответствующего собственного вектора смещения. Здесь и далее использованы обозначения из работы [Menke, 1979]. Искомое спектральное выражение для модуля отношения компонент, называемое эллиптичностью, имеет вид:

Функция ${{{{\chi }}}^{{\left( {syn} \right)}}}\left( {{\omega }} \right) = {{{{\chi }}}^{{\left( {syn} \right)}}}\left( {{{\omega }};\vec {m}} \right)$ параметрически зависит от многомерного вектора $\vec {m}.$ Его компоненты составлены из параметров модели: скоростей продольных и поперечных волн, плотности в каждом слое и толщин слоев.

Полученное выражение позволяет нам восстановить параметры слоистой модели под сейсмической станцией по записям сейсмических шумов. Математически, это означает, что нам необходимо найти набор параметров $\vec {m}$, который минимизирует невязку между экспериментально наблюдаемыми значениями спектрального отношения ${{{{\chi }}}^{{\left( {obs} \right)}}}$ и рассчитанными теоретически, в рамках слоистой модели, значениями эллиптичности фундаментальной моды поверхностной волны Рэлея ${{{{\chi }}}^{{\left( {syn} \right)}}}\left( {{{\omega }};\vec {m}} \right)$.

Целевая функция, подлежащая минимизации, определена на основе ${{L}_{2}}$-нормы:

Суммирование под знаком радикала выполняется по всем циклическим частотам, для которых проводились измерения. Явный вид вспомогательной функции $q\left( {\vec {m}} \right)$ будет выписан ниже, при описании процедуры инверсии.

Как показали результаты численных экспериментов, при изменении параметров модели значение эллиптичности ${{{{\chi }}}^{{\left( {syn} \right)}}}$ может принимать крайне большие значения. Это означает, что изменение параметров приводит к изменению направления движения частиц. По этой причине при вычислении значений функции ${{{{\chi }}}^{{\left( {syn} \right)}}}$, превышающие заданное значение ${{\chi }}_{M}^{{\left( {syn} \right)}}$, заменялись этой константой. Ее значение выбрано исходя из максимума измеренного спектра ${{\chi }}_{M}^{{\left( {obs} \right)}}$, а именно, ${{\chi }}_{M}^{{\left( {syn} \right)}} = 1.01~{{\chi }}_{M}^{{\left( {obs} \right)}}$.

СКОРОСТНАЯ МОДЕЛЬ СРЕДЫ С УЧЕТОМ H/V СПЕКТРА

Для достижения поставленной цели работы покажем, что модель ${{\vec {m}}_{0}}$, полученная при обращении приемной функции на сейсмической станции “Монаково”, не противоречит спектральным характеристикам сейсмического шума, измеренного на той же станции. Чтобы это продемонстрировать, необходимо сначала найти параметры модели ${{\vec {m}}_{1}}$, при которых целевой функционал (3) имеет минимум. Наконец, необходимо убедиться, что найденная модель ${{\vec {m}}_{1}}$ близка к модели ${{\vec {m}}_{0}}$ и не приводит к значительному изменению формы приемной функции. Математически, это означает, что минимумы этих двух функционалов достигаются на близких моделях. Это условие необходимо для совместной интерпретации приемной функции и H/V-спектра.

Аргументом в пользу существования модели ${{\vec {m}}_{1}}$, близкой к ${{\vec {m}}_{0}}$, может служить пунктирная кривая на рис. 2. Она является результатом расчета по формуле (2) спектра отношения H/V ${{\chi }^{{\left( {syn} \right)}}}\left( {\omega ;{{{\vec {m}}}_{0}}} \right)$ для скоростной модели ${{\vec {m}}_{0}}$, полученной по приемной функции. Верхняя часть этой модели приведена на рис. 1 также пунктиром. Видно, что синтетическая спектральная кривая близка к наблюдениям по форме, но немного смещена в высокочастотную область. Для поиска интересующей нас модели мы воспользовались симплекс-методом, адаптированным для решения задач нелинейной оптимизации [Lagarias et al., 1998], используя модель ${{\vec {m}}_{0}}$ в качестве начального приближения. Так как этот метод основан на линеаризации функционала вблизи заданной точки, то найденное решение соответствует ближайшему локальному экстремуму функционала (3).

Как уже отмечалось, плоскослоистая модель $\vec {m}$ определяется параметрами составляющих ее слоев – скоростью упругих продольных ${{V}_{P}}$ и поперечных ${{V}_{S}}$ волн, плотностью и толщиной. В предлагаемой работе, как и при инверсии приемной функции в работе [Гоев и др., 2018], варьировались только скорости поперечных волн. Остальные параметры – плотность и отношение ${{V}_{P}}$/${{V}_{S}}$ – были зафиксированы в соответствии с исходной моделью. Чтобы избежать решений, в которых значение ${{V}_{{S,i}}}$ в i-м слое выходит за физически допустимые пределы $\left[ {V_{{S,i}}^{{\left( {{\text{min}}} \right)}},V_{{S,i}}^{{\left( {{\text{max}}} \right)}}} \right]$, мы воспользовались следующим приемом. К исходному функционалу были добавлены слагаемые, равные нулю, если ${{V}_{{S,i}}}$ лежит внутри указанного интервала, и быстро растущие, если ${{V}_{{S,i}}}$ выходит за его пределы. Кроме того, для стабилизации решения в функционал было введено слагаемое, пропорциональное численному аналогу модуля градиента скорости поперечных волн $\left| {{{\partial {{V}_{S}}} \mathord{\left/ {\vphantom {{\partial {{V}_{S}}} {\partial z}}} \right. \kern-0em} {\partial z}}} \right|$. Таким образом, функция $q\left( {\vec {m}} \right)$ из выражения (3) принимает следующий вид:

где ${{\theta }}\left( x \right)$ означает ступенчатую функцию Хевисайда; a, b и с – эмпирические константы.

Модель, полученная в результате инверсии отношения H/V сейсмического шума, и рассчитанный по ней спектр эллиптичности волны Рэлея изображены толстыми сплошными линиями на рис. 1 и рис. 2 соответственно. Чтобы оценить надежность определения параметров, было выполнено численное моделирование. Измеренный спектр отношения H/V случайным образом модифицировался. Каждая спектральная амплитуда изменялась на случайное число, из равномерного распределения в интервале от –0.05 до + 0.05. Среднее изменение амплитуд составило 5%. Модифицированный спектр использовался в качестве входных данных в процедуре инверсии. Эта процедура была повторена около 1000 раз. Полученные таким образом модели использованы для расчета среднеквадратичных отклонений значений поперечных скоростей в каждом слое. На рис. 1 рассчитанные отклонения изображены горизонтальными отрезками.

Из графиков можно сделать вывод, что спектр, соответствующий оптимальной модели, мало отличается по форме от спектра модели стартовой. Однако максимумы измеренного и синтетического спектров теперь совпадают, и, в целом, спектральные кривые удовлетворительно согласуются друг с другом. Основную роль в минимизации играет резонансная часть спектра. В свою очередь, величина резонансной частоты f₀ определяется свойствами верхнего, приповерхностного слоя. Поэтому не удивительно, что в результате инверсии значимые изменения параметров произошли только в верхней части разреза. На средних частотах (0.5–1 Гц) имеется заметное несоответствие наблюдаемого и синтетического спектров. По-видимому, это несоответствие связано как с выбором параметризации модели, так и пренебрежением вклада в спектр волн Лява, высших рэлеевских мод и объемных волн.

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

ЗАКЛЮЧЕНИЕ

Подводя итог, можно сказать, что несмотря на значительные отличия в частотном составе сейсмического шума и обменных волн, оба типа данных могут быть объяснены близкими по форме моделями. Подчеркнем, что в методах используются принципиально разные исходные данные, которые разделяются не только по частотному составу, но и по типу волн. Совместная интерпретация этих данных имеет смысл, если непосредственно под станцией имеется низкоскоростной, например осадочный, слой со сложной структурой. В терминах данных наблюдений это означает, что в спектре шума должен присутствовать “резонанс”. Это, собственно говоря, и указывает на то, что под станцией находится тонкий (около 100 м) слой с пониженной скоростью поперечных волн ${{V}_{S}}$, так как в слое большей мощности резонансная частота лежит вне частотного диапазона наблюдения. В то же время, именно наличие такого слоя приводит к серьезным трудностям при интерпретации волновых форм приемных функций. Из приведенного анализа также следует, что спектральные данные сейсмического шума дают возможность ограничить параметры самого верхнего слоя. Этих данных, по-видимому, недостаточно для организации последовательной процедуры построения сейсмического разреза целиком, без частичной фиксации.