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

На правах рукописи

СИЛЬВЕСТРОВ ИЛЬЯ ЮРЬЕВИЧ

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ В ОБРАТНОЙ ДИНАМИЧЕСКОЙ ЗАДАЧЕ ВЕРТИКАЛЬНОГО СЕЙСМИЧЕСКОГО ПРОФИЛИРОВАНИЯ С ВЫНОСНЫМ ИСТОЧНИКОМ

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

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

АВТОРЕФЕРАТ

НОВОСИБИРСК 2008

003168893

Работа выполнена в Институте нефтегазовой геологии и геофизики им А А Трофимука Сибирского отделения Российской академии наук

Научный руководитель Официальные оппоненты

Ведущая организация

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

доктор физико-математических наук, профессор Крукиер Лев Абрамович

кандидат физико-математических наук, доцент Белоносов Андрей Сергеевич

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

Защита состоится 10 июня 2008 г в 17 часов на заседании диссертационного совета ДМ 003 046 01 по защите диссертаций на соискание ученой степени доктора наук при Институте вычислительных технологий СО РАН по адресу 630090, Новосибирск, проспект Академика М А Лаврентьева, 6 (с!$оуе1@1с1 пбс ги)

С диссертацией можно ознакомиться в специализированном читальном зале вычислительной математики и информатики ГПНТБ СО РАН

Автореферат разослан 8 мая 2008 г

И о ученого секретаря диссертационного совета доктор технических наук, профессор ' АД Рычков

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

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

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

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

ной функции Для решения рассматриваемой задачи в рамках двумерных уравнений изотропной теории упругости такой подход был впервые применен совсем недавно в работах М A Roberts (2006) и М A Roberts и В Е Hornby (2007), которые к сожалению, имеют исключительно практический характер в силу того, что метод был применен в первом случае для сложных синтетических данных, модель которых автор не приводит, а во втором случае - для реальных данных Поэтому оценить его эффективность можно только по косвенным признакам

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

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

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

Основные этапы исследования.

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

2 Выполнить численный анализ сингулярного разложения полученного линейного оператора и на его основе установить структуру ре-

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

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

4 Провести представительную серию численных экспериментов по решению обратной задачи на основе синтетических данных в средах различной степени сложности

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

Изучение свойства решения обратной задачи, получаемого методом Ньютона, происходило в рамках теории условно-корректных задач При этом существенно использовалось обобщение понятия г-решения, возникающее в теории вычислительной линейной алгебры на случай компактных операторов в гильбертовых пространствах, разработанное в работах В И Костина иВА Чеверды (1995, 1998)

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

Фактическим материалом для тестирования разработанного алгоритма обращения являлись синтетические данные для сред различной степени сложности (с одиночным рассеивателем, с горизонтальным слоем, с наклонным слоем, с несколькими наклонными слоями)

Результаты исследования свойств решений рассматриваемой обратной задачи, получаемые с использованием анализа сингулярного разложения линеаризованного оператора динамической теории упругости, сравнивались с результатами A Tarantola, основанными на исследовании диаграмм рассеивания от точечных объектов, а также с результатами D Lebrun и A Nicolao, полученными для других систем наблюдений

На защиту выносятся

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

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

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

Новизна работы.

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

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

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

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

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

только первого шага процесса Ньютона, то есть линейного обращения в поставленной задаче, - плотность оказывается "почти" ортогональной весьма большому числу старших векторов и поэтому не может быть определена при достигаемых на практике точностях 2 Исходя из установленной структуры численного решения обратной задачи разработан, программно реализован и протестирован оригинальный алгоритм определения местоположений и амплитуд разрывов (в линейном приближении) упругих импедансов среды на основе итерационного метода с использованием при моделировании волновых процессов конечно-разностной схемы Вирье на сдвинутых сетках с ограничением расчетной области при помощи идеально согласованных слоев (РМЬ) Личный вклад.

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

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

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

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

Научные результаты диссертации докладывались на Международной конференции, посвященной 75-летию академика М М Лаврентьева "Обратные и некорректные задачи математической физики" (Новосибирск, 2007), VII Международной конференции "Галь-перинские чтения-2007", (Москва, 2007), 69-ой Международной конференции европейской ассоциации геофизиков и инженеров "69th EAGE Conference & Exhibition" (Лондон, Великобритания, 2007), III Всероссийской конференции "Актуальные проблемы прикладной математики и механика" (Абрау-Дюрсо, 2006), Международной конференции, организованной ассоциациями геофизиков ОЕАГО, EAGE, SEG "Geosciences -То Discover and Develop" (Санкт - Петербург, 2006), V Международной научно-практической геолого-геофизической конференции-конкурсе молодых ученых и специалистов Геофизика-2005 (Санкт - Петербург, 2005), XI Всероссийской школе-семинаре "Современные проблемы математического моделирования" (Абрау-Дюрсо, 2005), семинарах в Институте вычислительной математики и математической геофизики СО РАН, Институте математики СО РАН, Институте вычислительных технологий СО РАН, Институте вычислительного моделирования СО РАН Публикации.

Результаты исследования по теме диссертации изложены в 8 опубликованных работах Из них - 2 статьи в ведущих российских рецензируемых журналах по Перечню ВАК, 1 работа - в трудах российской конференции, 2 работы - в расширенных тезисах докладов международных конференций, 3 работы - в тезисах российских и международных конференций

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

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

Автор выражает благодарность всем сотрудникам Лаборатории за полезные научные дискуссии в ходе выполнения работы Успешному выполнению работы во многом способствовали ценные и критические замечания В И Костина Хочется поблагодарить С Б Горшкалева и

Г В Решетову за то, что они взялись за труд ознакомиться с работой и высказать о ней свое мнение

Необходимо отметить неоценимую помощь, оказанную В И Самойловой при работе над рукописью

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

Объем н структура работы.

Диссертация состоит из введения, трех глав, заключения и двух приложений Общий объем диссертации составляет 90 страниц, включая 28 рисунков Список используемой литературы содержит 60 наименований

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

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

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

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

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

Ши + 2^1 - /-Г, А + ^ = ПфСх _ 5д

аГ ах1 ^ ах^) вх2 ^ дх2 ах1 ) ах1

Э2н, Э ( ,Эи. Эы,ч| Э [,, _ „ Эм, | . Э Г - — + Т-) "Г Мт + Г "

вг дх, ^ дх2 ох, ) ск, ( дх2) ах2

«,] =0, и, I = 0,

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

«и = 0<1<Т,х2е [хи, ,хт К < < хт },

о режиме колебаний, записанной в приемниках, расположенных в вертикальной скважине, восстановить параметры среды Л,/и,р ниже забоя скважины, то есть ниже линии, на которой располагается последний приемник Эту задачу можно рассматривать как нелинейное операторное уравнение

В <т>= йоЫ,

где т = {Л,р,р)Т - параметры среды, /7"'" - данные наблюдений (сейсмограммы), В - нелинейный оператор, неявно определяемый уравнениями теории упругости Для его решения предлагается использовать итерационный метод Ньютона

ОВ{тк) < тк+] -щ >=моМ -В{щ) (1)

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

Рис 1 Геометрия наблюдений, используемая для анализа сингулярного разложения

Для выяснения ответов на эти вопросы в работе используется подход к анализу и построению численных методов решения линейных операторных уравнений с компактным оператором в гильбертовых пространствах, основанный на применении анализа его сингулярного разложения (в дальнейшем 5УО-анализа) Следуя работам В И Костина, В А Чеверды и др (2004, 2005 гг) предлагается анализировать свойства решения с помощью усечения сингулярного разложения Эта процедура заключается в том, что решение х системы линеиных алгеораических уравнений Ах = у, полученное методом наименьших квадратов, приближается вектором * (г- решением), который является проекцией

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

■120 2000 4000

Номер сингулярлм о числа, п

В разделе 2 2 с учетом этих предположений выводится явный вид оператора ЭВ С этой целью используется приближение однократного рассеяния и явный вид функции Грина рассматриваемой задачи в однородной среде Оператор БВ оказывается матричным интегральным

2000

4000

6000

Рис 2 Сингулярные числа оператора ОВ в логарифмической

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

с

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

В разделе 2 3 выполняется анализ сингулярного разложения построенной матрицы для геометрии наблюдений, изображенной на Рис 1 Получающиеся в данном случае сингулярные числа представлены на Рис 2 Видно, что они действительно быстро стремятся к нулю, что обуславливается компактностью исходного оператора Считая допустимым число обусловленности равное 104 видно, что только проекция на 1500 старших сингулярных векторов может быть восстановлена с заданной точностью

Одним из наиболее принципиальных вопросов, возникающих при решении обратной динамической задачи сейсмики, является выбор оптимальной параметризации среды Наиболее распространенными наборами параметров для описания упругих сред являются, в дополнение к плост-ности (р), параметры Ламе (Я,//), скорости продольных и поперечных

волн (УР, У3) и упругие импедансы (1Р = рУр, /5 = рУ3) Известно, что эти

параметры не эквивалентны при решении обратной задачи

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

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

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

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

' I--'' / V

500 1 ООО 1500 2DOO 2500

Число старших сингулярных векторов. Р- импеданс

1500 2000 2500

Число старших сингулярных векторов. Плотность (в базисе импедансов)

500 1000 1500 2000 2500

Число старших сингулярных векторов.

500 1000 1500 2000 2500

Число старших сингулярных векторов.

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

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

Для решения линейного уравнения (1) применяется итерационный метод LSQR, предложенный в работе С.С. Paige (1982). Аналитически последовательность этого метода эквивалента последовательности сопряженных градиентов для уравнения, приведенного к нормальной форме. Как и для всех методов подобного типа, для реализации LSQR необходимо только знание действия операторов DB и DB1 на векторы

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

чае роль регуляризирующего алгоритма с параметром регуляризации, равным числу итераций.

, 0

Д)

100 200 300 -5

кг/(л?-с)* М'

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

Для моделирования волновых полей в работе используется модификация конечно-разностной схемы второго порядка на сдвинутых сетках, предложенной в работе ]. Утеих (1986). В качестве неотражающих условий для ограничения расчетной области применяются идеально согласованные поглощающие слои (РМЬ), предложенные в работе .1.Р. Вегег^ег (1994). Программный код, реализующий данную схему, был разработан в Лаборатории вычислительных методов геофизики ИНГГ СО РАН и модифицирован соискателем для использования в решении обратной задачи. Все программное обеспечение, созданное в ходе выполнения, работы написано на языке программирования С++. Раздел 3.3. посвящен тонкостям реализации алгоритма и разработанному программному обеспечению.

В разделе 3.4 приводятся результаты решения обратной задачи по синтетическим данным для сред различной степени сложности. В каче-

стве первого теста использовалась среда с одним точечным рассеивате-лем. Было показано, что с числом итерации метода Ь8<ЗГ£ локализуется в пространстве, а так же подавляются ошибки в решении (известные в сейсморазведке как артефакты миграции), связанные с ограниченностью системы наблюдений и с наличием в данных обменных волн. В качестве следующего тестового примера была взята среда с двумя горизонтальными отражающими границами. Результат решения обратной задачи приведен на Рис. 4. Видно, что местоположение (по времени) и амплитуды разрывов упругих импедансов восстановлены достаточно точно. Однако форма полученного решения осложнена формой импульса зондирующего сигнала. Это было предварительно предсказано с использованием 5УО-анализа.

Удаление от скважины, жфг-г) х Ю8 юЦлг с) х.10®

метры

Рис. 5. Результат решения обратной задачи для среды с наклонными ТОНКОСЛОИСТЫМИ пропластками а) истинная среда; б) восстановленные границы разрыва Р-импеданса; в) амплитуды истинного (сплошная линия) и восстановленного разрыва Р импеданса на первой границе вдоль вертикальной линии, соответствующей 225м по горизонтали; г) амплитуды истинного (сплошная линия) и восстановленного разрыва Р импеданса на второй границе вдоль вертикальной линии, соответствующей 300м по горизонтали.

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

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

ЗАКЛЮЧЕНИЕ

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

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

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

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

Путем решения обратной задачи на основе синтетических данных для различных сред была показана эффективность алгоритма и достоверность получаемых результатов как в простых средах (с точечным рассеивателем, с горизонтальным слоем, с одним наклонным слоем), так и в достаточно сложноустроенных средах (с несколькими наклонными слоями) При этом не накладываются ограничения ни на макроскорост-ную модель, ни на углы наклонов отражающих границ, так как используется конечно-разностное моделирование волновых полей Естественно, что это требует существенных временных затрат, но, как было показано, при использовании небольшого количества источников задачу удается решить на современных персональных компьютерах в реальные сроки Увеличение числа источников влечет необходимость продолжения разработки алгоритма, основным аспектом которой должно быть построение эффективного предобуславливателя для линеаризованного оператора прямой задачи Это позволит снизить количество итераций ЬБС^К, и, следовательно, существенно увеличит быстродействие всей процедуры обращения

СПИСОК ОСНОВНЫХ РАБОТ ПО ТЕМЕ ДИССЕРТАЦИИ Рецензируемые журналы по Перечню ВАК

1 Сильвестров И Ю Анализ сингулярного разложения линеаризованного оператора динамической теории упругости для случая вертикального сейсмического профилирования // Вычислительные технологии - 2007 -Т 12 -№6 - С 90-100

2 Сильвестров И Ю Прогнозирование строения среды ниже забоя скважины с помощью многокомпонентного обращения данных ВСП с выносным источником // Технологии сейсморазведки - 2007 -№3 - С 44-50

Труды конференций

3 Сильвестров И Ю Применение SVD-анализа для выбора оптимальной параметризации среды при решении обратной задачи вертикального сейсмического профилирования // Сб тр XI Всероссийской школы-семинара «Современные проблемы математического моделирования» - Ростов н/Д, Ростовский гос унив , 2005, с 452 - С 369-376

Расширенные тезисы международных конференций

4 Silvestrov I Full waveform inversion of multicomponent offset VSP data for recovery impedances bellow borehole bottom // Extended Abstracts 69th EAGE Annual Conference&Exhibition - 2007 -London UK - CD-ROM - P4 ISBN978-90-73781-54-2

5 Silvestrov I Full waveform inversion of VSP data for prediction of impedances below borehole bottom // Extended Abstracts EAGE Conference& Exhibition «Geosciences - To Discover and Develop» -2006 - SaintPetersburg - CD-ROM -P4 ISBN978-90-73781-64-1

Тезисы российских и международных конференций

6 Сильвестров И Ю Анализ сингулярного разложения линеаризованного оператора динамической теории упругости для задачи прогнозирования строения среды ниже забоя скважины по данным НВСП И Тез докл VII Ежегодной международной Конференции «Гальперинские чтения-2007» - М - 2007 - С 50-54

7 Сильвестров И Ю Итерационный метод решения обратной динамической задачи вертикального сейсмического профилирования // Тез докл Ш Всероссийской конференции «Актуальные проблемы прикладной математики и механика» -Абрау-Дюрсо -2006 -С 94-96

8 Сильвестров И Ю О выборе параметров среды при обращении данных вертикального сейсмического профилирования // Тез докл V Международной научно-практической геолого-геофизической конференции-конкурса молодых ученых и специалистов "Геофизика 2005" -Санкт-Петербург - С 262-264

_Технический редактор О М Вараксина_

Подписано к печати 25 04 2008 Формат 60x84/16 Бумага офсет №1 Гарнитура Тайме

_Печ л 0,9 Тираж 100 Заказ № 7_

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

ВВЕДЕНИЕ.

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

Глава 2. ТЕОРЕТИЧЕСКИЙ АНАЛИЗ СИНГУЛЯРНОГО

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

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

2.2. Построение линеаризованного оператора в однородной среде с использованием явного вида матрицы Грина.

2.3. Численный анализ сингулярного разложения.

Глава 3. РАЗРАБОТКА ЧИСЛЕННОГО АЛГОРИТМА РЕШЕНИЯ

ОБРАТНОЙ ЗАДАЧИ.

3.1. Описание алгоритма решения обратной задачи.

3.2. Численное построение линеаризованного оператора и сопряженного к нему.

3.3. Реализация алгоритма и разработанное программное обеспечение.

3.4. Численные эксперименты для различных сред (с точечным рассеивателем; с горизонтальным слоем; с наклонным слоем; с несколькими наклонными слоями).

Актуальность. Одним из подходов к численному решению обратных задач сейсморазведки в сложноустроенных средах является применение оптимизационных методой. При этом задача рассматривается как операторное уравнение, в правой части которого стоит зарегистрированное волновое поле, а неизвестными являются сейсмические параметры среды. Оператор задачи, отображающий эти параметры в данные наблюдений, определяется математической моделью, описывающей распространение волн в среде, и неявно задается, как правило, волновым уравнением либо уравнениями линейной теории упругости. Учитывая различные помехи, происходящие при записи данных, и неточности задания самой модели, решение при этом ищется в смысле наименьших квадратов. Для его поиска применяются, как правило, локальные итерационные методы, привлекающие производную возникающего нелинейного оператора. Это может быть либо метод Ньютона решения нелинейного операторного уравнения, либо градиентные методы минимизации функционала, представляющего собой среднеквадратичное уклонение зарегистрированного волнового поля от рассчитанного для '.'пробной" модели среды. В том и в другом случае свойства данной производной играют ключевую роль в сходимости итерационного процесса. В связи с этим их изучение является необходимым этапом при разработке численных методов решения обратной задачи и их применении. Исследованием данного вопроса для традиционной поверхностной системы наблюдений, когда источники и приемники упругих колебаний находятся па поверхности земли, занимались В.А. Чеверда, В.И. Костин, G. Chavent, F. Clement, D. Lebrun, D. Mace, A. Nicolao, A. Tarantola и др. Было показано, что попытки решения обратной задачи "в лоб"могут приводить к результатам с невысокой степенью достоверности. Так, например, в работе F. Assous и F. Collino (1990) проведено исследование явления связанности" упругих параметров. Оказывается, что при численном решении обратной задачи существуют такие тройки параметров, описывающие упругую изотропную среду, что истинная неоднородность в среде только по одному из них восстанавливается как возмущение по всем параметрам. Это означает, что решение задачи получается абсолютно неверным.

Для широко внедряемой в последнее время геометрии наблюдений вертикального сейсмического профилирования (ВСП) с выносным источником, при которой приемники колебаний располагаются в скважине, а источник находится на поверхности, на некотором удалении от нее, исследование свойств производной исходного нелинейного оператора для уравнений теории упругости проведено не было. В то же время с использованием именно такой системы наблюдений связывают в настоящее время возможность решения актуальной прикладной задачи определения упругих параметров среды ниже забоя скважины при бурении на нефть и газ в регионах со сложным геологическим строением. Считается, что с использованием этих параметров возможно прогнозирование наличия зон аномального давления, отсутствие информации о которых влечет аварийные ситуации вплоть до полного выхода скважины из строя. Впервые применение оптимизационного подхода для решения данной задачи в рамках двумерных уравнений динамической теории упругости было проделано совсем недавно в работах М.А. Roberts (2006) и М.А. Roberts и В.Е. Hornby (2007). Однако эти работы имеют исключительно практический характер в силу того, что подход был применен в первом случае для сложных синтетических данных, модель которых автор не приводит, а во втором случае - для реальных данных. Поэтому оценить достоверность полученных результатов можно только по косвенным признакам.

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

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

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

Основные этапы исследования.

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

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

3. С учетом этих свойств создать и программно реализовать алгоритм для численного решения рассматриваемой обратной динамической задачи теории упругости на основе итерационного метода LSQR решения систем линейных уравнений, с использованием для моделирования волнового процесса конечно-разностной схемы Вирьё на сдвинутых сетках и идеально согласованных слоев (PML).

4. Провести представительную серию численных экспериментов по решению обратной задачи на основе синтетических данных в средах различной степени сложности.

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

Основным методом исследования является математическое моделирование.

Изучение свойства решения обратной задачи, получаемого методом Ньютона, происходило в рамках теории условно-корректных задач. При этом существенно использовалось обобщение понятия г-решения, возникающее в теории вычислительной линейной алгебры, на случай компактных оператором в гильбертовых пространствах, разработанное в работах В.И. Костина и В.А. Чеверды (1995, 1998).

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

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

Результаты исследования свойств решений рассматриваемой обратной задачи, получаемые с использованием анализа сингулярного разложения линеаризованного оператора динамической теории упругости, сравнивались с результатами A. Tarantola, основанными на физических соображениях, а так же с результатами D. Lebrun и A. Nicolao полученными для других систем наблюдений.

На защиту выносятся

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

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

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

Новизна работы.

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

• опираясь на анализ поведения г-решения в зависимости от используемой параметризации среды обоснованно, что наименее связанными параметрами при численном решении обратной задачи с использованием метода Ньютона являются упругие импедансы;

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

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

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

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

• Исходя из полученной структуры численного решения обратной задачи разработан, программно реализован и протестирован оригинальный алгоритм определения местоположений и амплитуд разрывов (в линейном приближении) упругих импедансов среды на основе итерационного метода LSQR с использованием при моделировании волновых процессов конечноразпостной схемы Вирьё на сдвинутых сетках с ограничением расчетной области при помощи идеально согласованных слоев (PML).

Личный вклад.

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

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

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

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

Научные результаты диссертации докладывались на Международной конференции посвященной 75-летию академика М.М.Лаврентьева "Обратные и некорректные задачи математической физики" (Новосибирск, 2007); VII Международной конференции "Гальперинские чтения 2007" (Москва, 2007); 69-ой Международной конференции европейской ассоциации геофизиков и инженеров "69th EAGE Conference & Exhibition" (Лондон, Великобритания, 2007); III Всероссийской конференции "Актуальные проблемы прикладной математики и механика" (Абрау-Дюрсо, 2006); Международной конференции, организованной ассоциациями геофизиков ОЕА-ГО, EAGE, SEG "Geosciences - То Discover and Develop" (Санкт - Петербург, 2006); V Международной научно-практической геолого-геофизической конференции-конкурсе молодых ученых и специалистов Геофизика-2005 (Санкт - Петербург, 2005); XI Всероссийской школе-семинаре: "Современные проблемы математического моделирования" (Абрау-Дюрсо, 2005)

Публикации.

Результаты исследования по теме диссертации изложены в 8 опубликованных работах. Из них - 2 статьи в ведущих российских рецензируемых журналах по Перечню ВАК, 1 работа - в трудах российской конференции, 2 работы - в расширенных тезисах докладов международных конференций, 3 работы - в тезисах российских и международных конференций.

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

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

Автор выражает благодарность всем сотрудникам Лаборатории за полезные научные дискуссии в ходе выполнения работы. Успешному выполнению работы во многом способствовали ценные и критические замечания В.И. Костина. Хочется поблагодарить С.Б. Горшкалева и Г.В. Решетову за то, что они взялись за труд ознакомиться с работой и высказать о ней своё мнение.

Необходимо отметить неоценимую помощь, оказанную В.И. Самойловой при работе над рукописью.

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

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

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

Заключение

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

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

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

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

1. Алексеев А. С. Обратные динамические задачи сейсмики // Некоторые методы и алгоритмы интерпретации данных. - М. - 1967. - С. 9-84.

2. Алексеев А.С., Авдеев А.В., Фатьянов А.Г., Чеверда В.А. Волновые процессы в вертикально-неоднородных средах: прямые и обратные задачи. Новосибирск, 1994. - (Препринт / РАН. Сиб. отд-пие. ВЦ;1035).

3. Бакушинский А. Б., Гончарский А.В. Итерационные методы решения некорректных задач. М.: Наука. Гл. ред. физ.-мат. лит., 1989. - 128 с.

4. Бляс Э.А. Оптимизационное решение одномерной обратной динамической задачи (по данным ВСП) // Геология и геофизика. 1997. - Т. 38. - № 8. - С. 1398-1410.

5. Вишневский Д.М., Костин В.И., Чеверда В.А. Возбуждение сейсмических волн источником, расположенным в скважине, заполненной жидкостью // Физическая мезомеханика. 2003. - Т. 5 - № 5 - С. 85 - 92.

6. Владимиров B.C. Уравнения математической физики. М.: Наука, 1972. - 395 с.

7. Воронина Т.А., Чеверда В.А. Обращение полных волновых полей при обработке данных метода сейсмического профилирования // Доклады Академии Наук. -1994. Т. 35. - № 4. - С. 408-412.

8. Воронина Т.А., Чеверда В.А. Оптимизационный подход к обработке данных метода вертикального сейсмического профилирования // Геология и геофизика. 1994. -Т. 35. - № 5. - С. 127-139.

9. Гальперин Е.И. Вертикальное сейсмическое профилирование. М.: Недра, 1982. -344 с.

10. Гальперин Е.И. Вертикальное сейсмическое профилирование: опыт и результаты. М.: Наука, 1994. - 320 с.

11. Голуб Дж.,Ван Лоун Ч. Матричные вычисления. М.: Мир, 1999. - 548 с.

12. Годунов С.К., Антонов А.Г., Кирилюк О.П. , Костин В.И. Гарантированная точность решения систем линейных уравнений в евклидовых пространствах. Новосибирск, Наука, 1990. - 352 с.

13. Доледенок В.Г., Костип В.И., Чеверда В.А. r-решение уравнения первого рода с компактным оператором: существование и устойчивость. Новосибирск, 1994. -(Препринт / РАН. Сиб. отд-ние. ВЦ;1035).

14. Колмогоров А.Ф., Фомин С.В. Элементы теории функции и функционального анализа. М.: Наука, 1989. - 624 с.

15. Трехмерные задачи математической теории упругости и термоупругости / Под ред. В.Д. Купрадзе. М: Наука, 1976.

16. Романов В.Г. Обратные задачи математической физики. М.: Наука, 1984

17. Самарский А.А. Теория разностных схем. М.: Наука, 1983. - 616 с.

18. Сильвестров И.Ю. Прогнозирование строения среды ниже забоя скважины с помощью многокомпонентного обращения данных ВСП с выносным источником // Технологии сейсморазведки. 2007. - № 3. - С.44-50.

19. Сильвестров И. Ю. Анализ сингулярного разложения линеаризованного оператора динамической теории упругости для случая вертикального сейсмического профилирования // Вычислительные технологии. 2007. - Т. 12. - № 6. - С. 90-100.

20. Сильвестров И.Ю. О решении обратной задачи динамической теории вертикального сейсмического профилирования // Материалы XILV МНСК "Студент и научно-технический прогресс", Новосибирск, 2006. С. 160-161.

21. Сильвестров И.Ю. Итерационный метод решения обратной динамической задачи вертикального сейсмического профилирования // Тез.докл. III Всероссийской конференции "Актуальные проблемы прикладной математики и механика". Абрау-Дюрсо. - 2006. - С.94-96.

22. Тихонов А.Н., Самарский А.А. Уравнения математической физикию М.: Наука, 1972. - 736 с.

23. F. Assous, F.Collino A numerical method for the exploration of sensitivity: the case of the identification of the 2D stratified elastic medium // Inverse Problems. 1990. -Vol.6. - N. 4. - P. 487-513.

24. Bamberger A., Chavent G., Lailly P. About, the stability of the inverse problem in 1-D wave equations. Application to the interpretation of seismic profiles // Appl. Math. Optim. 1979. - Vol. 5. - P. 1-47.

25. Boschetti F., M. C. Dentith, R. D. List Inversion of seismic refraction data using genetic algorithms // Geophysics. 1996. - Vol. 61. - R 1715-1727.

26. J.P. Berenger. A perfectly matched layer for the absorption of electromagnetic waves // Journal of Сотр. Physics. 1994. - Vol. 114. - P. 185-200.

27. W. Chang and McMechan G. Reverse-time migration of offset vertical seismic profiling data using the exitation-time imaging condition. Geophysics. - 1986. - Vol. 51. - P. 67-84.

28. Cheverda V.A., Clement F., Khaidukov V.G., Kostin V.I. Linearized inversion of multi-offset data for vertically-inhomogeneous background // Journal of Inverse and Ill-Posed problems. 1998. - Vol. 6. - P. 455 - 477.

29. Claerbout, J. F. Toward a unified theory of reflector mapping // Geophysics. 1971. -Vol. 36 - P. 467-481.

30. Collino F., Tsogka C. Application of PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media // Geophysics. 2001. -Vol. 66. - P. 294 - 307.

31. Collino F. Perfectly matched layers for the paraxial equations //J. Comput. Phys. -1997. Vol. 131(1). - P. 164-180.

32. Esmersoy C. Inversion of P and SW waves from multicomponent offset vertical seismic profiles // Geophysics. 1990. - Vol. 55. - P. 39-50.

33. Gauthier O., Virieux J., Tarantola A. Two-dimensional nonlinear inversion of seismic waveforms. Numerical results // Geophysics. 1986. - Vol. 51. - P. 1387-1403.

34. Harlan W. Separation of signal and noise applied to vertical seismic profiles // Geophysics. 1988. - Vol. 53. - № 7. - P. 932 946.

35. Lebrun D., V. Richard, D. Mace, M. Cuer. SVD for multioffset linearized inversion: Resolution analysis in multicomponent acquisition // Geophysics. 2001. - Vol. 66. -P. 871-882

36. Kostin V.I., Tcheverda V.A., r-Pseudoinverse for compact operators in Hilbert space: existence and stability j j J.Inverse and Ill-Posed Problems. 1995. - Vol. 3. - P. 131148.

37. W. S. Leaney, В. E. Hornby, A. J. Campbell, S. Viceer, M. Albertin, A. Malinverno Sub-salt velocity prediction with a look-ahead AVO walkaway // Extended abstracts "74th Meeting SEG". 2004.

38. Mace D, P. Lailly Solution of the VSP one-dimensional inverse problem // Geophysical Prospecting. 1986. - Vol. 34. - P. 1002-1021.

39. Mittet R., Ketil Hokstad, Jan Helgesen, Guy Canadas Imaging of offset VSP data with an elastic iterative migration scheme // Geophysical Prospecting. 1997. - Vol. 45. -№ 2. - P. 247-267.

40. Mora P. Nonlinear two-dimensional elastic inversion of multioffset seismic data // Geophysics. 1987. - Vol. 52. - P. 1211-1228.

41. Mora P. Elastic wave-field inversion of reflection and transmission data // Geophysics.- 1988. Vol. 53. - P. 750-759.

42. Nemeth Т., Chengjun Wu, Schuster G.T. Least-squares migration of incomplete reflection data // Geophysics. 1999. - Vol. 64. - P. 208-221.

43. Nicolao A., Drufuca G., Rocca A. Eigenvalues and eigenvectors of linearized elastic inversion // Geophysics. 1993. - Vol. 58. - P. 670-679.

44. Paige С. С., M. A. Saunders LSQR: An algorithm for sparse linear equations and sparse least squares // TOMS. 1982. - Vol. 8. - № 1. -P. 43-71.

45. Roberts M. A. Elastic Parameter estimation from the 2d waveform inversion of a look-ahead walk-away VSP survey // Extended abstracts "68th Conference and Exhibition".- Viena. 2006.

46. Roberts M.A., B.E. Hornby Application of 2d full waveform inversion to walkaway VSP data for the estimation of sub-salt elastic parameters // Extended abstracts "69th Conference and Exhibition". London. - 2007.

47. Sambridge M. Monte Carlo methods in geophysical inverse problems // Reviews of Geophysics. 2002. - Vol. 40. - N. 3. - P. 1009

48. Silvestrov I. Full waveform inversion of multicomponent offset VSP data for recovery impedances bellow borehole bottom // Extended Abstracts 69th EAGE Annual Conference and Exhibition. 2007. - London. UK. - CD-ROM. - P.4. ISBN 978-9073781-54-2.

49. Snieder R., M. V. Xie, A. Pica and A. Tarantola. Retrieving both the impedance contrast and background velocity: A global strategy for the seismic reflection problem // Geophysics. -1989. Vol.54. - P. 991-1000.

50. Tal-Virsky B.B., Tabakov A. A. High-resolution prediction of acoustic impedances below bottom-of-hole // Geophysical Prospecting. 1983. - Vol. 31. - P .225-236.

51. A. Tarantola, B. Valette Generalized non linear inverse problems solved using the least-squares criterion // Rev. of Geophys.and Space Physics. 1982. - № 20. - P. 219-232.

52. A. Tarantola Inversion of seismic reflection data in the acoustic approximation // Geophysics. 1984. - Vol. 49. - P. 1259-1266.

53. A. Tarantola A strategy for nonlinear elastic inversion of seismic reflection data // Geophysics. 1986. - Vol. 51. - № 10. - P .1893-1903.59| Albert Tarantola Inverse problem theory and methods for model parameter estimation.- SIAM, 2005.

54. Virieux J. P-SV wave propagation in heterogeneous media: velocity-stress Finite difference method // Geophysics. 1986. - Vol. 51. - № 4. - P. 889-901.