\documentclass[11pt,reqno]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage[colorlinks]{hyperref}


\usepackage{titlesec}
\titlespacing{\section}{0pt}{2.5ex}{1.5ex}
\titlespacing{\subsection}{0pt}{1.5ex}{1ex}
\titlespacing{\subsubsection}{0pt}{1.5ex}{1ex}
\titleformat{\section}{\large\bfseries\centering}{\thesection}{1em}{}
\titleformat{\subsection}[runin]{\bfseries}{\thesubsection.}{0.5em}{}[.\mbox{\ }]
\titleformat{\subsubsection}[runin]{\bfseries}{\thesubsubsection.}{0.4em}{}[.\mbox{\ }]



%%%%%---------- ORCiD
\newcommand\orcidicon[1]{\href{https://orcid.org/#1}{\includegraphics[scale=0.02]{orcid.pdf}}}

\def\udcs{550.34} %Здесь автор определяет УДК своей работы
\def\mscs{86-08} %Здесь автор определяет классификаторы AMS своей работы
\setcounter{page}{144}

\newtheorem{lemma}{Лемма}
\newtheorem{theorem}{Теорема}
\newtheorem{definition}{Определение}
\newtheorem{corollary}{Следствие}
\newtheorem{proposition}{Предложение}


\definecolor{myblue}{rgb}{0.03, 0.27, 0.79}

\hypersetup{linkcolor=blue,filecolor=black,urlcolor=blue, citecolor=blue}

\def\logo{\textcolor{blue}{{\bf\Huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}}

\def\semrtop
     {
  \vbox{
     \noindent\logo \hspace{-3mm}  \parbox{12cm}{     \begin{center}
     {\LARGE \textsf{СИБИРСКИЕ ЭЛЕКТРОННЫЕ}} \\[2mm]
     {\LARGE \textsf{МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ}} \\[2mm]
     {\Large \textsf{Siberian Electronic Mathematical Reports}} \\[1mm]
     {\Large\tt{\textcolor{blue}{http://semr.math.nsc.ru}}}\\[1mm]
     \raisebox{1ex}{ISSN 1813-3304 }
     \vspace*{-3mm}
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
}
     \vspace{-1mm}
     \noindent
      \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa!} \\
     \hline\hline
     \end{tabular}

     \vspace{3mm}
     {\flushleft\it Том 20, № 2, стр. 144--145 (2023) \hspace{\fill}{\rm\small УДК \udcs \hspace{-4mm}}}
     \newline
     	{\rm\small \textcolor{blue}{https://doi.org/10.33048/semi.2023.20.???}}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa\;}{\rm\small MSC\ \mscs }
  }%\vbox
}



\begin{document}
\renewcommand{\refname}{References}


\thispagestyle{empty}
\title[\sc Синхронная инверсия с точным решением Цеппритца]{\large Синхронная сейсмическая инверсия на основе точного решения уравнений Цеппритца для восстановления упругих свойств среды}
\author[Д.К. ДМИТРАЧКОВ, М.И. ПРОТАСОВ]{{\bf \large \textsc{Д.К. Дмитрачков}}\protect\orcidicon{0009-0005-1059-2999}, {\bf \large \textsc{М.И. Протасов}}\protect\orcidicon{0000-0003-0449-1639}\\ \\
\textit{Представлено} \textsc{Д.К. Дмитрачковым}
%\textit{(заполняется редактором)}
}%
\address{Danil Konstantinovich Dmitrachkov  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} Institute of Petroleum Geology and Geophysics,
\newline\hphantom{iii} Koptyug St. 3,
\newline\hphantom{iii} 630090, Novosibirsk, Russia}%
\email{\textcolor{blue}{dmitrachkovdk@ipgg.sbras.ru}}%

\address{Maxim Igorevich Protasov  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} Institute of Petroleum Geology and Geophysics,
\newline\hphantom{iii} Koptyug St. 3,
\newline\hphantom{iii} 630090, Novosibirsk, Russia}%
\email{\textcolor{blue}{protasovmi@ipgg.sbras.ru}}%

\thanks{\sc Dmitrachkov, D.K., Protasov, M.I.,
Prestack seismic inversion based on the exact solution of Zoeppritz equations for medium elastic properties reconstruction}
\thanks{\copyright \ 2025 Дмитрачков Д.К., Протасов М.И.}
\thanks{\rm Работа поддержана РНФ (грант 21-71-20002-П)}
\thanks{\it Поступила 1 января 2025 г., опубликована 31 декабря 2025 г.}%


\semrtop \vspace{1cm}
\maketitle {\small
\begin{quote}
\noindent{\bf Abstract:} Seismic inversion allows to estimate the elastic medium properties. Simultaneous seismic inversion provides a detailed model of the elastic medium properties using the results of prestack seismic data processing in the time domain. The relationship of reflection coefficients with elastic properties is given by the Zoeppritz equations. Classic simultaneous inversion approaches are based on the approximations of the Zoeppritz equations exact solution, which are valid for small contrasts of properties in the medium and for small angles of incidence of the wave. In this paper we propose an algorithm for simultaneous inversion based on the exact solution of the Zoeppritz equations. The proposed algorithm is compared with classic simultaneous inversion algorithms using realistic synthetic models.\medskip

\noindent{\bf Keywords:} prestack inversion, Zoeppritz equations.
 \end{quote}
}

\bigskip

\section{Введение}
Сейсмическая инверсия представляет собой инструмент для определения коллекторских свойств горных пород. Эта процедура широко распространена на практике, так как она позволяет получить детальную модель упругих свойств исследуемого геологического объекта по результатам временной обработки сейсмических данных. Синхронная сейсмическая инверсия даёт возможность построить модель среды по результатам обработки сейсмических данных во временной области до суммирования \cite{Russell_1991}. Алгоритмы сейсмической инверсии основаны на одномерной свёрточной модели сейсмической трассы \cite{Russell_1988, Ampilov_2009}. В такой постановке предполагается, что сейсмическая трасса является результатом свёртки сейсмического импульса с трассой коэффициентов отражения. В случае синхронной инверсии зависимость коэффициента отражения от упругих свойств среды даётся системой уравнений Цёппритца. Чтобы упростить работу с решением этих уравнений, на практике используются различные аппроксимации точного решения \cite{Bortfeld_1961, Shuey_1985}. На основе данных аппроксимаций в настоящее время реализовано множество алгоритмов синхронной инверсии, нашедших своё применение в коммерческом программном обеспечении. Преимуществом аппроксимаций точного решения является возможность линеаризации оператора прямой задачи, что облегчает процесс оптимизации целевой функции при решении обратной задачи. К недостаткам аппроксимаций можно отнести их применимость для небольших углов падения волны (как правило, до 30°) и небольших контрастов упругих свойств, а также необходимость дополнительной априорной информации, такой как гладкий низкочастотный тренд упругих свойств среды или оценка величины отношения скорости продольных волн к скорости поперечных волн. Погрешность расчёта значений коэффициентов отражения в прямой задаче оказывает влияние на результат решения обратной задачи. Отметим, что даже при выполнении всех условий, обеспечивающих справедливость аппроксимаций, они по-прежнему остаются применимы для сравнительно небольших углов падения, однако большая доля информации о распределении упругих свойств среды содержится в сейсмических трассах, соответствующих большим углам падения \cite{Romanenko_2010}.

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


\section{Синхронная инверсия на основе точного решения системы уравнений Цеппритца}
Свёрточная модель сейсмической трассы может быть представлена уравнением вида
\begin{equation}    \label{eq:conv-1d}
    W(t, \theta_{ip}) * R_{pp}(m, \theta_{ip}) = S(t, \theta_{ip}),
\end{equation}
где $\theta_{ip}$~---угол падения плоской волны, $W$~--- сейсмический импульс, $R_{pp}$~--- коэффициент отражения продольной волны, $m$~--- модель упругой среды, $S$~--- сейсмическая трасса. В качестве модели упругих свойств среды $m$ обычно выступает набор $(v_p, v_s, \rho)$ скорости продольных волн, скорости поперечных волн и плотности соответственно. В синхронной инверсии имеется набор углов падения, каждому из которых соответствует свой сейсмический импульс и своя сейсмическая трасса. Требуется определить такую модель среды, которая бы обеспечивала выполнение равенств вида \eqref{eq:conv-1d} для всех имеющихся углов падения. От уравнений \eqref{eq:conv-1d} переходят к задаче минимизации целевого функционала вида
\begin{equation}    \label{eq:cost-function}
    f(m) = \sum \limits_{i=1}^{N_a}\|W_{\theta_{ip}^i} R_{pp}(m, \theta_{ip}^i) - S(\theta_{ip}^i)\|_2^2
    + \sum \limits_{k=1}^M \alpha_k \Omega_k(m),
\end{equation}
где $N_a$~--- количество углов падения, для которых имеются сейсмические трассы, $\alpha_k > 0$~--- параметры регуляризации, $\Omega_k(m)$~--- стабилизаторы задачи. Модель $m$, доставляющая минимум целевому функционалу \eqref{eq:cost-function}, будет являться решением обратной задачи. Стабилизаторы дают возможность учесть необходимую априорную информацию о решении и повысить устойчивость обращения.

Рассмотрим среду, состоящую из двух слоёв и горизонтальной границы между ними. Верхний слой описывается набором упругих свойств $(v_{p1}, v_{s1}, \rho)$, а нижний~--- набором $(v_{p2}, v_{s2}, \rho)$ соответственно. Пусть плоская продольная волна распространяется в верхнем слое и падает на границу между слоями под углом $\theta_{ip}$. В общем случае образуется четыре волны~--- продольная отражённая, продольная преломлённая, а также две обменные волны, то есть поперечная отражённая и поперечная преломлённая. Соответствующие коэффициенты отражения и преломления связаны с упругими свойствами слоёв системой уравнений Цёппритца, в параметризации $(v_p, v_s, \rho)$ имеющей вид
\begin{equation}    \label{eq:zoeppritz-system}
    Z(m) R = y(m),
\end{equation}
где
\begin{equation}
\begin{gathered}
     Z(m) = \left(
    \begin{array}{cccc}
    -\sin{\theta_{ip}}
    & -\cos{\theta_{rs}}
    &  \sin{\theta_{tp}}
    &  \cos{\theta_{ts}}\\
    
    \cos{\theta_{ip}}
    & -\sin{\theta_{rs}}
    &  \cos{\theta_{tp}}
    & -\sin{\theta_{ts}} \\
    
     -\rho_1 v_{p1} \cos{2\theta_{rs}}
    & \rho_1 v_{s1} \sin{2\theta_{rs}}
    & \rho_2 v_{p2} \cos{2\theta_{ts}}
    &-\rho_2 v_{s2} \sin{2\theta_{ts}}\\
    
      \dfrac{\rho_1 v_{s1}^2 \sin{2\theta_{ip}}}{v_{p1}}
    & \rho_1 v_{s1} \cos{2\theta_{rs}}
    & \dfrac{\rho_2 v_{s2}^2 \sin{2\theta_{tp}}}{v_{p2}}
    & \rho_2 v_{s2} \cos{2\theta_{ts}}
    \end{array}
    \right),
    \\
    R = \left(
    \begin{array}{c}
    R_{pp}\\
    R_{ps}\\
    T_{pp}\\
    T_{ps}
    \end{array}
    \right),
    \quad
    y(m) = \left(
    \begin{array}{c}
    \sin{\theta_{ip}}\\
    \cos{\theta_{ip}}\\
    \rho_1 v_{p1} \cos{2\theta_{rs}}\\
    \dfrac{\rho_1 v_{s1}^2 \sin{2\theta_{ip}}}{v_{p1}}
    \end{array}
    \right).
\end{gathered}
\end{equation}
Углы отражения и преломления связаны законом Снеллиуса в виде
\begin{equation}
    \frac{\sin{\theta_{ip}}}{v_{p1}}
    = \frac{\sin{\theta_{rs}}}{v_{s1}}
    = \frac{\sin{\theta_{tp}}}{v_{p2}}
    = \frac{\sin{\theta_{ts}}}{v_{s2}},
\end{equation}
где $\theta_{rs}$~--- угол отражения обменной поперечной волны, $\theta_{tp}$~--- угол преломления продольной волны, $\theta_{ts}$~--- угол преломления обменной поперечной волны. Таким образом, система уравнений \eqref{eq:zoeppritz-system} задаёт коэффициент отражения $R_{pp}$ как неявную функцию, зависящую от шести упругих параметров среды.

Разрешая систему уравнений \eqref{eq:zoeppritz-system}, можно получить явную формулу зависимости коэффициента отражения $R_{pp}$ от упругих параметров среды. Поскольку явная формула громоздка и сложна для дальнейшей работы, были разработаны различные аппроксимации для коэффициента отражения. Одной из самых распространённых аппроксимаций является аппроксимация Шуе, имеющая вид
\begin{equation}    \label{eq:shuey-approximation}
    R_{pp}(\theta_{ip}) = R_{pp, 0} + G \sin^2 \theta_{ip} + F (\tg^2 \theta_{ip} - \sin^2 \theta_{ip}),
\end{equation}
где
\begin{equation}
\begin{gathered}
    R_{pp, 0} = \frac{1}{2} \left( \frac{\Delta v_p}{v_p} + \frac{\Delta \rho}{\rho} \right),
    \\
    G = \frac{1}{2} \frac{\Delta v_p}{v_p}
    - 2\frac{v_s^2}{v_p^2} \left( \frac{\Delta \rho}{\rho} + 2\frac{\Delta v_s}{v_s} \right),
    \\
    F = \frac{1}{2} \frac{\Delta v_p}{v_p}.
\end{gathered}
\end{equation}
В данной аппроксимации двухслойная модель рассматривается как сумма гладкого тренда и контраста свойств на границе в виде
\begin{equation}    \label{eq:model-hf-lf}
\begin{gathered}
    v_p = \frac{v_{p2} + v_{p1}}{2},
    \quad
    v_s = \frac{v_{s2} + v_{s1}}{2},
    \quad
    \rho = \frac{\rho_2 + \rho_1}{2},
    \\
    \Delta v_p = v_{p2} - v_{p1},
    \quad
    \Delta v_s = v_{s2} - v_{s1},
    \quad
    \Delta \rho = \rho_2 - \rho_1.
\end{gathered}
\end{equation}
В рамках аппроксимаций точного решения системы уравнений Цёппритца предполагается, что угол падения $\theta_{ip}$ и контрасты упругих свойств среды между слоями малы.

Система уравнений Цёппритца может быть переписана в рамках иной параметризации, например, используя импедансы продольной и поперечной волн вместо скоростей. Помимо аппроксимации Шуе, существуют и другие, например, аппроксимация Бортфилда или аппроксимация Фатти в терминах импедансов и плотности. Система уравнений Цёппритца задаёт зависимость коэффициента отражения $R_{pp}$ от упругих параметров среды в виде неявной и нелинейной функции, однако точное решение системы Цёппритца справедливо для всех докритических углов падения волны и для любых контрастов свойств среды. Аппроксимации точного решения задают эту зависимость как явную и линейную, однако при этом для аппроксимаций требуется дополнительная информация, такая как гладкий тренд параметров среды, а также соблюдение предположений о малости угла падения и контрастов свойств. При нарушении этих предположений коэффициент отражения, вычисляемый согласно аппроксимациям точного решения, может содержать большую погрешность, то есть прямая задача свёрточного моделирования для расчёта синтетической сейсмической трассы будет решаться с ошибкой, что отрицательно сказывается на качестве результата решения обратной задачи.

Точное решение системы уравнений Цёппритца может быть использовано без дополнительной информации в виде гладких трендов параметров среды и применимо для больших углов падения и контрастов свойств, однако явная формула слишком громоздка и неудобна для работы. Более того, для поиска минимума функционала \eqref{eq:cost-function} хорошо подходят градиентные методы, но тогда возникает необходимость расчёта производных коэффициента отражения по упругим параметрам среды. В данной работе для расчёта производных предлагается использовать формулу производной неявной функции. Пусть $m_j$~--- один из шести упругих параметров, определяющих значение коэффициента отражения $R_{pp}(m)$. Дифференцируя каждое уравнение \eqref{eq:zoeppritz-system}, получаем уравнения
\begin{equation}
    \frac{\partial Z(m)}{\partial m_j}R + Z(m)\frac{\partial R}{\partial m_j} = \frac{\partial y(m)}{\partial m_j}.
\end{equation}
Выражая вектор производных коэффициентов отражения и преломления, получаем формулу
\begin{equation}    \label{eq:rpp-derivative}
    \frac{\partial R}{\partial m_j} = Z^{-1}(m) \left(\frac{\partial y(m)}{\partial m_j} - \frac{\partial Z(m)}{\partial m_j}R\right),
\end{equation}
где производные от матрицы $Z(m)$ и вектора правой части $y(m)$ понимаются поэлементно, а вектор коэффициентов отражения $R$ выражается как
\begin{equation}
    R = Z^{-1}(m) y(m).
\end{equation}
Таким образом, для вычисления производных коэффициента отражения $R_{pp}$ по переменным $m_j$, определяющим модель упругой среды, требуется обращать одну и ту же матрицу $Z(m)$ с различными правыми частями, поэтому при расчёте градиента целевого функционала выгодно использовать прямой метод обращения матрицы. Для обращения матриц малого размера существуют вычислительно эффективные формулы вычисления обратной матрицы напрямую. Численные эксперименты показывают, что для данной задачи этот подход является наиболее быстрым и численно устойчивым.

Для учёта априорной информации о решении можно использовать регуляризацию. Опишем, как выглядит целевой функционал синхронной инверсии на основе точного решения уравнений Цёппритца.

В предположении о том, что среда состоит из конечного числа слоёв с однородными свойствами, вводятся стабилизаторы, обеспечивающие разреженный вид трассы коэффициентов отражения и блочную структуру решения в виде
\begin{equation}    \label{eq:stab-tv}
\begin{gathered}
    \Omega_{r}(m) = \sum \limits_{i=1}^{N_a}\|R_{pp}(m, \theta_{ip}^i)\|_1,
    \\
    \Omega_{\nabla}(m) = \|\nabla v_p\|_1 + \|\nabla v_s\|_1 + \|\nabla \rho\|_1.
\end{gathered}
\end{equation}
При выборе нормы $L_1$ в результате обращения будет получаться наиболее разреженная трасса коэффициентов отражения и наиболее разреженные трассы градиентов модельных параметров, то есть модель, получающаяся в результате обращения, будет иметь блочную структуру. В численной реализации используется приближение
\begin{equation}
    |x| \approx \sqrt{x^2 + \delta^2},
\end{equation}
где $\delta \approx 10^{-6}$~--- малый параметр. Это обеспечивает гладкость целевой функции и позволяет сохранить блочную структуру модели.

В сейсмических трассах отсутствует информация о гладкой низкочастотной составляющей модели, и для учёта тренда модельных параметров в процессе инверсии используется стабилизатор вида
\begin{equation}
    \Omega_{m_0}(m) = \|v_p - v_{p,0}\|_2^2 + \|v_s - v_{s,0}\|_2^2 + \|\rho - \rho_0\|_2^2,
\end{equation}
где $(v_{p,0}, v_{s,0}, \rho_0)$~--- заданные гладкие тренды упругих параметров. Выбор нормы $L_2$ в этом стабилизаторе позволяет итоговому решению отклоняться от заданного тренда, обеспечивая достаточно свободы для учёта контрастных границ, вовлекая при этом заданный тренд модели в решение.

Поскольку в основе сейсмической инверсии лежит одномерная свёрточная модель сейсмической трассы, то при обращении необходимо учитывать информацию о соседних трассах для обеспечения пространственной связности решения. Для этого рассматривается стабилизатор вида
\begin{equation}
    \Omega_{n}(m) = \|v_p - \hat{v}_p\|_2^2 + \|v_s - \hat{v}_s\|_2^2 + \|\rho - \hat{\rho}\|_2^2,
\end{equation}
где $(\hat{v}_p, \hat{v}_s, \hat{\rho})$~--- решение, полученное для соседней трассы модели.

Таким образом, итоговый целевой функционал синхронной инверсии на основе точного решения уравнений Цёппритца имеет вид
\begin{equation}
\begin{split}
    f(m) &= \sum \limits_{i=1}^{N_a}\|W_{\theta_{ip}^i} R_{pp}(m, \theta_{ip}^i) - S(\theta_{ip}^i)\|_2^2
    \\
    &+ \alpha_1 \Omega_{r}(m) + \alpha_2 \Omega_{\nabla}(m) + \alpha_3 \Omega_{m_0}(m) + \alpha_4 \Omega_{n}(m).
\end{split}
\end{equation}
В качестве параметризации модели $m$ выбран набор параметров $(v_p, v_s, \rho)$. Для использования другой параметризации необходимо соответствующим образом модифицировать вид стабилизаторов, а способ расчёта производных коэффициента отражения по параметрам модели остаётся прежним. Для численной оптимизации используется квазиньютоновский метод BFGS, который требует вычисления градиента целевой функции.


Упругая среда может быть описана различными наборами параметров, и все эти наборы связаны друг с другом. Так, например, популярна параметризация с помощью скоростей и плотности $(v_p, v_s, \rho)$. Наряду с ней можно рассмотреть параметризацию через импедансы и плотность в виде $(I_p, I_s, \rho)$. Одни параметры могут быть пересчитаны в другие по известным формулам. Импедансы, как известно из литературы, восстанавливаются наиболее устойчиво \cite{Silvestrov_2011}. Синхронная инверсия на основе точного решения уравнений Цёппритца может быть сформулирована в любой параметризации. Например, для параметризации $(I_p, I_s, \rho)$ система уравнений Цёппритца примет вид
\begin{equation}
\begin{gathered}
    Z(m) = \left(
    \begin{array}{cccc}
    -\sin{\theta_{ip}}
    & -\cos{\theta_{rs}}
    &  \sin{\theta_{tp}}
    &  \cos{\theta_{ts}}\\
    
    \cos{\theta_{ip}}
    & -\sin{\theta_{rs}}
    &  \cos{\theta_{tp}}
    & -\sin{\theta_{ts}} \\
    
     -I_{p1} \cos{2\theta_{rs}}
    & I_{s1} \sin{2\theta_{rs}}
    & I_{p2} \cos{2\theta_{ts}}
    &-I_{s2} \sin{2\theta_{ts}}\\
    
      \dfrac{I_{s1}^2 \sin{2\theta_{ip}}}{I_{p1}}
    & I_{s1} \cos{2\theta_{rs}}
    & \dfrac{I_{s2}^2 \sin{2\theta_{tp}}}{I_{p2}}
    & I_{s2} \cos{2\theta_{ts}}
    \end{array}
    \right),
    \\
    y(m) = \left(
    \begin{array}{c}
    \sin{\theta_{ip}}\\
    \cos{\theta_{ip}}\\
    I_{p1} \cos{2\theta_{rs}}\\
    \dfrac{I_{s1}^2 \sin{2\theta_{ip}}}{I_{p1}}
    \end{array}
    \right),
\end{gathered}
\end{equation}
а закон Снеллиуса запишется в виде
\begin{equation}
    \frac{\rho_{1} \sin{\theta_{ip}}}{I_{p1}}
    = \frac{\rho_{1} \sin{\theta_{rs}}}{I_{s1}}
    = \frac{\rho_{2} \sin{\theta_{tp}}}{I_{p2}}
    = \frac{\rho_{2} \sin{\theta_{ts}}}{I_{s2}}.
\end{equation}
Для вычисления градиента целевого функционала производная коэффициента отражения по параметрам $(I_p, I_s, \rho)$ может быть вычислена по формуле \eqref{eq:rpp-derivative}, как и для любой другой параметризации. 


\section{Сравнение с синхронной инверсией на основе аппроксимаций}
Как правило, в синхронной инверсии в качестве оператора прямой задачи выступают аппроксимации точного решения уравнений Цёппритца, представляющие собой линейные операторы относительно контрастов свойств модели на границах. Рассмотрим синхронную сейсмическую инверсию на основе аппроксимации Шуе. Расчёт коэффициентов отражения по формуле \eqref{eq:shuey-approximation} можно представить в матрично-векторном виде как
\begin{equation}
    R_{pp}(m, \theta_{ip}) = G_{\theta_{ip}} m,
\end{equation}
где $G_{\theta_{ip}}$~--- матрица, учитывающая в себе гладкий тренд модели и угол падения волны, а $m = (\Delta v_p, \Delta v_s, \Delta \rho)$~--- параметризация модели. Для обеспечения блочной структуры модели используются стабилизаторы, аналогичные стабилизаторам \eqref{eq:stab-tv}. После того, как в результате инверсии определены контрасты свойств $(\Delta v_p, \Delta v_s, \Delta \rho)$, по формулам \eqref{eq:model-hf-lf} восстанавливается модель упругих свойств среды $(v_p, v_s, \rho)$.

Рассмотрим на примере, как ведут себя аппроксимации и точное решение для различных сред. На рисунке~\ref{fig:rpp-approx-errors} представлено поведение коэффициента отражения, рассчитанного согласно точному решению системы уравнений Цёппритца и рассчитанного согласно аппроксимациям Шуе и Бортфилда. На рисунках~\ref{fig:rpp-approx-errors}a и~\ref{fig:rpp-approx-errors}б контраст упругих свойств между двумя слоями мал, и аппроксимации справедливы на достаточно широком диапазоне углов. Однако на рисунках~\ref{fig:rpp-approx-errors}в и ~\ref{fig:rpp-approx-errors}г контраст свойств среды выбран достаточно большим, и аппроксимации отличаются от точного решения даже для небольших углов падения волны.
\begin{figure}[ht]
    \includegraphics[width=1.0\linewidth]{figures/rpp_approximations_errors.pdf}
    \caption{Коэффициент отражения по точному решению системы уравнений Цёппритца и по аппроксимациям Шуе и Бортфилда для слабоконтрастных (а, б) и контрастных (в, г) сред}
    \label{fig:rpp-approx-errors}
\end{figure}
При работе с моделями, в которых имеется большой контраст упругих свойств, алгоритмы синхронной сейсмической инверсии, основанные на аппроксимациях, будут проигрывать в точности в силу того, что приближённый оператор прямой задачи при больших контрастах упругих свойств может значительно отличаться от истинного. Также для аппроксимаций необходим гладкий тренд упругих параметров среды, который будет влиять на итоговое решение и на конкретный вид линейного оператора, соответствующего аппроксимации коэффициента отражения.

Сравним результаты синхронной инверсии на основе аппроксимаций и на основе точного решения на простых трёхслойных моделях. В качестве данных рассмотрим синтетические сейсмические трассы, полученные свёрткой трассы коэффициентов отражения, рассчитанных согласно системе уравнений Цёппритца, с импульсом Риккера с главной частотой 30~Гц. На рисунке~\ref{fig:low-contrast-1d} показаны результаты инверсии для слабоконтрастной трёхслойной среды. Алгоритмы инверсии на основе аппроксимаций и на основе точного решения позволили восстановить корректные значения упругих свойств среды. Этого удаётся достичь за счёт того, что при малых контрастах упругих свойств среды линейный оператор прямой задачи на основе аппроксимаций достаточно точно описывает исходный нелинейный оператор на основе точного решения уравнений Цёппритца.
\begin{figure}[ht]
     \includegraphics[width=1.0\linewidth]{figures/zoeppritz_shuey_low_contrast.pdf}
     \caption{Результат инверсии с аппроксимацией Шуе и с точным решением уравнений Цёппритца для слабоконтрастной среды}
     \label{fig:low-contrast-1d}
\end{figure}

Рассматривая среду со значительными контрастами упругих свойств, можем заметить, что инверсия на основе аппроксимаций уступает инверсии на основе точного решения, как это показано на рисунке~\ref{fig:high-contrast-1d}. Ошибки в прямом моделировании, возникающие в силу линеаризации исходной задачи, сказываются на результате решения обратной задачи. Однако синхронная инверсия на основе точного решения в случае контрастной среды позволяет получить достаточно точное решение, в то время как инверсия на основе аппроксимаций выходит за рамки применимости.
\begin{figure}[ht]
     \includegraphics[width=1.0\linewidth]{figures/zoeppritz_shuey_high_contrast.pdf}
     \caption{Результат инверсии с аппроксимацией Шуе и с точным решением уравнений Цёппритца для контрастной среды}
     \label{fig:high-contrast-1d}
\end{figure}

Таким образом, синхронная инверсия на основе точного решения уравнений Цёппритца способна корректно восстанавливать модель упругих свойств среды при наличии больших контрастов свойств, в то время как инверсия на основе аппроксимаций хорошо справляется со слабоконтрастными средами. Вычислительные затраты при инверсии на основе точного решения могут превышать таковые для инверсии на основе аппроксимаций в 2-3 раза, что является приемлемым и достигается за счёт использования вычислительно эффективных формул обращения матриц малого размера. Область применимости инверсии на основе точного решения шире, чем у инверсии на основе аппроксимаций, поскольку она способна корректно восстанавливать упругие свойства при наличии больших контрастов между слоями.


\section{Восстановление распределения упругих свойств среды}

Рассмотрим трассу реалистичной модели Marmousi, в которой имеется гладкий тренд упругих параметров, а также значительные контрасты свойств. Результат применения синхронной инверсии на основе точного решения уравнений Цёппритца в параметризации $(I_p, I_s, \rho)$ показан на рисунке~\ref{fig:marmousi-1d-impedance}. Импеданс продольной волны $I_p$ восстановлен с достаточной точностью, однако благодаря тому, что синхронная инверсия работает с данными до суммирования, в сейсмических трассах, соответствующих ненулевым углам падения волны, содержится информация об импедансе поперечной волны $I_s$, который также восстанавливается с хорошей точностью. Известно, что при синхронной инверсии в классической параметризации $(v_p, v_s, \rho)$ наиболее трудно определяемым параметром является плотность среды $\rho$. Для параметризации через импедансы и плотность этот эффект становится ещё более выраженным. В параметризации $(I_p, I_s, \rho)$ практически вся информация об отражающих свойствах среды заложена в импедансы $I_p$ и $I_s$, и о надёжном определении плотности среды говорить не приходится. 
\begin{figure}[ht]
     \includegraphics[width=1.0\linewidth]{figures/marmousi_1d_impedance.pdf}
     \caption{Результат инверсии с точным решением уравнений Цёппритца в параметризации через импедансы $(I_p, I_s, \rho)$ для трассы модели Marmousi}
     \label{fig:marmousi-1d-impedance}
\end{figure}

Синхронная инверсия на основе точного решения системы уравнений Цёппритца позволяет извлечь из сейсмических трасс информацию не только об импедансе продольной волны $I_p$, но и об импедансе поперечной волны $I_s$. Следовательно, при синхронной инверсии появляется возможность оценить такой параметр, как соотношение скоростей $v_p / v_s$. На рисунке~\ref{fig:marmousi-2d-impedance}а показан импеданс продольной волны, полученный в результате акустической инверсии, а на рисунке~\ref{fig:marmousi-2d-impedance}б~--- импеданс продольной волны, полученный в результате синхронной инверсии на основе точного решения уравнений Цёппритца. Отличия между импедансами $I_p$ могут быть обусловлены тем, что в программной реализации синхронной инверсии не использовалась нормализация данных, и сейсмические трассы, соответствующие нулевому углу падения волны, имеют меньшую амплитуду, чем остальные, и, следовательно, меньший вклад в общую невязку по всем углам. Тем не менее, при синхронной инверсии дополнительно может быть определён импеданс поперечной волны, как это показано на рисунке~\ref{fig:marmousi-2d-impedance}в, а используя импедансы $I_p$ и $I_s$ появляется возможность оценить соотношение скоростей $v_p / v_s$, изображённое на рисунке~\ref{fig:marmousi-2d-impedance}г. При анализе соотношения скоростей выделяется аномалия, расположенная около $x=3$~км по латерали и около $t=1500$~мс по времени. Таким образом, синхронная инверсия позволяет получать больше информации о геологической среде, используя данные до суммирования. Применение точного решения уравнений Цёппритца при синхронной инверсии даёт возможность успешного её применения для сложных и контрастных сред.
\begin{figure}[ht]
     \includegraphics[width=1.0\linewidth]{figures/marmousi_2d_impedance.pdf}
     \caption{а) Импеданс $I_p$, восстановленный акустической инверсией, б) импеданс $I_p$, восстановленный синхронной инверсией, в) импеданс $I_s$, восстановленный синхронной инверсией, г) отношение скоростей $v_p / v_s$, синхронная инверсия}
     \label{fig:marmousi-2d-impedance}
\end{figure}

Рассмотрим пример применения синхронной инверсии на основе точного решения уравнений Цёппритца в классической параметризации через скорости и плотность к реалистичной модели Marmousi. На рисунке~\ref{fig:marmousi-2d-velocity} показаны начальная, истинная и восстановленная скоростные модели, а также одна трасса при $x=3$~км. В результате синхронной инверсии получается детальная модель скорости продольной волны $v_p$, на которой хорошо прослеживаются высокие контрасты свойств и границы слоёв истинной модели. Использование акустической инверсии в данном примере позволило бы определить детальную модель акустического импеданса, однако синхронная инверсия позволяет определить скорость продольной волны $v_p$ отдельно от плотности. Возможность реконструкции скоростной модели среды позволяет использовать синхронную инверсию в качестве инструмента уточнения скоростной модели с целью дальнейшего применения процедур миграции.
\begin{figure}[ht]
     \includegraphics[width=1.0\linewidth]{figures/marmousi_2d_velocity.pdf}
     \caption{а) Начальная скоростная модель, б) трасса при $x=3$~км, в) истинная скоростная модель, г) восстановленная синхронной инверсией на основе точного решения уравнений Цёппритца скоростная модель}
     \label{fig:marmousi-2d-velocity}
\end{figure}


\section{Заключение}
Предложенный алгоритм синхронной инверсии на основе точного решения системы уравнений Цёппритца позволяет расширить область применимости синхронной инверсии. Классические методы синхронной инверсии на основе аппроксимаций точного решения уравнений Цёппритца способны корректно восстанавливать упругие свойства в случае слабоконтрастных сред, в то время как предложенный алгоритм синхронной инверсии на основе точного решения применим к средам с большим контрастом упругих свойств.

Записывая систему уравнений Цёппритца в параметризации через импедансы, использование данных до суммирования в синхронной инверсии позволяет извлечь из сейсмических трасс не только информацию об импедансе продольной волны $I_p$, но и об импедансе поперечной волны $I_s$. Таким образом, при синхронной инверсии имеется возможность оценить отношение скоростей $v_p/v_s$, которое является важным параметром для дальнейшей интерпретации.

Сейсмическая инверсия может использоваться как инструмент для уточнения скоростной модели среды при совместном применении с глубинной миграцией. В отличие от акустической инверсии, которая работает с импедансом продольной волны $I_p$, синхронная инверсия позволяет отделить скорость продольной волны $v_p$ от плотности $\rho$. Реконструкция скоростной модели среды вместо модели импедансов даёт возможность применения синхронной инверсии в качестве инструмента уточнения скоростной модели с целью дальнейшего применения процедур миграции.


\bigskip
\begin{thebibliography}{1} % Ссылки на литературу приводятся только на английском языке (даже если у источника нет англоязычной версии).

\bibitem{Russell_1991} B.H. Russell, D.P. Hampson,
\href{10.1190/1.1888870}{\it Comparison of poststack seismic inversion methods}, SEG Annual Meeting, Expanded Abstracts, Houston, 1991, 876--878.

\bibitem{Russell_1988} B.H. Russell,
\href{https://doi.org/10.1190/1.9781560802303}{\it Introduction to Seismic Inversion Methods}, Course Notes Series, Society of
Exploration Geophysicists, S. N. Domenico, 1988, 80--101.

\bibitem{Ampilov_2009} Yu.P. Ampilov, A.Yu. Barkov, I.V. Yakovlev, K.E. Filippova, I.I. Priezzhev, 
{\it  Almost everything of seismic inversion. Part 1}, Seismic technologies, 2009, {\bf 4}, 3--16.

\bibitem{Bortfeld_1961} R. Bortfeld,
\href{https://doi.org/10.1111/j.1365-2478.1961.tb01670.x}{\it Approximation to the reflection and transmission coefficients of plane longitudinal and transverse waves}, Geophysical Prospecting, 1961, {\bf 9}, 485--502.

\bibitem{Romanenko_2010} M.Yu. Romanenko, I.N. Kerusov, D.E. Miroshnichenko, Yu.V. Masalkin,
{\it Estimation of the simultaneous seismic data inversion method effectiveness in relation to the low-contrast reservoirs models}, Seismic technologies, 2010, {\bf 2}, 55--61.

\bibitem{Silvestrov_2011} I. Silvestrov, V. Tcheverda,
\href{https://doi.org/10.1088/1742-6596/290/1/012014}{\it SVD analysis in application to full waveform inversion of multicomponent seismic data}, Journal of Physics: Conference Series, 2011, {\bf 290} (1), 12--14.

\bibitem{Shuey_1985} R. Shuey,
\href{https://doi.org/10.1190/1.1441936}{\it A simplification of the zoeppritz equations}, Geophysics, 1985, {\bf 50} (4), 609--614.

\bibitem{Loris_2007} I. Loris, G. Nolet, I. Daubechies, F.A. Dahlen,
\href{https://doi.org/10.48550/arXiv.physics/0608094}{\it Tomographic inversion using $l_1$-norm regularization of wavelet coefficients}, Geophysical Journal International, 2007, {\bf 170} (1), 359--370.

\bibitem{Aggrawal_2021} H.O. Aggrawal, J. Modersitzki,
\href{https://doi.org/10.48550/arXiv.2103.10010}{\it Hessian initialization strategies for L-BFGS solving non-linear inverse problems}, SSVM, 2021, 216--228.

\bibitem{Zhang_2011} R. Zhang, J. Castagna,
\href{https://doi.org/10.1190/geo2011-0103.1}{\it Seismic sparse-layer reflectivity inversion using basis pursuit decomposition}, Geophysics, 2011, {\bf 76}, R147--R158.

\bibitem{Zhang_2014} F. Zhang, R. Dai, H. Liu,
\href{https://doi.org/10.1016/j.jappgeo.2014.07.024}{\it Seismic inversion based on L1-norm misfit function and total variation regularization}, Journal of Applied Geophysics, 2014, {\bf 109}, 111--118.

\bibitem{Belge_2002} M. Belge, M.E. Kilmer, E.L. Miller,
\href{https://doi.org/10.1088/0266-5611/18/4/314}{\it Efficient determination of multiple regularization parameters in a generalized L-curve framework}, Inverse Problems, 2002, {\bf 18} (4), 1161--1183.

\bibitem{Tarantola_1986} A. Tarantola,
\href{https://doi.org/10.1190/1.1442046}{\it A strategy for nonlinear elastic inversion of seismic reflection data}, Geophysics, 1986, {\bf 51} (10), 1893--1903.

\bibitem{Misra_2008} S. Misra, M.D. Sacchi,
\href{https://doi.org/10.1190/1.2958008}{\it Global optimization with model-space preconditioning: Application to AVO inversion}, Geophysics, 2008, {\bf 73} (5), R71--R82.

\bibitem{Tran_2012} K.T. Tran, D.R. Hiltunen,
\href{https://doi.org/10.2113/JEEG17.4.197}{\it One-dimensional inversion of full waveforms using a genetic algorithm}, Journal of Environmental and Engineering Geophysics, 2012, {\bf 17} (4), 197--213.

\end{thebibliography}
\end{document}
