Доклады Российской академии наук. Науки о Земле, 2021, T. 496, № 1, стр. 63-66

Построение решения обратной задачи по ансамблю моделей на примере инверсии приемных функций

И. М. Алёшин 1*

1 Институт физики Земли им. О.Ю. Шмидта Российской академии наук
Москва, Россия

* E-mail: ima@ifz.ru

Поступила в редакцию 16.09.2020
После доработки 19.10.2020
Принята к публикации 19.10.2020

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

Аннотация

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

Ключевые слова: обратная задача, приемная функция

Метод приемных функций широко применяется для изучения свойств коры и верхней мантии [1]. Неотъемлемой частью метода является процедура инверсии. Она состоит в подборе таких значений сейсмических скоростей и плотности среды под станцией, которые позволят достаточно хорошо описать наблюдения. Обычно предполагают, что в окрестности станции среда однородна по латерали, и ее параметры зависят только от глубины. Тогда инверсию наблюденной волновой формы S0(t) можно выполнить в рамках слоистой модели: пачки упругих, однородных слоев на упругом, однородном полупространстве. В этом случае для расчета синтетической волновой формы S(L; t) применимы эффективные матричные методы. Здесь символом L = {li} обозначена совокупность параметров слоистой модели; li означает характеристики i-го слоя: ${{v}_{{Pi}}}$, ${{v}_{{Si}}}$, ρi, hi – скорости продольных и поперечных волн, плотность и толщина соответственно. Обратная задача может быть сформулирована как поиск модели L, которая минимизирует расстояние

(1)
$W\left( L \right) = \left\| {{{S}_{0}}\left( t \right),S\left( {L,t} \right)} \right\|.$

В качестве метрики обычно выбирают среднеквадратичное отклонение, однако возможен и другой выбор [3]. При вероятностном подходе решением задачи минимизации функционала (1) служит не одна модель, а функциональная зависимость f = f(L; S0(t)). Последняя имеет смысл распределения вероятности того, что для имеющихся измерений S0(t) модель L имеет заданные параметры. Суть таких методов сводится к следующему. На основе датчика случайных чисел по тем или иным правилам генерируется ансамбль из N моделей LN. На его основе строится плотность распределения вероятности f(L; S0(t)). Оптимальная модель обычно определяется как совокупность модельных параметров, вычисленных по распределению f. Это может быть математическое ожидание, медиана или мода распределения f(L; S0(t)).

Отличия разных вероятностных методов обусловлены способом создания ансамбля моделей, на основе которых строится искомое распределение. Для этого используют случайное однородное распределение моделей, правило Гиббса, цепи Маркова и др. [13]. Любая из этих процедур требует большого объема вычислений, экспоненциально растущего с увеличением числа слоев в составе модели. Поэтому при моделировании коры и даже коры и верхней мантии приходится ограничиваться небольшим (порядка 10) количеством слоев [13]. Это может оказаться недостаточным для адекватного описания среды, например, привести к появлению в модели фиктивных скачков сейсмических скоростей, обусловленных на самом деле выбором параметризации. Кроме того, малое количество слоев приводит к трудностям интерпретации результатов инверсии и при последующем использовании модели.

Описанный ниже метод позволяет избежать этих недостатков и вообще существенно снизить зависимость решения задачи от параметризации. Мы исходим их того, что в геофизических приложениях интерес представляют не сами параметры слоев, а некоторые их комбинации. В нашем случае речь идет о зависимости сейсмических скоростей среды и плотности от глубины D: ${{v}_{S}} = {{v}_{S}}(D)$, ${{v}_{P}} = {{v}_{P}}(D)$ и т.п. Поэтому вместо расчета плотности f(L; S0(t)) нужно строить плотность вероятности непосредственно для этих величин. Ниже для таких распределений мы будем использовать обозначение q = qn(D), где q означает одну из величин ${{v}_{P}}$, ${{v}_{S}}$, ${{v}_{P}}{\text{/}}{{v}_{S}}$ и т.п.

Искомое распределение представляет собой полевую величину Φq(Q, D), зависящую от двух переменных. Зависимость q = qn(D), составленная для n-й модели ансамбля LN, является однозначной функцией аргумента и представляет собой линию на плоскости (Q, D). Чтобы перейти к полевым величинам, введем сингулярное распределение

${{\varphi }_{q}}\left( {Q,D} \right) = \left( {1{\text{/}}N} \right)\mathop \sum \limits_{n = 1}^N \,\delta \left( {Q - {{q}_{n}}\left( D \right)} \right),$
где δ(x) – дельта-функция Дирака. Искомый результат получается после усреднения φq(Q, D) по интервалу ΔQ по формуле:

${{\Phi }_{q}}\left( {Q,D} \right) = \left( {1{\text{/}}\Delta Q} \right)\mathop \smallint \limits_{\Delta Q}^{} \,{{\varphi }_{q}}\left( {Q + \xi ,D} \right)~d\xi .$

Отметим, что определенная таким образом величина действительно имеет смысл плотности вероятности. В частности, интеграл по всем допустимым значениям переменной Q равен 1:

$\int {\Phi \left( {Q,D} \right)dQ} = 1,\quad \forall D.$

Как следует из определения qn(D), эта нормировка имеет место при любом значении глубины. Распределение Φq(Q, D) позволяет непосредственно рассчитать усредненную зависимость от глубины интересующих нас параметров среды. Например, для математического ожидания имеем:

$\left\langle q \right\rangle \left( D \right) = \int {Q~{{\Phi }_{q}}\left( {Q,D} \right)dQ} .$

Плотность Φq(Q, D) позволяют вычислить зависимость от глубины всех статистических моментов соответствующей величины. Поэтому совокупность плотностей Φq(Q, D), например, для $q = \{ {{v}_{P}},{{v}_{S}},r\} $ полностью определяет решение обратной задачи минимизации функционала (1).

Существенным моментом предлагаемого метода является метод оценки качества решения. В качестве критерия мы использовали усредненную модель 〈L〉: совокупность зависимостей ${{v}_{P}}(D)$, ${{v}_{S}}(D)$ и ρ(D), полученных в результате усреднения по соответствующим распределениям. Качество подбора можно оценить по среднеквадратичному отклонению от наблюдений синтетических сейсмограмм, рассчитанных по усредненной модели. Фактически, мы возвращаемся к задаче минимизации функционала (1), где под L теперь следует понимать усредненную модель 〈L〉.

Для иллюстрации в качестве исходных данных мы взяли приемные функции P- и S-волн, рассчитанные по модели, изображенной на рис. 1 пунктиром. Для ее построения в стандартную модель IASP91 [4] были добавлены два слоя с пониженными значениями скорости поперечных волн ${{v}_{S}}$ и высокими значениями параметра ${{v}_{P}}{\text{/}}{{v}_{S}}$. Первый слой мощностью приблизительно 3 км расположен непосредственно у поверхности. Второй, 50-километровый слой находится в верхней мантии, на глубине около 100 км. Отметим, что согласованное изменение ${{v}_{S}}$ и ${{v}_{P}}{\text{/}}{{v}_{S}}$ приводят к существенно меньшей модификации скоростей продольных волн в обоих слоях, так что в мантийном слое снижение значений ${{v}_{P}}$ с глубиной и вовсе отсутствует. Наконец, преобразованная модель была сглажена гауссовым фильтром шириной 5 км.

Рис. 1.

Плотности распределения сейсмических параметров среды в зависимости от глубины. Средние значения распределений изображены сплошной линией. Цвет линий – черный или белый – определяется фоном рисунка. В нижней части панелей приведены среднеквадратичные отклонения $\delta {{v}_{S}}$, $\delta ({{v}_{P}}{\text{/}}{{v}_{S}})$ и $\delta {{v}_{P}}$ усредненной модели от исходной. “Истинная” модель изображена пунктиром.

Ансамбль моделей был составлен на основе модели из 15 слоев переменной толщины, аналогично технике из статьи [3]. Мощность каждого слоя и скорости ${{v}_{P}}$, ${{v}_{S}}$ в нем выбирались случайно, на основе равномерного распределения. Полученная модель служила начальным приближением при локальной квазилинейной минимизации функционала (1). Затем процедура многократно повторялась. После 20 000 испытаний было отобрано около 1000 моделей, которым соответствуют меньшие значения целевого функционала (1). На их основе были получены плотности распределения вероятности для ${{v}_{P}}$, ${{v}_{S}}$ и ${{v}_{P}}{\text{/}}{{v}_{S}}$.

На рис. 1 видно, что выполненная инверсия позволила восстановить основные геофизически значимые особенности модели: слой пониженной скорости у поверхности и повышенное значение отношения ${{v}_{P}}{\text{/}}{{v}_{S}}$ в этом слое, плавное возрастание ${{v}_{S}}$ в коре, глубину границы Мохо и, наконец, слой пониженной скорости в мантии на глубине 100–130 км.

Убедимся, что усредненная модель $\left\langle L \right\rangle $ является приближенным решением задачи (1). Для этого следует рассчитать для этой модели синтетические приемные функции и сравнить их с наблюдениями. Отметим прежде, что усредненные сейсмические параметры среды непрерывно зависят от глубины. Поэтому для вычисления приемных функций мы аппроксимировали усредненную модель слоистой, толщина слоев равна 1 км. Результат расчета представлен на рис. 2. Видно хорошее соответствие результатов инверсии и исходных данных. Среднеквадратичное отклонение между ними составило 0.008 для обоих типов волн.

Рис. 2.

Зависимость приемных функций S волны (вверху) и P волны (внизу) от времени. Пунктиром изображены сейсмограммы, соответствующие “истинной” модели. Сплошные линии – волновые формы, рассчитанные по усредненной модели. Отсчет времени и амплитуд ведется по соответствующим значениям основной фазы.

Как и следовало ожидать, надежнее всего определяется ${{v}_{S}}$, значения ${{v}_{P}}$ восстанавливаются хуже. Можно надеяться, что качество подбора улучшится, если включить в инверсию дополнительные данные – невязки времен пробега телесейсмических фаз (см. [1, 5]), дисперсию поверхностных волн (см. [2, 3, 5] и пр.).

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

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

  1. Винник Л.П. Сейсмология приемных функций // Физика Земли. 2019. № 1. С. 16–27.

  2. Bodin T., Sambridge M., Tkalčić H., Arroucau P., Gallagher K., Rawlinson N. Transdimensional Inversion of Receiver Functions and Surface Wave Dispersion // J. Geophys. Res. Solid Earth. 2012. V. 117. № B2. P. 1978–2012.

  3. Hedjazian N., Bodin T., Métivier L. An Optimal Transport Approach to Linearized Inversion of Receiver Functions // Geophys. J. Int. 2019. V. 216. C. 130–147.

  4. Kennett B.L., Engdahl E.R. Traveltimes for Global Earthquake Location and Phase Identification // Geophysical J. Int. 1991. №. 2. V. 105. P. 429–465.

  5. Kozlovskaya E., Kosarev G.L., Aleshin I.M., Riznichenko O.Yu., Sanina I.A. Structure and Composition of the Crust and Upper Mantle of the Archean-Proterozoic Boundary in the Fennoscandian Shield Obtained by Joint Inversion of Receiver Function and Surface Wave Phase Velocity of Recording of the SVEKALAPKO Array // Geophys. J. Int. 2008. V. 175. P. 135–152.

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