Динамический подход к обращению временных невязок

Краткое описание

ОГЛАВЛЕНИЕ

ВВЕДЕНИЕ...3

ГЛАВА 1. ИЗУЧЕННОСТЬ ВОПРОСА...8

ГЛАВА 2. ДИНАМИЧЕСКИ ОБОСНОВАННЫЙ КИНЕМАТИЧЕСКИЙ

ОПЕРАТОР (КО)...14

2.1. Построение динамически обоснованного кинематического оператора...14

2.2. Сопряженный КО. Миграция временных невязок...28

2.3. КО в случае вертикально-неоднородной вмещающей среды ...41

ГЛАВА 3. СОВМЕСТНОЕ ОБРАЩЕНИЕ КИНЕМАТИЧЕСКОГО И

ДИНАМИЧЕСКОГО ОПЕРАТОРОВ В СЛУЧАЕ ВЕРТИКАЛЬНО -НЕОДНОРОДНОЙ ВМЕЩАЮЩЕЙ СРЕДЫ...52

3.1. Динамический оператор для вертикально-неоднородной вмещающей среды...52

3.2. Составной оператор для вертикально - неоднородной

вмещающей среды...'...56

3.3. Численная реализация совместного обращения...58

ЗАКЛЮЧЕНИЕ...,...65

ЛИТЕРАТУРА...67
Введение

ВВЕДЕНИЕ

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

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

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

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

Основные этапы исследований:

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

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

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

Методы исследований и фактический материал

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

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

Основным фактическим материалом являются результаты обработки данных численного моделирования процессов распространения сейсмических волн. Для опробования алгоритмов и ключевых процедур использовались как синтетические данные, так и данные реальных сейсмических наблюдений по методу вертикального сейсмического профилирования, полученные в Западной Сибири (предоставлены СБ. Горшкалевым в 2001 г.).

Сформулированы и защищаются следующие научные результаты

1. Линеаризованный кинематический оператор

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

- сводится к линейному интегральному уравнению Фредгольма первого рода относительно возмущения скорости.

2. Решение обратной кинематической задачи

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

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

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

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

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

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

Новизна работы. Личный вклад.

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

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

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

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

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

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

Апробация работы и публикации

Основные положения и результаты докладывались на международной конференции молодых ученых, специалистов и студентов "Геофизика - 2001" (Новосибирск, 2001), международной конференции "Математические методы в геофизике" (Новосибирск, 2003).

Результаты исследований по теме диссертации изложены в 4 опубликованных работах.

Диссертация выполнена в Лаборатории динамических проблем сейсмики Института геофизики СО РАН.

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

ГЛАВА 1 ИЗУЧЕННОСТЬ ВОПРОСА

Предлагаемый подход к решению линеаризованной обратной задачи является естественным продолжением методов сейсмической томографии, в случае, когда возникает необходимость учитывать спектральную ограниченность сейсмических сигналов. Традиционная кинематическая томография основывается на использовании лучевой теории для расчета времен пробега между источником и приемником. Основным преимуществом такого подхода к вычислению времен пробега является тот факт, что лучи вычисляются очень быстро (по сравнению с расчетом полного волнового поля) и, следовательно, открывается возможность изучения строения трехмерных структур на современной вычислительной технике. Кроме того, теория лучевого метода весьма существенно разработана (см., например, Ccrveny V., et. al. [28]) и имеется широкий набор алгоритмов и программных комплексов. Одним из основных постулатов лучевого метода является использование сигналов с бесконечно высокой частотой. Именно для них имеет смысл понятие луча и именно для них справедливо выражение для времени распространения волны вдоль луча - интеграл от величины, обратной скорости, вычисленный вдоль луча, соединяющего источник с приемником.

Различным разновидностям сейсмической томографии (межскважинной, отражательной, телесейсмической) посвящено большое количество различных публикаций. Подробную библиографию иностранных источников можно найти в [17]. Из отечественных публикаций необходимо отметить теоретические работы СВ. Гольдина [6,7,8], в которых строится континуальная теория лучевой сейсмической томографии, рассматриваются условия, обеспечивающие единственность решения в межскважинной и отражательной томографии. Работы СМ. Зеркаля и др. [9,14] посвящены обратной кинематической задаче с неполными данными. Показано, что она сводится к угловому преобразованию Радона и может решаться по известным явным формулам только в случае, когда используется фоновая модель с линейным нарастанием скорости по глубине.

8

Особый интерес для нас в дальнейшем будет представлять томография на отраженных волнах. Результаты изучения обратных кинематических задач отражательной сейсмики для ряда частных моделей можно найти в монографии СВ. Гольдина [5].

Принципиальный момент в томографии на отраженных волнах заключается в том, что неизвестным является не только распределение скоростей в среде, но и глубина до отражающей границы. В связи с этим выделяются три основных подхода (Williamson, P.R [62]). 1). Предполагается, что глубина до отражающей границы и ее форма известны, и необходимо найти только скорость в покрывающей среде (Neumann, G. [46]). Этот подход применим в достаточно простых геологических ситуациях. 2). Обращение невязок времен производится таким образом, чтобы одновременно определить и глубину до отражателя, и скорость в среде. Для этого вводится некоторая подходящая единообразная параметризация отражающих границ и скоростей (например, с помощью сплайнов (Мао, W., Stuart G.W. [41]) ) см. также Bishop, Т. ct. al. [24], Farra, V., Madariaga, R. [35]. Однако, как показано в Lines, L [39], Ticman, H.J. [57], невозможно однозначно определить и параметры границы (форма и положение), и скорость в среде без априорной информации о положении границы. 3). Итеративный алгоритм, который заключается в последовательном применении миграционной процедуры на основе некоторой текущей скоростной модели для определения положения отражающей границы и вычислении поправки к скоростной модели с помощью томографии, обращая невязки времен пробега волны, отраженной от найденной при миграции границы Mora, P. [44], Stork, С, Clayton, R.W. [55]. Недостатком данного подхода является то, что миграция в латерально - неоднородную среду (особенно 3-D) требует значительных затрат машинного времени, тем более, что здесь приходится применять ее многократно.

Во всех вышеупомянутых подходах совершенно не учитывается тот факт, что реальные сейсмические сигналы имеют ограниченный спектр, и что на времена распространения волны между заданной парой источник - приемник будет влиять скоростное строение в некоторой окрестности луча. Чтобы определить эту область среды, в работе Ю.Л. Кравцова, Ю.И. Орлова [13]

вводится понятие фрснелевского объема (или "физического" луча, Ccrvcny, V., Soares J.E. [29]) . По сути, фрепелевский объем - совокупность первых зон Френеля, рассчитанных в каждой точке заданного луча. С их помощью была глубже изучена разрешающая способность сейсмических методов, определены условия применимости лучевого метода (Клем-Мусатов К.Д. [11], Ю.А. Кравцов, Ю.И. Орлов [13], Pant, D.R., Grccnhalgh, S.A. [49], Williamson, P.R., Worthington М.Н. [63]). В работе [29] предложен алгоритм быстрого вычисления фреиелевских объемов, сопоставимый по затратам машинного времени с алгоритмами трассировки лучей. Это позволяет построить на их основе эффективную процедуру обращения временных невязок, учитывающую влияние возмущений скорости, расположенных внутри френелевского объема. (Vasco, D.W. et. al.[60] (рассматривается на примере межскважииного просвечивания)). Следует отметить, что фрепелевские объемы точно определяются лишь для монохроматических сигналов (или сигналов с узким спектром). Существует другой способ для того, чтобы учесть влияние конечного объема среды на времена пробега, который заключается в расчете так называемых функций чувствительности (ядер чувствительности) (Li, X.-D., Tanimoto, T. [38], Stark, P., Nikolayev, D. [54], Vasco, D.W., Majer E.L. [59], Woodward, M. J. [64]). В этих работах времена пробега определяются как функционалы от параметров среды по некоторому объему. Далее проводится процедура линеаризации. Невязки времен пробега определяются как аргумент максимума функции взаимной корреляции волн в референтной и реальной средах. Воспользовавшись борцовским приближением, получают линейное соотношение между возмущением среды и невязками времен. Ядра этих интегральных уравнений и называются функциями чувствительности (фактически, функции чувствительности представляют собой производные Фреше). Функции чувствительности для трехмерных сред строятся в работах Dahlen, F.A. et. al. [32,33], Marquering, H. et. al. [42,43], Zhao, L. [66]. Отметим интересный факт: как показано в [42,43] возмущения скорости вдоль луча не вызывают задержки (используется частотный диапазон 0.01-0.06 Hz.). В работе Vasco, D.W., Majer E.L.[59] предложен алгоритм обращения временных невязок на основе

10

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

Функции чувствительности более правильно описывают зависимость невязок времен от скоростных возмущений, расположенных в окрестности луча, однако их вычисление занимает гораздо больше машинного времени, чем вычисление френелевских объемов (необходимо два раза решать прямую динамическую задачу). Как показано на примере задач межскважинного просвечивания в Vasco, D.W. ct. al. [60], при определенных условиях результаты обращения невязок времен, основанные на использовании френелевских объемов и функций чувствительности достаточно похожи, что делает первый подход предпочтительнее (если используются достаточно высокочастотные сигналы).

В работе Luo, Y., Schuster, G.T. [40] предложен альтернативный подход к сейсмической томографии, совмещающий обращение невязок времен и обращения полного волнового поля. Он заключается в минимизации градиентным методом целевой функции, которая представляет собой сумму квадратов временных невязок по всем положениям источников и приемников. Невязки определяются через функцию взаимной корреляции и могут быть вычислены непосредственно по сейсмограммам. В работе показывается, что при стремлении частоты зондирующего сигнала к бесконечности все сводится к традиционной лучевой томографии. В случае, когда начальная модель близка к точной, предлагаемый подход эквивалентен линеаризованному обращению полного волнового поля. Однако в общем случае он позволяет получать лишь сглаженное изображение среды.

Наилучшей разрешающей способностью обладают методы обращения полного волнового поля, т.е. методы, основанные на решении обратной динамической задачи. Этой тематике посвящена достаточно обширная литература (см. например [1,20,26,27,30,31,34,53,65]). Общим недостатком всех

11

линеаризованных подходов, основанных на борновском приближении является то, что они мало применимы в ситуациях с резким перепадом модельных параметров. Кроме того, как это было установлено на практике, а затем обосновано в Clement F. ct. al. [31] при анализе сингулярного спектра динамического оператора, в результате обращения полного волнового поля принципиально невозможно устойчиво восстановить треидовую составляющую скоростного разреза.

Нелинейные методы обращения (см. например Tarantola, A. [56]) в большинстве своем представляют собой последовательность обращений, на основе борновского приближения - проблема пропагатора актуальна и для них. В этом случае существует нелинейная зависимость целевой функции (как правило, квадрат нормы (Z,2) невязок волновых полей) от параметров модели (или, по мнению Clement F. et. al. [31], Mora, P. [44] от возмущений пропагатора). Это говорит о том, что сходимость итерационной процедуры к удовлетворительному решению возможна лишь в том случае, когда начальная модель достаточно близка к реальной.

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

В работе Neele, F. ct. al. [45] для повышения детальности изображений, полученных при обращении временных невязок, используются амплитуды волн, которые чувствительны к границам раздела и представляют помимо времен пробега еще один независимый набор данных. Амплитуды рассчитываются модифицированным лучевым методом. Приводится пример применения такого подхода для томографических исследований верхней мантии.

Развитие идеи, предложенной в [40], приводится в работе Zhou, С. et. al. [67]. Здесь изображение, полученное при обращении невязок времен используется как референтная модель для обращения волнового поля, причем обе процедуры очень схожи и требуют двукратного решения прямой

12

динамической задачи. В результате строится итеративный процесс, который позволяет получать весьма подробные изображения среды, что показано на ряде примеров для задачи межскважинного просвечивания. Работа Zhou, С. ct. al. [68] обобщает этот подход для уравнений упругости.

Предлагаемый в данной диссертации подход к обращению временных невязок предназначен прежде всего для обработки данных поверхностных систем многократного перекрытия, хотя, конечно, может применяться и для любых других сейсмических методов. Его оригинальность заключается в том, что в случае поверхностной системы наблюдений, когда используется вертикально-неоднородная вмещающая модель, исходное интегральное уравнение сводится к распадающейся системе одномерных интегральных уравнений. При этом становится возможным совместное обращение полученного кинематического оператора с динамическим, рассматриваемым в работах [2,31]. В данной работе рассматривается 2-D постановка. Переход к полной 3-D постановке не вносит существенных изменений, кроме более громоздких технических выкладок и значительно более ресурсоемких вычислений.

13

ГЛАВА 2

ДИНАМИЧЕСКИ ОБОСНОВАННЫЙ КИНЕМАТИЧЕСКИЙ

ОПЕРАТОР

2.1. Построение динамически обоснованного кинематического оператора

Динамический подход к определению невязок времен пробега волн.

Пусть полуплоскость z?0 заполнена средой с неизвестной скоростью v{x,z), которая представляется в виде:

v(x,z) = vQ(x,z) + Sv(x,z), | 5' /

juo(rg;rs;t)---------^-------Л

Л В таком виде выражение для невязок времен пробега получено в [43].

Построение кинематического оператора

Для последующих выкладок нам будет необходимо переписать (2.1.11) в спектральной области. Но прежде сделаем несколько предположений. Если используется волна, отраженная от некоторой известной границы, будем считать, что vo(r)- гладкая функция, не имеющая разрывов вплоть до этой границы. В случае, когда используется прямая волна (наблюдения в скважине) будем считать v0(г) гладкой. Далее, полагаем, что носитель Sv(r) содержится в

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

оо

jiw ¦ SU(rg ;rs;w)-U0 (rg \rs;w) dw

16

Формула (2.1.12) устанавливает связь между невязкой времени пробега волны из точки rs в точку rg с возмущением волнового поля в точке rg. Чтобы

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

1

(2.1.13)

и

ди

obs

obs\,-0

dt

= 0, и

obs\z=o

= 0

/=0

Л волновое поле мо(г,г,;О, рассчитанное для фоновой модели vo(r)

dt2 диг

(2.1.13а)

ur

'0|