Содержание к диссертации
Введение
1 Теоретические основы методов обращения полного волнового поля 15
1.1 Обратная динамическая задача сейсмики 16
1.1.1 Явные динамические методы 16
1.1.2 Динамические методы решения линеаризованной обратной задачи 18
1.1.3 Связь линеаризованного обращения с методами сейсмической миграции 24
1.2 Метод обращения полных волновых полей 26
1.2.1 Вычисление градиента функционала невязки 27
1.2.2 Определение направления поиска 29
1.3 Потенциальные возможности метода и практические затруднения 29
1.3.1 Проблема определения трендовой составляющей 29
1.3.2 Модификации метода обращения полных волновых полей 30
2 Теория метода спектральных чувствительностей 33
2.1 Введение 33
2.2 Построение аналитических решений 36
2.2.1 О временах пробега рефрагированных волн и волн шепчущей галереи 41
2.3 Спектральные чувствительности 44
2.4 О связи спектров рассеивающей неоднородности и рассеянных волновых полей 48
3 Метод спектральных чувствительностей, применение 55
3.1 Спектральные чувствительности для среды из двух полупространств 55
3.1.1 Коэффициенты отражения 55
3.1.2 Аналитическое разделение волновых полей 63
3.2 Выводы 64
3.3 Спектральные чувствительности для среды в виде слоя над полупространством и роль кратных волн и волн-спутников
3.3.1 О решении прямой задачи 66
3.3.2 Осредненные или интегральные чувствительности 68
3.3.3 Выводы 75
3.4 Спектральные чувствительности для среды с постоянным градиентом скорости и роль рефрагированных волн в восстановлении скоростных неоднородностей 75
3.4.1 Квази-плоские волны 76
3.4.2 Среда с постоянным вертикальным градиентом скорости 77
3.4.3 Асимптотика функции Макдональда 79
3.4.4 О решении обратной задачи 80
3.4.5 Спектральные чувствительности 83
3.4.6 Выводы 84
4 Псевдоспектральное обращение полных волновых полей 86
4.1 Псевдоспектральное моделирование 86
4.2 Псевдоспектральное обращение полных волновых полей в отсутствии низких временных частот в наблюденных данных
4.2.1 Теория 90
4.2.2 Пример работы алгоритма на синтетической модели Мармузи 91
4.2.3 Фильтр масштабирующий частоту
4.3 Псевдоспектральное обращение полных волновых полей в присутствии шума в наблюденных данных 97
4.4 Выводы 101
Заключение 102
Благодарности 103
Литература 104
- Связь линеаризованного обращения с методами сейсмической миграции
- О временах пробега рефрагированных волн и волн шепчущей галереи
- Спектральные чувствительности для среды в виде слоя над полупространством и роль кратных волн и волн-спутников
- Псевдоспектральное обращение полных волновых полей в отсутствии низких временных частот в наблюденных данных
Введение к работе
Актуальность темы исследований. Сейсмическая разведка является одним из наиболее достоверных источников информации о недрах Земли при поисках полезных ископаемых и мониторинге их месторождений. Современные высокоточные технологии бурения скважин позволяют увеличить добычу углеводородов, однако для этого также требуются более точные представления о геологическом строении месторождений. Эти представления можно получить из данных о распределении физических параметров недр (в первую очередь распределения скорости продольных волн), содержащихся в сейсмограммах, решая обратную задачу сейсмики. Классические методы сейсмической томографии предполагают лучевое описание распространения волновых фронтов в среде, что является кинематической обратной задачей сейсмики. Динамическая обратная задача сейсмики предполагает работу не только с годографами зарегистрированных волн, но и с волновыми полями. Ближе всего к кинематическим методам по сложности вычислений является обратная динамическая задача для волнового уравнения, описывающего распространение волн в акустической среде постоянной плотности. В рамках такой физической модели и предполагается проводить инверсию в настоящей работе.
Использование обращения полных волновых полей при обработке данных сейсмической разведки позволяет получить распределение физических параметров среды в более высоком, в сравнении с лучевыми методами, разрешении. Обращение полных волновых полей также не требует пикирования первых вступлений на сейсмограммах. Основными недостатками метода, препятствующими его внедрению на практике, являются необходимость наличия низких частот в наблюденных данных и относительно высокие требования к вычислительным ресурсам. Моделирование полных волновых полей во временной или частотной области - наиболее ресурсоемкий процесс при применении данного метода. Как правило, для моделирования сейсмических полей используются конечно-разностные методы. В работе предлагается реализация метода обращения полных волновых полей с использованием псевдоспектрального моделирования, что позволяет существенно сократить количество необходимых для описания сейсмической среды параметров за счет использования более разреженных пространственных сеток. Разрешение сетки в данном случае оказывается в полном соответствии с разрешающей способностью метода обращения полных волновых полей. Для преодоления известной проблемы обращения полных волновых полей, состоящей в отсутствии низких частот в наблюдениях, в работе предложен способ регуляризации, включающий новый метод нестационарной фильтрации мигрированных
изображений, используемых в ходе решения обратной задачи сейсмики. Тестирование метода на акустической модели Мармузи показало, что в идеальных условиях наличия
низких частот в наблюденных данных псевдоспектральное обращение полных волновых полей по классической многомасштабной схеме (сначала обращаются низкие частоты, затем высокие) позволяет полностью восстановить строение данной модели. В отсутствии низких частот классический метод обращения полного волнового поля не дает удовлетворительных результатов решения обратной задачи сейсмики, однако применение предложенного в работе алгоритма регуляризации позволяет ослабить требования к наличию в данных низких частот при условии регистрации данных для больших выносов.
Целью данной работы является усовершенствование (исследование возможностей и последующая модернизация) метода обращения полных сейсмических волновых полей.
Для достижения поставленной цели было необходимо решить следующие задачи:
1. Исследовать потенциальные возможности метода обращения полного волнового поля,
что в рамках настоящей работы предполагает:
исследование зависимости разрешающей способности и скорости сходимости метода от особенностей трендовой (макро) составляющей опорной скоростной модели;
определение вкладов различных типов сейсмических волн в волновое поле, рассеянное на различных скоростных аномалиях, помещаемых в опорные модели;
исследование зависимости пространственного спектра рассеянного поля от пространственного спектра рассеивающей скоростной аномалии;
определение чувствительностей различных типов акустических волн (прямых,
отраженных, рефрагированных, кратных волн и волн шепчущей галереи) к
изменениям в пространственном спектре исследуемой скоростной аномалии.
2. Улучшить устойчивость и скорость сходимости метода обращения полного волнового
поля, а именно:
улучшить сходимость метода в отсутствии наблюдений для низких временных
частот;
- предложить способ улучшения изображения в областях неравномерной и недостаточной освещенности;
сократить число параметров, используемых для описания сейсмической среды
при моделировании, с сохранением точности решения прямых и обратных задач сейсмики; реализовать предложенные алгоритмы в виде программного комплекса, работающего по улучшенному алгоритму метода обращения полных волновых полей.
Основные положения, выносимые на защиту:
-
Предложенный в работе метод спектральных чувствительностей позволяет исследовать вклады различных типов акустических волн в решение обратной динамической задачи с применением метода обращения полных волновых полей.
-
Псевдоспектральное моделирование улучшает эффективность метода обращения полных волновых полей при восстановлении гладких моделей за счет использования меньшего количества параметров для описания среды по сравнению с конечно-разностным моделированием.
-
Предложенный в работе алгоритм регуляризации метода обращения полных волновых полей с использованием нестационарной фильтрации, позволяет сместить диапазон наблюденных данных, необходимых для успешного применения метода обращения полных волновых полей, в область более высоких временных частот.
Научная новизна:
-
Предложен метод спектральных чувствительностей.
-
Предложен алгоритм псевдоспектрального обращения полных волновых полей.
-
Предложен алгоритм нестационарной фильтрации "масштабирующий частоту".
Научная и практическая значимость. Предложено обобщение метода сфер Эваль-да в применении к задачам сейсморазведки. Продемонстрировано, что метод спектральных чувствительностей позволяет проводить подробный анализ освещенности исследуемых неоднородностей. С помощью метода спектральных чувствительностей можно оценивать восстановление неоднородностей в условиях неравномерной освещенности. Исследована роль различных типов акустических волн в решении обратной задачи сейсмики по методу обращения полных волновых полей. Для решения проблемы отсутствия низких частот в наблюдениях был предложен оригинальный метод регуляризации. Разработан комплекс алгоритмов и программ, позволяющий проводить псевдоспектральное обращение полных
волновых полей в приближении акустических моделей постоянной плотности.
Апробация работы. Результаты работы обсуждались и докладывались на научных семинарах: кафедры физики Земли физического факультета СПбГУ, компании Шелл (Shell GSI&BV, Netherlands), Колорадского университета горного дела (Colorado School of Mines), Университета им. короля Абдуллы (KAUST), института физики Земли им. О.Ю. Шмидта Российской Академии Наук. Основные результаты исследований докладывались на следующих конференциях:
Научно-практическая конференция "Сейсмические Технологии - 2015 Москва, 2015 // Казей В. В., Каштан Б. М., Троян В. Н, Tessmer E. Псевдоспектральное обращение полных волновых полей.
74th Annual EAGE Conference and Exhibition 2012, Copenhagen, Denmark // Kazei V.V., Ponomarenko A.V., Kashtan B.M., Troyan V.N., Mulder W.A. On the Contribution of Head Waves to Full Waveform Inversion. — 2012.
75th Annual EAGE Conference and Exhibition 2013, London, UK // Kazei V.V., Kashtan B.M., Troyan V.N., Mulder W.A. Spectral Sensitivity Analysis of FWI in a Constant-gradient Background Velocity Model.
SEG workshop on full waveform inversion: From near surface to deep 2013, Muscat, Oman // Kazei V.V., Kashtan B.M., Troyan V.N., Mulder W.A. FWI sensitivity analysis in the presence of free-surface multiples.
International workshop on Multi-scale Waveform Modeling and Inversion, March, 2015, KAUST, Saudi-Arabia // Kazei V.V., Tessmer E., FWI without low frequencies, non-stationary gradient filtering approach.
77th Annual EAGE Conference and Exhibition 2015, Madrid, Spain // Kazei V.V., Kashtan B.M., Troyan V.N., Mulder W.A. Free-surface Multiples and full-waveform inversion spectral resolution.
SEG International Exposition and 85th Annual Meeting in New Orleans, Louisiana, 2015 // Kazei V.V., Kashtan B.M., Troyan V.N., Mulder W.A. FWI spectral sensitivity analysis in the presence of a free surface.
Личный вклад.
Автором самостоятельно были разработаны:
метод спектральных чувствительностей для исследования возможностей метода обращения полных волновых полей в зависимости от полученной из априорной информации опорной модели и типов акустических волн зарегистрированных на сейсмограммах;
программный код в среде MATLAB, реализующий метод спектральных чувствительностей для ряда опорных скоростных моделей и типов волн;
программный код на языке С для псевдоспектрального обращения полных волновых полей с использованием, предложенной автором, методики регуляризации фильтром "масштабирующим" частоту.
Публикации. Основные результаты по теме диссертации изложены в 10 печатных работах, 3 из которых изданы в журналах, рекомендованных ВАК РФ для опубликования научных результатов диссертаций на соискание ученой степени кандидата наук, 7 — в тезисах докладов.
Объем и структура работы. Диссертация состоит из введения, четырех глав и заключения. Полный объем диссертации составляет 118 страниц с 50 рисунками. Список литературы содержит 168 наименований, из них 126 - иностранных.
Связь линеаризованного обращения с методами сейсмической миграции
Скоростную модель среды условно можно разделить на две составляющие: плавно меняющуюся трендовую составляющую, которая характеризует время распространения волн между удаленными точками среды, и составляющую, которая включает все локальные возмущения параметров среды, не оказывающие значительного влияния на времена пробега волн между удаленными точками среды. Зачастую обращение волновых полей рассматривается в рамках теории обратного рассеяния, то есть в предположении, что трендовая составляющая модели известна, а неизвестны лишь локальные неоднородности – возмущения физических параметров среды. Именно наличие таких возмущений приводит к возникновению регистрируемых отраженных, дифрагированных и рассеянных волн. В предположении заданной опорной модели можно линеаризовать задачу, воспользовавшись приближением Борна, и тем самым записать систему линейных уравнений относительно модели локальных возмущений. Линеаризованная обратная задача для волнового уравнения в случае вертикально-неоднородной модели среды впервые исследована в работах [147; 155], где было показано, что задача восстановления локальных неоднородностей скорости распространения волн сводится к одной из множества задач дифференциальной геометрии. Решение одномерной обратной задачи в линеаризованной постановке, по-видимому, было впервые получено в статье [36] для однородной акустической вмещающей среды, после чего развито в работах [17; 37]. Впоследствии метод получил название “борновское обращение” (Born inversion) и был обобщен на случай вертикально-неоднородной вмещающей среды [15]. Однако впервые неоднородность опорной модели среды была учтена в [33], где за основу взято лучевое представление (приближение Венцеля-Краммерса-Бриллюэна) для функции Грина в акустической среде. Структура возникающего при этом интегрального оператора и особенности численной реализации подхода исследованы в работе [38]. Аспекты, связанные с численной реализацией линеаризованного обращения на основе лучевого представления для функции Грина, представлены в работе [18].
Использование параксиального лучевого метода для приближенного вычисления тензора Грина в случае трехмерной упругой среды обсуждается в статье [11]. Линеаризованное обращение в совокупности с лучевой теорией, используемой для построения аппроксимации волновых полей в опорной среде, привело к появлению термина “лучевое борновское обращение” (Ray-Born inversion) [11; 46]. В рамках динамической теории упругости линеаризованное обращение данных многократных перекрытий для вертикально-неоднородной опорной среды рассмотрено в работе [154], где Романов В.Г. получил и проанализировал систему интегральных уравнений и доказал единственность определения с ее помощью упругих параметров Ламе и плотности. Впоследствии похожий подход к обратной задаче динамической теории упругости применен в работе [60] в спектральной области, позже обобщив его на случай анизотропной среды [61]. Белькин [12] предложил метод решения линеаризованной обратной задачи для акустического уравнения, основанный на обобщенном обратном преобразовании Радона и впоследствии расширенный на случай упругой среды в работе [27]. Совершенствование регистрирующей аппаратуры и вычислительной мощности позволило перейти от метода лучевой сейсмической томографии к более полной и потому более сложной методике дифракционной томографии (см.: [42; 158—160; 165; 166]). Метод дифракционной томографии, как и все ранее перечисленные подходы, основан на теории обратного рассеяния в предположении, что известное распределение скорости распространения сейсмических волн содержит неизвестные малые локальные скоростные аномалии. Термин “дифракционная томография” применительно к сейсморазведке введен в работе [42] с целью обозначить принципиальное отличие метода от ЛСТ. В общем случае искомыми параметрами могут являться произвольные параметры упругой среды — плотность, скорости продольных и поперечных волн например (см.: [165]). Поскольку термин “томография” исторически связан с определением именно скорости распространения волн, то и термин “дифракционная томография” в настоящее время используется в основном для задач восстановления скоростных аномалий в опорной среде (см., напр.: [123]). Если же ставится задача нахождения локальных аномалий упругих параметров среды, то используется термин “линеаризованное обращение” (linearized inversion; см., напр.: [56]). В работе [108] линеаризованная обратная динамическая задача сейсмики трактуется с точки зрения теории оптимизации, что позволяет находить решение путем последовательных итераций, связанных с минимизацией среднеквадратического функционала невязки между наблюденным и модельным волновыми полями. Метод получил название “итеративное линеаризованное обращение” (iterative linearized inversion). Подобные статистические модели интерпретации сейсмических данных предложены ранее [138; 166]. Итеративный метод асимптотического линеаризованного обращения на базе метода наименьших квадратов обсуждается также в работе [63]. К реальным сейсмическим данным метод был, по-видимому впервые, применен в работе [76]. Два эффективных алгоритма линеаризованного обращения волновых полей (во временной и частотной областях) в случае упругой среды представлены в статье [49]. Рассмотрим основные теоретические аспекты линеаризованного обращения волновых полей.
О временах пробега рефрагированных волн и волн шепчущей галереи
Времена пробега для рефрагированных, головных волн и волн шепчущей галереи первого порядка представлены на Рис. 2.2. Очевидно, что для не слишком больших выносов и не слишком больших частот времена пробега всех этих волн практически совпадают и волны складываются в высокоамплитудные волны-суперпозиции.
В этом параграфе вводится определение спектральных чувствительностей, при этом многие формальности опущены для краткости изложения, поскольку общих представлений о спектральных чувствительно стях может быть достаточно для понимания конкретных примеров следующей главы. В следующем параграфе приводится более детальное теоретическое обоснование метода спектральных чувствительностей.
Приближение Борна используется далее для аппроксимации рассеянного поля в частотной области для источника и приемника с горизонтальными координатами xs и хг соответственно, расположенными на одной глубине zs = zr: малое возмущение квадрата медленности в опорной среде. G(xs,r) и G(xg,r) — функции Грина от источника и приемника в точку среды г, соответственно. Мы предполагаем, что функция источника совпадает для всех подобных экспериментов. Поскольку далее рассматривается монохроматический случай, этот множитель отбрасывается. После преобразований Фурье по горизонтальным координатам источника и приемника, возмущение волновых полей выглядит следующим образом:
Теперь, G - функция Грина для опорной среды в спектральной области, или преобразование Фурье от стандартной функции Грина в частотной области G. В параграфе 2.2.1 выводятся явные аналитические выражения для G для моделей типов 0 и 3. Используя явные аналитические выражения для опорной среды, можно получить линейное соотношение между спектром возмущения полных волновых полей и пространственным спектром рассеивающей неоднородности:
Фурье возмущения квадрата медленности 8W(r), где К = (KX,KZ)T — вектор в спектральной области. С{к3) и С{кд) - коэффициенты отражения от границы между полупространствами для монохроматических плоских волн от источников и приемников соответственно. Первый член в уравнении (2.16) соответствует волнам, которые идут от источников, рассеиваются неоднородностью и затем возвращаются к линии приемников. Обращение этого члена соответствует миграции в однородной опорной среде. Аналогично, обращение второго и третьего членов соответствует томографии на отраженных волнах, а четвертый член миграции на основе рассеянных волн, отраженных от рефлектора [85]. Подробный вывод формулы (2.16) приведен в параграфе 2.4. Векторы источника s и приемника g определяются следующим образом:
Здесь 7S и 7ff квадратные корни: Эти корни предполагаются неотрицательными на всей области определения. Следует отметить, что поскольку вертикальная ось направлена вниз, положительное направление в пространстве Фурье образов также вниз, поэтому, когда векторы s+ и g+ направлены вверх, их вертикальные компоненты отрицательны. Рис. 2.8 показывает обыкновенные направления векторов источника и приемника вдоль путей распространения волн через центр неоднородности. источники приемники f s+\ /8+ / рефлектор
Уравнение (2.16) линейно зависит от пространственного преобразования Фурье неоднородности квадрата медленности. Следут подчеркнуть, что эта линейная зависимость специального вида выполняется только для локально-однородной вмещающей среды, позволяющей разложить волновые поля по плоским монохроматическим волнам. Если неоднородность погружена в область с изменяющейся скоростью это простое соотношение между спектрами неоднородности и волновых полей в точном виде исчезает, тем не менее, если асимптотически разложить поля можно использовать МСЧ и в среде с постоянным вертикальным градиентом скорости [64], как показано далее. В соответствии с уравнением (2.16) можно построить набор точек K в спектре неоднородности, которые оказывают влияние на амплитуду спектра рассеянного поля в произвольной заданной точке (ks,kg). Эти точки определяются выражениями:
Формула (2.18) определяет компоненты неоднородности, которые восстанавливаются с помощью миграции, в то время как уравнения (2.19) характеризуют компоненты, которые могут быть восстановлены томографией на отраженных волнах. Несмотря на то, что одной точке в спектре рассеянного поля соответствуют четыре точки в спектре неоднородности, обратное преобразование однозначно. Для вычислений параграфа 3.1, необходима информация об обратном преобразовании или выражение (к8,кд) через К = (KX,KZ) Поскольку правая часть уравнения (2.16) симметрична относительно перестановки (к3,кд)— (кд,к3), данные пары эквивалентны и различать их нет необходимости. В параграфе 2.3 показано, что существует единственная пара (к3, кд) для каждого вектора К, (К —), с точностью до перестановки к3 и ка.
Соотношение между спектром рассеянного волнового поля и спектром рассеивающей неоднородности удобно характеризовать c помощью спектральных чувствительностей. однородна, обратная задача хорошо обусловлена, если же чувствительность сильно неоднородна, решение обратной линеаризованной задачи затруднительно, что можно пояснить следующим образом: Если функционал невязки определен, как L2 норма в пространстве функций волновых чисел (KX,KZ), то квадрат чувствительности, S2(KX,KZ) пропорционален матрице Гессе для данного функционала. Гессиан диагонален в спектральной области, если неоднородность находится в локально однородной области опорной модели. Предполагается, что все данные имеют равный вес в области (KX,KZ), поэтому другое взвешивание может улучшить чувствительности и следовательно улучшить сходимость решения обратной задачи по методу МОП. Определенная уравнением (3.45) чувствительность характеризует изменение спектра наблюденных полей при изменении спектра рассеивающей неоднородности. Замечательным свойством чувствительности является отсутствие зависимости от помещаемых в модель неоднородностей, что позволяет выделить и исследовать влияние опорной среды на работу метода обращения полных волновых полей. В следующем параграфе теоретически обосновывается это свойство чувствительности, а также дается более подробный алгоритм получения диаграмм следующей главы. Если читатель более заинтересован в изучении результатов практического применения МСЧ, то следующий параграф может быть пропущен.
Спектральные чувствительности для среды в виде слоя над полупространством и роль кратных волн и волн-спутников
В этом параграфе МСЧ обобщается на случай среды с постоянным вертикальным градиентом скорости, которая зачастую используется в качестве стартовой модели для МОП ([2; 104]). В среде такого типа единственный наблюдаемый тип волн — рефрагированные волны ("diving waves"), лучи которых представляют из себя дуги окружностей. Спектральные чувствительности в данном случае зависят как от частоты зондирующего сигнала, так и от глубины залегания исследуемой неоднородности, поскольку модель является вертикально неоднородной в области неоднородности. Сравнивая диаграммы чувствительностей, можно прийти к выводу, что низкие временные частоты в наблюденных данных крайне полезны при первых итерациях МОП, однако, высокие временные частоты также содержат информацию о низких волновых числах в возмущении скорости, что дает возможность для МОП в отсутствии низких частот. При исследованиях также было выяснено, что низкие волновые числа лучше освещены в более глубоких частях модели в случае наличия неограниченных выносов в наблюденных данных.
Снова запишем волновое уравнение, описывающее давление p(r,t) как функцию координат г и времени t в акустической среде постоянной плотности здесь c(z) - скорость в среде, а /(г — го,) - функция источника. В предположении горизонтально однородной среды преобразования Фурье по времени и горизонтальной координате сводят задачу к одномерному уравнению, где p(kx,z,u) функция распределения давления в новых координатах. В случае вертикально неоднородной среды решение может быть формально записано в виде интеграла - суперпозиции квази-плоских волн [22]:
Для того чтобы обобщить МСЧ на случай неоднородной в области исследуемого возмущения скорости опорной среды предлагается сначала обратиться к обыкновенному однородному случаю. Решение в локально однородной среде может быть представлено в виде линейной комбинации двух экспонент, exp ±.iz {ш/с)2 — к2., соответствующих восходящим (+) и нисходящим (—) волновым полям. Отдельная квази-плоская волна в данном случае вырождается в обыкновенную плоскую волну:
Локально, решение уравнения 3.28 можно представить в виде A(z)exp(izB(z)), где А — амплитуда, которая изменяется медленно с изменением глубины z, а В — также меняется медленно и положительно для восходящих волн. Потому, можно определить вектор направления распространения волны на глубине z как (kx,B(z)), по аналогии с однородными плоскими волнами. Такую технику предлагается применять лишь только в ограниченной области, где располагается исследуемая неоднородность.
Точное решение в области временных частот для бесконечной акустической среды с постоянным вертикальным градиентом скорости было найдено в работе [93] в аналитическом виде для трехмерного случая и в работе [74] для двумерного случая в виде функции Лежандра. В работе [120] исследуются асимптотические свойства аналитического решения в данном типе среды. Для записи функции Грина предлагается воспользоваться обозначениями из работы [74], G(kx,Z ,OJ) = exp(—ikxxs)yz z Ily(kxz )Kly(kxz ), где z = z + со/а (3.31) для источника в точке (xs,zs) и приемника на глубине zr. Здесь, Iv и Kv — модифицированные функции Бесселя чисто мнимого порядка v = iy {ш/а) — 1/4. Скорость изменяется по закону c(z) = со + az для произвольного значения глубины z, а также введены обозначения z =mm(zs,zr) (3.32) z =max(zs,zr). (3.33)
Далее предполагается рассматривать рассеяние неоднородностью на глубине zp za, где а = s или г, а потому z = za и z zp. Асимптотические разложения модифицированных функций бесселя при больших по модулю чисто мнимых значениях и, соответствующие высоким временным частотам для источников стали темой нескольких исследовательских работ во второй половине ХХ века (см. например [43; 91]). В работе [91] выведены формулы в виде асимптотических рядов для Iv и Kv для высоких частот. Если ввести сокращенное обозначение = 2 +2, формулы Олвера
В работе [54] эти асимптотические разложения использованы для расчета коэффициента отражения от линейного переходного слоя и показано, что они совпадают в первом порядке с обыкновенным лучевым методом. К сожалению, Гупта [54] использовал формулы Олвера на границе их области применимости, что ставит заключения, сделанные в его работе, под сомнение. Разложения в по возрастающим степеням Z для модифицированных функций Бесселя широко известны [1]. Эти разложения, справедливые для произвольных чисто мнимых значений v, имеют бесконечный радиус сходимости в плоскости переменной Z и не имеют особенностей при положительных значениях Z:
Из формулы 3.31, можно заметить, что при малых значениях кх, что соответствует волнам, распространяющимся практически вертикально, аргумент Z также становится малым, а потому первый член в разложении Iv является хорошей аппроксимацией всего ряда. Потому при не слишком больших глубинах и малых углах падения можно получить: Iv(Z) (Z/2)1 = expl i/ln (Z/2) ) expl і/ln (ZQ/2) lexpl UZ/ZQ ) . (3.36) Здесь предполагается, что мы рассматриваем значения функций Бесселя в районе небольшого, но отличного от нуля значения глубины ZQ. Записывая Z = ZQ + Z, можно выделить часть решения в виде локально плоской волны, используя ряд Тейлора для степени экспоненты в выражении 3.36.
Мнимая часть v при этом определяет направление распространения волн, соответствующих решениям одномерного уравнения пропорциональным Iv. Таким образом, если Im (и) 0, то волна распространяется вверх, иначе — вниз. Первый член разложения определяет направление распространения волн для небольших значений Z, но для того чтобы проанализировать решения вплоть до точки поворота следует рассмотреть первые несколько членов ряда:
Данное приближение хорошо работает и при достаточно низких частотах для всех рассматриваемых в настоящей работе значений Z (экспоненциально затухающие волны за точкой поворота не представляют интереса), что наглядно иллюстрирует Рис. 3.16. Наиболее важным отличием формулы (3.39) от приближения методом ВКБ или асимптотики Олвера является отсутствие особенностей вблизи точки поворота Z = -. В окрестности этой точки приближения, предложенные ранее, расходятся. Другим преимуществом формулы (3.39) является сравнительная легкость анализа решения по сравнению с полной суммой ряда для модифицированной функции Бесселя Iv. К сожалению асимптотическая формула (3.39) не обладает высокой точностью при сравнительно больших значениях Z и -, что очевидно из Рис. 3.17, и Рис. 3.18.
Псевдоспектральное обращение полных волновых полей в отсутствии низких временных частот в наблюденных данных
Поскольку на каждом шаге работы метода обращения полных волновых полей необходимо несколько раз сопоставлять смоделированные и наблюденные поля, моделирование волновых полей во временной или частотной области — наиболее ресурсоемкий процесс при применении МОП. Как правило, для моделирования сейсмических полей используются конечно-разностные методы. В настоящей диссертации предлагается реализация метода обращения полных волновых полей с применением простейшего псевдоспектрального моделирования, что позволяет сократить количество необходимых для описания сейсмической среды параметров за счет использования более разреженных пространственных сеток. Теория псевдоспектрального метода моделирования рассматривается на примере акустического уравнения в среде с постоянной плотностью
В уравнениях (4.6), (4.7) не использована функция источника, но ее можно было бы добавить, вычислив свертку с функцией Грина. Функция Грина может быть найдена, как решение задачи Коши с нулевым значением поля и производной, равной дельта-функции в нулевой момент времени G(x,xo,t) = VL sin(VLt)o(x — хо). (4.8) Для нахождения функции Грина необходимо построение решения задачи для источника в виде дельта-функции в пространстве и времени, что требует некоторой аппроксимации для нее и ее производных в виде обыкновенных функций по времени, что существенно ухудшило бы точность решения. Альтернативным способом является добавление источника в уравнение (4.1) непосредственно на каждом шаге моделирования по времени, с выделением при этом аппроксимации для второй производной по времени. Решение для уравнения с источником может быть записано как сумма решения однородного уравнения с истинными начальными условиями и решения уравнения с источником с нулевыми начальными условиями. Решение уравнения с нулевыми начальными условиями может быть получено сверткой с функцией Грина [94; 109]
Формула (4.12) позволяет использовать конечные разности 4-го порядка по времени, оперируя всего тремя временными срезами поля (снэпшотами). Метод, известный как схема Лакса-Вендроффа [19; 39; 77; 80; 107] сокращает требования к используемой компьютерной памяти во время моделирования. Последнее делает данный метод предпочтительным при реализации МОП с использованием графических процессоров. Псевдоспектральный подход в данном случае сводится к конечно-разностному по времени, в то время как для пространственных производных используются производные, вычисленные с помощью преобразования Фурье. Классический конечно-разностный подход [145] требует частоты дискретизации не менее четырех точек сетки на минимальную длину волны для продолжения полей на расстояния порядка 10 длин волн [47; 48], что приводит к значительным требованиям к объему компьютерной памяти, необходимому для сохранения всех данных модели, а также накладывает ограничение Куранта-Фридрихса-Леви на длину шага по времени (CFL). Использование преобразования Фурье для вычисления пространственных производных позволяет избежать численной дисперсии при решении прямой задачи сейсмики, поскольку действие оператора Лапласа вычисляется точно для всех волновых чисел, представляющих модель. Псевдоспектральный метод можно сравнить с конечно-разностным методом бесконечного порядка или порядка минимального линейного размера сетки параметров. Однако за счет применения быстрых преобразований Фурье псевдоспектральное моделирование работает значительно быстрее, нежели конечно разностное порядка сравнимого с размером сетки параметров (соответствует по эффективности дискретным преобразованиям Фурье, но для несколько меньших массивов). Дисперсии при распространении полей не возникает вплоть до дискретизации в две точки сетки на наименьшую длину волны, что также позволяет сделать шаг дискретизации по времени более крупным. Псевдоспектральное моделирование успешно применяется для решения прямой задачи сейсмики и процедур миграции в обратном времени, однако до сих пор не было использовано для метода обращения полных волновых полей. Число необходимых для описания модели среды параметров сокращается в 36 раз в двумерных задачах и в 216 раз в трехмерных, что делает возможной реализацию метода на графических процессорах в ряде случаев без декомпозиции модели. Метод обращения полных волновых полей имеет разрешающую способность в половину локальной длины волны [113], которая в низкоскоростных частях моделей совпадает с разрешением сетки при псевдоспектральном моделировании. Данное обстоятельство приводит к полному восстановлению низкоскоростной части разреза в численных экспериментах с синтетическими данными. Отметим, что для точного описания скоростных моделей сейсмических сред иногда может потребоваться большее разрешение, поэтому псевдоспектральный подход идеально подходит для гладких сред. В случае же сред с разрывами параметров точность описания скорости определяется максимальным волновым числом, которое может быть описано данной сеткой параметров в соответствии с теоремой отсчетов. Для описания распространения волн в упругих сейсмических средах со свободными границами сложной топологии, затуханием или анизотропией могут потребоваться более сложные точные и современные методы моделирования [149; 151; 152], однако рассматриваемая в данном случае однопараметрическая акустическая инверсия позволяет использовать простое и эффективное псевдоспектральное моделирование
Классическая многомасштабная инверсия. Для улучшения сходимости инверсии предлагаются различные методы. Наиболее распространенным на данный момент является метод последовательной обработки данных на разных временных частотах, от низких до высоких [25]. На низких частотах целевой функционал имеет более плавный характер, а потому вероятность спуска в локальный минимум вместо глобального значительно ниже. Авторы указанной работы демонстрируют это на примере запаздывания сигнала в опорной модели относительно реального зарегистрированного. Расстояние между локальными минимумами функционала невязки соответствует кинематической ошибке в период регистрируемого сигнала. Чем ниже минимальная частота используемых данных, тем большие кинематические ошибки допустимы. Несмотря на многочисленные модификации, для практического применения метода обращения полных волновых полей используется следующий подход (см. [117]): сначала трендовая составляющая находится с помощью скоростного анализа, использующего только кинематические атрибуты волновых полей (времена пробега), далее распределение локальных неоднородностей определяется с помощью локальных методов минимизации целевого функционала. Основными требованиями к зарегистрированным данным для успешного применения метода в его современной форме остаются наличие низких временных частот в спектре наблюденных данных во избежание эффекта пропускания цикла, наличие сейсмических данных для значительных выносов, гарантирующих полноценное “просвечивание” сейсмическими волнами локальных неоднородностей на больших глубинах [64; 66]. Отсутствие низких временных частот в зарегистрированных трассах в некоторых случаях может быть скомпенсировано выбором функционала невязки в виде весового интеграла от корреляции смоделированного и наблюденного полей [112] или 2 нормы разности их огибающих [21].