Разработка помехоустойчивых алгоритмов динамической инверсии сейсмических данных тема диссертации и автореферата по ВАК РФ 25.00.10, кандидат наук Ли Цян

• Ли Цян
• кандидат наук
• 2018

Оглавление диссертации кандидат наук Ли Цян

СОДЕРЖАНИЕ

ВВЕДЕНИЕ

ГЛАВА 1. ОБЗОР МЕТОДОВ СЕЙСМИЧЕСКОЙ АМПЛИТУДНОЙ ИНВЕРСИИ

1.1 Акустическая инверсия

1.1.1 Рекурсивная инверсия

1.1.2 Цветная инверсия

1.1.3 Инверсия редких импульсов

1.1.4 Инверсия, основанная на модели

1.2 Упругая инверсия

1.2.1 Расширенная упругая инверсия

1.2.2 Синхронная инверсия

1.3 Стохастическая инверсия

1.4 Методы регуляризации и оптимизации для решения обратной задачи

1.4.1 Метод регуляризации

1.4.1.1 Метод на основе сингулярного разложения

1.4.1.2 Метод Тихонова

1.4.1.3 Итерационный метод

1.4.1.4 Байесовский подход

1.4.2 Метод выбора параметра регуляризации

1.4.3 Метод оптимизации

1.5 Метод оценивания импульса

1.5.1 Статистический метод

1.5.2 Детерминистический метод

1.5.2.1 Метод оценивания импульса с постоянным фазовым спектром

1.5.2.2 Полноволновый метод

1.5.2.3 Алгоритм Рой-Уайт

1.5.2.4 Байесовский метод извлечения импульсов

1.5.3 Анализ влияния импульса на результат сейсмической амплитудной

инверсии

1.6 Анализ устойчивости амплитудной сейсмической инверсии на

ПРИМЕРЕ СТАТИСТИЧЕСКИХ МОДЕЛЬНЫХ ДАННЫХ

1.6.1 Создание набора статистических моделей

1.7 Выводы

ГЛАВА 2. ОПТИМИЗАЦИОННЫЙ АЛГОРИТМ ОЦЕНИВАНИЯ ИМПУЛЬСА

2.1 Математические основы оптимизационного алгоритма оценивания импульса

2.2 Блок-схема алгоритма

2.3 ПРИМЕНЕНИЕ АЛГОРИТМА НА МОДЕЛЬНЫХ ДАННЫХ

2.3.1 Применение на наборе статистических моделей

2.3.2 Применение на синтетических данных для оценки АУО/АУА импульсов

2.4 Выводы

ГЛАВА 3. РАЗРАБОТКА АЛГОРИТМОВ ПОМЕХОУСТОЙЧИВЫХ МЕТОДОВ АКУСТИЧЕСКОЙ И СИНХРОННОЙ ИНВЕРСИИ СЕЙСМИЧЕСКИХ ДАННЫХ

3.1 Анализ применимости различных методов регуляризации в задачах

ДИНАМИЧЕСКОЙ ИНВЕРСИИ СЕЙСМИЧЕСКИХ ДАННЫХ

3.1.1 Алгоритм сейсмической инверсии, основанный на модели

3.1.2 Метод регуляризации

3.1.3 Создание набора статистических моделей

3.1.4 Результат тестирования стабилизаторов

3.2 Автоматический подбор оптимального параметра регуляризации

3.2.1 Алгоритм подбора параметра регуляризации

3.2.2 Постановка задачи

3.2.3 Результат вычислительных экспериментов

3.3 Синхронная инверсия на основе разработанных алгоритмов

3.3.1 Математические основы алгоритма синхронной инверсии с регуляризации минимальным носителем градиента

3.3.2 Применение алгоритма на модельных данных

3.5 Выводы

ГЛАВА 4. ПРИМЕНЕНИЕ РАЗРАБОТАННЫХ АЛГОРИТМОВ НА РЕАЛЬНЫХ ДАННЫХ

4.1 Эталон 1. Карбонаты осинского горизонта нижнего кембрия. Восточная Сибирь

4.1.1 Оценивание импульсов и увязка ГИС и сейсмики

4.1.2 Применение метода акустической детерминистической инверсии с регуляризацией MGS

4.1.2.1 Построение низкочастотной модели

4.1.2.2 Результаты акустической инверсии

4.2 Эталон 2. Тонкий пласт песчаника ярактинского горизонта вендского возраста

4.2.1 Петроупругое обоснование методик исследования

4.2.2 Применение метода акустической инверсии с регуляризацией MGS

4.2.2.1 Оценивание импульсов и привязка скважин к сейсмическим данным

4.2.2.2 Результат акустической инверсии

4.2.3 Применение синхронной инверсии с регуляризацией MGS

4.2.3.1 Подготовка сейсмических данных

4.2.3.2 Извлечение импульсов

4.2.3.3 Результаты синхронной инверсии

4.2.4 Интерпретация результатов инверсии сейсмических данных

4.3 Выводы

ЗАКЛЮЧЕНИЕ

СПИСОК ЛИТЕРАТУРЫ

140

ВВЕДЕНИЕ

Актуальность темы

Количественная интерпретация сейсмических данных на данный момент является самой современной технологией для поиска и разведки залежей углеводородов в неструктурных ловушках. В основу данного подхода положен комплексный анализ разномасштабных данных (рисунок 1.1): литологические исследования, интерпретация данных ГИС, петроупруое моделирование, сейсмическая инверсия и другие. При этом, сейсмическая амплитудная инверсия является ключевой технологией в комплексе количественной интерпретации [18, 19, 20]. В настоящее время сейсморазведка применяется не только для изучения структурного строения нефтяных и газовых месторождений, но и для прогноза их геологических свойств (литология, пористость и насыщение углеводородами) [3, 7, 21]. Этот успех связан с технологией сейсмической амплитудной инверсии -совокупности методов оценки упругих свойств геологической среды по сейсмическим данным на основе решения обратной динамической задачи [10]. Данная работа посвящена разработке инструментов сейсмической амплитудной инверсии для улучшения качества прогноза залежей углеводородов.

Сейсмическая инверсия преобразует параметры отражений от границ в свойства пласта и позволяет получить упругие параметры - акустический импеданс Zp, сдвиговый импеданс Zs, отношение скоростей Vp/Vs и плотность р [2, 36]. Эти упругие параметры могут пересчитываться в геологические параметры (пористость, эффективные толщины коллекторов, литотип, тип флюида и другие) на основе петроупругого моделирования (rock physics).

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

Рисунок 1.1. Количественная интерпретация сейсмических данных. По материалам кафедры разведочной геофизики РГУ нефти и газа им. И.М.Губкина

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

2. Надёжность результата сейсмической инверсии в различных геологических условиях зависит от устойчивости алгоритма к параметрам слоистой среды (толщин пластов, контрастности среды) и уровня шума сейсмических данных. Исследование этого вопроса позволит повысить эффективность метода в различных геологических условиях.

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

4. Выбор оптимального параметра регуляризации оказывает существенное влияние на точность получаемых регуляризированных решений. Необходимо разработать алгоритм автоматического подбора параметра регуляризации для задачи динамической инверсии сейсмических данных независящий от субъективного фактора.

Введение диссертации (часть автореферата) на тему «Разработка помехоустойчивых алгоритмов динамической инверсии сейсмических данных»

Цель работы

Разработка инструментов сейсмической амплитудной инверсии для улучшения

качества прогноза залежей углеводородов.

Задачи исследований

1. Разработка оптимизационного алгоритма оценивания сейсмического импульса.

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

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

4. Развитие метода синхронной инверсии на основе разработанных автором алгоритмов.

5. Апробация предложенных инструментов в реальных условиях.

Защищаемые положения

1. Разработанный автором оптимизационный алгоритм извлечения импульса позволяет уменьшить фазовую ошибку.

2. На основе исследования различных методов регуляризации в задачах динамической инверсии сейсмических данных выбран оптимальный стабилизирующий функционал.

3. Разработанный автором алгоритм позволяет автоматизировать подбор оптимального параметра регуляризации для сейсмической инверсии.

4. Авторские алгоритмы, реализованные в методах акустической и синхронной инверсий, позволили получить более точную оценку упругих свойств среды.

Научная новизна

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

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

3. Автором предложен алгоритм автоматического подбора оптимального параметра регуляризации для задач динамической сейсмической инверсии.

4. В результате выполненных автором исследований усовершенствованы методы акустической и синхронной инверсий сейсмических данных.

Личный вклад

Все представленные в работе идеи и выводы получены и программно реализованы автором лично на кафедре разведочной геофизики и компьютерных систем РГУ нефти и газа (НИУ) имени И.М. Губкина. Практическая значимость работы

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

Основные тезисы и результаты диссертационной работы были представлены на следующих международных конференциях: 87-я Ежегодная конференция и выставка общества геофизиков (SEG), г. Хьюстон, США, 2017; 70-ая международная молодежная научная конференция «Нефть и газ - 2016», г. Москва, 18-20 апреля 2016г.; 71-ая международная молодежная научная конференция «Нефть и газ - 2017», г. Москва, 23-26 апреля 2017г.; XII Всероссийская конференция молодых ученых, специалистов и студентов «Новые технологии в газовой промышленности» (газ, нефть, энергетика). г. Москва, 24-27 октября 2017г.

По теме работы опубликованы 2 статьи в журналах из перечня ведущих рецензируемых изданий, рекомендованных ВАК. Объем и структура диссертации

Диссертация состоит из 4-х глав, введения и заключения, содержит 147 страниц, в том числе 97 рисунка и 5 таблиц. Список литературы включает 93 наименований. Благодарности

Автор глубоко признателен заведующему кафедрой, профессору В.И. Рыжкову, который является научным руководителем и координатором работы, а также ассистенту Данько Д.А. и ведущему инженеру Филимоненко С.В. за содействие в разработке алгоритмов. Хочу выразить искреннюю благодарность и признательность учителям и коллегам, которые помогали мне во время обучения на кафедре разведочной геофизики и компьютерных систем РГУ нефти и газа (НИР) имени И.М. Губкина. Среди них проф. Воскресенский Ю.Н., доц. Шубин А.В. доц. Белоусов А.В., доц. Варов Е.Б., асс. Варов Ю.Е., доц. Жуков А.М, доц. Карапетов Г. А., асп. Гаркин А.С., вед.инж. Локшина Л.Б., вед.инж. Зяйнетдинова Р.М., проф. Ермолкин В.И.

ГЛАВА 1. ОБЗОР МЕТОДОВ СЕЙСМИЧЕСКОЙ АМПЛИТУДНОЙ ИНВЕРСИИ

Сейсмическая инверсия представляет собой совокупность методов оценки упругих свойств слоистой среды по сейсмическим данным на основе решения обратной динамической задачи с учетом тонкослоистого характера геологической среды [10].

Основная цель инверсии сейсмических данных - восстановление упругих свойств горных пород, которые связаны с геологическими свойствам (литология, пористость, насыщение углеводородами) или с физическими условиями (давление, температура) при которых они находятся [34, 35]. Сейсмическая инверсия успешно применяется для картирования коллекторов. Различные модификации инверсии позволяют восстанавливать акустический импеданс, сдвиговый импеданс и плотность и прогнозировать литологию и емкостные свойства пород [12, 13, 30].

Сейсмическую инверсию можно разделить на детерминистическую и стохастическую по методу оптимизации, так же ее можно подразделить на инверсию после суммирования и до суммирования сейсмических данных. Кроме того, ее делят на акустическую инверсию, служащую для оценки атрибутов - импедансов слоистой среды, получаемых по суммированным сейсмическим данным и упругую инверсию - для оценки атрибутов в виде упругих свойств (продольных, поперечных импедансов и плотностей), получаемых по сейсмическим трассам для различных удалений от источника

[9].

В первых 4-х разделах главы проведен обзор различных видов сейсмической инверсии (1.1-1.4). Далее рассмотрен метод регуляризации и оптимизации для решения обратной проблемы (1.5), метод извлечения импульса для сейсмической инверсии (1.6). В конце приведен анализ устойчивости амплитудной сейсмической инверсии на примере статистических модельных данных (1.7).

1.1 Акустическая инверсия

Суть амплитудной инверсии сейсмических данных - оценка акустического импеданса горных пород. Акустическая инверсия для получения импеданса впервые была выполнена Lindseth в 1979 г. [63]. В последующие годы были разработаны различные подходы для инверсии суммарных сейсмических данных, в том числе инверсия с ограниченной полосой пропускания, инверсия редких импульсов, инверсия, основанная на модели и др.

Отражательная способность среды может быть извлечена из сейсмических данных на основе сверточной модели сейсмической трассы:

5 = Ш*Я + Ы, (1.1)

где 5 - сейсмическая трасса однократно отраженных волн, Ж - сейсмический импульс, Я - коэффициент отражения, N - шум и символ * обозначает свертку. Шум предполагается случайным и некоррелированным.

Коэффициент отражения Я определяется как контраст импедансов между двумя соседними пластами:

^ ' (1.2) где 2 - акустический импеданс в пласте, являющийся произведением скорости и плотности.

Исходными данными для акустической инверсии служат сейсмические кубы или профили после процедуры миграции и скважинные данные (акустический каротаж, ВСП).

1.1.1 Рекурсивная инверсия

Рекурсивная инверсия или псевдоакустический каротаж (ПАК) - это наиболее ранний и простой метод акустической инверсии [63]. Преобразования сейсмических записей по методике ПАК сводятся к трем основным операциям:

1) преобразованию сейсмических записей в трассы, близкие к импульсным трассам (псевдоимпульсные трассы). Это можно выполнить путем деконволюции;

2) масштабированию амплитуд псевдоимпульсных трасс к значениям коэффициентов отражения;

3) переходу от последовательности коэффициентов отражения к распределению импедансов в слоистой среде по рекуррентной формуле

/п = /сП?=11^, (1.3)

где /0 - акустический импеданс верхнего слоя; rt - коэффициент отражения от подошвы -ого слоя; п - номер слоя в слоистой среде.

Из формулы (1.3) только можно рассчитать последовательность относительных значений импеднсов. Возможно восстановление абсолютных значений акустических импеднсов с добавлением низкочастотной фоновой модели. Основным недостатком метода рекурсивной инверсии является возможность накопления ошибок со временем.

1.1.2 Цветная инверсия

Цветная инверсия - хороший ускоренный метод инверсии сейсмических данных. Процесс инверсии может быть аппроксимирован фильтрацией с простым оператором [61]. Фаза такого оператора должна быть постоянной -90o, что соответствует представлению о трансформации нуль-фазового отражения в скачок акустического импеданса. Кроме того, Walden and Hosken [86] показали, что спектр акустического импеданса, полученный из данных ГИС, может быть представлен трендом , где Д - константа, f - частота. Это означает, что для выполнения инверсии можно использовать один сверточный оператор. Этот подход предполагает нуль-фазовый импульс, поэтому не нужно оценивать амплитудный спектр импульса. Амплитудный спектр оператора инверсии получается объединением амплитудного сейсмического спектра и спектра акустического импеданса в скважинах. После того, как был получен оператор цветной инверсии, его можно просто применить к сейсмическим данным при интерпретации как «определенный фильтр».

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

1.1.3 Инверсия редких импульсов

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

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

Существуют несколько подходов к инверсии редких импульсов, среди которых основными являются: 1) инверсия максимального правдоподобия [74]; 2) инверсия с ограничением нормы-£1 [65]; 3) инверсия редких импульсов с про странственными ограничениями.

При инверсии максимального правдоподобия, итоговая целевая функция для инверсии редких импульсов может быть записана следующим образом:

/(Г) = Й^^Т + - 2ш • 1п(Л) - 2(1 - т) 1п(1 - X), (1.4)

где

г(Ь) - последовательность коэффициентов отражения; т - количество коэффициентов отражения; Ь - общее количество временных отсчетов; Я - стандартное отклонение коэффициента отражения; а - стандартное отклонение шума; п(1) - отсчеты шума;

X - правдоподобие того, что данный отсчет соответствует отражению.

При инверсии редких импульсов с ограничением (норма-£1), итоговая целевая функция может быть записана следующим образом:

где \\Wr — dobs\\ — невязка между сейсмическими (dobs) и синтетическими (Wr) данными; W —матрица импульса; г—вектор коэффициента отражения; X^—i^l - ограничение минимума суммы модулей искомых коэффициентов отражения; lrtl — отсчеты модулей коэффициента отражения.

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

Программное обеспечение Jason (компания CCG) предоставляет современный алгоритм инверсии редких импульсов с пространственными ограничениями (Constrained Sparse-Spike Inversion - CSSI). При необходимости к целевой функции добавляется ограничение и тогда большие группы трасс инвертируются одновременно.

1.1.4 Инверсия, основанная на модели

Идея инверсии, основанной на модели, состоит в подборе такой акустической модели слоистой среды (акустические импедансы), чтобы соответствующие ей синтетические трассы совпали с реальными сейсмическими трассами. Решение начинается с задания начальной модели, которая затем постепенно корректируется, пока не будет достигнуто совпадение синтетических и реальных сейсмических данных. Считается, что достижение совпадения свидетельствует об адекватности подобранной модели строению реальной среды.

В алгоритме инверсии, основанной на модели, используется низкоконтрастное приближение [56, 75]:

J(r)= wwr — d^w+xl^lni,

(1.5)

Ri^-Mn(Zi) = -[ln(Zi+i) — ln(Zi)]l

i

i

(1.6)

2

2

где I представляет собой границу между слоями I и I + 1. Если мы рассмотрим N коэффициентов отражения, уравнение (1.6) можно записать в матричной форме как:

№1

_ 1

= 2

Яы.

-1 1 0 -0 -1 1 ••. 0 0 -1 1

¿1 ¿2

(1.7)

где = 1п^).

Далее, если мы представим сейсмические трассы как свертку сейсмического импульса с коэффициентами отражения для ^слоистой модели, мы можем записать результат в матричном виде как:

0 0 -

Г511

=

w1

w2 0

w1

*2

Я

N J

(1.8)

- отсчеты извлеченного

5 = -ШОЬ,

2

где - отсчеты сейсмической трассы, сейсмического импульса, ^ - последовательность коэффициентов отражения.

Объединение уравнений (1.6) и (1.8) дает нам окончательное выражение, которое связывает сейсмические трассы с логарифмом продольного импеданса:

(1.9)

где Ж - матрица сейсмического импульса, определенная в уравнении (1.8), О -производная матрица в уравнении (1.7). Если уравнение (1.9) инвертируется с использованием стандартной методики обратной матрицы, существуют две проблемы. Во-первых, вычисление обратных матриц ресурсоемкая и потенциально нестабильная операция. Что еще более важно, матрица инверсии не будет восстанавливать низкочастотную составляющую импеданса. Альтернативной стратегией является построение начальной приближенной модели импеданса, а затем итеративный поиск решения с использованием метода сопряженных градиентов. Инверсия на основе модели использует линейный алгоритм инверсии, который предполагает, что сейсмические трассы (5) и сейсмический импульс (Ж) известны. Итерации, изменяющие первоначальную

модель, производятся до тех пор, пока полученные синтетические трассы не будут соответствовать реальным сейсмическим трассам. Начальная (фоновая) модель строится с помощью интерполяции сглаженных данных ГИС вдоль прослеженных горизонтов.

1.2 Упругая инверсия

P. Connolly (1999) [44] ввел понятие упругого импеданса, как обобщение акустического импеданса для ненулевых углов падения. Упругий импеданс является функцией скорости продольной Vp, поперечной волны Vs, плотности р и может быть получен из формулы Аки-Ричардса [37], которая является линеаризацией системы уравнений Цеппритца [93]:

т =1№ + Ï£) + +-^sin2etan2e, (1.10)

W 2\Vp р ) \2Vp V* Vs V* Р ) 2 Vp > \ ;

_ 1

где AV = V(ti) - V(ti-1), V = - [V(ti) + V(ti-1)], Ар = p(tt) — p(ti-1),

n (a+ß) n

в = , a — угол падения, ß — угол отражения.

Упругий импеданс El является функцией, которая обладает свойствами, аналогичными акустическому импедансу, зависит от угла падения, скоростей и плотности. Запишем через EI(t) формулу для коэффициентов отражения (1.2) для переменного угла падения:

(MD

Подставляем (1.10) в (1.11). После преобразований получим выражение для упругого импеданса:

EJ — y(1+tan2e)y(-8Ksin2e)p(1-4Ksin26) (1.12)

vs

где К = 77*2.

V2 VP

Из уравнения (1.11) следует, что упругий импеданс Е1 является фактическим обобщением понятия акустического импеданса А1 для случая наклонного

падения волны на границу, так как при 0 = 0 он равен акустическому импедансу Al.

Размерность упругого импеданса из выражения (1.12) зависит от угла. Это вносит проблемы при практическом применении. Whitcombe [88] предложил нормировать выражение (1.12) для EI:

,v \1+tan2d/ ^-8Ksin26 .!-4Ksin2e

= fé) Ш , а»)

где VPo> VSo> p0 — константы, которые получают усреднением соответствующих кривых ГИС.

Основные этапы упругой инверсии сейсмических данных сводятся к следующему:

1) Проведению частичного углового суммирования.

2) Определению индивидуального импульса для каждой угловой суммы.

3) Выполнению инверсии согласно любому алгоритму акустической инверсии.

1.2.1 Расширенная упругая инверсия

Двухчленная аппроксимация для коэффициентов отражения из выражения (1.10):

+ + № — sin26. (1.14)

W 2\Vp р ) \2Vp V2 Vs V2 р ) V J

Заменив sin2в на tgx, а затем умножив на cosx [89] получим:

R(6) = А • cosx + В • sinx, (1.15)

. l(bVv , Др\ (ДУр Av2 AVS „Vs2Ap\

где А = -(^ + -г), В = (^г- — 4——— 2——), отсюда следует уравнение 2 \Vp р / \2Vp Vp Vs Vp P /

нормированного упругого импеданса (1.13):

/v \c0sx+sinx ( \ -QKsin^ cosx-4Ksinx

™M = VPoPO(%) © (f) , (1.16)

EEI(x) — это расширенный упругий импеданс.

Ниже приведены характерные особенности EEI инверсии:

1) Коэффициент отражения на ненормальном падении определяется по формуле (1.2), и имеет ту же форму, что и при нормальном падении.

2) Алгоритмы пост-стек инверсии могут быть применены к частичной угловой сумме для получения ЕЕ1.

3) Процедура аналогична акустической инверсии. Импульс определяется индивидуально для различных угловых сумм. ЕЕ1 является основой для привязки сейсмических данных, а также их инверсии на дальних удалениях. Стоимость ЕЕ1 инверсии очень низкая.

4) Коэффициент Пуассона и параметры Ламе могут быть получены из модели ЕЕ1. Кроссплот различных полученных параметров может быть использован в динамической интерпретации для определения литологии и флюидонасыщения [4, 38, 89].

5) Этот подход достаточно хорошо оценивает продольный и поперечный импедансы.

1.2.2 Синхронная инверсия

Для изучения литологического состава и УВ насыщения горных пород важно иметь оценки отношений скоростей продольных и поперечных волн. Дополнительная информация о поперечных волнах может быть извлечена из материалов, полученных при ненулевых углах падения волны на границу. Оценивая изменение амплитуды отраженных продольных волн в зависимости от угла их падения на границу, можно получить информацию о свойствах, образовавшихся на границе поперечных волн. Получить информацию о поперечных волнах можно только из обработанных сейсмограмм, которые несут информацию об углах падения волны на границу [43, 52, 55, 71, 72, 73]. Синхронная инверсия позволяет восстановить модель продольных скоростей Ур, поперечных скоростей У3 и плотностей р одновременно.

Применяется линейная аппроксимация системы уравнений Цеппритца в записи Фатти [37, 48, 93]:

Ярр(в) = (1 + 1ап2в)ЯР + (-8у25т2в)Яп + (4у2^т2в - 1ап2в)Яв, (1.17)

^ 2 Ур р^ ъ 2\У5 р' 2 р у '

где ур , У8 , и ~р - средняя продольная скорость, поперечная скорость и плотность над отражающей границы, АУР, АУБ и Ар - соответствующие контрасты, и в - угол падения. у - отношение продольной скорости к поперечной скорости.

При низкоконтрастнном приближении:

^ = + (120)

+ (121)

где 2Р — продольный импеданс; 15 — поперечный импеданс; ЬР = 1п(ЕР) , ЬБ = 1п(75) и Ьр = 1п(р) — соответствующие натуральные логарифмы. t — время. Для У-слоистой модели, мы можем записать результат в матричном виде как:

11 1

Яр = , = 10Ь5, и Яр = \Ыр , (1.22)

где Б - производная матрица:

—1 1 0

Б =

0 —1 1 \ 0 0 —1 1

(1.23)

Далее, если мы представим сейсмические трассы как свертку сейсмического импульса (зависит от угла падения) с коэффициентами отражения, то уравнение (1.17) можно записать в матричной форме как:

Т(в) = 0.5с1ш(6)0ЬР + 0.5с2ш(6)0Ь5 + ш(в)с3ОЬр , (1.24)

где с± = 1 + 1ап2в, с2 = —8у21ап2в, с3 = —0.51ап26 + 2у2Бт2в, У = ^ .

Существует линейная связь между логарифмами продольного импеданса 2р, поперечного импеданса 23 и плотностью р [42, 51,56].

ln(Zs) = kln(Zp) + kc+АLs, (1.25)

1п(1р) = тЫ(1Р) + тс + АLp . (1.26)

Объединяя уравнения (1.24), (1.25) и (1.26), получим:

Т(в) = с^ш(в)тР + с2^(в)ЭМ5 + ш(в)с3ОМр , (1.27)

где с1 = 0.5^ + 0.5кс2 + тс3, с2 = 0.5с2.

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

-Т(в1у

Т(в2) =

-Т(вм).

c^(e1)W(e1)D с2(61)W(в1) D с3(61)W(в1) D ct( e2)W(e1)D С2( e2)W(e1)D сз( e2)W(e1)D

£i(eM)W(вм) D c2(eM)W(вм) D сз(вм)W(eM) D

Lp

ALS

ALp

(1.28)

Это уравнение решается итерационно методом сопряженных градиентов при использовании начального приближения.

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

1.3 Стохастическая инверсия

Результат работы процедуры детерминистической инверсии относительно гладкий (или блочный). Стохастическая инверсия сейсмических данных выдает набор различных разнородных реализаций атрибутов сейсмических данных. Стохастические методы могут давать существенное увеличение разрешающей способности, хорошо фиксируя детали разреза вне сейсмической полосы частот [49, 50]. Любой из методов инверсии может быть реализован в стохастическом варианте для получения значений акустического и упругого импеданса.

Существует множество различных алгоритмов для создания геостатистических реализаций, удовлетворяющих критериям измеренных выборок, гистограммам, и пространственной корреляции. Самый легкий алгоритм геостатистического моделирования - это последовательное гауссово моделирование (Sequential Gaussian Simulation - SGS). Метод последовательного

гауссового моделирования комбинирует моделирование и метод кригинга (Кп§т§). ЭиЬги1е О. [46] предложил алгоритм стохастической инверсии (геостатистической инверсии), используя последовательное гауссово моделирование. Процедура выполнения стохастической инверсии представляет следующее:

1) Определяется случайный путь через все узлы (х, у), в которых производится стохастическое моделирование.

2) Для каждого узла генерируется множество реализаций локальных трасс акустического импеданса (акустическая инверсия) или продольного и поперечного импедансов и плотности (синхронная инверсия);

3) Производится свертка последовательности коэффициентов отражения с импульсом;

4) Сравниваются результаты моделирования с наблюденной сейсмической трассой;

5) Выбирается лучшая модель и заносится в глобальную реализацию.

6) Производится переход к следующему узлу. И повторяются шаги 2 - 5, до заполнения всех узлов.

После заполнения всего пространства эти трассы рассматриваются как одна глобальная реализация. Производится много таких реализаций. Осреднение всех возможных реализаций стохастической инверсии является детерминистической инверсией.

Сейсмические и скважинные данные представляются в качестве функции плотности вероятности, которая дает геостатистическое описание, основанное на гистограммах и вариограммах [16]. В геостатистической инверсии предпочитают использовать вариограммы, оценивающие пространственную корреляцию данных, получамые согласно выражению:

А _

У^^^^-г+к)2, (1.29)

где 21 и 21+И - пары величин, отстоящих на расстояние И, а N - число значений вариограммы, попадающих в малые интервалы ДИ.

Существуют горизонтальные и вертикальные вариограммы. Моделью вертикальной вариограммы, которая подбирается с использованием данных ГИС, определяются высокочастотные вариации. Горизонтальная вариограмма контролирует вариации импеданса (или других упругих свойств среды) внутри пластов и определяется из сейсмических данных. Стохастическая инверсия может быть выполнена для любой заданной частоты, так же как данные ГИС. Это замечание не подразумевает что стохастическая инверсия, имеет более высокую разрешающую способность, так как разрешающая способность определяется частотным составом и полосой частот сейсмических данных. Наконец заметим, что метод стохастической сейсмической инверсии может дать вероятность реализации. Эта полезная информация помогает количественно оценить неопределенность результата инверсии с целью определения возможностей резервуара в рамках имеющихся ограничений.

1.4 Методы регуляризации и оптимизации для решения обратной задачи

Традиционный подход к анализу геофизических полей заключается в построении геологических моделей и сравнении теоретических геофизических данных, вычисленных для этих моделей, с наблюдёнными данными. Численное моделирование геофизических полей для параметров заданной модели обычно называют прямой задачей. Обычно мы описываем реальную геологию относительно простой моделью и пытаемся определить параметры модели по геофизическим данным. Такую задачу мы называем обратной задачей [1, 14]. Сейсмическая инверсия является обратной задачей.

При решении любой обратной задачи встают три важных вопроса:

1) Существует ли решение?

2) Является ли оно единственным?

3) Устойчиво ли оно?

Если все три вопроса, поставленные выше, имеют положительный ответ, эта математическая задача сформулирована корректно. В соответствии с терминологией Адамара [54] задача является некорректно поставленной, если

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

К сожалению, оказалось, что задачи сейсмической инверсии являются некорректно поставленными. В 60-е годы прошлого века российский математик Андрей Николаевич Тихонов [31, 32, 33] разработал основы теории решения некорректно поставленных задач. Он ввёл в решение обратной задачи метод регуляризации, который был основан на приближении некорректно поставленной задачи некоторой последовательностью корректно поставленных задач. В данном разделе главы рассмотрены методы регуляризации и оптимизации решения некорректно поставленных задач.

1.4.1 Метод регуляризации

Формальное решение некорректно поставленной обратной задачи может привести к неустойчивым, нереалистичным моделям. Метод регуляризации - это метод добавления некоторой дополнительной информации к условию с целью решить некорректно поставленную задачу [27, 28, 29]. Эта информация часто имеет вид штрафа за сложность модели. Например, это могут быть ограничения гладкости результирующей функции или ограничения по норме векторного пространства.

Список литературы диссертационного исследования кандидат наук Ли Цян, 2018 год

СПИСОК ЛИТЕРАТУРЫ

1. Алексеев А. С., Обратные динамические задачи сейсмики: Некоторые методы и алгоритмы интерпретации геофизических данных: М., Наука, 1967. С. 9-84.

2. Ампилов Ю. П., Барков А. Ю., Яковлев И. В., Филиппова К. Е., Приезжев И. И. Почти все о сейсмической инверсии. Часть 1 // Технологии сейсморазведки, 2009. № 4. С. 3-16.

3. Ампилов Ю.П. От сейсмической интерпретации к моделированию и оценке месторождений нефти и газа. М.: Геоинформмарк, 2008. 384 с.

4. Ампилов Ю.П., Барков А.Ю., Яковлев И.В., Роль сейсмической инверсии в геологическом моделировании морских газовых месторождений // Газовая промышленность. 2011. № 12. С. 69-74.

5. Багринцева К.И. Условия формирования и свойства карбонатных коллекторов нефти и газа. М.: РГГУ, 1999 (II). 285 с.

6. Беркинблит М.Б. нейронные сети: Учебное пособие. М.:МИРОС и ВЗМШ РАО, 1993. 96 с.

7. Боганик Г.Н., Гурвич И.И. Сейсморазведка. Учебник для вузов. Тверь: Издательство АИС, 2006. 744 с.

8. Василенко Г. И., Тараторин А. М. Восстановление изображений. М: Радио и связь, 1986. 304 с.

9. Воскресенский Ю.Н. Изучение изменений амплитуд сейсмических отражений для поисков и разведки залежей углеводородов. М.: РГУ нефти и газа, 2001. 68 с.

10. Воскресенский Ю.Н. Полевая геофизика. М.: Недра. 2010. 480 с.

11. Гилл Ф., Мюррей У, Райт М. Практическая оптимизация. М.: «Мир», 1985. 509 с.

12. Данько Д.А. Сравнение методов детерминистической акустической инверсии для выделения акустически контрастных объектов по сейсмическим данным // Геофизика, 2016. № 1. С. 2-11.

13. Данько Д.А., Рыжков В.И., Филимоненко С.В. Изучение перспективных акустически контрастных объектов Восточной Сибири методом пластовой амплитудной инверсии // EAGE GeoBaikal, Expanded Abstracts, 2016.

14. Жданов М.С. Теория обратных задач и регуляризации в геофизике. М.: Научный мир, 2007. 712 с.

15. Жиглявский А. А., Жилинкас А. Г. Методы поиска глобального экстремума. М.: Наука, Физматлит, 1991. 248 с.

16. Задорина Е.А. Исследование параметров геостатистической инверсии для прогноза коллекторских свойств по данным сейсморазведки: дис. ... канд. тех. наук : 25.00.10. М.: МГУ, 2015. 117 с.

17. Кащеев Д.Е., Кирнос Д.Г. Использование имитационного аннилинга для инверсии данных сейсморазведки // Геофизика, 2002. № Специальный выпуск. С. 75-79.

18. Кондратьев И.К., Бондаренко М.Т., Камнев С.П. Динамическая интерпретация данных сейсморазведки при решении задач нефтегазовой геологии // Геофизика, 1996. №5. С. 41-47.

19. Кондратьев И.К., Лисицын П.А., Киссин Ю.М. Детальность и точность решений в задаче сейсмической волновой инверсии // Геофизика, 2005. №3. С. 19-25.

20. Кондратьев И.К., Рыжков В.И., Бондаренко М.Т., Лапина Е.В. Эффективность прогнозирования коллекторов способами динамической интерпретации в Восточной Сибири // Технологии сейсморазведки, 2010. №4. С. 26-34.

21. Кондратьев И.К., Рыжков В.И., Киссин Ю.М., Шубин А.В. Способы реализации и оценка эффективности сейсмической инверсии. Учебное пособие. М.: РГУ нефти и газа, 2011. 63 с.

22. Кузнецов В.М., Жуков А.П., Шнеерсон М.Б. Введение в сейсмическую анизотропию: теория и практика. Тверь: ООО Изд-во ГЕРС, 2006. 160 с.

23. Ли Цян, Анализ устойчивости амплитудной сейсмической инверсии на примере статистических модельных данных. // Труды РГУ нефти и газа (НИУ)

имени И.М. Губкина. 2017. №3. С.15-28.

24. Ли Цян, Рыжков В. Анализ применимости различных методов регуляризации в задачах динамической инверсии сейсмических данных. // Технологии сейсморазведки. 2017. №2. С. 15-28..

25. Моисеев Н. Н., Иванилов Ю. П., Столярова Е. М. Методы оптимизации. М: Наука, 1978. 352 с.

26. Морозов В. А., О регуляризации некорректно поставленных задач и выборе параметра регуляризации, Ж. вычисл. матем. и матем. физ., 1966, том 6, номер 1, С.170-175.

27. Морозов В.А. Регулярные методы решения некорректно поставленных задач. М.: Наука., 1987. 240 с.

28. Морозов В. А. Об оптимальных методах решения некорректных задач // Вычислительные методы и программирование: Новые вычислительные технологии (Электронный научный журнал). 2006. Т. 7. С. 105-107.

29. Морозов В. А. Применение метода регуляризации к решению одной некорректной задачи // Вестник Московского университета. Серия 1: Математика. Механика. 1965. № 4. С. 13-21.

30. Приезжев И.И., Шмарьян Л.Е., Солоха Е.В. Методика сейсмической инверсии с помощью генетического алгоритма с последующим использованием результатов инверсии при моделировании коллекторских свойств резервуара // Технологии сейсморазведки, 2009. № 2. С. 18-23.

31. Тихонов А.Н. К математическому обоснованию теории электромагнитных зондирований // Журн. Вычислит, матем. и матем. физики. 1965. 5. 3. С.545-547.

32. Тихонов А.Н. О некорректных задачах линейной алгебры и устойчивом методе их решения // ДАН СССР, 1965, т. 163, № 3, с. 591-594.

33. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. Учебное пособие для вузов. М.: Наука, 1986. 287 с.

34. Шубин А.В. Рыжков В.И. Изучение эффекта засолонения порового пространства терригенного коллектора по сейсмическим данным // Геофизика. 2013. №5. С. 17-25.

35. Шубин А.В. Методика изучения сложно построенных природных резервуаров на основе петроупругого моделирования и инверсии сейсмических данных: дис. ... канд.техн. наук: 25.00.10. М.: РГУ нефти и газа им. И.М. Губкина, 2014. 146 с.

36. Яковлев И. В., Ампилов Ю. П., Филиппова К. Е. Почти все о сейсмической инверсии. Часть 2 // Технологии сейсморазведки, 2011. № 1. С. 5-15.

37. Aki L., Richards P.G. Quantitative seismology, 2nd edition. University Science Books, 2009. 700 p.

38. Arsalan S., Yadav A. Application of extended elastic impedance: A case study from Krishna-Godavari Basin, India: The Leading Edge, 2009. № 28. Р. 1204-1209.

39. Blaschek R., Hördt A., and Kemna A. A new sensitivity-controlled focusing regularization scheme for the inversion of induced polarization data based on the minimum gradient support. GEOPHYSICS, 2008,V0L. 73, NO. 2 MARCH-APRIL; P. F45-F54.

40. Bogert B. P., Healy M. J. R., and Tukey J. W.: "The Quefrency Alanysis [sic] of Time Series for Echoes: Cepstrum, Pseudo Autocovariance, Cross-Cepstrum and Saphe Cracking". Proceedings of the Symposium on Time Series Analysis (M. Rosenblatt, Ed). New York: Wiley, 1963. Chapter 15. P. 209-243.

41. Buland, A. Bayesian seismic AVO inversion: Thesis for the degree of doctor, Norwegian University of Science and Technology. 2002. 159p.

42. Castagna, J.P., Batzle, M.L., and Eastwood, R.L. Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks: Geophysics. 1985. 50. P. 571-581.

43. Chopra Satinder and Castagna John, AVO. Society of Exploration Geophysicists. 2014. 288p.

44. Connolly P. Elastic impedance // The Leading Edge. 1999.18. P. 438-452.

45. Cooke, D.A., Schneider W.A. Generalized linear inversion of reflection seismic data. Geophysics. 1983. 48. P. 665-676.

46. Dubrule O. Geostatistics for seismic data integration in earth models. SEG-EAGE Publications, 2003. 273 p.

47. Edgar, J.A. and van der Baan, M. How reliable is statistical wavelet estimation? Geophysics. 2011. 76, P. 59-68.

48. Fatti, J., G. Smith, P. Vail, P. Strauss, and P. Levitt. Detection of gas in sandstone reservoirs using AVO analysis: a 3D seismic case history using the Geostack technique: Geophysics. 1994. 59. P. 1362-1376.

49. Francis, A. Understanding stochastic inversion: Part 1: First Break. 2006a. 24. P. 69-77.

50. Francis, A. Understanding stochastic inversion: Part 2: First Break. 2006b. 24. P. 79-84.

51. Gardner, G.H.F., Gardner, L.W. and Gregory, A.R. Formation velocity and density -The diagnostic basics for stratigraphic traps: Geophysics. 1974. 50. P. 2085-2095.

52. Goodway, W., J. Varsek, and C. Abaco. Practical applications of P-wave AVO for unconventional gas resource plays, part 2: Detection of fracture prone zones with azimuthal AVO and coherence discontinuity: CSEG Recorder, 31. 2006. № 4. P. 52-65.

53. Gunning J. and Glinsky M. Wavelet Extractor: A Bayesian well-tie and wavelet extraction program. Computers and Geosciences. 2006. 32. P. 681-695.

54. Hadamard J. Sur les problemes aux derivees partielles et leur signification physique // Princeton Univ. Bull. 1902. 13. P. 49-52: reprinted in his Oeuvres, Vol. HI, Centre Nat. Recherche Sci., Paris, 1968. P. 1099-1105.

55. Hampson D. AVO inversion, theory and practice: The Leading Edge. 1991.№ 10. P. 39-42.

56. Hampson D.P., Russell B.H., Bankhead B. Simultaneous inversion of pre-stack seismic data // SEG Annual Meeting, Expanded Abstracts. 2005. Pp. 1633-1637.

57. Hansen P.C. Rank-deficient and discrete ill-posed problems. Numerical aspects of linear inversion. Lyngby: Department of mathematical modeling. Technical University of Denmark, 1998. 247 p.

58. Hansen, P. C. The truncated SVD as a method for regularization. BIT. 1987. 27. P. 534-553.

59. Hansen. P. C. Discrete Inverse Problems: Insight and Algorithms. SIAM. 2010. 213p.

60. Ikelle, L.T., Roberts, G., and Weglein, A.B. Source signature estimation based on the removal of first-order multiples. Geophysics. 1997. 62. P. 1904-1920.

61. Lancaster, S. and Whitcombe, D. Fast track "coloured" inversion. Expanded abstracts, 70th SEG Annual Meeting, Calgary. 2000. 1572-1575.

62. Li Qiang, Yang Pei, and Ryzhkov Valeriy Ivanovich. Focusing AVO inversion based on the minimum gradient support regularization. SEG Technical Program Expanded Abstracts 2017: P. 682-686.

63. Lindseth, R.O. Synthetic sonic logs - A process for stratigraphic interpretation: Geophysics.1979. 44. P. 3-26.

64. Ma Hong-Da, Roy E. White, and Hu Tian-Yue. Wavelet estimation by matching well-log, VSP, and surface seismic Data. APPLIED GEOPHYSICS. 2010.Vol.7, No.4. P. 384 - 391.

65. Oldenburg D. W., T. Scheuer, and S. Levy. Recovery of the acoustic impedance from reflection seismograms: Geophysics. 1983. 48. P.1318-1337.

66. Oppenheim A. V., Schafer R. W., Homomorphic analysis of speech, IEEE Trans. Audio Electroacoust. AU-16. 1968. P. 221-226.

67. Portniaguine O., and Castagna J. P. Inverse spectral decomposition: 74th Annual International Meeting, SEG, Expanded Abstracts. 2004. P.1786-1789.

68. Portniaguine O., and Zhdanov M. S. Focusing geophysical inversion images: Geophysics. 1999. 64. P. 874-887.

69. Portniaguine O., and Zhdanov M. S. 3-D magnetic inversion with data compression and image focusing: Geophysics. 2002. 67. P.1532-1541.

70. Rudin L.I., Osher S., Fatemi E. Nonlinear total variation based noise removal algorithms // Physica D. 1992. 60. P. 259-268.

71. Ruger A. P-wave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry: Geophysics. 1997. 62. P. 713-722.

72. Ruger A. Variation of P-wave reflectivity with offset and azimuth in anisotropic media: Geophysics. 1998. 63. P. 935-947.

73. Ruger A. Reflections coefficients and azimuthal AVO analysis in anisotropic media. Geophysical Monograph Series. SEG Publication P. 2002. P. 39-62.

74. Russell B. H. Introduction to Seismic Inversion Methods. Course Notes Series, Society of Exploration Geophysicists. 1988. P. 80-101.

75. Russell B.H. and Hampson D.P. Comparison of poststack seismic inversion methods // SEG Annual Meeting, Expanded Abstracts. 1991. P. 876-878.

76. Sen Mrinal K., Stoffa Paul L. Global optimization methods in geophysical inversion // Cambridge University Press - Second edition. 2013. 289p.

77. Silvia M. T. and Robinson E.A. Use of the cepstrum in signal analysis. Geoexploration. 1978. vol.16. P. 55-73.

78. Smith R.T., Zoltani C.K., Klem G.J., Coleman M.W. Reconstruction of the tomographic images from sparse data sets by a new finite element maximum entropy approach // Appl. Opt. 1991. 30. P. 573-582.

79. Tarantola A. Inverse Problem Theory. Amsterdam, Oxford, New York, Tokyo: Elsevier, 1987. 613 p.

80. Thomsen L. Weak elastic anisotropy // Geophysics. 1986. № 51. P. 1954-1966.

81. Tria Mohamed, Mirko Van Der Baan, Larue Anthony, Mars Jerome. Wavelet estimation in homomorphic domain by spectral averaging, for deconvolution of seismic data. Fifth International Symposium on Physics in Signal and Image Processing, PSIP2007 Mulhouse, France. 2007. P. 2-6.

82. Vogel C.R., Oman M.E. Fast total variation based reconstruction of noisy, blurred images // IEEE T. Image Process. 1998. 7. P. 813-824.

83. Wahba, G. Spline Models for Observational Data, CBMS-NSF Regional Conference Series in Applied Mathematics 59, SIAM, Philadelphia, 1990. 180p.

84. Walden, A. T., and R. E. White. On errors of fit and accuracy in matching synthetic seismograms and seismic traces: Geophysical Prospecting. 1984. 32. P. 871-891.

85. Walden, A. T., and White R. Seismic wavelet estimation: a frequency domain solution to a geophysical noisy input-output problem. IEEE Transactions on Geoscience and Remote Sensing 36(1), 1998. P. 287-297.

86. Walden, A.T. and Hosken, J.W.J. An investigation of the spectral properties of primary reflection coefficients. Geophysical Prospecting. 1985. 33. 400-435.

87. Wernecke S.J., D'Addario L.R. Maximum entropy image reconstruction // IEEE T. Comput. 1977. 26. P. 351-364.

88. Whitcombe, D. N. Elastic impedance normalization: Geophysics. 2002. 67. P. 60-62.

89. Whitcombe, D. N., Extended elastic impedance for fluid lithology prediction: Geophysics. 2002. 67. P. 63-67.

90. White R.E. Partial coherence matching of synthetic seismograms with seismic traces. Geophysical Prospecting. 1980. 28. P. 333-358.

91. White R.E., Simm, R. Tutorial: Good practice in well ties. First Break. 2003. 21. P. 75-83.

92. Zhdanov M. S. Geophysical Inverse Theory and Regularization Problems. Elsevier - Second edition. 2015. 704p.

93. Zoeppritz Karl. Erdbebenwellen VII. VIIb. Über Reflexion und Durchgang seismischer Wellen durch Unstetigkeitsflächen. Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen // Mathematisch-physikalische Klasse. 1919, 66-84.