Миграция сейсмических данных в истинных амплитудах на основе метода параксиальной аппроксимации тема диссертации и автореферата по ВАК РФ 25.00.10, кандидат физико-математических наук Киященко, Денис Александрович

• Киященко, Денис Александрович
• кандидат физико-математических наук
• 2006

Оглавление диссертации кандидат физико-математических наук Киященко, Денис Александрович

Введение.1

Глава 1. Решение волнового уравнения в параксиальном приближении . 14

1.1 Параксиальная аппроксимация решения уравнения Гельмгольца. Уравнения квазипараболического типа.14

1.2 Решение неоднородного уравнения Гельмгольца. Алгоритм продолжения волнового поля .18

1.3 Поглощающие граничные условия.21

1.4 Моделирование волновых полей в трехмерном случае.22

1.5 Корректность определения амплитуд методом параксиальной аппроксимации в однородной среде.25

Глава 2. Метод многоступенчатой параксиальной аппроксимации.28

2.1 Уравнение Гельмгольца и уравнения квазипараболического типа.28

2.2 Приближенное решение уравнения Гельмгольца методом многоступенчатой параксиальной аппроксимации.30

2.3 Граничные условия для уравнений метода многоступенчатой параксиальной аппроксимации.31

2.4 Решение уравнений квазипараболического типа с ненулевой правой частью . .34

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

3.1 Вертикально-неоднородные модели.38

3.2 Горизонтально-неоднородные модели.38

3.3 Модели со скачкообразным изменением скорости.47

3.4 Модель соляного тела БЕС/БАСЕ.53

3.5 Трехмерная модель.59

Глава 4. Метод многоступенчатой параксиальной аппроксимации и другие подходы к решению проблемы корректной оценки амплитуды волнового поля методами параксиальной аппроксимации.63

4.1 Метод параксиальной аппроксимации с лучевой поправкой за неоднородность среды.63

4.2 Модифицированный метод многоступенчатой параксиальной аппроксимации .66

4.3 Связь метода локально-многоступенчатой параксиальной аппроксимации и метода параксиальной аппроксимации с лучевой поправкой за неоднородность среды .68

4.4 Вопросы практической реализации алгоритма продолжения волново- го поля в нижнее полупространство.70

Глава 5. Модифицированный принцип построения изображения. Связь с методом пссвдоинверсии . 74

5.1 Модифицированный принцип построения изображения.74

5.2 Пример применения модифицированного принципа построения изображения для миграции.77

5.3 Метод псевдоинверсии.80

5.4 Высокочастотная асимптотика модифицированного принципа построения изображения .82

5.5 Эффекты, связанные с ограниченностью апертуры системы наблюдений . 84

5.6 Примеры применения миграции на основе модифицированного принципа построения изображения: синтетические данные.89

5.7 Результат миграции сейсмических данных Мексиканского залива . 97

5.8 Алгоритм миграции в трехмерном случае.99

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

Обзор работ по теме диссертации. На сегодня методы миграции, основанные на решении волнового уравнения в параксиальном приближении (см. [1], [2]), широко используются в сейсморазведке для построения изображений среды. Однако существенным недостатком большей части разработанных методов является невозможность корректного определения параметров среды по данным изображениям. Задачей данной работы являлась разработка метода миграции сейсмических данных, позволяющего восстанавливать параметры среды на основе методов параксиальной аппроксимации. Исследования в данной работе были ограничены акустическим случаем с постоянной плотностью, при котором процессы распространения волн описывались волновым уравнением.

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

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

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

Идея параксиальной аппроксимации была впервые предложена в работах [5] и [б], где рассматривалось решение волнового уравнения, отвечающее волнам, распространяющимся в направлениях, близких к выделенному. Для решения задач геофизики этот подход был впервые использован в работах Клербо (см. [1], [2]). Методы параксиальной аппроксимации позволяют корректно оценивать время прихода волн при условии, что угол между направлением их распространения и выделенным направлением (вертикалью) не превышает некоторого известного значения. Так, в работе [2] были предложены два типа параксиальной аппроксимации: 15-градусная и 45-градусная, для которых вышеупомянутое значение максимально допустимого угла составляет соответственно 15 и 45 градусов. С течением времени появились схемы, для которых значение данного угла составляет 60 и более градусов (см. [7|, [8], |9]). Вышеупомянутые реализации параксиальной аппроксимации сводятся к конечно-разностному решению уравнений параболического (в случае 15-градусной аппроксимации) или квазипараболического (в случае аппроксимации высшего порядка) типа. Уравнением квазипараболического типа здесь и в дальнейшем названо уравнение, символьная запись которого асимптотически совпадает с уравнением параболы при малых углах между вертикалью и направлением распространения волны. Процедура решения данных уравнений представляет собой продолжение волнового поля в частотной области вдоль выделенного пространственного направления. Так как эта группа методов использована в настоящей работе, их описание дано в главе 2.

Альтернативным подходом является продолжение волновых полей в области частота-волновое число. Этот подход впервые был предложен в работах [10], [11] (метод фазового сдвига) и со временем породил группу методов GSP (Generalized Screen Propagators, см. [12-18]). Методы GSP основаны на учете вклада двух составляющих волнового поля на каждом шаге его продолжения через тонкую прослойку среды. Первая составляющая расчитывается методом фазового сдвига для однородной опорной модели слоя. Вторая составляющая учитывает вклад неоднородности прослойки в приближении Борна.

Методы параксиальной аппроксимации позволяют корректно моделировать время прихода волн в некотором диапазоне направлений их распространения. Однако, эти методы обладают следующим недостатком: амплитуда волн, полученная на их основе, верна лишь в однородной среде. Таким образом, классические методы параксиальной аппроксимации позволяют получить лишь структурное изображение среды, но не информацию о ее параметрах.

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

В связи с проблемой оценки амплитуды волновых полей с использованием методов параксиальной аппроксимации следует упомянуть работы [21], [22], в которых предлагается интересный подход к описанию взаимосвязи волнового уравнения и уравнений квазипараболического типа, описывающих распространение волн в противоположных направлениях. Волновое уравнение разбивается на два зависимых уравнения квазипараболического типа, которые становятся независимыми в случае вертикально-однородной среды. Такое представление позволяет учитывать влияние вертикального градиента скорости в среде на волновое поле. Однако, полученные уравнения содержат оператор квадратного корня из оператора Гельмгольца, который не может быть выписан явным образом. Аппроксимация данного оператора ведет к ошибкам в оценке ампрлитуд волновых полей. В работе [21] предлагается выразить оператор квадратного корня через собственные числа и собственные функции (вектора в дискретном случае) оператора Гельмгольца. Но практическое использование такого подхода затруднено в силу чрезвычайных вычислительных затрат, требуемых для его реализации в трехмерном случае. В работе [22] предлагается использовать однородное асимптотическое разложение для получения функции Грина уравнения квазииараболического типа. Однако метод вычислений, следующий из построений работы [22], является слишком громоздким для практического использования.

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

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

Первый подход предложен в работах Клербо [1, 2]. Данный метод позволяет восстанавливать коэффициент отражения от границ и основан на обработке данных, отсортированных по пункту взрыва. Изображение среды строится как результат кросс-корреляции волнового поля, созданного источником, с волновым полем, полученным в результате обращенного продолжения ноля, зарегистрированного приемниками. В отечественной литературе вопросы построения обращенно-продолженного волнового поля обсуждались в работах [23, 24, 25]. В работе [23] было показано, что восстановление положения отражающих границ определяется кинематическими характеристиками волновых полей. Для того, чтобы восстановить коэффициент отражения, результат кросс-корреляции нормируется на мощность волнового поля, созданного источником. Этот принцип построения изображения среды будет далее в работе называться классическим. При вышеописанной постановке задачи не делается никаких предположений относительно метода моделирования волнового поля. Таким образом, принцип построения изображения Клербо может применятся совместно с лучевыми и конечно-разностными методами моделирования.

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

В восьмедесятых годах задача миграции сейсмических данных была сформулирована как обратная задача в терминах метода наименьших квадратов, см. [26], или в терминах обобщенного преобразования Радона, см. [27]. В отечественной литературе постановка задачи восстановления параметров среды как обратной геофизической задачи, решаемой методом наименьших квадратов, обсуждалась в работах [28, 29].

Использование лучевого метода для представления волновых полей позволило разработать методы псевдоинверсии, см. [27], [30-33]. Идея методов псевдоинвсрсии заключается в подборе такого оператора отображения из пространства сейсмических данных в пространство параметров среды, копмозиция которого с высокочастотным представлением оператора моделирования (отображения из пространства параметров в пространство данных) давала бы оператор тождественного равенства в высокочастотном приближении и при неограниченной апертуре сейсмических наблюдений. Такой подход напоминает метод Бейкуса-Гильберта, см. [34]. К этой-же группе методов можно отнести методы миграции Киргхоффа, см. [35-37]. Миграция Киргхоффа основана на взвешенном диффракционном преобразовании сейсмических данных. Данное преобразование определяет отображение из пространства данных в пространство параметров. Весовая функция подбирается так, чтобы параметры восстанавливались корректно. Такой подход позволил также разработать метод восстановления параметров анизотропных сред, см. [38]. Миграция Киргхоффа тесно связана с методами псевдоинверсии, см. [37], и классическим изображающим принципом, см. [39].

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

Теоретически, наиболее общий подход к решению обратных задач сейсмической разведки основан па методе наименьших квадратов, см. [26, 29]. При постановке задачи на основе метода наименьших квадратов имеется возможность работы с любой конфигурацией системы наблюдений, в том числе и системы с неплотным покрытием участков профилей (площадей) приемниками (источниками). Также подход на основе метода наименьших квадратов позволяет работать с любыми методами моделирования волновых полей. Использование регуляризации позволяет учитывать априорную информацию при восстановлении параметров среды.

В работах [26, 29, 41] было показано, что градиент функционала невязки между модельными и измеренными данными, взятый но параметрам среды, имеет структуру оператора миграции на основе классического изображающего принципа (т. е. формула для градиента представляет из себя результат кросс-корреляции между полем, порожденным источником, и полем, полученным в результате обратного продолжения волновых полей, зарегистрированных приемниками). Но сам по себе градиент функционала невязки не дает истинных значений параметров среды. Набор параметров среды, доставляющий минимум функционалу невязки, получается в результате умножения матрицы, обратной Гессиану (матрицы вторых производных функционала невязки по параметрам), на градиент функционала (с обратным знаком). На сегодняшний день ввиду больших размеров матрицы Гессиана не представляется возможным рассчитать обратную матрицу. Для того чтобы аппроксимировать обратную матрицу, был предложен ряд диагональных аппроксимаций матрицы Гессиана, см. [42, 43, 44]. Однако, на практике матрица Гессиана не является диагональной, поэтому подходы, основанные на диагональной аппроксимации, могут дать правильную оценку параметров в ограниченном числе случаев. Этот вопрос будет рассмотрен в главе 7 настоящей работы. Тем не менее, для того, чтобы улучшить оценку параметров среды, в рамках метода наименьших квадратов существует возможность итерационного подхода. Итерационный алгоритм восстановления параметров среды может быть также использован в случае недостаточного количества источников (приемников) на поверхности для построения качественного изображения среды в результате миграции (первого шага итерационного алгоритма). Однако, применение итерационных алгоритмов требует значительных вычислительных затрат.

Таким образом, рассмотрены основные методы восстановления параметров среды: метод, основанный на классическом принципе построения изображений, метод псевдо-инверсии и метод наименьших квадратов. Метод псевдо-инвсрсии основан на лучевом методе расчета волновых полей. Однако, когда модель среды является достаточно сложной, возникают зоны многолучевости и каустик. Это приводит к практическим трудностям вычисления волновых полей лучевым методом в таких зонах, а также невозможности корректной оценки амплитуды волновых полей в зонах каустик. В таких случаях более предпочтительно рассчитывать волновые поля конечно-разностными методами, в том числе и основанными на параксиальной аппроксимации. Использование метода псевдо-инверсии совместно с конечно-разностными методами расчета волновых полей не представляется практически возможным в силу излишней вычислительной громоздкости такого подхода. Следовательно, в данной работе метод псевдо-инверсии в его классической формулировке не используется для практических расчетов. Для того чтобы построить алгоритм восстановления параметров среды с использованием конечно-разностных методов моделирования, нужно либо адаптировать метод псевдо-инверсии для использования совместно с конечно-разностным методом, либо модифицировать классический принцип построения изображения, либо использовать метод наименьших квадратов, основанный на диагональной аппроксимации матрицы Гессиана, см. [43, 44].

В этом отношении следует упомянуть метод восстановления коэффициентов отражения волн от границ раздела в среде, предложенный в работе [45]. Метод является по сути коррекцией аналогичного метода, предложенного в работе [46]. Данный метод предназначен для обработки данных многократного перекрытия и основан на I продолжении волновых полей в нижнее полупространство методами параксиальной аппроксимации. Первым шагом такого алгоритма является снятие эффекта среды выше некоторого заданного глубинного уровня. Такое снятие достигается 'перемещением' источников и приемников на заданный уровень и осуществляется с использованием теоремы взаимности и алгоритмов продолжения волновых полей. Вторым шагом данного алгоритма является вычисление коэффициента отражения в результате преобразования типа Радона сейсмических данных, 'перемещенных' на заданный глубинный уровень. Вышеописанный метод требует достаточно плотного покрытия профилей приемниками и источниками. Также общей проблемой технического характера, касающейся методов этого типа, является трудность хранения огромных объемов информации, требующихся для вычисления коэффициентов отражения для ряда углов падения в каждой точке среды.

Таким образом, на пути восстановления параметров на основе методов параксиальной аппроксимации стоит две проблемы: проблема некорректной оценки амплитуды в рамках классической параксиальной аппроксимации и проблема разработки метода восстановления параметров среды по данным многократного перекрытия. Цели работы

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

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

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

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

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

4. В высокочастотном приближении модифицированный принцип построения изображения тесно связан с методом псевдо-инверсии. На основе изучения этой связи предлагается формулировка модифицированного изображающего принципа, позволяющая восстановить относительное возмущение медленности как функцию пространственных координат. Также взаимосвязь между модифицированным изображающим принципом и методом псевдо-инверсии позволяет показать, что суммирование изображений среды, построенных по данным с общим путктом взрыва (common shot data), с изображениями, построенными по данным с общим пунктом приема (common receiver data), дает в высокочастотном приближении тот-же результат, что и суммирование изображений, построенных по данным с общим удалением приемника (common-oiFset data) методом псевдо-инверсии, которые являются на практике более информативными, чем изображения на основе данных с общим пунктом взрыва и общим пунктом приема.

5. Показано, что в высокочастотном приближении миграция на основе метода наименьших квадратов приводит к тому-же результату, что и метод псевдо-инверсии, при условии, что волновые поля, рассеянные на неоднородностях среды, могут быть описаны в приближении Борна, а также при условии неограниченности апертуры сети приемников. Учитывая связь между модифицированным изображающим принципом и методом псевдо-инверсии, этот результат позволяет утверждать, что модифицированный принцип построения изображения может быть использован в итерационном алгоритме уточнения параметров среды. Следует отметить, что метод исевдо-инверсии использовался для этой цели в работах [47, 48]. Полученный в данной работе результат формально оправдывает эти подходы. Практическая значимость результатов

Предложена методология миграции сейсмических данных, основанная на методе многоступенчатой параксиальной аппроксимации и модифицированном принципе построения изображения, см. [64-70]. Метод многоступенчатой параксиальной аппроксимации является возможным решением проблемы некорректной оценки амплитуд волновых полей классическими методами параксиальной аппроксимации. По построению реализация метода многоступенчатой параксиальной аппроксимации основана на классическом методе параксиальной аппроксимации. В данной работе применялся метод параксиальной аппроксимации, который был предложен в работах [8] и [9], основанный на рациональной аппроксимации оператора квадратного корня и конечно-разностной схеме продолжения волнового поля внутрь исследуемой среды. Но метод многоступенчатой параксиальной аппроксимации может применяться совместно с любым другим алгоритмом продолжения волнового поля на основе параксиальной аппроксимации, который дает верную оценку времен прихода волн в некотором диапазоне направлений их распространения, но не дает корректной оценки их амплитуд, например, с методом, основанном на методах вЭР (см. [12], [17]).

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

Основные положения, выносимые на защиту

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

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

3. Предлагаемая методология миграции опробована на ряде полевых и синтетических сейсмичеких данных, в том числе и созданных для модели соляного тела, принятой в международных геофизических обществах БЕС и ЕАСЕ в качестве эталона для тестирования алгоритмов миграции.

Краткое содержание работы

В главе 1 дано описание методов параксиальной аппроксимации, предложенных в работе [9], которые были использованы автором для расчетов волновых полей.

В разделе 1.1 введено понятие параксиальной аппроксимации решения однородного уравнения Гельмгольца, т.е. решения, описывающего распространение волн в направлениях, близких к выделенному. Вводится понятие уравнений квазипараболического типа: и(х,г,ш) = 1Ааи(кх,г,ш), (1) где

К = + (2) 1 п - целое положительное число, Ла - псепдо-дифференциальный оператор, который является аппроксимацией оператора Л = ^к2 + ^, к - волновое число. Решение уравнений вида (1) в случае однородной среды корректно описывает распространение воли в направлениях, близких к выделенному. Практическое использование этих уравнений обусловлено тем, что в неоднородной среде их решение корректно воспроизводит времена прихода волн. Следовательно, миграция на основе уравнений квазипараболического типа дает корректные структурные изображения среды.

В разделе 1.2 описывается способ решения неоднородного уравнения Гельмгольца с правой частью, отвечающей точечному источнику внешних сил. Приведено описание алгоритма продолжнения волнового поля вглубь исследуемой среды. В разделе 1.3 дано описание поглощающих граничных областей, позволяющих подавить отражения от границ модели. В разделе 1.4 дано описание алгоритма продолжения волнового поля в трехмерном случае. В разделе 1.5 показано, что даже в случае однородной среды метод параксиальной аппроксимации дает корректную оценку амплитуды волнового поля лишь в ограниченной ее области.

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

В разделе 2.1 однородное уравнение Гельмгольца представлено в терминах операторов продолжения волновых нолей, используемых в методе параксиальной аппроксимации:

§-г+гЛа)(£-гК)и = Еи, (3) где

Е = + (АЦ - А2). (4) является оператором расхождения, содержащим информацию о неоднородности среды.

В разделе 2.2 предлагается решение уравнения (3), основанное на методе последовательных приближений, которое сводится к последовательному решению трех уравнений квазипараболического типа: = гАаи0, (5) дг гАа)щ/2 = Ещ, (6) - г'Ла)"1 = «1/2. (7) щ является решением уравнения (3) в нулевом приближении и отвечает методу классической параксиальной аппроксимации, щ является поправкой, учитывающей влияние неоднородности среды и ошибки аппроксимации оператора Л0. ио + щ является решением (3) в первом приближении. Последовательное решение уравнений (5), (б), (7) названо методом многоступенчатой параксиальной аппроксимации.

В разделе 2.3 получены граничные условия для уравнений метода многоступенчатой параксиальной аппроксимации: ио|г=+0 = (8) где к(х) является решением уравнения:

9) с поглощающими граничными условиями; 0; (ю) где является решением уравнения:

Дх + ^^оу)Л1=«1/2ио- (12)

Таким образом, практическая реализация метода многоступенчатой параксиальной аппроксимации сводится к трем шагам:

1. Решение уравнения (5) с граничным условием (8) на верхней границе модели (г = 0). На данном этапе вычисляется значение оператора расхождения Ещ на поле и0.

2. Решение уравнения (б) с правой частью Ещ, вычисленной на предыдущем шаге, и граничным условием (10) на нижней границе модели (г — гтах).

3. Решение уравнения (7) с правой частью 111/2 и условием (11) на верхней границе модели. В результате, вычисляется поправка г^ к решению уравнения Гельмгольца в нулевом приближении. Итоговое решение дается суммой щ + щ.

В разделе 2.4 приведено формальное решение неоднородных уравнений квазипараболического типа (6) и (7) с граничными условиями (10), (И) и указан способ его практической реализации.

В главе 3 приведены примеры опробования метода многоступенчатой параксиальной аппроксимации на различных моделях сред. Для сравнения амплитуд волновых полей, полученных в рамках классической и многоступенчатой параксиальной аппроксимаций, была рассчитана относительная разность г = (А — Н)/Н между амплитудами А, даваемыми параксиальной аппроксимацией (классической или многоступенчатой) и амплитудой II результата конечно-разностного решения уравнения Гельмгольца.

В разделах 3.1 и 3.2 приведены результаты сравнения на вертикально-неоднородных и горизонтально-неоднородных скоростных моделях. В разделе 3.3 приведены результаты сравнения на моделях сред, содержащих разрывы скорости распространения воли (отражающие границы) и показано, что метод многоступенчатой параксиальной аппроксимации позволяет моделировать сейсмограммы отраженных волн. В разделе 3.4 метод многоступенчатой параксиальной аппроксимации опробован на сглаженной модели соляного тела БЕС/ЕАСЕ. В разделе 3.5 приведен пример применения метода многоступенчатой параксиальной аппроксимации в трехмерном случае.

В главе 4 предложен метод локально-многоступенчатой параксиальной аппроксимации, позволяющий скорректировать амплитуду волн, распространяющихся в определенном направлении, но не воспроизводящий отраженных волн. Необходимость такой модификации вызвана тем, что наличие отраженных волн при миграции сейсмических данных может приводить к артефактам в изображениях среды. Проведено сравнение метода локально-многоступенчатой параксиальной аппроксимации с методом параксиальной аппроксимации с лучевой поправкой за неоднородность среды (предложенным в [20]).

В разделе 4.1 приведено описание метода параксиальной аппроксимации с лучевой поправкой за неоднородность среды (см. [20]). Данный метод основан на решении модифицированного уравнения квазипараболического типа: - гЛ'и -ви = 0. (13) аг

Лд является интегральной аппроксимацией оператора квадратного корня, задаваемой выражением к{и = к{и + - С ¿5\/1 - з'21тХ{з,х,г,ш){у^х)2и), (14)

7Г J-l где г(5, х, 2, и) = и2 + з2(г;У1)2, (15) г;Ух)2 = Уг^у-^, а результат действия оператора Ь?1 на функцию / (д = Ь^1}) определяется решением уравнения: и2 + з2(ь^х)2)д = /. (16)

Уравнение переноса для амплитуды высокочастотного решения уравнения (13) совпадает с уравнением переноса, полученного для уравнения Гельмгольца. Однако путь практической реализации решения уравнения (13) не вполне ясен.

В настоящей работе предлагается метод локально-многоступенчатой параксиальной аппроксимации в разделе 4.2. Этот метод основан на решении уравнений метода многоступенчатой параксиальной аппроксимации (5), (6) и (7) не во всей области модели, от 2г = 0 до г = гтах, а локально в каждом интервале глубин (г\,г\ 4- Дг). Таким образом можно осуществить продолжение волнового поля с верхней границы {г\) интервала глубин на нижнюю его границу (21 + Лг) и построить волновое поле во всей среде. Предельный переход при Аг —* 0 позволяет получить уравнение для волнового поля в рамках метода локально-многоступенчатой параксиальной аппроксимации: ди - ih.aU -Ти = 0, (17) ах где

Ти~шЕи. (18)

Оказывается, что если Л„ = А1а (где Л„ задано выражением (14)), уравнения переноса для уравнений (13) и (17) совпадают. Это показано в разделе 4.3. В разделе 4.4 обсуждаются вопросы практической реализации алгоритмов продолжения волновых полей, описанных в данной главе.

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

В разделе 5.1 обсуждаются недостатки классического принципа построения изображения (который предназначен для восстановления коэффициента отражения В.(9), зависящего от угла падения в по данным общего пункта взрыва) и предлагается модифицированный принцип построения изображения среды, предназначенный для восстановления относительного возмущения медленности среды, определяемого следующим образом: г(х) = 2!<М, (и) где <Т( - истинная медленность среды, зависящая от координат, и сг0 - медленность в опорной модели среды, используемой для миграции.

Модифицированный принцип построения изображения среды формулируется в виде: 1 дщ(х„х, и) дЩ(х*, х, и>) duj(xs,x,u>) дщ(хя,х,ш) ох ох cfz oz

20)

Результат применения данного принципа /¿(х3, х) дает оценку функции г+(х), которая совпадает с г(х) на верхней кромке отражающих границ. Изображение 1^(хя,х) выгодно отличается от результата применения классического принципа построения изображения тем, что почти не зависит от положения пункта взрыва. Поэтому при суммировании изображений, построенных по данным разных пуктов взрыва, информация о параметрах среды сохраняется.

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

Сравнение модифицированного принципа построения изображения в высокочастотном приближении с методом нсевдо-инверсии, предложенным в [30], позволяет сформулировать модифицированный принцип построения изображения таким образом, чтобы он позволял восстанавливать функцию г(х): г(х) » Im(*s,х) = Л>еп duisgn(u)(iid(x„ х,ш)щ(хя, х,ш)+ ,21)

1 rl диы(х,,х,ш) 1 9щ(х,,х,и>) , 1 дил(х,,х,ш) 1 дЩ(х,,х,ш)-1\ к2 И дх i дх ^ 1 dz i dz J/'

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

В разделе 5.6 приведены примеры опробования разработанной методологии на примере синтетических данных морского типа. Один из примеров посвящен обработке данных, созданных для модели соляного тела, которая считается в международных геофизических обществах ЭЕС и ЕАСЕ одним из эталонов для тестирования методов миграции.

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

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

В разделе 6.1 приводится интегральное уравнение для определения возмущения параметров среды. В матричной записи ему аналогично уравнение

ВтВг = Втс1, (22) где г - искомое возмущение параметров среды, й - данные, измеренные на дневной поверхности, В - оператор, сопоставляющий малому возмущению параметров среды сейсмические данные на дневной поверхности. Оператор ВТВ называется Гессианом. Оператор Вт является оператором миграции сейсмических данных.

В разделе 6.2 приводится решение данного уравнения в высокочастотном приближении. Оказывается, что результат решения полностью аналогичен результату применения метода псевдо-инверсии, и, следовательно, результату применения модифицированного принципа построения изображения среды в высокочастотном приближении.

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

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

Практическая реализация миграции па основе метода наименьших квадратов (см. [44]) основана на диагональной аппроксимации Гессиана: ВТВ ~ На = (Иад{ВтВ). Оценка параметров среды является взвешенным результатом миграции данных: г ~ УГН~1ВТ<1, (23) где IV - диагональная весовая матрица (см. раздел 7.1).

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

Список литературы диссертационного исследования кандидат физико-математических наук Киященко, Денис Александрович, 2006 год

1. J.F. Claerbout, Toward a unified theory of reflector mapping, Geophysics, V. 36, N. 3, pp. 467-481, 1971.

2. J.F. Claerbout, Imaging the earth's interior, Blackwell Scientific Publication со., 1983.

3. W.A. Mulder and R.-E. Plessix, A comparision between one-way and two-way wave-equation migration, Geophysics, V. 69 N. 6, pp. 1491-1504, 2004.

4. M. A. Leontovich and V.A. Fock, Solution of the problem of propagation of electromagnetic waves along the earth's surface by the method of parabolic equation, J. Exp. Theoret. Phys., 16, pp. 557-573, 1946.

5. F. Tappert, The parabolic approximation method, in: J.B. Keller, J.S. Papadakis (Eds.), Wave Propagation and Underwater Acoustics, American Institute of Physics, v. 70, Springer, New York, pp. 224-280, 1977.

6. D. Lee and A.D. Pierce, Parabolic equation development in recent decade, J. of Сотр. Acoustics, 3(2), pp. 95-173, 1995.

7. D. Ristow and T. Ruhl, 3-D implicit finite-difference migration by multi-way splitting, Geophysics, V. 62 N. 2, pp. 554-567, 1997.

8. F. Collino and P. Joly, Splitting of operators, alternate directions, and paraxial approximations for three-dimensional wave equation, SIAM J. Sci. Comput., Vol. 16, No. 5, pp. 1019-1048, 1999.

9. Gazdag J., Wave equation migration with the phase shift method, Geophysics, 43, pp. 1342-1351, 1978.

10. J. Gazdag and P. Squazerro, Migration of seismic data by phase shift plus interpolation, Geophysics, 48, pp. 124-131, 1984.

11. R.-S. Wu, Wide-angle elastic wave one-way propagation in heterogeneous media and an elastic wave complex-screen method, Journal of Geophysical Research, V. 99, No Bl, pp. 751-766, January 10, 1994.

12. R.-S. Wu and L.-J. Huang, Reflected wave modeling in heterogeneous acoustic media using De Wolf approximation, Mathematical Methods in Geophycal Imaging, edited by Hassenzadeh, SPIE, pp. 176-186, 1995.

13. R.-S. Wu, Accuracy analysis of screen propagators for wave extrapolation using a thin-slab model, SEG 67th Annual Meeting, 1997.

14. X.-B. Xie and R.-S. Wu, 3D elastic wave modeling using the complex screen method, SEG 66th Annual Meeting, 1996.

15. X.-B. Xie and R.-S. Wu, Improve the wide-angle accuracy of screen method under large contrast, SEG 68th Annual Meeting, 1998.

16. S. Jin and R.-S. Wu, Prestack depth migration using a hybrid pseudo-screen propagator, SEG 69th Annual Meeting, 1999.

17. S. Jin, Experimenting with the hybrid pseudo-screen migration, SEG 70th Annual Meeting, 2000.

18. A. Bamberger, B. Engquist, L. Halpern and P. Joly, Parabolic wave equation approximations in heterogeneous media, Siam J. Appl. Math., Vol. 48, No. 1, pp. 99128, 1991.

19. Yu. Zhang, G. Zhang and N. Bleistein, Wave equation migration arising from true amplitude one-way equations, Problems Inverses, 19, pp. 1113-1138, 2003.

20. C. Wapenaar and J. Grimbergen, Reciprocity theorems for one-way wavefields, Geophys. J. Int., 127, pp. 169-177, 1996.

21. M. De Hoop and A. Gautesen, Uniform asymptotic expansion of the square root Helmholtz operator and the one-way wave propagator, Siam. J. Appl. Math., V. 63, No. 3, pp. 777-800, 2003.

22. Петрашень Г.И., Нахамкин С.И. Продолжение волновых полей в задачах сейсморазведки. Л.: Наука, 1973, 170 с.

23. Гольдин С.В. Интегральные продолжения волновых полей. Геология и геофизика, 1985, No. 5, с. 103-113.

24. Гольдин С.В. К общей теории продолжения волновых полей. Новосибирск, 1986. 45 с. (Препринт No. 8, Институт геологии и геофизики СО АН СССР).

25. Tarantola, А., 1984. Inversion of seismic reflection data in the acoustic approximation, Geophysics, V. 49, N. 8, pp. 1259-1266.

26. Beylkin, G., 1985. Imaging of discontinuities in the inverse scattering problem by inversion of a causal generalized Randon transform, J. math. Phys., 26, 99-108.

27. Рыжиков Г.А., Троян B.H. Применение методов вычислительной томографии при решении интерпретационных сейсмических задач // Вестник Ленинградского университета, 1985, N 25.

28. Рыжиков Г.А., Троян В.Н. Определение волнового поля и скорости в среде с использованием метода теории возмущений // Вестник Ленинградского университета, 1986, N 25.

29. Bleistein, N., 1987, On the imaging of the reflectors in the earth, Geophysics, 52, 931-942.

30. Bleistein, N., J. Cohen and J. Stockwell, 2001, Mathematics of multi-dimensional seismic imaging and inversion, Springer, New York, 511 p.

31. ЬатЬагё, G., Virieux, J., Madariaga, R., and Jin, S., 1992. Iterative asymptotic inversion in the acoustic approximation, Geophysics, 57, 1138-1154.

32. Thierry, P., Operto, S. and Lambar6, G., 1999. Fast 2D ray+Born migration/inversion in complex media, Geophysics, 34, 162-182.

33. Backus G., Gilbert F. Numerical applications of a formalism for geophysical inverse problems // Geophys. J. R. Astr. Soc. 1967, Vol. 13.

34. Goldin, S. V., 1987, Dynamic analysis of images in seismics, Geology and Geophysics, Allerton Press, 28, no. 2, 84-93 (репринт, Геология и Геофизика, 28, N 2, с. 90-98)

35. Keho, Т.Н., and Beydoun, В., 1988. Paraxial ray Kirchhoff migration Geophysics, 53, 12, 1540-1546.

36. Schleicher, J., Filpo, E., Hanitzsch, C. Hubral, P. and Tygel, M., 1993, Amplitude-preserving migration using diffraction stacks, 38th. Ann. Int. Mtg., SPIE, Proceedings, Mathematical Methods in Geophysical Imaging, vol. 2003, 97-108.

37. Kiyashchenko, D., Kashtan, В., Plessix, R.-E., 2004, Anisotropic migration weight for amplitude-preserving miration and sensitivity analysis, Geophys. J. Int., V. 157, 753763.

38. Docherty, P., 1991. A brief comparison of some Kirchhoff integral formula for migration and inversion, Geophysics, 56, 1164-1169.

39. A. Ziolkowski, P. Hannsen, R. Gatliff, H. Jacubowicz, A. Dobson, G. Hampson, X.Y. Li and E. Liu, Use of low frequencies for sub-basalt imaging, Geophysical Prospecting, 51, pp. 169-182, 2003.

40. Lailly, P., 1983. The seismic inverse problem as a sequence of before stack migration, in Proc. Conf. on Inverse Scattering, Theory and Application, SIAM, Philadelphia.

41. Chavent, G. & Plessix, R.-E., 1999. An optimal true amplitude least-squares prestack depth migration, Geophysics, 64, 508-515.

42. Shin, C., Jang, S. k. Min, D., 2001. Improved amplitude preservation for prestack depth migration by inverse scattering theory, Geophysical Prospecting, 49, 592-606.

43. Plessix, R.E. & Mulder, W.A., 2004. Frequency-domain finite-difference amplitude-preserving migration, Geophys. J. Int., 157, 975-987.

44. F. Joncour, J. Svay-Lucas, G. Lambare and B. Duquet, True amplitude wave equation migration, Expanded Abstracts of 67th EAGE International Conference and Exhibition, June 2005, Madrid.

45. Sava, P., Fomel, S., Amplitude-preserved common gathers by wavefield continuation methods, Geophysics, Vol. 68, N. 3, pp. 1065-1074, 2003.

46. Lambar<5, G., Virieux, J., Madariaga, R., & Jin, S., 1992. Iterative asymptotic inversion in the acoustic approximation, Geophysics, 57, 1138-1154.

47. Sevink, A., Asymptotic seismic inversion, Ph.D thesis, Delft University of Technology, The Netherlands, 1996.

48. George, A. and Liu, J.W, Computer solution of large sparse positivie definite systems, Prentice-Hall, Englewood Cliffs, New-York, 1981.

49. Erlanga, Y.A., Vuik, C., and Oosterlee, C., 2004, On a class of preconditioned for the Helmholtz equation, Applied Num. Math, 50, 409-425.

50. J.P. Berenger. A perfectly matched layer for the absorption of electromagnetic waves. J. of Сотр. Phys, 114:185-200, 1994.

51. H.H. Калиткин, Численные методы, M., Наука, 1978.

52. George, A. and Liu, J.W, Computer solution of large sparse positivie definite systems, Prentice-Hall, Englewood ClifTs, New-York, 1981.

53. Robert W. Clayton and Bjorn Engquist, Absorbing boundary conditions for wave equation migration, Geophysics, V. 45, N. 5, p. 895-904.

54. Aminzadeh, F., Brae, J. and Kunz, Т., 1997. 3D Salt and Overthrust Model, Soc. of Expl. Geophysicists.

55. F. Collino. Numerical analysis of mathematical models for wave propagation. Report, PSI Consortium, IFP, Rueil-Malmaison, France, 1993.

56. A. ten Kroode, Smit, D. &c Verdel, A., 1994. A microlocal analysis of the migration, Wave Motion, 28, 149-172.

57. Treves, F., 1980. Introduction to Pseudo-differential and Fourier Integral Operators, Plenum Press, New-York.

58. В. Троян, Ю. Киселев, Статистические методы обработки и интерпретации геофизических данных, Издательство СПбГУ, 2000, 578 с.

59. Lambar6, G., Virieux, J., Madariaga, R., & Jin, S., 1992. Iterative asymptotic inversion in the acoustic approximation, Geophysics, 57, 1138-1154.

60. Sevink, A., Asymptotic seismic inversion, Ph.D thesis, Delft University of Technology, The Netherlands.

61. K. Aki and P.G. Richard, Quantitative Seismology, Theory and Methods, Freeman abd company, San Franscico, 1980.

62. Федорюк M.B., Метод перевала, M., Наука, 1977, 368 с.

63. Kiyashchcnko, D., Kashtan, В., 2004: Multi-one-way modeling method, Proceedings of 5-th international conference "Problems of Geocosmos", May 24-28, St. Petersburg, Russia, p. 255

64. Kiyashchenko, D., Plessix, R.-E., Kashtan, B., 2004: Improved amplitude multistep one-way modelling, 66-th EAGE Conference and Exhibition, 7-10 June 2004, Paris, France, p. 4

65. Kiyashchenko, D., Plessix, R.-E., Kashtan, B., Multi-one-way modeling method, Proceedings of 7-th SEGJ International Symposium, 24-26 November 2004, Sendai, Japan, p. 162.

66. Kiyashchenko, D., Plessix, R.-E., Kashtan, B., Troyan, V., Improved amplitude multi-one-way wave equation migration, Expanded Abstracts of 67th EAGE Conference and Exhibition, 13-16 June 2005, Madrid, Spain, p. 4.

67. Kiyashchenko, D., Plessix, R.-E., True amplitude multi-one-way modeling method: comparison of different imaging principles, Expanded Abstracts of 75th SEG Annual Conference and Exhibition, 6-11 November, Houston, USA, p. 4.

68. Kiyashchenko, D., Plessix, R.-E., Kashtan, B., Troyan, V., True amplitude thin slab multi-one-way inversion application to marine-type seismic data, Expanded Abstracts of 68th EAGE Conference and Exhibition, 12-15 June 2006, Vienna, Austria, p. 4.

69. Kiyashchenko, D., Plessix, R.-E., Kashtan, B., Troyan, 2005, Improved amplitude multi-one-way modeling method, Wave Motion, 43, p. 99-115.