Математические и вычислительные подходы к повышению качества сейсмических изображений на основе моделирования упругих волновых полей тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат физико-математических наук Байдин, Василий Григорьевич

• Байдин, Василий Григорьевич
• кандидат физико-математических наук
• 2013

Оглавление диссертации кандидат физико-математических наук Байдин, Василий Григорьевич

Содержание

Введение

Глава 1. Обзор технологий и методов сейсмической миграции

1.1. Понятие сейсмической миграции

1.1.1. Геометрическая сейсмика

1.1.2. Построение эффективной модели и миграция — две стадии реконструкции среды

1.2. Многообразие методов миграции

1.2.1. Структура среды

1.2.2. Параметризация модели среды

1.2.3. Модели распространения колебаний

1.2.4. Типы систем наблюдения

1.3. Миграция в обратном времени

1.4. Постобработка

1.4.1. Нормализация

1.4.2. Высокочастотная фильтрация High-pass

1.5. Финальное описание алгоритма миграции в обратном времени с использованием экстраполяционного подхода

1.6. Итоги главы

Глава 2. Математические модели для задачи миграции

2.1. Уравнения распространения сейсмических колебаний

2.1.1. Общий вид линейно-упругой среды

2.1.2. Параметризация Фойгта

2.1.3. Акустическая модель среды

2.1.4. Изотропная упругая модель среды

2.1.5. Оценка направления движения волны. Вектор Умова-Пойнтинга

2.2. Математическая модель системы наблюдения

2.3. Разделение упругих волновых полей на продольные и поперечные волны

2.3.1. Разделения поля с помощью дифференциальных операторов

2.3.2. Разделение поля с помощью вектора Умова-Пойнтинга

Глава 3. Оптимизационный подход к построению метода RTM

3.1. Акустическая миграция

3.1.1. Постановка задачи акустической миграции

3.1.2. Обратная динамическая задача сейсмики

3.1.3. Оценка градиента функционала

3.2. Упругая миграция в обратном времени

3.2.1. Прямая задача

3.2.2. Постановка обратной задачи

3.2.3. Оценка градиента функционала

3.2.4. Выражение градиента для различных параметризаций модели среды

3.3. Построение изображений на основе градиента

3.3.1. От градиента по параметру — к изображению

3.4. Описание алгоритма миграции в обратном времени

3.5. Итоги главы

Глава 4. Классификация артефактов миграции

4.1. Артефакты многократных отражений

4.2. Апертурные артефакты (улыбки миграции)

4.3. Артефакты симметрии

4.4. Артефакты конвертации

4.5. Артефакты двойственности

4.6. Итоги главы

Глава 5. Повышение качества сейсмических изображений

5.1. Миграция на основе дивергентно-роторного функционала

5.1.1. Оценка градиентов дивергентного и роторного функционала

5.1.2. Применение в практических задачах

5.2. Влияние системы наблюдения на качество изображения

5.3. Подавление артефактов двойственности

5.3.1. Подавление с помощью дивергентно-роторного функционала

5.3.2. Другие способы подавления артефактов двойственности

5.4. Фильтр по углу наклона границы (Dip-фильтр)

5.4.1. Акустический случай

5.4.2. Изотропный упругий случай

5.4.3. Алгоритм Dip-фильтра

5.5. Итоги главы

Глава 6. Модельные расчёты

6.1. Акустическая модель EAGE ВР2004 Benchmark

6.2. Наклонная акустическая модель

6.3. Упругая модель с плоской отражающей границей

6.4. Трёхслойная упругая модель

6.5. Трёхслойная модель с включением

6.6. Модель «Покрывающая толща»

6.6.1. Расчёт изображений

6.7. Модель Marmousi-mod

6.7.1. Расчёт изображений

6.8. Итоги главы

Заключение

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

Приложение А. Экстраполяция с помощью граничных условий

Приложение Б. Эквивалентность экстраполяционного подхода и оптимизационного подхода с некоторым функционалом

Приложение В. Программная реализация

В.1. Аппаратные и программные требования

В.2. Код Ыау1ег-5

В.З. Алгоритм упругой миграции

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

Введение

Актуальность темы исследования

В задачах сейсморазведки математическое моделирование играет определяющую роль, поскольку позволяет «заглянуть» вглубь земной толщи, представляющей собой «чёрный ящик», рассеивающий и отражающий приходящие от искусственных источников волны (частоты таких колебаний обычно находятся в диапазоне 0.1Гц - 100Гц). Численное моделирование распространения упругих волн в среде позволяет сформулировать обратную задачу: восстановление упругих параметров среды по отклику сейсмических колебаний, записанных на сейсмоприёмниках. Исторически обратная задача сейсморазведки включает в себя две основные подзадачи: 1) определение макроскоростной модели среды — низкочастотной по пространству «трендовой» составляющей и 2) восстановление резких перепадов упругих параметров среды — интерфейсов, обычно на границах слоёв. Для обратной задачи в привычном математическом смысле используется термин полноволновое обращение. Подзадача 2) носит название сейсмическая миграция или построение сейсмических изображений отражающих границ.

Качество изображений сейсмической миграции зависит от сложности структуры исследуемой среды и от степени детализации приближённой модели земной толщи, полученной после решения подзадачи 1). Классические методы миграции, основанные на приближении геометрической сейсмики, такие как ОГТ (общая глубинная точка), Кирхгофа, лучевая миграция, не способны адекватно работать на моделях со сложным геологическим строением. Благодаря существенному росту вычислительных мощностей стало возможным решать задачи миграции на основе прямого моделирования волновых полей. В данных методах структура модели среды может быть уже достаточно сложной и более близкой к реальности. В диссертации рассмотрен метод миграции в обратном времени (RTM — reverse time migration), основанный на неполной итерации обратной динамической задачи упругости. Этот метод даёт возможность построения корректных изображений в средах со сложной структурой и, как следствие, является наиболее ресурсоёмким из известных методов миграции.

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

Методология и методы исследования. Математическое моделирование, числен-

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

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

2. Сам по себе метод ЯТМ основан на решении большого числа прямых задач, моделирующих распространение волн в среде. Его создание связано с передовыми вычислительными технологиями.

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

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

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

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

3. Разработка и анализ ряда фильтров, подавляющих различные типы артефактов в процессе построения изображения.

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

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

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

1. Произведена модификация задачи упругой миграции в обратном времени для скважинной сейсмики путём введения дивергентно-роторного функционала.

2. Разработан алгоритм решения задачи упругой миграции на основе предложенного дивергентно-роторного функционала.

3. Проведен анализ причин возникновения ряда известных и новых типов артефактов, дана их классификация.

4. Предложен dip-фильтр по углу наклона границы для уничтожения ряда артефактов (результат получен в соавторстве и на него подана и опубликована заявка на патент [80]).

5. Создан комплекс программ, позволяющий осуществлять миграцию и получать различные изображения на основе РР, PS, SP, SS событий с гибкой настройкой фильтрации, нормализации, суммирования по источникам и т.п.

Теоретическая и практическая значимость.

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

Реализованный метод акустической и упругой миграции в двумерной постановке входит в комплекс программ миграции в обратном времени (Reverse Time Migration -RTM), разработанный совместно с Л.Е. Довгиловичем и И.Л. Софроновым. Комплекс применяется в Московском научно-исследовательском центре «Шлюмберже». Положения, выносимые на защиту:

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

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

3. Комплекс программ по акустической и упругой миграции в обратном времени.

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

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

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

Апробация результатов

Основные результаты диссертации докладывались на следующих научных конференциях и семинарах:

1. 53-я, 54-я, 55-я конференции МФТИ, Долгопрудный, 2010 - 2012;

2. XIV молодёжная конференция-школа с международным участием «Современные проблемы математического моделирования», Абрау-Дюрсо, 2011;

3. 74-я конференция и выставка EAGE (European Association of Geoscientists & Engineers), Копенгаген, 2012;

4. XIV конференция EAGE «Геомодель-2012», Геленджик, 2012;

5. XII международная конференция «Гальперинские чтения-2012», Москва, 2012;

6. Семинар кафедры сейсмометрии и геоакустики МГУ, Москва, 2013;

7. Семинар Института физики Земли им. О.Ю.Шмидта РАН, Москва, 2013;

8. Семинар кафедры информатики МФТИ, Долгопрудный, 2013;

9. Научные семинары Московского научно-исследовательского центра «Шлюмбер-же», Москва, 2010 - 2013;

10. Объединённый семинар сейсмических лабораторий Института нефтегазовой геологии и геофизики им. А.А.Трофимука СО РАН, Новосибирск, 2013.

Публикации и личный вклад

Результаты диссертации опубликованы в 8 печатных работах, из них одна [53] в издании из списка, рекомендованного ВАК РФ, подана и опубликована заявка на патент [80].

Все научные результаты, вынесенные на защиту, получены лично автором.

Заключение диссертации по теме «Математическое моделирование, численные методы и комплексы программ», Байдин, Василий Григорьевич

Основные результаты работы

1. Проведено обобщение различных интерпретаций метода 11ТМ и его конкретизация для уравнений акустики и упругости на основе экстраполяционного и оптимизационного подходов.

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

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

4. На основе вектора Умова-Пойнтинга предложен и реализован сНр-фильтр, работающий как в акустическом, так и в изотропном упругом приближении.

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

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

Благодарности

Автор выражает искреннюю благодарность своим старшим коллегам: Софронову Ивану Львовичу за научное руководство, постоянную поддержку, многочисленные советы и помощь в подготовке защиты диссертации;

Шевченко Алексею Александровичу за плодотворное обсуждение работы и консультации;

Бородину Игорю Вячеславовичу за возможность пройти стажировку и обмен опытом;

Чеверде Владимиру Альбертовичу за поддержку и содействие, возможность выступить на семинаре и плодотворные замечания;

Байдину Григорию Васильевичу за обсуждение и конструктивную критику; а также Довгиловичу Леониду Евгеньевичу и Сотрудникам Московского научно-исследовательского центра Шлюмберже за помощь в работе и поддержку.

Заключение

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

Глава 3 посвящена описанию возникающих артефактов и способам их подавления, как общеизвестных, так и не рассматривающихся в применении к скважинной миграции в обратном времени.

В главе 6 рассматриваются различные примеры работы алгоритма миграции и производится анализ разных режимов работы алгоритма. Рассмотрены следующие модели:

1. EAGE ВР2004 Benchmark. Общепринятая реалистичная акустическая модель для отладки алгоритма миграции и проверке программного кода на работу в параллельном режиме.

2. Наклонная акустическая модель. Тестовый расчёт для проверки работоспособности фильтра по углу наклона границы. Результаты опубликованы на ежегодной международной конференции EAGE [31].

3. Двуслойная упругая модель с горизонтальной границей. Простейший пример для иллюстрации образования и фильтрации артефактов двойственности.

4. Трёхслойная упругая модель. Очевидный пример преимущества упругой миграции перед акустической. Поскольку источник рассматривался типа центра расширения, все поперечные волны, возникающие в среде, — конвертированные. Упругая миграция позволяет корректно построить изображения конвертированных волн, на что акустической миграции требуется решить несколько прямых задач, со скоростью поперечной волны в некоторых участках модели.

5. Трёхслойная модель со включением. Пример, в котором дивергентный и роторный функционалы используются по отдельности для «улавливания» отражений определённого типа.

6. Модель «Покрывающая толща». Моделирование той стадии реконструкции скоростной модели среды, когда сложная структура над скважиной уже определена и необходимо найти положение отражающей границы, находящейся глубже. Расчёт изображений данной модели показал, что у четырёх типов изображений (РР, PS, SP,

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

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

Список литературы диссертационного исследования кандидат физико-математических наук Байдин, Василий Григорьевич, 2013 год

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

1. Aarre Victor. Globally consistent dip estimation// 2010 SEG Annual Meeting. - 2010.

2. Baysal Edip, Kosloff Dan D., Sherwood John W. C. Reverse time migration//Geophysics. — 1983. - November. - Vol. 48, no. 11. - P. 1514-1524.

3. Billette Frederic, Brandsberg-Dahl Sverre. The 2004 BP velocity benchmark. // 67th Annual Internat. Mtg., EAGE, Expanded Abstracts. - EAGE, 2005. - P. B035.

4. Chang Wen-Fong, McMechan George A. Reverse-time migration of offset vertical seismic profiling data using the excitation-time imaging condition // Geophysics. — 1986. — January. — Vol. 51, no. 1. — P. 67.

5. Chapman C. Fundamentals of Seismic Wave Propagation. — Cambridge University Press, 2004.

6. Clapp Robert G. Reverse time migration with random boundaries // SEG Technical Program Expanded Abstracts 2009. - P. 2809-2813.

7. Edme P., Muyzert E. Rotational Data Measurement // 75th EAGE Conference & Exhibition. — 2013. — June.

8. Etgen John T.Prestacked migration of P and Sv-waves//1988 SEG Annual Meeting. — 1988.

9. Fichtner A, Bunge H-P, Igel H. The adjoint method in seismology: I. Theory // Physics of the Earth and Planetary Interiors. - 2006. - Vol. 157, no. 1. - P. 86-104.

10. Fletcher Robin P, Du Xiang, Fowler Paul J. Reverse time migration in tilted transversely isotropic (TTI) media// Geophysics. - 2009. - Vol. 74, no. 6. - P. WCA179-WCA187.

11. Gray Samuel H. True-amplitude seismic migration: A comparison of three approaches // Geophysics. - 1997. - Vol. 62, no. 3. - P. 929-936.

12. Gray Samuel H, Etgen John, Dellinger Joe, Whitmore Dan. Seismic migration problems and solutions //Geophysics. - 2001. - Vol. 66, no. 5. - P. 1622-1640.

13. Gray Samuel H, May William P. Kirchhoff migration using eikonal equation traveltimes // Geophysics. - 1994. - Vol. 59, no. 5. - P. 810-817.

14. Greenwood Andrew, Dupuis Christian J, Urosevic Milovan, Kepic Anton. Hydrophone VSP surveys in hard rock// Geophysics. - 2012. - Vol. 77, no. 5. - P. WC223-WC234.

15. Horn Berthold KP, Schunck Brian G. Determining optical flow // Artificial intelligence. — 1981. - Vol. 17, no. 1. - P. 185-203.

16. Kaelin Bruno, Guitton Antoine. Imaging condition for reverse time migration // 2006 SEG Annual Meeting. - 2006.

17. Klockner Andreas, Pinto Nicolas, Lee Yunsup et al. PyCUDA and PyOpenCL: A Scripting-Based Approach to GPU Run-Time Code Generation // Parallel Computing. — 2012. — Vol. 38, no. 3.- P. 157-174.

18. Liu Hong-Wei, Li Bo, Liu Hong et al. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation // The Chinese Journal of Geophysics. — 2010.- P. 1725-1733.

19. Loewenthal D, Lu L, Roberson R, Sherwood J. The Wave Equation Applied to MIGRATION // Geophysical Prospecting. - 1976. - Vol. 24, no. 2. - P. 380-399.

20. Lu R., Traynin P., Anderson J.E., Dickens T. Impact of TTI Anisotropy on Elastic and Acoustic Reverse Time Migration// 72nd EAGE Conference & Exhibition. — 2010. — June.

21. Lu Rongrong, Traynin Peter, Anderson John E. Comparison of elastic and acoustic reverse-time migration on the synthetic elastic Marmousi-II OBC dataset // 79th Annual International Meeting, SEG, Expanded Abstracts. - 2009. - P. 2799-2803.

22. Lu Rongrong, Willis Mark E., Mateeva Albena et al. Imaging a salt dome flank by directional redatuming of a field 3-D VSP survey// SEG Technical Program Expanded Abstracts 2009. — P. 4085-4089.

23. Ma Zhixia, Sun Sam Zandong, Bai Haijun, Sun Wenbo. 3C-3D VSP PP and PSV imaging in the Lungu 38 area, Tarim Basin // 2009 SEG Annual Meeting. - 2009.

24. Martin Gary S, Wiley Robert, Marfurt Kurt J. Marmousi2: An elastic upgrade for Marmousi // The Leading Edge. - 2006. - Vol. 25, no. 2. - P. 156-166.

25. Mittet Rune, Helgesen Jan. Elastic iterative migration of offset VSP data//1992 SEG Annual Meeting. - 1992. - P. 1-4.

26. Mittet Rune, Hokstad Ketil, Helgesen Jan, Canadas Guy. Imaging of offset VSP data with an elastic iterative migration scheme // Geophysical prospecting. — 1997. — Vol. 45, no. 2. — P. 247-267.

27. Mora Peter. Elastic wave-field inversion of reflection and transmission data // Geophysics. - 1988. - Vol. 53, no. 6. - P. 750-759.

28. Mulder WA, Plessix RE. One-way and two-way wave-equation migration // 2003 SEG Annual Meeting. - 2003.

29. Pavloske Ken, Beyoglu Jana, Eshete Tefera, Baldock Simon. Pre-processing Improvements For Enhanced PSTM And TTI PSDM Images In the Santos Basin: A Case History // 2011 SEG Annual Meeting. - 2011.

30. Plessix R-E. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications // Geophysical Journal International. — 2006. — Vol. 167, no. 2. - P. 495-503.

31. Sofronov I.L., Baidin V., Borodin I. et al. Model Dip-guided Imaging Condition for Reversetime Migration in VSP Scenario // 74th EAGE Conference & Exhibition. - 2012. -June. — P. 1-5.— URL: http://www.earthdoc.org/publication/publicationdetails/ ?publication=58908.

32. Song Xiaolei. Application of Fourier Finite Differences and Lowrank Approximation Method for Seismic Modeling and Subsalt Imaging : Ph. D. thesis / Xiaolei Song ; THE UNIVERSITY OF TEXAS AT AUSTIN. - 2012.

33. Sun J, Zhang Y. Practical issues of reverse time migration: true amplitude gathers, noise removal and harmonic-source encoding // SEG Extended Abstracts. — 2009. — Vol. 2009, no. 1,- P. 1-5.

34. Sun Robert, McMechan George A. Scalar reverse-time depth migration of prestack elastic seismic data//Geophysics. - 2001. - Vol. 66, no. 5. - P. 1519-1527.

35. Sun Robert, McMechan George A. Pre-stack reverse-time migration for elastic waves with application to synthetic offset vertical seismic profiles // Proceedings of the IEEE. — 1986. - Vol. 74, no. 3. - P. 457-465.

36. Sun Robert, McMechan George A, Hsiao Hsu-Hong, Chow Jinder. Separating P-and S-waves in prestack 3D elastic seismograms using divergence and curl // Geophysics. — 2004. — Vol. 69, no. 1,- P. 286-297.

37. Symes William W. Reverse time migration with optimal checkpointing // Geophysics. — 2007. - Vol. 72, no. 5. - P. SM213-SM221.

38. Tarantola Albert. Theoretical background for the inversion of seismic waveforms including elasticity and attenuation // Pure and Applied Geophysics. — 1988. — March. — Vol. 128, no. 1-2,- P. 365-399.

39. ThomsenL. Weak elastic anisotropy//GEOPHYSICS. - 1986. - Vol. 51, no. 10. - P. 19541966.

40. Verschuur Dirk J, Berkhout AJ, Wapenaar CPA. Adaptive surface-related multiple elimination//Geophysics. - 1992. - Vol. 57, no. 9. - P. 1166-1177.

41. Whitmore ND. Iterative depth migration by backward time propagation//1983 SEG Annual Meeting. - 1983.

42. Xiao Xiang, Leaney W Scott. Local vertical seismic profiling (VSP) elastic reverse-time migration and migration resolution: Salt-flank imaging with transmitted P-to-S waves // Geophysics. - 2010. - Vol. 75, no. 2. - P. S35-S49.

43. Yan Jia, Sava Paul. Elastic wave-mode separation for VTI media // Geophysics. — 2009. — Vol. 74, no. 5. - P. WB19-WB32.

44. Yan Jia, Sava Paul. Isotropic angle-domain elastic reverse-time migration // Geophysics. — 2008. - Vol. 73, no. 6. - P. S229-S239.

45. Yan Jia, Sava Paul. Elastic wavefield imaging with scalar and vector potentials // 2007 SEG Annual Meeting. - 2007.

46. Yoon Kwangjin, Marfurt K.J. Reverse-time migration using the Poynting vector // Exploration Geophysics. - 2006. - Vol. 37, no. 1. - P. 102-107.

47. Zhdanov M.S. Geophysical Inverse Theory and Regularization Problems. Methods in Geochemistry and Geophysics. - Elsevier Science, 2002. - ISBN: 9780080532509.

48. Zhe Jingping, Greenhalgh Stewart A. Prestack multicomponent migration // Geophysics. — 1997. - Vol. 62, no. 2. - P. 598-613.

49. Zhu Jinming, Lines Larry, Gray Sam. Smiles and frowns in migration/velocity analysis // Geophysics. - 1998. - Vol. 63, no. 4. - P. 1200-1209.

50. А.Г. Куликовский, Е.И. Свешникова. Нелинейные волны в упругих средах. — Калининград, Янтарный сказ, 1998.

51. Авербух А.Г. Изучение состава и свойств горных пород при сейсморазведке. — М., Недра, 1982.

52. Аки К., Ричарде П.Г. Количественная сейсмология: теория и методы. — М., Мир, 1983.

53. Байдин В.Г. Построение изображения сейсмического разреза по модели ВР2004 Benchmark//Труды МФТИ. - 2013. - Т. 5, № 2. - С. 150-159.

54. Байдин В. Г. Фильтрация изображений с помощью вектора Умова-Пойнтинга для задач сейсмической миграции в обратном времени // Труды XIV конференции EAGE «Геомодель-2012». — 2012.— С. 1-5.— URL: http://www.earthdoc.org/publication/ publicationdetails/?publication=62397.

55. Байдин В. Г. Применение вектора Умова-Пойнтинга для разделения полей в задачах упругой миграции в обратном времени // Труды 55-й научной конференции МФТИ. Аэрофизика и космические исследования. — Т. 2. — М.: МФТИ, 2012. — С. 83-84.

56. Байдин В. Г. Проблемы качества сейсмической миграции в обратном времени на примере модели ВР2004 Benchmark // Труды XIV молодёжной конференции-школы с международным участием «Современные проблемы математического моделирования». — Ростов-на-Дону: изд-во Южного федерального университета, 2011. — С. 5155.

57. Байдин В. Г. Проблемы качества сейсмической миграции в обратном времени на примере модели ВР2004 // Труды 54-й научной конференции МФТИ. Управление и прикладная математика. — Т. 2. - М.: МФТИ, 2011. - С. 16-18.

58. Байдин В. Г. Методы улучшения сейсмических изображений. Разработка и реализация на высокопроизводительных вычислительных системах // Труды 53-й научной конференции МФТИ. Часть VII. Управление и прикладная математика. — М.: МФТИ, 2010. - С. 7-10.

59. Владимиров B.C. Уравнения математической физики. — М. Наука, 1971.

60. Гурвич И.И, Боганик Г.Н. Сейсмическая разведка. — М., Недра, 1980.

61. Гурвич Н.И., Номоконов В.П. Сейсморазведка. Справочник геофизика. — М., Недра, 1981.

62. Довгилович JI.E. Анализ производительности алгоритма решении обратной задачи сейсморазведки методом миграции в обратном времени с применением компактных схем // Труды 53-й научной конференции МФТИ. — 2010.

63. Довгилович JI.E, Софронов И.Л. Анализ явных и неявных центрально-разностных операторов для вычисления второй производной на равномерных сетках // Труды МФТИ. - 2013. - Т. 5, № 2. - С. 175-182.

64. Е.А. Козлов. Миграционные преобразования в сейсморазведке. — М., Недра, 1986.

65. Клаербоут Д.Ф. Сейсмическое изображение земных недр. — М., Недра, 1989.

66. Кочин Н.Е. Векторное исчисление и начала тензорного исчисления. — М., Наука, 1965.

67. Кравцов Ю.А., Орлов Ю.И. Геометрическая оптика в неоднородных средах. — М., Наука, 1980.

68. Кутовенко М. П., Протасов М. И., Чеверда В. А. Использование Гауссовских пучков для построения сейсмических изображений в истинных амплитудах по многокомпонентным данным//Технологии сейсморазведки. — 2010. — по. 4.

69. Лехницкий С.Г. Теория упругости анизотропного тела. — М., Наука, 1977.

70. Ляв А. Математическая теория упругости. — М., НКТП СССР, 1935.

71. Неклюдов Д.А., Бородин И. Миграция данных непродольного ВСП для построения глубинных изображений среды с неизвестной верхней частью разреза в районах со сложным геологическим строением // Seismic Technology. — 2008. — Vol. 5, no. 4. — P. 6-14.

72. Сильвестров И. Ю., Неклюдов Д. А. Многокомпонентная миграция данных НВСП по методу наименьших квадратов с подавлением артефактов//Технологии сейсморазведки. - 2008. - по. 4. - Р. 15-24.

73. Умов H.A. Уравнения движения энергии в телах. — Одесса, 1874.

74. Урупов А.К. Основы трехмерной сейсморазведки. — М., ФГУП Изд-во «Нефть и газ», 2004.

75. Федоренко Р.П. Приближенное решение задач оптимального управления. — М., Наука, 1978.

76. Черняк В. Г., Суетин П.Е. Механика сплошных сред. - М., ФИЗМАТЛИТ, 2006.

77. Шевченко А. А. Сейсмические исследования в скважинах. Учебное пособие. — М.: Издательство МГУ, 2007.

78. Шерифф Р., Гелдарт Л. Сейсморазведка. История, теория и получение данных. — М., Мир, 1987.- Т. 1.

79. Шерифф Р., Гелдарт Л. Сейсморазведка. Обработка и интерпретация данных. — М., Мир, 1987. - Т. 2.

80. Baydin Vasily, Dovgilovich Leonid, Lin Kui et al. Reverse time migration model dip-guided imaging. Appl. No.: 13/725,154; Pub. No.: US 2013/0182538 Al; заявл. 21.12.2012; опубл. 18.07.2013, 14c.

Приложение А. Экстраполяция с помощью граничных

условий

Рисунок А. 1. «Опоясывающая» система наблюдения. £2 — исследуемая область. Л — поверхность приёмников, I) — внутренность поверхности Я, Я — область локализации источника.

Рассмотрим «опоясывающую» систему наблюдения (рис. А.1), когда сейсмоприём-ники равномерно и плотно расположены на поверхности Я, охватывающей исследуемую

■-у

область О,. При этом задача ставится на всём пространстве К . Внутренность поверхности Я обозначим как I). Скорость звука в истинной среде обозначим V, в приближенной эффективной среде скорость звука обозначим Vа. Будем считать, что вне исследуемой области £1 скорость V известна точно, то есть:

¿V ={уг-Уа) =0.

о\п п\п

Истинные колебания в среде задаются следующим уравнением:

(угу

х е к2, г е [о, г],

рг (х,г = 0) = о,

(А.1)

где / (х, /) — локализованный в пространстве (х е 5) источник сейсмических колебаний, действующий в течение времени ? (Е ^0, Г1^, то есть:

/М

/М

= 0,

?е[г5/г]

= 0.

Запишем значение давления на сейсмоприёмниках:

Предположим, что на момент времени I = Т все колебания рг(х, Т) ушли за пределы области £>, то есть:

= 0.

/ (х,Г)

Сформулируем следующую краевую задачу:

хея

1 д

- Арг = 0,

(у) "

х е £>+, г е [0, г], рГ (х,г = г) = о,

(А.2)

/ М

= 5й (х,г).

хея

Заметим, что вместо начальных условий здесь конечные для ? = Т. Очевидно, что решение рг данной задачи и решение рг задачи (А.1) совпадают на области: х е / е

(г5, г).

Теперь вместо истинной скорости звука V возьмём приближённую скорость Vа, отличающуюся от истинной только в области £1. Аналогично задаче (А.2) сформулируем

следующую задачу:

—Ц дира - Ара = О,

(V)2 пи и

х£])+ле [О, г],

ра (х,Г = Г) = О,

(А.З)

Ра М

= 9> (х,/) .

хел

~г

Рассмотрим отклонение истинного поля рг от приближённого ра. Запишем

рг = ра + 8р,

Подставляем данные соотношения в уравнение из задачи(А.2):

( + 3" ^ + 5Р)~А ^ + 6Р) =

и вычитаем из полученного выражения уравнение из задачи (А.З), получаем:

(у;)2 д„8р-Ь5р = -6±д„ (рг).

Таким образом, отклонение 8р подчиняется волновому уравнению с некоторой правой частью, которая не равна нулю только в области 8Ф 0 (то есть, в области £2) и, соответственно, только когда рг в области £2 не равен нулю. Получается, что до того момента, как волны, сгенерированные на границе вследствие граничного условия в задаче (А.2), не придут в область П, они эквивалентны волнам, выпущенным в задаче (А.З) с такими же граничными условиями.

Приложение Б. Эквивалентность экстраполяционного подхода и оптимизационного подхода с некоторым

функционалом

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

Дана «опоясывающая» система наблюдения (см. приложение А). Линия приёмников Я охватывает всю исследуемую область, источник тоже находится внутри. Сейсмограммы представляют из себя следующее поле:

41 (х,г) = иг (х,г)

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

Ьаи; (х,г) =0, хед/е [о, т],

= (иг-иа) (х,г)

IV

я

(Б.1)

хек

и; (г = г) = 0, и>, (г = т) = 0.

Сформулируем следующие требования, которые нам понадобятся для математического вывода:

1. Поле, записанное на приёмниках, представляет собой волны, идущие изнутри области Б и не возвращающиеся обратно.

2. В момент времени Т и далее значения поля и> и его временной производной и>{ равны нулю.

3. Подстановка в качестве граничных условий в уравнение (Б.1) означает экстраполяцию волн с приёмников обратно в среду I) в том смысле, что существует некоторый момент времени 0 < < Т, когда носитель функции IV лежит внутри области Б, то есть

V* е ЯЫ < 7*,и;(х,г) = 0.

При соблюдении данных условий можно сформулировать эквивалентную уравнению

(Б.1) задачу Коши:

Ьаи>* (х,г) = О,

[г*, -Ноо] ,

(х,1 = г*) = IV [х, г*),

(Б.2)

IV* (х,Г = г*) = н;, (х,?*),

тогда в области В и для Г е поля и; тл и?* будут равны. Доопределим поле ш*

нулём при ? < полученное поле назовём м;^. Тогда поле и^ удовлетворяет следующему уравнению в неограниченной области (см. [59]):

Ьаи>* (х, г)= р(и> (х, Г) 8' (г - Г) + 8 (г - Г) и>{ (х, г*)),

(Б.З)

Найдём такую функцию / (х,?), сосредоточенную на линии сейсмоприёмников Я, то есть Ух £ .КУ/ : / (х, = 0, что поле ¿V является решением следующего уравнения в

'у

неограниченном пространстве К :

Ьай; (х,г) = / (х, г), хеи2,! е [-оо, Г],

й> (х,Г = т) = О, й>! (х, Г = т) = О,

(Б.4)

причём и?* = й; при х £ I), г £ (?*, Т). Для этого воспользуемся представлением матрицы Грина оператора Ь, и свойством двойственности матрицы Грина (здесь и в дальнейшем для простоты мы опускаем индекс •а у оператора Ь)\

ьх,г°и (*> т> х, £) = т, х, £) = <5 (Г - т) <5 (х - £) 81у

Тогда можно записать:

м)((х,г)= У J (£,т) х,£)

(Б.5)

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

ций ф, х выполняется:

I I (фЬХ - хЩ ахй1 =

лес г=г*

/ дх дф4

хеи

¿х + 11 сик1 [хАФк - ФАхк]

? г=1* л

Воспользовавшись этим тождеством, мы можем переписать выражение (Б.5):

+ I Лм^.т)

дви(^т,х,£) ди>и(£,т) --Оч (г,т,х,£)---

дт

гео

+ I I С]тк1

т=г*

°л (*>т'х' £) д1Сдк г) - (I.т) д^ь (г, т, х,

¿¿■„¿т. (Б. 6)

Первый интеграл равен нулю, поскольку = 0 внутри области В. Во-вторых, в момент времени Т значение поля й>, как и его временной производной, равно нулю. В-третьих, в подынтегральных выражениях в правой части функция Ни берется внутри области х е Б, X е где Ни = и>* = и>. Таким образом, равенство (Б.6) упрощается:

(о, (х,г) = - у р

т

ш

дв (/,т,х,£)

^-К-

дш (£,т)

-оч (/, £,*,£)—-—

г=?

т=1

+

1]тк1

°л (*>т' I) д1юк (£> г) ~ и)] (%> Т) дрь ('> Т> х> £)

й8тйх.

(0

Подействовав оператором Ьх { на данное выражение, мы получим:

Ыи1 (х, ?) =

-I

гей

ш

(М

дЬх ,Сч (г, г,х,£)

дт

т=1*

\

дг т=Г,

+

lJmkl

т=1*

/

дwJ (£,т)

с1Бт(1г =

дт

Г

+ I I CJmkl т=1*

с18тс1т =

8 (Г- т)5(х-£) <5„] д,шк (£, г) - ^ (£, г) д, [5 (г-т)5(х- «5, = „ ^ (х, Г*) (/ - Г) + 5 (г - Г*) +

£&Я

То есть, поле й?1 удовлетворяет волновому оператору Навье с такой правой частью. А поскольку и> (ух Е I), V? > Т) = Оиш( (\/х Е > Т) = 0, то решение задачи Коши (Б.4) действительно совпадает с«;в области х Е I), г Е (/*,Т).

Первое слагаемое в этой сумме — источник, действующий в момент времени г*, а так как мы рассматриваем задачу с конца по времени, то он не влияет на поле при I Е (/*, Т). Таким образом, чтобы промоделировать уравнение (Б.1) в неограниченной области, используя правую часть вместо граничных условий, необходимо в качестве источника в обратном времени использовать следующую функцию:

А = } тк1 {я (* - £) д1шк м 5и - (£.0 д,8 (х -

или, в терминах точечных источников, в каждой точке £ Е Я находится точечный источник:

Р1(х,г)= с1]к1д1ьик-п] <5 (х - £) - С,тии;,(£^)-пт д/(х - , (Б.7)

где п — нормаль к Я в данной точке. Данные источники приведены для уравнения упру-

гости в форме Навье. Переведём их в уравнение (2.5):

рш, - = /„

°1] ~ Сик1д1Шк =

при этом

Сцк1д1Шк ■ П]

б(х-£)= 5ип/ {х - £),

где эквивалентен тензору упругости, записанному в сейсмоприёмниках, для источника напряжений:

ёи = -счк1Щпк5 (х - £).

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

Фх= I I

(Б.8)

х€Л 1=0

где п, — нормаль к контуру Я. Вариация данного функционала будет следующей:

5ФХ= 11 („; - „;) „, + - ^) 5и«„,

хет=о

(1х(1{ =

= / /[к («?-■;)».к».

хев. ?=0

¿хй1 -

1

= / / + (<7'-/) 5и°П'

хея г=о

где и = иа — иг, с = оа - ог. Дальнейшие выводы аналогичны проведённым в параграфе 3.2.3, в результате уравнение в сопряжённой задача для данного функционала будет следующей:

Ьаьи, =

<Х,уИ,<5 (х - £) - Сци^П&б (х - £)

{ЕЯ

Заметим, что правая часть в данном уравнении эквивалентна правой части (Б.7). Таким образом, экстраполяционный подход эквивалентен оптимизационному с функционалом (Б.8). Естественно, не стоит забывать, что мы показали это с учётом серьёзных

ограничений на систему наблюдения и модель среды.

Приложение В. Программная реализация

В данной главе описан программный комплекс, созданный для решения задачи упругой миграции.

В.1. Аппаратные и программные требования

Задачи миграции требуют колоссальных вычислительных ресурсов. Типичные размеры областей для проведения сейсморазведки — сотни квадратных километров по поверхности и несколько километров вглубь. Типичная задача ВСП чуть меньше — десятки квадратных километров на поверхности. А поскольку для проведения миграции нам необходимо хранить в памяти всё поле смещений на всём отрезке времени, очевидно, «бутылочным горлышком» вычислительной системы становится файловое хранилище.

Оценим размер занимаемой памяти для типичной трёхмерной задачи упругой миграции ВСП [22]. Пусть размер области равен 12х 12x4 км , предполагаемое время расчёта 10 секунд, центральная частота импульса сигнала — 40 Гц. Исходя из предположения, что необходимо сохранять не менее 5 точек на период, количество записанных шагов по времени будет равно 10 • 40 • 5 = 2000. Типичную скорость поперечной волны в среде, где проводится миграция, оценим 1500 м/с. Для разрешения 5 точек на длину волны шаг сетки должен быть равен 1500/40/5=7.5 м. Тогда количество точек в пространстве и вре-

о 10

мени будет равняться 12000 • 12000 • 4000/7.5 • 2000 и 2.7 • 10 . Для хранения смещений поля на каждую точку нам необходимо содержать три числа с плавающей точкой. Минимальный объём памяти получается равным 2.7 • 1012 • 3 • 8 = 7.2 • 1013, то есть около 75 Т Байт. Можно существенно уменьшить задействованный объём дисковой памяти, используя различные приёмы, описанные, например, в статьях [37], [6]. В любом случае приходится делать ставку на распределённые вычислительные системы, обладающие существенной локальной памятью.

В первую очередь требования к аппаратным ресурсам диктуются кодом для расчёта волнового уравнения упругости. В данной работе использовался код NAVIER-S (см. следующий параграф).

Для тестирования предложенных в диссертации подходов миграции был написан программный код на языке Python, с использованием библиотек PyCuda [17], NumPy, SciPy, Matplotlib, PyOt.

B.2. KoANavier-S

Для расчёта распространения упругих колебаний согласно уравнению Навье (2.6) использовался код, созданный Л. Е. Довгиловичем в Московском исследовательском центре «Шлюмберже».

Его основные характеристики:

• Расчёт уравнения Навье в общем случае (произвольный вид тензора жёсткости)

• Конечно-разностный подход с использованием несимметричных шаблонов.

• Возможность использовать криволинейные сетки, изоморфные прямоугольной сетке (в данной работе не использовалась).

• Восьмой порядок аппроксимации по пространству, второй порядок — по времени.

• Выполняется параллельно на CPU с использованием технологии MPI, либо на графических ускорителях с использованием CUDA.

• Сейсмоисточники реализованы трёх типов: сосредоточенной силы, центра расширения и роторные.

• В качестве сигнатуры сейсмоисточника берётся либо некоторая зависимость от времени (из файла), либо импульс Рикера:

где vq — центральная частота. В.З. Алгоритм упругой миграции

Рисунок В.1. Схема работы алгоритма упругой миграции.

Схема работы алгоритма миграции показана на рисунке (В.1). После чтения глобальных параметров задачи выполняется цикл по всем источникам, результатом на каждой итерации является построение «локального» изображения. Для его построения сначала решается задача прямого моделирования (1) с помощью кода Кау1ег-5, соответствующая конечно-разностному решению дифференциального уравнения (3.24). Затем по рассчитанному полю на сейсмоприёмниках и по «истинным» сейсмограммам вычисляется правая часть для построения «обратного» поля, согласно формулам (5.6), (5.7).

Построение источников для «обратного» поля

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

Для построения источников, входящих в правую часть «обратного» поля, мы разделяем сейсмограммы на продольные и поперечные волны и на их основе оцениваем значение дивергенции и ротора поля смещения. В случае «приближённых» сейсмограмм, мы можем рассчитать дивергенцию и ротор численно. Далее разница роторов и дивергенций подаётся на вход программы Нау1ег-5, и решается задача по расчёту «обратного» поля (рисунок В.1, 3). После этого по сохранённым «прямому» и «обратному» полям производится построение изображения.

Рисунок В.З. Схема процедуры построения изображения.

Процедура построения изображения

Для построения изображений «прямое» и «обратное» поле разделяются на продольные и поперечные волны, как показано в параграфе 2.3. Далее в каждой точке среды в каждый момент времени оценивается направление распространения волн и рассчитывается предполагаемый угол наклона границы, на основе которого строится фильтр по углу наклона границы (см. главу 5.4). Затем из разделенных полей собираются 4 типа изображений (РР, РБ, БР, ББ) с использованием соответствующих 01р-фильтров.