Реконструкция глубинной скоростной модели упругой среды методом обращения полного волнового поля

Гадыльшин Кирилл Геннадьевич; Чеверда Владимир Альбертович

УДК 550.34.01

DOI: 10.183 03/2618-981X-2018-4-194-203

РЕКОНСТРУКЦИЯ ГЛУБИННОЙ СКОРОСТНОЙ МОДЕЛИ УПРУГОЙ СРЕДЫ МЕТОДОМ ОБРАЩЕНИЯ ПОЛНОГО ВОЛНОВОГО ПОЛЯ

Кирилл Геннадьевич Гадыльшин

Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, кандидат физико-математических наук, старший научный сотрудник, тел. (383)330-27-96, e-mail: gadylshinkg@ipgg.sbras.ru

Владимир Альбертович Чеверда

Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, доктор физико-математических наук, профессор, зав. лабораторией, тел. (383)333-00-54, e-mail: gadylshinkg@ipgg.sbras.ru

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

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

RECOVERY OF ELASTIC DEPTH MACROVELOCITY MODEL BY FULL WAVEFORM INVERSION

Kirill G. Gadylshin

Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 3, Prospect Аkademik Koptyug St., Novosibirsk, 630090, Russia, Ph. D., Senior Researcher, тел. (383)330-27-96, e-mail: gadylshinkg@ipgg.sbras.ru

Vladimir A. Tcheverda

Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 3, Prospect Аkademik Koptyug St., Novosibirsk, 630090, Russia, D. Sc., Professor, Head of Laboratory, тел. (383)333-00-54, e-mail: cheverdava@ipgg.sbras.ru

Depth velocity model is necessary for correct placement of migration images of local geological objects. Reconstruction of such kind of model by seismic observations is the most laborintensive procedures in seismic data processing and claims huge computer resources. We propose original numerical technique providing stable recovery of macrovelocity model for both P- and S-waves. This technique is based on decomposition of full velocity models onto two constituents: Smooth slowly variated in space propagators and Sharp reflectors.

Key words: seismic waves, macrovelocity model, full waveform inversion, singular spectrum.

Введение

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

Весьма многообещающим для развития методов реконструкции макроско-ростной модели представлялся подход, основанный на обращении полного волнового поля. Его исходная формулировка сводилась к нелинейному методу наименьших квадратов [5], где в роль целевого функционала выполняло среднеквадратичное уклонение зарегистрированного волнового поля м0 от синтетического, рассчитанного для текущей модели среды:

Ф(т) = |^(т) - ^> (!)

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

Однако весьма быстро выяснилось, что реконструкция макроскоростной модели даже в простейшем случае скалярного волнового уравнения по-прежнему остается чрезвычайно неустойчивой процедурой [8] и требует либо привлечения достаточно близкого к истинной модели начального приближения, либо наличия нереально низких временных частот в спектре зондирующего сигнала, либо нереально больших выносов источник-приемник [4]. Проведенное в работах [1, 10] исследование этого оператора путем изучения сингулярного спектра его производной Фреше, показало, что только при выполнении упомянутых жестких ограничений устойчивые подпространства сингулярных векторов содержат плавно изменяющиеся в пространстве функции. Следовательно, никакие ухищрения в организации итерационного процесса не позволят извлечь информацию о поведении макроскоростной составляющей без модификации исходного целевого функционала (1).

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

Пусть полупространство Х3 > 0 заполнено неоднородной упругой средой с параметрами Ламэ Х(хьХ2,Х3), хьХ2,Х3) и плотностью р(хьХ2,Х3). Рас-сматрим упругое волновое уравнение в области временных частот (формулировка в скоростях-напряжениях):

1 Здесь В - пространство данных.

• 1 ^ \Гт ъди г

тМ{хъхъ)и-Р---= /

дхл

дх?

где вектор неизвестных и = (и\,и3,<3\ 1,а33,а|з) , а матрицы определяются следующими соотношениями:

м = Р-!2x2 0 ^ , Р = ( 0 А 1 , е= Г 0 в1

10 53х3 1А 0 , 0)

(1 0 01 (0 0 11

А = , в =

у0 0 1) V0 1 0)

5 =

а - Ь 0Ч -Ь а 0

Ч 0 0 С ;

а =

Х + 2ц 4ц(Х + ц)

Ь=

X

4ц(Х + ц)

1

с = —. ц

Правая часть в уравнении (2) задает воздействие источника. В наших численных экспериментах использовался источник типа центра расширения. При этом волновое поле предполагалось известным для некоторой

конфигурации приемников хг = (х[,х'^ ) и источников х5 =(х(,хз).

Модификация целевого функционала

Предлагаемая нами модификация целевого функционала основывается на декомпозиции пространства полных скоростных моделей на два подпространства [2]:

• Подпространство медленно изменяющихся пропагаторов, предписывающих времена распространения волн, но не меняющих их направление;

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

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

Ур,5 = (РУ )Р5 + (гУ )р,

5 =

(3)

где р означает пропагатор, т. е. макроскоростная модель, а г - пространствен-

Л

ный рефлектор (см. рис. 1).

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

Напомним, что при известной макроскоростной модели расположение рефлектора корректно определяется в результате выполнения миграционного преобразования зарегистрированных данных. Более того, существует ряд методов, позволяющих строить изображения рефлекторов в «истинных» амплитудах, т. е. получения изображений, интенсивность которых связана с отражающей способностью соответствующих локальных геологических объектов, напрямую определяемой контрастом в акустических/упругих параметрах [3, 6]. Исходя из этого, искомая полная скоростная модель при точно известном пропа-гаторе может быть представлена следующим образом:

Ур,£ = (рУ)р,£ + ЪГШЫР8 (рУР8 ) < * > (4)

где Мр £(рУ) < * > есть результат применения миграционного преобразования

для продольных/поперечных волн (в дальнейшем просто «миграционный оператор») к так называемому «временному» рефлектору, который, в свою очередь, есть результат продольной/поперечной «демиграции» [9] пространственного продольного/поперечного рефлектора г. Оператор Ж обеспечивает корректировку получаемых при этом амплитуд изображения [3, 6], а проектор Пг обеспечивает погружение полученного результата в подпространство рефлекторов.

скорость, м/с скорость, м/с „Д.

Рис. 1. Разложение полной скорости на пропагатор и рефлектор (слева направо: полная скоростная модель, макроскоростная составляющая,

пространственный рефлектор)

Таким образом, вводятся два новых подпространства в пространстве искомых моделей:

• Подпространство пропагаторов р¥. Пропагатор строится как для продольных, так и для поперечных волн;

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

Такая декомпозиция позволяет выполнить модификацию целевого функционала:

2

Ф(pVP,pVs,s)= F(pVPS+UrWMPS(pVPS)<s>)-u0 т . (5)

Ы D )'

и рассматривать следующую формулировку обратной задачи:

Ф (pVP,pVs,s) -> min , (6)

PVp,PVs>s

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

rP S = UrWMP S(pVP s) < s >. (7)

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

Организация итерационного процесса минимизации модифицированного функционала

Поиск точки минимума модифицированного функционала (5) реализуется путем организации двух последовательных процессов:

1. Квадратичная минимизация по временному рефлектору s.

Временной рефлектор лежит в пространстве данных и соответствует сейсмическим волнам, возвращающимся на систему наблюдений после взаимодействия падающей волны с пространственным рефлектором. Следовательно, временной рефлектор должен обладать высокочастотной составляющей, которая определяется спектром зондирующего сигнала. Это означает, что дискретизация временного рефлектора должна совпадать с дискретизацией зарегистрированных волновых полей. Но тогда и дискретизация пространственного рефлектора должна быть соответствующей и составлять доли доминирующей длины волны. Отметим, что это не вызывает сколько-нибудь значительных затруднений, так как применение миграционной процедуры Mp s(pVp s) < s > сводится к применению сопряженного линеаризованного оператора решения прямой задачи и не требует итераций.

2. Нелинейная минимизация по пропагаторам pVp и pVs.

Зависимость модифицированного функционала от пропагаторов существенно нелинейная и вычисление его градиента по этим направлениям требует вычисления второй производной от оператора F(m) (см. соотношение (1)). Однако, в силу того, что пропагаторы, суть, плавно меняющиеся функции, для их

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

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

Для численных экспериментов была выбрана представленная на рис. 2 синтетическая модель упругой среды, описывающая строение месторождения Gullfaks в Северном море [11].

1500 2000 2500 3000 3500 4000

Расстояние(м)

Рис. 2. Распределение скоростей продольных волн для синтетической модели Gullfaks, использованной для выполнения численных экспериментов

Для решения обратной задачи итерационным методом, изложенным в [7], были рассчитаны синтетические монохроматические волновые поля, рассчитанные для 21 временной частоты, равномерно распределенных в интервале 5-25 Гц. При этом использовалась следующая расстановка источников и приемников:

• 39 равномерно расположенных с шагом 100 м на глубине 5 м ниже сво-

-5

бодной поверхности источников типа центра расширения ;

3 Самый верхний слой в этой модели - вода со скоростью распространения продольных волн 1500 м/сек и нулевой скоростью распространения поперечных волн.

• 198 равномерно расположенных с шагом 20 м на той же глубине двух-компонентных приемников.

В качестве начального приближения использовалась вертикально-неоднородная среда, взятая в самом левом вертикальном сечении модели (см. рис. 2). Полученные при этом скорости продольных и поперечных волн были заглажены с использованием Гауссова ядра (см. рис. 3).

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

ооэ

Рис. 3. Начальное приближение (слева направо: скорость продольных волн, скорость поперечных волн, отношение у)

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Рис. 4. Слева заданная модель распределения продольной скорости (справа - результат восстановления с использованным начальным приближением (рис. 3))

Модель Ув Восстановленная Уб

О 1000 2000 3000 4000 0 1000 2000 3000 4000

Расстояние, м Расстояние, ы

Рис. 5. Слева заданная модель распределения поперечной скорости (справа - результат восстановления)

Модель УрА/э Восстановленное УрЛ/э

0 1000 2000 3000 4000 0 1000 2000 3000 4000

Расстояние, м Расстояние, м

Рис. 6. Слева соотношение у в исходной модели (справа - в восстановленной)

Заключение

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

• Плавно меняющуюся в пространстве макроскоростную компоненту;

• Резко изменчивый пространственный рефлектор.

Работа выполнена при поддержке Российского Научного Фонда, проект 17-17-01128. Расчеты проведены на вычислительных мощностях Сибирского Суперкомпьютерного Центра. Мы благодарны компании StatoilHydro за разрешение использовать в численных экспериментах упругую модель Gullfaks.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Алексеев А. С., Костин В. И., Хайдуков В. Г., Чеверда В. А. Восстановление двумерных возмущений скорости вертикально-неоднородной акустической среды по данным мно-

гократного перекрытия (линеаризованная постановка) // Геология и геофизика. - 1997. -Т. 38 (№ 2). - С. 1980-1992.

2. Гадыльшин К. Г., Чеверда В. А. Реконструкция глубинной макроскоростной модели путем обращения полного волнового поля // Доклады Академии наук. - 2017. - Т. 476(6). -С. 693-697.

3. Протасов М. И., Чеверда В. А. Построение сейсмических изображений в истинных амплитудах // Доклады РАН. - 2006. - Т. 407 (4). - С. 528-532.

4. Alekseev A. S., Avdeev A. V., Fatianov A. G., Cheverda V. A. Wave processes in vertically-inhomogeneous media: a new strategy for a velocity inversion // Inverse Problems. -1993. - V. 9(3). - P. 367-390.

5. Bamberger A., Chavent G., Hemon C., Lailly P. Inversion of normal incidence seismograms. // Geophysics - 1982. - V. 47(5). - P. 757-770.

6. Beylkin G. Imaging of discontinuities in the inverse scattering problem by inversion of a causal generalized Radon transform // Journal of mathematical physics. - 1985. - V. 26. - P. 99108.

7. Belonosov M., Tcheverda V., Neklyudov D., Kostin V. and Dmitriev M. 3D elastic frequency-domain iterative solver for full waveform inversion // SEG International Exposition and 86th Annual Meeting. - 2016. - Expanded abstracts. - P. 3825-3829.

8. Gauthier O., Virieux J., Tarantola A. Two-dimensional nonlinear inversion of seismic waveforms: Numerical results // Geophysics. - 1986. - V. 51(7). - 1387-1403.

9. Iversen E., Tygel M., Ursin B., de Hoop M. Kinematic time migration and demigration of reflections in pre-stack seismic data // Geophysical Journal International. - 2012. - V. 189. -p. 1635-1666.

10. Kostov C., Neklyudov D., Silvestrov I., Tcheverda V. Full waveform inversion for macrovelocity model reconstruction in look-ahead offset VSP: Numerical SVD-based analysis // Geophysical Prospecting. - 2013. - V. 61(6). - P. 1099-1113.

11. Stovas A., Landro M., Arnsten B. A sensitivity study based on 2D synthetic data from the Gullfaks field using PP- and PS time-lapse stacks for fluid pressure discrimination // Journal of Geophysics and Engineering. - 2006. - V. 3. - P. 314-328.

REFERENCES

1. Alekseev A. S., Kostin V. I., Hajdukov V. G., Cheverda V. A. Vosstanovlenie dvumernyh vozmushchenij skorosti vertikal'no-neodnorodnoj akusticheskoj sredy po dannym mnogokratnogo perekrytiya (linearizovannaya postanovka) // Geologiya i geofizika. - 1997. - T. 38 (№ 2). -S.1980-1992.

2. Gadyl'shin K. G., Cheverda V. A. Rekonstrukciya glubinnoj makroskorostnoj modeli putyom obrashcheniya polnogo volnovogo polya // Doklady Akademii nauk. - 2017. - T. 476(6). -S. 693-697.

3. Protasov M. I., Cheverda V. A. Postroenie sejsmicheskih izobrazhenij v istinnyh amplitudah // Doklady RAN. - 2006. - T. 407 (4). - S. 528-532.