Построение изображения сейсмического разреза по модели bp2004 Benchmark Текст научной статьи по специальности «Математика»

УДК 519.688

В. Г. Байдин

Московский физико-технический институт (государственный университет)

Московский исследовательский центр Шлюмберже

Построение изображения сейсмического разреза по модели ВР2004 Benchmark

Задачи сейморазведки всегда требовали серьёзных вычислительных мощностей.

В данной работе рассматривается метод миграции в обратном времени (Reverse time migration — RTM). Автором проделано тестирование программного кода [9] на модели EAGE ВР2004 Benchmark. Модель представляет собой двумерный геологический разрез с заданным распределением скоростей (рис. 1). Модель создана по инициативе геофизического общества EAGE и предназначена для сравнительного тестирования программных разработок.

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

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

В работе также освещены вопросы фильтрации изображений. Рассмотрены различные методы нормализации и высокочастотная фильтрация.

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

1. Введение

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

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

С ростом вычислительных мощностей усложнялись и методы построения изображений. Метод миграции в обратном времени (I! ГМ) — один из наиболее продвинутых методов построения изображений [11], [8]. Построение изображения идёт в три стадии. Сначала моделируется волновое поле от каждого источника, затем моделируется поле от сейсмограмм, запущенных в обратном времени и взятых в качестве источников. Далее выполняется корреляция двух полученных полей и последующее суммирование по всем источникам. Полученное поле и является изображением среды.

Теория метода I! ГМ описана в работах [7], [10]. Отталкиваясь от этих результатов, приведём математическое обоснование данного метода и выведем основные формулы.

1.1. Математическая постановка задачи

В данной работе рассматривается акустическое приближение, описываемое волновым уравнением. Среда О С Кп, п = 2, 3 с границей Г задаётся одним скалярным параметром

— скоростью звука с, зависящей от пространственных координат и неизменной во времени.

Удобно ввести параметр к = ^2- Выпишем уравнения относительно давления и:

яиц — Аи = / (х,ї) ,х Є 0,ї Є [0, Т] , (1)

Ъки = /, (2)

где / (х, і) — функция сейсмоисточника. В качестве начальных условий обычно берётся состояние покоя:

(и (х,г = 0) = 0,

(х, і = 0) = 0.

Г

Неймана, так и более сложные граничные условия, имитирующие неотражение волн:

и(х,ї) х& = 0, (4)

ди

дп

ди ди дп + С ді

ХЄГ = 0, (5)

жег = 0. (6)

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

Задача миграции в обратном времени ставится следующим образом. В некоторых точках среды расположены СеЙСМОИСТОЧНИКИ 5 = {$г} и сейсмоприём ники К = {гг}. Известно значение волнового поля (г%, ¿) при 0 < Ь < Т в точке сейсмоприёмника г% от источника Sj. Также известна некоторая приближённая модель среды к. Задача: по имеющимся данным восстановить местоположение отражающих границ, т.е. поверхностей резких градиентов

скорости звука Со; соответственно к0 = А-

с0

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

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

КГМ как первый шаг полного волнового обращения

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

Обратная задача ставится как минимизация функционала ошибки:

т

Е = /* ^ 2 (и3 (гі,1) — ио (гі, ¿)) М,

о ^

где и,3 (гг,Ь) — значение сейсмограммы на г-м приёмнике от _/-го источника, и,3 (п) — значение решения уравнения

Ьки3 = 83 (8)

на сейсмоприёмнике Гг.

В терминах обобщённых функций Е можно представить в виде

т

Е = 1 [ {У'3 (Х,^) — ио (Х,^ ^ ^ (Х — гг) ЛЫх, (9)

3 ° о г

где и0 (х,Ь) — гладкое продолжение точечной функции, определённой в точках г* и имеющей в них значение и0 (г*, ¿).

Метод градиентного спуска

Разложим функционал Е по формуле Тейлора по к:

Е (к + 5к) = Е(к) + УкЕ ■ 5к + О (\\5к\\2) . (10)

Будем считать, что вариация модели 5к достаточно мала, чтобы пренебречь порядками,

начиная со второго.

Согласно методу градиентного спуска первый шаг итерации равен:

5к = —аУкЕ, (11)

где а > 0 — некоторое число.

Для нахождения шага необходимо оценить градиент функционала.

Нахождение градиента функционала

Запишем уравнение (1) в вариациях:

д 2

(к + 5к) —г (и + 5и) — Л (и + 5и) = /. иЪ2

Раскроем скобки, отбросив малые величины порядка больше первого:

кии + к5и,и + ёкиц — Ли — Л5и = /.

Используя (1), получаем

к5и,и + ёкиц — Л5и = 0, или, по определению оператора Ьк (2):

Б к&и = —5 кии. (-^)

Проварьируем функционал Е (9) относительно и:

Е + 5Е = Е (и + 6и) =

^ I I (и3 (х, ¿) + 5и,3 (х, Ь) — и0 (х, ¿)^ ^ 5 (х — Гг) (М(1х.

Тогда

о

т

6Е = 1 ^ J ! | (и3 + 5и3 — и0^ — (и3 — и0^ | ^ 6 (х — гг) dtdx

о

^11 I 6и3 ^ — иоо) ^ 5 (х — гг) 1 Мйх.

3 В 0 I г )

Запишем тождество Лагранжа. Для любых гладких функций и и ф выполняется

(13)

т

т

J ! [иЬкф — фЬки] йЫх = J J [кифи — кфии — иАф + фАи] йЫх =

= у [ифг — фщ]

В

т

т

(1х +

иУф — фУи

■ йБ<И.

о дВ

Теперь пусть функция ф3 является решением следующей задачи:

Ьк Ф3

' = 'и3 — и1

,) Е* 5 (Х — П)

ф3 (х, £ = Т) = 0, Ф{ (х, г = т) = о, ф3(х,Ь) =0,

жег

(14)

(15)

а функция и3 — решение задачи (1), (3), (4), / = в3, см. также (8). Тогда правая часть тождества Лагранжа (14) обратится в ноль. Согласно (13) и (12)

т

т

5Е =

£

5и3 Ьк ф3сМс!х

т

т

/ / ф3Ьк5и3dtdx = ^ / / ф3и3и5кМйх

В

0

В

0

Используя нулевые граничные условия функции Ф1, имеем

т

т

(И ■ Ьк.

В соответствии с (10) получаем следующую оценку:

т

м.

0

(16)

Вычисление изображения

Градиет функционала Е (16) даёт нам возможность строить изображения, то есть проводить миграцию. Действительно, если предположить, что приближённая скоростная модель с достаточно хорошо описывает истинную кинематику распространения волн, можно ожидать, что локально максимальные амплитуды УкЕ (16) будут соответствовать интерфейсам, то есть поверхностям разрыва скорости звука. Именно в этих местах сглаженная скоростная модель сильно отличается от истинной со- Поэтому на основе (16) мы определяем изображения, где Б3 (х,Ь) — решения волнового уравнения (1), а К3 (х,Ь) — решения уравнения (15), то есть волнового уравнения с сейсмоприёмниками в правой части, запущенными назад во времени. Полученное выражение имеет вид

о

о

о

h (х) = J St (x,t)Rt (x,t)dt, (17)

На практике, помимо формулы (17), также применяется другая, дающая схожий результат, например:

h (х) = J S (x,t) R (x,t)dt. (18)

1.3. Фильтрация изображений

Изображения, полученные с помощью формул (17, 18) имеют ряд существенных недостатков [5], которые ухудшают качество изображения.

Неравномерная освещённость

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

освещённости и вследствие невысокого динамического диапазона финального изображе-

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

^ (х) = JS• (1Э)

1ап (х)= f R2 (x?t)it ’ (20)

I (х')

Isan (х) f |S (х, I • IR (х, f) l^f.

Высокочастотный фильтр

В полученных избражениях также зачастую присутствует низкочастотный шум [6]. Чтобы его убрать применяют фильтрацию, так называемый High-Pass Filter:

Г 1 (z—)2

Ihpf (х) =І(х) - —r=I (0 e -2 , (22)

J ay 2n

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

2. Миграция для модели ВР2004 Benchmark

2.1. Конфигурация расчётной системы

Для тестирования методов расчёта волновых полей и реализованных подходов улучшения качества изображения была использована модель ВР 2004 Benchmark размером 82, 5 х 12 км2. Скоростная модель представлена на рис. 1. По данной модели рассчитываются волновые поля, причём схема расположения сейсмических источников и приёмников является такой же, как принято использовать в практических сейсмических наблюдениях. Скоростная модель задана на регулярной сетке с шагом 6, 25 х 6, 25 м, размеры модели — 13201 х 1921 точек. Всего имеется 1348 сейсмограмм длительностью 12 секунд. Для каж-

12, 5

м.

0.0 km 11.8 km 23.6 km 35.4 km 47.1 km 58.9 km 70.7 km

Рис. 1. Скоростная модель среды

2.2. Код для расчёта миграции и параметры дискретной модели

Для вычислений использовался конечно-разностный код, написанный J1.E. Довгиловичем, на основе схем высокого порядка аппроксимации по пространству и второго порядка по времени. Код позволяет использовать возможности архитектуры с общей памятью на основе POSIX Threads и работать распределённо, используя интерфейс MPI. Расчёты проводились на кластере МФТИ-60. Одна подзадача (из 1348) имела размеры 13201 х 1921 точек по пространству, 30000 шагов по времени, занимала 16 четырёхъядерных узлов и работала в среднем 45 минут. Одновременно на кластере выполнялось до пяти независимых подзадач. Полное время расчёта составило около двух недель. Суммирование и фильтрация происходили на отдельной машине с помощью программ, написанных на языке Python.

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

3. Обработка результатов расчётов

В результате работы программы было получено следующее изображение, фрагмент которого показан на рисунке 2.

Рис. 2. Часть изображения, посчитанного по исходной формуле (18)

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

Как было сказано ранее, полученное изображение необходимо нормализовать. Были проделаны различные нормализации (19) - (21) и показано, что нормализация по формуле

(20) даёт наилучший результат, см. рис. 3.

0.0 km

2.4 km

4.8 km

7.2 km

9.6 km

37.5 km 39.6 km 41.7 km 43.7 km 45.8 km

Рис. 3. Нормализованное изображение по формуле (20)

Применение High-Pass Filter (22) после нормализации делает границы более чёткими.

0.0 km

2.4 km

4.8 km

7.2 km

9.6 km

37.5 km 39.6 km 41.7 km 43.7 km 45.8 km

Рис. 4. Отфильтрованное HPF-фильтром изображение (ф-ла (22))

4. Сравнение компактных схем с центральными разностями

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

1) На сетке 6, 25 х 6, 25 м с использованием центральных разностей восьмого порядка (1348 источников).

2) На сетке 12, 5 х 12, 5 м с использованием компактных схем 12-го порядка (1348 источников).

3) На сетке 12, 5 х 12, 5 м с использованием центральных разностей 20-го порядка (100 источников).

4.1. Сравнение качества

Были получены изображения для первого и второго случаев:

Рис. 5. Изображение, полученное на сетке 6, 25 х 6, 25 метров с использованием центральных разностей четвёртого порядка

Рис. 6. Изображение, полученное на сетке 12, 5 х 12, 5 метров с использованием компактной схемы 12-го порядка

Сравнение позволяет сделать вывод о том, что оба изображения имеют одинаковое качество, то есть компактные схемы 12-го порядка на сетке 12, 5 х 12, 5 м работают не хуже центральных разностей 8-го порядка на сетке 6, 25 х 6, 25 м.

0.0 3.0 6.0

Length, km

Рис. 7. Часть изображения, полученного на сетке 12, 5 х 12, 5 метров с использованием центральных разностей 20-го порядка

0.0 3.0 6.0

Length, km

Рис. 8. Часть изображения, полученного на сетке 12, 5 х 12, 5 метров с использованием компактной схемы 12-го порядка

Также мы пришли к выводу, что качество изображения, полученного на компактной схеме 12-го порядка, приблизительно равно качеству изображения, полученного на центральной разности 20-го порядка. Для этого мы рассмотрели изображение от 100 источников, см. рис. 7, 8. Отметим, что этот вывод совпал с теоретической оценкой, сделанной для гладких одномерных функций [12].

4.2. Сравнение производительности

Для тестирования производительности различных режимов работы алгоритма использовалась рабочая станция с восьмиядерным процессором Intex Xeon 5560. Для минимизации влияния других факторов на скорость вычислений запись полей на диск была выключена, таким образом, сравнивалась именно скорость вычислений. Были получены следующие средние значения длительности одного шага по времени (рис. 9) — 4 : 1, 8:1 для ЦРС8 на сетке 12, 5 х 12, 5 м; ЦРС20 на 6, 25 х 6, 25 м и КС на 6, 25 х 6, 25 м.

Рис. 9. Время выполнения одного шага по времени

5. Выводы и дальнейшие исследования

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

цельная программа aRTM, получающая изображения сейсмического разреза в 2D и 3D случаях.

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

Литература

1. Shin С., Yoon К., Marfurt K.J. Efficient calculation of a partial-derivative wavefield using reciprocity for seismic imaging and inversion // Geophysics. — 2001. — V. 66, N 6. -P. V1856-V1863.

2. Pratt G., Shin C., Hicks Gauss-Newton and full Newton methods in frequencv-space seismic waveform inversion // Geophvs. J. International. — 1998. — V. 133, N 2. — P. V341-V362.

3. Kaelin B., Guitton A. Imaging condition for reverse time migration // SEG Technical Program Expanded Abstracts. - 2006. — V. 25, N 1. — P. V2594-V2598.

4. Billette F., Brandsberg-Dahl S. The 2004 BP velocity benchmark // 67th Annual Internat. Mtg., EAGE, Expanded Abstracts. — 2005. — P. B035.

5. Yoon K., Marfurt K. Challenges in reverse-time migration // 2004 SEG Annual Meeting. — 2004. - N October.

6. Kaelin B.,Carvajal C. Eliminating imaging artifacts in RTM using pre-stack gathers // Imaging. — 2011. — P. \ 3125 \ 3120.

7. Baysal E., Kosloff D., Sherwood J. Reverse time migration // Geophysics. — 1983. — V. 48,

N Ii. - P. V1514-V1524.

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

9. Довгилович Л.Е. О трехуровневом параллельном алгоритме обратной миграции во временной области (RTM) на основе компактных схем // XVIII конференция «Теоретические основы и конструирование численных алгоритмов решения задач математической физики». — 2010.

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

11. Курин Е.А. Сейсморазведка и суперкомпьютеры // Вычислительные методы и программирование. — 2011. — Т. 12. — С. \ 3 1 \ 30.

12. Довгилович Л.Е., Софронов И.Л. О применении компактных схем для решения волнового уравнения: препринт 1111 \ I им. Келдыша РАН. — М., 2008. — № 84.

Поступим в редакцию 23.04-2012.