автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное обращение времён первых вступлений для скважинных систем наблюдения в трансверсально-изотропных средах

ЧИСЛЕННОЕ ОБРАЩЕНИЕ ВРЕМЁН ПЕРВЫХ ВСТУПЛЕНИЙ ДЛЯ СКВАЖИННЫХ СИСТЕМ НАБЛЮДЕНИЯ В ТРАНСВЕРСАЛЬНО-ИЗОТРОПНЫХ СРЕДАХ

05.13.18 - математическое моделирование, численные методы и комплексы программ

2 6 Я Н В 2012

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

НОВОСИБИРСК - 2012

005009646

Работа выполнена в Федеральном государственном бюджетном образовательном учреждении высшего профессионального образования Новосибирском национальном исследовательском государственном университете и в Учреждении Российской академии наук Институте нефтегазовой геологии и геофизики им. A.A. Трофимука Сибирского отделения РАН.

Научный руководитель: доктор физико-математических наук, доцент

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

Официальные оппоненты: доктор физико-математических наук,

профессор Кабанихин Сергей Игоревич

кандидат физико-математических наук, доцент Дучков Антон Альбертович

Ведущая организация: Учреждение Российской академии наук

Институт вычислительного моделирования Сибирского отделения РАН (ИВМ СО РАН, г. Красноярск)

Защита состоится 21 февраля 2012 г. в 15 часов на заседании диссертационного совета Д 003.061.02 при Учреждении Российской академии наук Институте вычислительной математики и математической геофизики Сибирского отделения РАН (ИВМиМГ СО РАН), по адресу: пр-т Ак. Лаврентьева, 6, Новосибирск, 630090

С диссертацией можно ознакомиться в библиотеке Учреждения Российской академии наук Институте вычислительной математики и математической геофизики Сибирского отделения РАН (ИВМиМГ СО РАН).

Автореферат разослан 18 января 2012 г.

Ученый секретарь диссертационного совета доктор физико-математических наук V У 9 С.Б. Сорокин

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

Актуальность. Одним из основных этапов решения обратной кинематической задачи для уравнений теории упругости является построение лучей и вычисление времен пробега между приемниками и источниками. Известные методы двухточечного трассирования подбором начального направления луча (метод «пристрелки») или подбором траектории по минимизации времени пробега волны в сложных трехмерных высококонтрастных средах работают неустойчиво, требуют значительного времени вычисления.

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

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

Цель исследования - развитие метода и повышение достоверности результатов численного решения обратной кинематической задачи для изотропной и трансверсально-изотропной сред путём использования устойчивого двухточечного трассирования на

основе конечно-разностного решения уравнения эйконала и усечённого сингулярного разложения для томографической системы уравнений.

Научная задача - обращение времён первых вступлений для скважинных систем наблюдения в случае изотропной и трансверсально-изотропной неоднородных сред.

Этапы исследования

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

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

3. Разработка помехоустойчивого алгоритма обращения кинематических данных. Создание научно-исследовательской версии программного обеспечения.

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

Фактический материал. Научные методы исследования

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

- конечно-разностные методы решения нелинейных уравнений в частных производных;

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

- современные достижения численных методов линейной алгебры.

Для численного решения плохообусловленной системы алгебраических уравнений использовался метод усечённого сингулярного разложения, вычисляемого с помощью процедуры Арнольди, а также метод ШЬЗ-минимизации соответствующего целевого функционала в пространстве /, 5.

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

Защищаемые научные результаты.

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

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

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

Научная новизна. Личный вклад.

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

- решение двухточечной задачи трассирования сейсмических лучей с использованием конечно-разностных схем решения уравнения эйконала;

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

- на основе комбинации усеченного сингулярного разложения и обращения с использованием /|5 нормы разработан помехоустойчивый алгоритм решения обратной кинематической задачи.

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

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

На основе численного исследования пространств старших правых сингулярных векторов определён оптимальный набор

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

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

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

Результаты диссертации известны научной общественности: они докладывались и получили одобрение специалистов на международных научных конференциях «Обратные и некорректные задачи математической физики» (Новосибирск, 2007), «Математические методы в геофизике» (Новосибирск, 2008), «Days on diffraction» (Санкт-Петербург, 2009), «Inverse problems» (Китай, Ухань, 2010), международной студенческой конференции «Студент и научно-технический прогресс» (Новосибирск, 2010), на Второй молодежной международной научной школе-конференции «Теория и численные методы решения обратных и некорректных задач» (Новосибирск, 2010), международном научном конгрессе «ГЕО-Сибирь-2010» (Новосибирск, 2010), 73-й конференции EAGE (Австрия, Вена, 2011).

Полученные научные результаты полностью изложены в 12 публикациях, из которых 2 статьи - в ведущих рецензируемых научных журналах, определенных Высшей аттестационной комиссией (Технологии сейсморазведки, Физико-технические проблемы разработки полезных ископаемых), 1 статья в иностранном научном журнале (Journal of Mining Science), 9 - материалы российских и международных конференций и конгрессов.

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

Автор выражает искреннюю признательность научному руководителю д.ф.-м.н. В.А. Чеверде за всестороннюю поддержку, д.г. -м.н. В .Д. Суворову, к.т.н. С.Б. Горшкалеву, к.ф.-м.н. A.A. Дучкову, д.г,-м.н. И.Р. Оболенцевой, к.ф.-м.н. М.И. Протасову, к.ф.-м.н. Д.А. Неклюдову, к.ф.-м.н. A.M. Айзенбергу и другим коллегам по работе за обсуждения результатов, В.И.Самойловой за помощь в подготовке диссертации.

Объём и структура работы

Диссертация состоит из введения, 4 глав, заключения и списка литературы из 86 наименований. Работа содержит 133 страницы основного текста и 48 рисунков.

СОДЕРЖАНИЕ РАБОТЫ

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

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

Во второй главе изложен алгоритм двухточечного трассирования сейсмических лучей для неоднородных изотропных и трансверсально-изотропных сред на основе конечно-разностных схем решения уравнения эйконала; приведены примеры расчётов лучей в различных скоростных моделях в двумерном случае.

В разделе 2.1 даны постановка двухточечной задачи трассирования сейсмических лучей и алгоритм её решения, идея которого состоит в использовании информации, полученной с помощью конечно-разностного решения уравнения эйконала, при интегрировании системы уравнений луча:

ат

йрк = 1 Э V

V дхк '

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

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

В разделе 2.2 приведен алгоритм расчёта времён пробега сейсмических волн, разработанный с использованием конечно-разностной схемы \УЕ1ЧО-]ЖЗ. Даны примеры решения двухточечной задачи сейсмического трассирования для реалистичной сейсмогеологической модели одного из районов Северного моря в изотропном случае и при наличии двух трансверсально-изотропных слоев (рис. 1).

6=0.05 £=0.05

лучи в изотропной модели

Лучи в модели с двумя анизотропными слоями Влияние анизотропии

1000

ось Ъ

2000

6=0.15 £=0.20

— Изотроп

1000

ось X

2000,м

1000 2000, м

Рис. 1. Лучи, рассчитанные с помощью предлагаемого в работе алгоритма на основе схемы БЖ) \VENO-RK3: слева - в изотропной модели; в центре - в анизотропная модели; справа - сравнение лучей в изотропной и анизотропной модели.

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

попали в источник из-за эффекта накопления ошибок при интегрировании методом Рунге-Кутта. Приведенные примеры показывают преимущество алгоритма трассирования, предлагаемого в работе.

ось X

Рис. 2. Лучи, вычисленные с помощью алгоритма, предлагаемого в диссертации (показаны чёрным цветом), и лучи, вычисленные методом пристрелки (показаны красным).

В разделе 2.3 дан алгоритм расчёта времён пробега сейсмических волн в изотропных высококонтрастных средах с использованием конечно-разностной схемы FAST MARCHING. На рис. 3 приведены лучи в среде с высокоскоростным включением, рассчитанные с помощью предлагаемого в работе алгоритма. Скорость в соляном теле в несколько раз выше, чем во вмещающей среде. Широкий веер лучей из приёмников сходится в узкую лучевую трубку. Здесь наблюдается явление рефракции. Такие лучи затруднительно найти традиционными методами трассирования.

Рис. 3. Лучи, вычисленные с использованием алгоритма двухточечного трассирования для различного положения источника и приёмников для модели среды с высокоскоростным включением

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

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

Для изотропной вмещающей среды получено соотношение:

= ~ 1 \к К)4 + «)4 + 2«з + 2с1с'м)«г;1 )2

ось 1А

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

Приведены результаты исследования зависимости сингулярного спектра от числа N на основе сравнения между собой 500 старших сингулярных чисел и правых сингулярных векторов при увеличении N в два раза. Для чистоты и точности эксперимента рассмотрена изотропная однородная вмещающая среда, в которой трассирование и вычисление длины лучей в каждой ячейке выполнены по явным формулам и не вносят дополнительных ошибок в вычисления. Сингулярное разложение рассчитаны с помощью пакета АЯРЛСК.

Приведены результаты парного сравнения пространственной структуры правых старших сингулярных векторов и углы между соответствующими подпространствам для рассматриваемой последовательности параметров дискретизации № N=100, 200, 400, 800. По результатам сравнения определена достоверно вычисляемая часть сингулярного спектра, которую следует использовать при обращении кинематических данных для геометрии наблюдений ВСП и межскважинного просвечивания для рассмотренных параметров N.

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

Четвёртая глава посвящена обращению кинематических данных. Решение обратной кинематической задачи в линеаризованной постановке сводится к решению системы линейных алгебраических

уравнений большой размерности с использованием сингулярного разложения и минимизации в норме /, 5. Приведены примеры обращения

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

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

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

В разделе 4.2 дан пример обращения синтетических кинематических данных с использованием сингулярного разложения для реалистичной сейсмогеологической модели одного из месторождений Северного моря в случае межскважинного просвечивания. На первом этапе проводится обращение с использованием минимизации в норме /, 5 методом IRLS, что позволяет

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

Для верификации метода обращения кинематических данных был проведен численный эксперимент по обращению кинематических данных межскважинного просвечивания части сейсмогеологической модели одного из районов Северного моря (см. рис. 4).

осьХ

Рис. 4. Лучевые траектории в модели «Северное море», соответствующие синтетическим кинематическим данным, используемым в численном эксперименте.

В качестве начального приближения была взята однородная среда со скоростью 2700 м/с, это среднее значение скорости вдоль скважин. Для получения синтетических данных времен пробега в рассматриваемой части скоростной модели были рассчитаны сейсмограммы и годографы с использованием конечно-разностной схемы. Пример расчета волнового поля для одного из источников показан на рис. 5. Сейсмограммы и годограф для того же источника (источник № 50) показаны на рис. 6.

На рис.7 приведён результат обращения синтетических данных (слева) и точное решение (справа). Видно, что удаётся восстановить разлом и крупные слои.

10 м

ОСЬ Z ™

,М о 50 100 150 ZOO 250 100 350 400

ось X

Рис. 5. Слева - скоростная модель, источник №50 и приемники, справа -амплитуда волнового поля (возбужденного источником №50) в момент времени t=0,1524 с

0,0508

0,1016

0.2540 t, с.

Рис. 6. Лучевые траектории в модели «Северное море», соответствующие синтетическим кинематическим данным, используемым в численном эксперименте.

о

50 100 150 200

ОСЬ 2 250

300

350

400

450

500 ,М

Рис. 7. Результат решения обратной кинематической задачи: слева - точное решение; справа - результат обращения синтетических данных.

\

'О 50 100 150 200 250 300

ось X х т

™0СЬХ т 253 х\т 300

Рис. 8. Анизотропное возмущение (а) угольного пласта и результат его восстановления (б).

а)

е

N 6!

ось 2

40

ось 2

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

Рассматривается применение многократного просвечивания угольного пласта для обнаружения нестабильных участков анизотропного возмущения без детального изучения остальной части пласта. Приведены результаты решения модельной задачи (рис. 8)

ЗАКЛЮЧЕНИЕ

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

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

В работе рассматриваются два варианта алгоритма, основанные на разных конечно-разностных схемах решения. Первый вариант, на основе схемы WENO-RK3, позволяет решать двухточечную задачу в трансверсально-изотропной среде. Второй вариант, использующий схему FAST MARCHING, позволяет решать двухточечную задачу в высококонтрастных скоростных моделях.

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

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

Одними из направлений дальнейшего развития предложенных алгоритмов и созданного программного обеспечения является

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

ПУБЛИКАЦИИ ПО ТЕМЕ ДИССЕРТАЦИИ

1. Протасов М.И., Сердюков A.C., Чеверда В.А. Оптимальная параметризация трансверсально-изотропной среды для обращения времён первых вступлений для системы наблюдений вертикального сейсмического профилирования с выносными источниками» // Технологии сейсморазведки. 2010. №3. С.25-31.

2. Курленя М.В., Сердюков A.C., Сердюков C.B., Чеверда В.А. Локация очагов аккумулирования метана в угольном пласте сейсмическим методом // Физ.-тех. проб, разраб. полез, ископ. -2010. №6. С.39-49.

3. Kurlenya M.V., Serdyukov A.S., Serdyukov S.V. and Cheverda V.A.. Seismic approach to location of methane accumulation zones in a coal seam //Journal of Mining Science. 2010. 46, №6, p. 621-629.

4. Serdyukov A.S. Optimal parameterization in the inverse kinematic problem in anisotropic media / Abstacts of International conference "DAYS ON DIFRACTION'2009" (Saint Petersburg, May 26-29, 2009). Saint Petersburg: Universitas petropolitana MDCCXXIV. 2009, p.80.

5. Сердюков A.C. Оптимальная параметризация в обратной кинематической задаче теории упругости в случае транверсально-изотропной среды / Тезисы докладов Молодежной международ, науч. школы-конференции «Теория и численные методы решения обратных и некорректных задач» (Новосибирск, 10-20 августа 2009 г.). Новосибирск: ИМ СО РАН, 2009, с. 91.

6. Сердюков A.C. Обратная кинематическая задача в трансверсально-симметричной упругой среде / Материалы XLVIII Международной научной студенческой конференции «Студент и научно-технический прогресс» (Новосибирск, 10-14 апреля 2010 г.): Математика. Новосибирск: Новосиб. гос. ун-т, 2010, с. 221.

7. Сердюков A.C. SVD анализ линеаризованного оператора обратной задачи динамической теории упругости для анизотропных сред // Сб. мат-лов VI Международного научного конгресса «ГЕО-Сибирь-

2010» (19-29 апреля 2010г., Новосибирск). Том 2. ч. 2 «Недропользование. Горное дело. Новые направления и технологии поиска, разведки и разработки полезных ископаемых». Новосибирск: СГГА, 2010, с. 23-27.

8. Serdyukov A.S. The Linearized Traveltime Tomography in Tranversal Isotropic Elastic Media / Materials of International Conference on Inverse Problems (April 25-29, 2010, Wuhan University, Wuhan, P.R. China). Wuhan, Wuhan University, 2010, p. 21.

9. Сердюков A.C. Алгоритмы численного решения уравнения эйконала и двухточечная трассировка сейсмических лучей в сложных средах / Вторая молодежная международная научная школа-конференция «Теория и численные методы решения обратных и некорректных задач» (Новосибирск, 21-29 сентября 2010 г.). http.y/math.nsc.ru/conference/onzlO/thesis/abstracts.pdf.

10. Serdyukov A.S, Protasov M.I. Two-point Ray Tracing Algorithm for Isotropic/Anisotropic Media with Complicately Shaped High-contrast Interfaces / in «Unconventional Resources and the Role of Technology» -Proceedings of the 73rd EAGE Conference & Exhibition incorporating SPE EUROPEC 2011 (Wien, Austria, 23-26 May 2011). -Houten, Nederlalands: Co Production - CD-ROM. (paper #P048).

11. Serdyukov S.V., Lisitsa V.V, Serdyukov A.S. & Tcheverda V.A. Seismic Monitoring of the Fracture Zones of a Coal Seam / in «Unconventional Resources and the Role of Technology» - Proceedings of the 73rd EAGE Conference & Exhibition incorporating SPE EUROPEC 2011 (Wien, Austria, 23-26 May 2011). Houten, Nederlalands: Co Production - CD-ROM. (paper #P326).

12. Serdyukov S.V., Kurlenya M.V., Serdyukov A.S. and Tcheverda V.A. Location of Methane Accumulation Sites in a Coal Seam by Seismic Method / Proceedings of 45th US Rock Mechanics / Geomechanics Symposium (San Francisco, USA, June 26-29, 2011) p. 123-130.

_Технический редактор Е.В. Бекренева_

Подписано к печати 20.12.2011 Формат 60x84/16. Бумага офсет № 1. Гарнитура Тайме. _Печ. л. 0.9. Тираж 110. Заказ № 68_

ИНГГ СО РАН, ОИТ, 630090, Новосибирск, просп. Акад. Коптюга, 3

61 42-1/603

УЧРЕЖДЕНИЕ РОССИЙСКОЙ АКАДЕМИИ НАУК ИНСТИТУТ НЕФТЕГАЗОВОЙ ГЕОЛОГИИ И ГЕОФИЗИКИ

им. A.A. ТРОФИМУКА СИБИРСКОГО ОТДЕЛЕНИЯ РАН

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ НОВОСИБИРСКИЙ НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

ЧИСЛЕННОЕ ОБРАЩЕНИЕ ВРЕМЕН ПЕРВЫХ ВСТУПЛЕНИЙ ДЛЯ СКВАЖИННЫХ СИСТЕМ НАБЛЮДЕНИЯ В ТРАНСВЕРСАЛЬНО-ИЗОТРОПНЫХ СРЕДАХ

СЕРДЮКОВ Александр Сергеевич

Специальность 05.13.18 - математическое моделирование, численные методы и комплексы программ

ДИССЕРТАЦИЯ

на соискание ученой степени кандидата физико-математических наук

Научный руководитель:

доктор физико-математических наук, доцент Чеверда В. А.

НОВОСИБИРСК-2011

Оглавление

ВВЕДЕНИЕ...............................................................................4

Глава 1. ИЗУЧЕННОСТЬ ЧИСЛЕННЫХ МЕТОДОВ РЕШЕНИЯ

ОБРАТНОЙ КИНЕМАТИЧЕСКОЙ ЗАДАЧИ...................10

1.1 Краткий анализ основных работ в области сейсмической томографии...............................................................10

1.2 Анализ современных численных методов решения прямой и обратной кинематических задач......................................13

Глава 2. АЛГОРИТМ ДВУХТОЧЕЧНОГО ТРАССИРОВАНИЯ СЕЙСМИЧЕСКИХ ЛУЧЕЙ НА ОСНОВЕ КОНЕЧНО-РАЗНОСТНОГО РЕШЕНИЯ УРАВНЕНИЯ ЭЙКОНАЛА.....................................................................................31

2.1 Алгоритм решения двухточечной задачи...........................31

2.2 Трассировка в трансверсально-изотропной неоднородной среде на основе конечно-разностной схемы WENO-RK3.......35

2.3 Трассировка в высококонтрастной среде на основе схемы

Fast Marching..............................................................48

Глава 3. ИССЛЕДОВАНИЕ СВОЙСТВ ЛИНЕАРИЗОВАННОГО ТОМОГРАФИЧЕСКОГО ОПЕРАТОРА С ИСПОЛЬЗОВАНИЕМ ЕГО СИНГУЛЯРНОГО РАЗЛОЖЕНИЯ............................................................55

3.1 Вывод интегральных соотношений линеаризованной

задачи сейсмической томографии...............................................55

3.2 Исследование сходимости матричной аппроксимации томографического оператора.........................................60

3.3 Оптимальный набор восстанавливаемых параметров при решении обратной кинематической задачи для трансверсально-изотропной среды...................................74

Глава 4. АЛГОРИТМ УСТОЙЧИВОГО ОБРАЩЕНИЯ КИНЕМАТИЧЕСКИХ ДАННЫХ И ЕГО ПРИМЕНЕНИЕ.............................................................87

4.1 Устойчивое обращение кинематических данных в норме

/15 методом IRLS.........................................................87

4.2 Пример обращения синтетических данных в случае межскважинного просвечивания......................................98

4.3 Локация зон аккумулирования метана в угольном пласте сейсмотомографическим методом...................................106

ЗАКЛЮЧЕНИЕ.........................................................................................121

СПИСОК ЛИТЕРАТУРЫ..........................................................................123

ВВЕДЕНИЕ

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

Актуальность. Одним из основных этапов решения обратной кинематической задачи для уравнений теории упругости является построение лучей и вычисление времен пробега между приемниками и источниками. Известные методы двухточечного трассирования подбором начального направления луча (метод «пристрелки») или подбором траектории по минимизации времени пробега волны в сложных трехмерных высококонтрастных средах работают неустойчиво, требуют значительного времени вычисления.

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

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

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

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

Научная задача - обращение времён первых вступлений для скважинных систем наблюдения в случае изотропной и трансверсально-изотропной неоднородных сред.

Этапы исследования

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

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

3. Разработка помехоустойчивого алгоритма обращения кинематических данных. Создание научно-исследовательской версии программного обеспечения.

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

Фактический материал. Научные методы исследования

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

конечно-разностные методы решения нелинейных уравнений в частных производных;

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

- современные достижения численных методов линейной алгебры.

Для численного решения плохообусловленной системы алгебраических уравнений использовался метод усечённого сингулярного разложения, вычисляемого с помощью процедуры Арнольди, а также метод ШЬ8-минимизации соответствующего целевого функционала в пространстве /, 5.

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

Защищаемые научные результаты.

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

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

влиянием параметров при численном обращении времён первых вступлений.

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

Научная новизна. Личный вклад.

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

- решение двухточечной задачи трассирования сейсмических лучей с использованием конечно-разностных схем решения уравнения эйконала;

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

- на основе комбинации усеченного сингулярного разложения и обращения с использованием /]5 нормы разработан помехоустойчивый алгоритм решения обратной кинематической задачи.

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

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

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

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

Результаты диссертации известны научной общественности: они докладывались и получили одобрение специалистов на международных научных конференциях «Обратные и некорректные задачи математической физики» (Новосибирск, 2007), «Математические методы в геофизике» (Новосибирск, 2008), «Days on diffraction» (Санкт-Петербург, 2009), «Inverse problems» (Китай, Ухань, 2010), международной студенческой конференции «Студент и научно-технический прогресс» (Новосибирск, 2010), на Второй молодежной международной научной школе-конференции «Теория и численные методы решения обратных и некорректных задач» (Новосибирск, 2010), международном научном конгрессе «ГЕО-Сибирь-2010» (Новосибирск, 2010), 73-й конференции EAGE (Австрия, Вена, 2011).

Полученные научные результаты полностью изложены в 12

публикациях, из которых 2 статьи - в ведущих рецензируемых научных

журналах, определенных Высшей аттестационной комиссией (Технологии

сейсморазведки, Физико-технические проблемы разработки полезных

8

ископаемых), 1 статья в иностранном научном журнале (Journal of Mining Science), 9 - материалы российских и международных конференций и конгрессов.

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

Автор выражает искреннюю признательность научному руководителю д.ф.-м.н. В.А. Чеверде за всестороннюю поддержку, д.г. -м.н. В.Д. Суворову, к.т.н. С.Б. Горшкалеву, к.ф.-м.н. A.A. Дучкову, д.г.-м.н. И.Р. Оболенцевой, к.ф.-м.н. М.И. Протасову, к.ф.-м.н. Д.А. Неклюдову, к.ф.-м.н. A.M. Айзенбергу и другим коллегам по работе за обсуждения результатов, В.И.Самойловой за помощь в подготовке диссертации.

Глава 1

ИЗУЧЕННОСТЬ ЧИСЛЕННЫХ МЕТОДОВ РЕШЕНИЯ ОБРАТНОЙ КИНЕМАТИЧЕСКОЙ ЗАДАЧИ

1.1 Краткий анализ основных работ в области сейсмической томографии

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

Кинематический подход исторически является более ранним. Исходными данными для обратной кинематической задачи являются времена пробега упругих волн из источников в приёмники, которые могут быть определены с высокой точностью. Первые результаты восстановления скоростного строения по временам пробега упругих волн от землетрясений получены в начале двадцатого века Герглотцем и Вихертом [Herglotz, 1907; Wiechert, Zoeppritz, 1907]. Обратная кинематическая задача решалась аналитическими методами. Так, в работе Шлихтера [Schlichter, 1932] рассмотрена задача восстановления скоростного строения вертикально- неоднородной среды.

Дальнейшие усилия долгое время были направлены на развитие методов решения одномерной задачи. В частности, характер неоднозначности решения обратной кинематической задачи изучен M.JI. Гервером и В.М. Марушкевичем [Гервер, Марушкевич, 1965]. Одними из первых являются результаты для двумерной постановки обратной кинематической задачи, изложенные М.М.Лаврентьевым и В.Г.Романовым [Лаврентьев, Романов, 1966]. Основополагающими в области обратных кинематических задач являются работы А.С.Алексеева, М.М. Лаврентьева, В.Г.Романова и др., посвященные исследованию кинематической задачи в линеаризованной постановке [Алексеев и др., 1969; Алексеев и др., 1979].

В современной научной литературе обратную кинематическую задачу часто называют задачей сеисмическои томографии. Компьютерной томографией называется численное восстановление функций по их криволинейным или поверхностным интегралам, которое находит применение в различных областях науки и техники [Натеррер, 1990]. Впервые о компьютерной томографии заговорили в начале 70-х годов прошлого столетия в связи с рентгенодиагностикой в медицине [Натеррер, 1990; Ganzkörper - Computer - Tomographie, 1977].

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

В современной геофизике обратная задача кинематики решается численно. Применяется стандартная техника поиска минимума целевого функционала, представляющего собой норму разности измеренных данных и данных, рассчитанных с помощью решения прямой задачи для текущей модели среды. Первые работы в этом направлении появились в конце 70-х - начале 80-х годов [Dines, Lytle, 1979; Pederson et al., 1985].

Известно, что многие горные породы обладают анизотропными свойствами, и скорость распространения сейсмических волн зависит от направления [Thomsen, 1986]. Закон Гука, уравнения, описывающие распространение волн, законы преломления, трассировка лучей в анизотропной среде существенно отличаются от изотропного случая. В частности, широко распространённые глинистые сланцы обладают анизотропными свойствами вследствие своей микроструктуры, состоящей из тонких слоев [Sayers, 2005]. Распространение сейсмических волн в сланцах может быть описано в рамках модели трансверсально-симметричной среды с вертикальной осью симметрии. Методика построения эффективных трансверсально-изотропных моделей, описывающих тонкослоистые среды, представлена в работах Backus [1962] и Berryman [1979]. Модели для глин и сланцев рассматриваются в [Sayers, 2005] и [Xu, White, 1995]. Методы введения анизотропии для трещиноватых сред разработаны Hudson [1981] и Thomsen [1995].

Традиционно обработка сейсмических данных проводится в рамках модели изотропной среды. Для небольших выносов (меньше 30°) и незначительной анизотропии "изотропные" методы работают достаточно хорошо, но в ряде случаев неучёт анизотропии приводит к значительным ошибкам. В работе Levin [1990], Гурвич [1940] проведены модельные исследования ошибки, возникающей при неучёте анизотропии в методе отражённых волн. Задача томографии в анизотропной среде изучена слабо

и недостаточно представлена в литературе.

12

1.2 Анализ современных численных методов решения прямой и обратной кинематических задач

Известны три основных способа численного решения задачи томографии: методы матричного обращения, методы преобразования Фурье, алгебраические про�