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

Карчевский, Андрей Леонидович

Бесплатный автореферат и диссертация по наукам о земле на тему
Оптимизационный подход к решению обратной динамической задачи сейсмики для горизонтально-слоистых анизотропных сред на основе явного аналитического решения прямой задачи в частотной области
ВАК РФ 25.00.10, Геофизика, геофизические методы поисков полезных ископаемых

Автореферат диссертации по теме "Оптимизационный подход к решению обратной динамической задачи сейсмики для горизонтально-слоистых анизотропных сред на основе явного аналитического решения прямой задачи в частотной области"

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

КАРЧЕВСКИЙ Андрей Леонидович

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

25.00.10 — Геофизика, геофизические методы

поисков полезных ископаемых, 05.13.18 — Математическое моделирование, численные методы и комплексы программ

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

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

Работа выполнена в Лаборатории волновых процессов Института математики им. С.Л.Соболева СО РАН

Официальные оппоненты:

доктор физико-математических наук, член-корреспондент РАН Михайленко Борис Григорьевич доктор физико-математических наук, профессор Волчков Юрий Матвеевич

доктор физико-математических наук, профессор Каштан Борис Маркович

Ведущая организация: Институт геофизики УрО РАН

(ИГФ УрО РАН, г.Екатеринбург)

Защита состоится 9 декабря 2005 года в 10 часов на заседании диссертационного совета Д 003.050.05 при Объединенном институте геологии, геофизики и минералогии им. А.А.Трофимука СО РАН, в конференц-зале. """"

Адрес: 630090, Новосибирск, пр-т Ак. Коптюга, 3.

Факс: (383) 333 27 92

С диссертацией можно ознакомиться в библиотеке Объединенного института геологии, геофизики и минералогии им. А.А.Трофимука СО РАН.

Автореферат разослан и • I о_2005 года.

Ученый секретарь диссертационного совета

д.ф.-м.н. Ю.А.Дашевский

Общая характеристика работы

Актуальность. Как правило, использование математического моделирования в практической сейсморазведке ограничено. При этом в широко распространенных обрабатывающих и интерпретационных комплексах программ оно ограничивается акустическими моделями. Это связано с тем, что до настоящего времени основной целью стандартной обработки является получение временных разрезов, которые в некотором приближении могут соответствовать реакции акустической модели среды на нормальное падение плоской продольной волны. На данном представлении основаны многие интерпретационные подходы, например, псевдоакустический каротаж. Однако, при реализации этого подхода происходит значительная потеря информации о среде, содержащаяся в исходных сейсмограммах. Желание не терять данную информацию, а использовать ее при интерпретации, приводит к необходимости использования способов обработки сейсмических данных, в которых параметры среды определяются по полному волновому полю. В этом случае акустическое приближение становится не адекватным и должно быть заменено моделью упругой среды с возможным учетом анизотропии и поглощения в зависимости о ситуации. Естественно, возникает целесообразность использования не только вертикальной, но и горизонтальных компонент вектора смещений, в частности вовлечения в обработку волн РБ (хотя возможность определения упругих характеристик среды по лишь одной из компонент остается интересной и в этом случае).

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

1 библиотека I

1! ¿"39а8// \

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

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

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

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

В пользу такого подхода говорят следующие обстоятельства. Во-первых, существует метод решения прямой задачи для СДУ теории упругости в частотной области для частного вида горизонтально-слоистых однородных сред (Аккуратов Г.В. и Дмитриев В.И. (1979, 1984), Фатьянов А.Г. и Ми-хайленко Б.Г. (1988), Фатьянов А.Г. (1989)). Численный алгоритм на основе данного метода дает решение за малое время. Во-вторых, для построения функционала невязки могут быть использованы функции, значения которых известны только для ограниченного набора параметров интегральных преобразований, что приведет к существенному сокращению времени вычисления функционала невязки и числа решений прямой задачи.

При решении обратной задачи в частотной области встают две основные проблемы.

Первая проблема — переход к образам интегральных преобразований от реальных сейсмограмм.

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

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

Вторая проблема — создание метода решения прямой задачи для СДУ теории упругости в частотной области для горизонтально-слоистой однородной среды, способного дать решение за короткое время, создание устойчивого вычислительного алгоритма для получения решения данной задачи.

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

Для функционала невязки, во-первых, необходимо проведение исследований

• на наличие локальных минимумов и максимумов,

• на выявление характера овражности,

• на выяснение возможности достижения точки глобального минимума,

• на установление характера зависимости функционала невязки от его параметров.

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

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

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

• Какая часть спектральных данных пригодна для численного решения обратной задачи?

• Какова разрешающая способность обратной динамической задачи сей-смики в частотной области? (То есть необходимо установить степень влияния отдельных упругих параметров на поведение решения прямой задачи).

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

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

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

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

• создания устойчивого вычислительного алгоритма для решения данной задачи.

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

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

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

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

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

Р{1)Ч5{х1,х2,хг-х\). (1)

Источник вида (1) есть центр давления, который является моделью взрыва.

Для перехода в частотную область к СДУ теории упругости применено преобразование Фурье по пространственным переменным и х^ и преобразование Лапласа по временной переменной t. В случае изотропной среды СДУ теории упругости может быть записана в цилиндрической системе координат, тогда для перехода в частотную область по пространственной переменной г применено преобразование Фурье-Бесселя, а по временной переменной t — преобразование Лапласа.

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

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

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

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

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

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

3. Разработан метод построения явных аналитических формул для градиента функционала невязки на основе предложенного метода решения СДУ теории упругости в частотной области.

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

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

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

1. Прямая задача:

• построено явное аналитическое решение СДУ теории упругости в частотной области для горизонтально-слоистых однородных сред любого вида анизотропии в любой точке среды, не имеющее ограничений на толщину слоев (модель среды может содержать как очень толстые, так и очень тонкие слои), в основе которого лежит редукция СДУ второго порядка к ДМУР, для которого также построено решение в явной аналитической форме в каждом слое среды;

• предложен устойчивый алгоритм численного нахождения решения СДУ теории упругости, в основе которого лежит знание корней характеристического уравнения данной системы и их свойств;

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

2. Функционал невязки:

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

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

3. Градиент функционала невязки:

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

- введение ступенчатой функции специального вида,

- получение постановок прямых задач для приращений решения прямой задачи для СДУ теории упругости при вариации различных упругих постоянных и плотности,

- объединение в единую постановку прямой задачи для решения СДУ теории упругости и его приращений,

- применение к полученной СДУ схемы исследований, предложенной для решения прямой задачи для СДУ теории упругости.

4. Обратная задача: На основе знания аналитической структуры решения прямой задачи и численного исследования зависимости обратной задачи и функционала певязки от своих параметров

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

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

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

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

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

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

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

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

Результаты работы докладывались на следующих научных конференциях:

• II Международная конференция по математическому моделированию, Якутск, 28 июня - 2 июля, 1997.

• Всероссийская научная конференции "Алгоритмический анализ некорректных задач", Екатеринбург, 2-6 февраля, 1998.

• III Сибирский конгресс по прикладной и индустриальной математике, Новосибирск, 22-25 июня, 1998.

• Международная конференция "Обратные задачи математической физики", Новосибирск, 21-25 сентября, 1998.

• VI конференция "Обратные и некорректно поставленные задачи", Москва, 20-21 июня, 2000.

• Международная конференция "Некорректные и обратные задачи физики", Новосибирск, 5-9 августа, 2002.

• Международная конференция по Математическим Методам в Геофизике (ММГ-2003), Новосибирск, 8-12 октября, 2003.

• International Conference "Inverse Problems: Modeling and Simulation", Fethiye, Turkey, June 07-12, 2004.

• Международная конференция по вычислительной математике (МКВМ-2004), Новосибирск, 21-25 июня, 2004.

Результаты работы вошли в отчеты по следующим грантам:

• Гранты РФФИ: 05-01-00171, 05-01-00559, 02-01-00818, 02-01-00809, 99-01-00563, 96-01-01937, 93-01-01739.

• Грантам научных школ: НШ-1172.2003.1, 00-15-96183, 96-15-96284.

• Грантам Министерства образования РФ: УР.04.01.026, N-00-1.0-73.

Результаты работы докладывались на научных семинарах:

• Семинар лаборатории волновых процессов Института математики им. С.Л.Соболева СО РАН (руководитель семинара член-корр. РАН В.Г.Романов).

• Семинар лаборатории обратных задач математической физики Института математики им. С.Л.Соболева СО РАН (руководитель семинара д.ф.-м.н. Ю.Е.Аниконов).

• Семинар лаборатории численных методов решения обратных задач Института математики им. С.Л.Соболева СО РАН (руководитель семинара д.ф.-м.н. А.Л.Бухгейм).

• Семинар отдела механики деформируемого твердого тела Института гидродинамики им. М.А.Лаврентьева СО РАН (руководитель семинара член-корр. РАН Б.Д.Аннин)

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

Структура и объем диссертации. Диссертационная работа состоит из Введения, трех глав, Заключения, библиографического списка, включающего 391 наименование. Работа изложена на 198 страницах, содержит 37 рисунков и 12 таблиц.

Содержание работы

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

В Главе 1 приведен метод решения прямой динамической задачи сейсмики в частотной области, в основе которого лежит СДУ уравнений упругости.

Пусть среда - Л^-слойная структура с границами раздела х§, к = О, N1,

х\ = 0; ш-ый слой - интервал [я™-1, я™], N¡+1 (подстилающий) слой - [х£,оо).

Физические свойства каждого слоя характеризуются величинами модулей упругости Ст]Ы и плотности р, то есть Ст]ы и р- кусочно-постоянные функции ОТ Хз, 0 < Хз < 00.

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

*3'

д2и„

],к,1=1 ^

дщ

г М2

т = 1,3, Начальные условия:

дхк

Рис. 1: Слоистая среда.

х = (хь х2, х3) б!2х1+, г € К. Ит|г<о = 0, т = 1,3.

Краевые условия

к,1=1

дщ

дхк

= 0, 7 = 1,3.

х3=0

Условия склейки: з

к,1=1 ОХк

= О, [и,ь = о,

1,3.

х3

(2)

(3)

(4)

(5)

Нашей целью является определение величин

ит{их,и2,х3,р), тп = 1,3,

где функции •¡¿„(уъ и2,х3,р) есть образ функций ит(ж1, ж3, £), р——а+ш - параметр преобразования Лапласа по временной переменной V] и 1/2 - параметры преобразования Фурье по переменным и х2 соответственно.

Величины ит{и1,и2,х3,р) удовлетворяют задаче:

д_

дх3

А^-и + Ши) + гВ-~и - Ш = Р[р, х3, х*3),

= 0,

А-г—и + 1Ви ох з

х3=0

= о, [и]хз=хь = 0,

к = 1,п,

и —» 0 (хз —► оо). Здесь введены следующие обозначения:

С15 С56 С55 С56 С25 С45

Си С46 С45 + /9!/2 С46 С24 С44

.С13 С36 С35. .Сзб С23 С34.

"й5 Сц С1з" С56 С46 С36

С56 С46 С36 + С25 С24 С23

.С55 С45 С35. .С45 С44 С34.

В = рих

В = ррх

Б = рр2Е + ри\

2С16 С12 + С66 С14 + С56 +рь>11>2 С12 + С66 2С26 С25 + С46 С14 + С56 С25 + С46 2С45

Р(р,х3,х*3) = -Р{р)(Иг6{х3 - + Ы\хъ - Х%)),

здесь счр = СЦр!р — приведенные модули упругости.

"сц С16 С15 Сбб С26 С46

С16 Сбб С56 + С26 С22 С24

.С15 С56 С55. .С46 С24 С44.

(6) (7)

(8) (9)

"0" С55 С45 С35

и = и2 , к = V2 , к = 0 А — р С45 С44 С34

0 1 ,С35 С34 сзз.

Случай изотропной среды

Смещения продольных и поперечных волн могут быть определены из

д2щ д ( 0иЛ Р-,~ = — I м— I +

ЭР дх3 \ дх3)

д_ дг

- Г(р)6"(гЩх3-®5),

, <Эи2 д .

92Ц2 _ 8 7 дг2 ~ дх3

1 а

+"а~ 1 г гог

дщ д . .

-^(рМгМяз-хЗ),

Целью является определение величин

Щп(у,х3,р), т = Т72,

где функции ит(и,хз,р) есть образ функций ит(г, хз, £), р = —сН-го/ - параметр преобразования Лапласа по временной переменной - параметр преобразования Фурье-Бесселя по переменной г. Величины ип(1/,хз,р) удовлетворяют задаче:

А—V + иВи 9х3

Здесь введены обозначения:

= 0, и 0 (х3-»оо),

1з=0

= 0, [С/Ь = 0, * =

(11) (12)

«1

Р{р,и,хз,Хз) = ^(р) (I/

¿(х3 -х5) -

¿'(х3-х^) ,

Гг>2 0 ' 5 = р 0

5 0 < у2р- 2У] 0

£> = рр Е + ри2

у2р о о

Ранее известные результаты по разработке метода решения прямой динамической задачи сейсмики в частотной области для горизонтально-слоистых сред принадлежат Аккуратову Г.В. и Дмитриеву В.И. (1979, 1984) для изотропной среды, Фатьянову А.Г. и Михайленко Б.Г. (1988) для изотропной среды с поглощением, Фатьянову А.Г. (1989) для трансверсально-изотропной среды (с поглощением) с осю симметрии бесконечного порядка, сопадающей с Охз.

В нашем случае схема построения решения отличается от схем в вышеперечисленных работах.

Для прямой задачи (6)-(9) введем квадратную матрицу X порядка 3 следующим соотношением:

А~~и + гВИ = XII. (13)

дх3

ДМУР, которому удовлетворяет матрица X:

4-Х + (X + (В)А~1 (X - Ш) = £>. (14)

ах з

Схема решения прямой задачи:

• имеется постановка прямой задачи для СДУ теории упругости в горизонтально-слоистой среде;

• для СДУ теории упругости получаем ДМУР (14);

в каждом слое строим решение ДМУР в аналитическом виде; условия склейки [Х]хь = 0 в точках разрыва среды позволяют вести рекуррентный пересчет со слоя на слой;

• вычислив по формулам, представленным в аналитическом виде, необходимые величины, получаем [/ в любой точке £3.

Решение ДМУР (14) на интервале [я™-1, х™} получено в следующем виде: X = Х?+) + 1Г\ (15)

Ь = -

ХЦ, = АС 4- Ш, Х!1) = дА-1В, С - АСА^ + гВА"1.

Для того чтобы можно было использовать решение ДМУР в виде (15), были решены проблемы:

• построения матриц С, С и

• вычисления выражений типа второго слагаемого в (16);

• вычисления корней характеристического уравнения СДУ:

¿е^АХ2 + ¿В А - £>] = 0.

Таким образом, в Главе 1 с помощью явного аналитического решения ДМУР (15) построено явное аналитическое решение прямой задачи (б)-(9) и(щ,и2,хз,р) для горизонтально-слоистой однородной среды любого вида анизотропии, решение может быть вычислено в любой точке хз, не имеет ограничений на толщину слоев: модель среды может содержать как очень толстые, так и очень тонкие слои. Каждый этап алгоритма численного решения прямой задачи является устойчивым и содержит конечное количество шагов вычислений. Единственый этап — вычисление корней характеристического уравнения СДУ (6) — производится при помощи итерационного метода, поскольку характеристическое уравнение является алгебраическим уравнением шестого порядка. Для его решения предложен численный алгоритм, в основе которого лежит использование известной априорной информации о корнях уравнения и метод выделения квадратного множителя Хичкока. Результат тестирования алгоритма: при всех численных экспериментах для всех вычисленных корней уравнения ошибка вычислений присутствовала не ближе 8-12 знака после запятой. Для метода Хичкока требовалось 4-6 итераций.

Здесь в качестве примера приведем только один пример, относящийся к явлению квазианизотропии и согласующийся с результатами монографии Сибирякова Б.П., Максимова Л.А., Татарникова М.А. (1980).

Для численного эксперимента была выбрана среда, представленная на Рисунке 2. На глубине 1400 м был расположен слой мощностью 100 м. Средние значения продольной и поперечной скоростей были выбраны следующими: ьр = 2,8 км/с, = 1,5 км/с.

Доминантная частота источника — 30 Гц, что соответствует при скорости ьр = 2,8 км/с длине волны 93 м.

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

отвечающие различным значениям Avp и Ava. Для Avp последовательно были выбраны значения {0.4, 0.8, 1.2, 1.6}, в соответствии с ними Av3 принимало значения {0.2, 0.4, 0.6, 0.8}. Для определености на Рисунке 3 приведены только значения Avp.

Как видно на Рисунке 3 при разбиении 100-метрового слоя на 2-а слоя мощностью 50 метров хорошо видны отражения от трех границ: кровли 100-метрового слоя, его середины и подошвы (Рисунок 3 (а)). При разбиении 100-метрового слоя на 10-ть слоев мощностью 10 метров видно отражение от кровли слоя, но далее идет сложная картина накладывающихся

. v I iv|/>iv^ri пwnrnn iiv^vwiHriwunrivifi, ni/w/ibMi/iiwu^vri

друг на друга отражений (Ри- СВОйства однородной анизотропной среды, сунок 3 (Ь)). При разбиении

100-метрового слоя на 50 слоев мощностью 2 метра среда проявляет свойства квазианизотропии (Рисунок 3 (с)). В данном случае 100 метровый слой, составленный из 50-ти изотропных слоев, начинает себя вести как однородный трансверсально-изотропный слой с осью симметрии перпендикулярной напластованию: видны отражения от кровли и подошвы слоя.

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

Задержку по времени вк (к = 1,3) для толстого слоя, разбитого на 50 тонких изотропных слоев, нетрудно получить из кинематики, вычисляя время прохождения волны по вертикали от источника до подошвы толстого слоя

Vp-Avt>

Velocity, km/s

Vs.

Vp+Avp Vs-Avs Vs+Avs

Рис. 2: Модель горизонтально-слоистой изотропной спелы с прпиоличргким пеоеелаиванием. пооявляюшей

/ а \

4 > 2-а слоя

10-ть слоев

Рис. 3: Вертикальные компоненты смещений для случаев, когда слой мощностью 100 м разбит на (а) 2-а, (Ь) 10-ть и (с) 50-т слоев равной мощности. Сейсмотрассы приведены для различных значений Лур и Ли, (для удобства указаны только значения Лур). Удаление от источника 100 м.

Рис. 4: Вертикальные компоненты смещений для случаев, когда слой мощностью 100 м разбит на (а) 50-т и (Ь) 2-а слоя равной мощности. На рисунке представлены введенные обозначения для временных промежутков.

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

Результаты измерений, вычислений и сравнений временных интервалов сведены в Таблице 1. Значения величин без крышек измерены, с крышками — вычислены.

Таблица 1: Измеренные и вычисленные временные интервалы, (с)

# Тк вк вк Тк-вк ]вк ~ вк\ тк тк \Тк-Тк\

1 0.0085 0.0049 0.00486 0.0036 0.00004 0.0731 0.07365 0.0006

2 0.0244 0.0144 0.01458 0.0100 0.00014 0.0821 0.08089 0.0012

3 0.0517 0.0332 0.03314 0.0180 0.00006 0.0974 0.09520 0.0022

4 0.1243 0.12215 0.0021

Таким образом, очевиден факт, что в случае разбиеня 100-метрового слоя на два слоя смещение отражения от подошвы слоя с ростом Avp и Ava подчиняется кинематическим соотношениям, а в случае разбиения 100-метрового слоя на 50-т слоев с ростом Avp и Avs сказывается эффект квазианизотропии. Данный результат численного эксперимента вполне согласуется с данными лабораторных экспериментов, представленных в монографии Сибирякова Б.П., Максимова JI.A., Татарникова М.А. (1980).

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

Результаты исследований поведения функционала невязки можно найти в работе Чеверды В.А. (1990) (функционал невязки был вида

«2

I \n(u)\2du,

где п(и>) - невязка, и> - круговая частота), в работах Кабанихина С.И. (1995) и Карчевского А.Л. (1997, 1999) (исследовался функционал невязки вида

JN2(t)dt,

где N(t) - невязка, t - переменная времени).

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

В диссертационной работе функционал невязки строится следующим образом.

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

и(р,и,р) = и0(и,р), и0 =

«2,0

(17)

где Щ{у,р) — известная вектор-функция, а € Оа = {а | с*ь..., адга}, ш € = {а; | и е = {и \ и,..., 1/лгЛ и ЛГа, Л^ и —

конечные числа.

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

Обратная задача: Найти продольную ур и поперечную у3 скорости в каждом слое, где они неизвестны, если относительно прямой задачи известна дополнительная информация (17).

Функционал невязки может быть записан в следующей форме:

ЩУА = X! (18)

Такой вид функционала невязки отвечает практической постановке обратной задачи.

Обратная задача и, следовательно, функционал невязки (18) зависят от следующих параметров:

• значений пространственной и временной частот V и / (и; = 27г/);

• значений параметра затухания а, (р = —а + ш),

• числа пространственных и временнйх частот ЛГ„ и ЛГШ;

• интервалов, из которых они взяты.

На численных примерах показано, что функционал невязки (18) может иметь локальные минимумы и максимумы.

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

Было проведено численное исследование поведения функционала невязки в зависимости от выше перечисленных параметров. Выводы из численных экспериментов:

1. Рост параметра а оказывает сглаживающий эффект на {/.

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

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

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

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

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

Вышеперечисленные выводы 1-6 позволяют выбрать следующую общую стратегию минимизации функционала невязки:

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

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

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

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

• таким образом, восстановление каждого слоя будет проходить постепенно от верхих слоёв до нижних.

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

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

Вторая часть Главы 2 посвящена созданию метода получения аналитических формул для градиента функционала невязки вида (18) в случае горизонтально-слоистой среды.

Общий способ (использование сопряженной задачи) получения формул функционала невязки для СДУ теории упругости не дает аналитических формул.

Вводится ступенчая функция

и прямые

задачи Ц = 2, N1, I = 1,3, к = 1, N1):

дх3

(А^+.ву/)

дх3

1з=0

= - (л^—и+^и) , V? -> О (®з - 00),

13=0

где введены обозначения

А\ = ре3

М = ре3

О О

0 1

1 О О О

В{ = ре, В{ = -ре-з

0 О

1 О

О 1 2 О

Ь[ = ри2е.

Ц =

1 О О О

О О О 1

Аз — ез

Ч2 О О

, В]3 — е.

О -V:

УР~2 О

,21

, дз = е1

уур о О уУ5

г2 = Р_ +

V 9 ' '

р «р

^ = +

Иегр > О, 11егя > О.

Вектор-функция V* является линейной частью приращения для решения прямой задачи при вариации определенного параметра среды, индекс ] отвечает, что варьируемый параметр принадлежит слою с номером у, индекс I указывает, какой параметр варьируется: 1 — ьр, 2 — V2, 3 — р.

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

УФ[т£ VI р\ = ..., VI..., ..., , (19)

= Ие ^(ЛГь,, V/), Ыкп = и(0,ик,рп) - и0{ик,Рп),

к,п

1 = 1,3, ¿ = 2,м. Вводятся вектор-функции

(20)

которые удовлетворяют следующим прямым задачам: 1 дхз ■ / ■ с/х3

-С

+ - - Щ^х = Пр)*», (21)

= 0, И? -» 0 (ягз — оо), (22)

1з=0

= 0, Юхь = 0, к = Т^. (23) 23

К постановкам прямых задач (21)-(23) полностью применима схема исследования, предложенная для решения прямой динамической задачи сейсмики в Главе 1.

Таким образом, можем получить аналитические выражения для 1У/ и, следовательно, для V* и для градиента функционала невязки (18).

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

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

Пусть г будет первым слоем, а г+й — последним слоем этой пачки (см. Рис. 5).

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

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

• А: изотропной;

• В: орторомбической;

• С: трансверсально-изотропной с осью симметрии бесконечного порядка, лежащей в плоскости 0х\х2-

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

А: В: С:

ьр Сзз си

С55, С44 С\2, Сц

С13, с2з С\3

Си, С22, С\2 + Сбб сзз

С12, Сбб

х° х3

л3

х3 х"'

Рис 5: Модель среды. Серым отмечена пачка тонких слоев.

плотность р и точки разрыва среды

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

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

Показано, что улучшить точность восстановления элементов матриц В, В и D можно только за счет увеличения значений пространственных частот

V\ И V2-

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

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

Для орторомбической среды обратная задача по определению девяти приведенных модулей упругости сц, С12, С13, С22, С23, С33, С44, С55 и Сбб может быть сведена к решению серии обратных задач:

• Определяем С33.

• Определяем С55, С13 и сц, используя известный С33.

• Определяем С44, С23 и С22, используя известный С33.

• Определяем Сбб и С12, используя известные С33, С55, с44, С13, с2з, Сц, с22-

Для трансверсально-изотропной среды по определению пяти приведенных модулей упругости сц, С12, С13, С33, С44 обратная задача может быть сведена к решению серии обратных задач:

• Определяем сц.

• Определяем С12, используя известный сц.

• Определяем C13, C33 и С44, используя известные сц и Ci2.

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

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

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

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

• Неточное задание границ слоев в пачке. Численные эксперименты проводились для пачек с мощностью тонких слоев 7-12 метров с контрастностью 10-20%. Показано, что вариации значений в пределах ±2 метра оказывают слабое влияние на функционал невязки.

• Неточное задание плотности р. Для тех же пачек неточное задание плотности в пределах 10% может сказаться двояко па восстановлении скоростных параметров среды. Если значения плотности завышены или занижены, то скоростные параметры слоев могут быть определены с удовлетворительной точностью. Если все значения плотности заданы вразнобой с точным значением (и завышены, и занижены), то в этом случае возможно неудовлетворительное восстановление скоростей г>р и г>„.

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

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

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

Заключение

В работе предлагается решать обратную динамическую задачу сейсми-ки в частотной области методом минимизации функционала невязки. Известно, что скорость решения обратной задачи методом минимизации функционала невязки напрямую зависит от скорости решения прямой задачи, поэтому Глава 1 данной работы посвящена разработке метода решения прямой задачи. Было уделено максимальное внимание вычислительным вопросам, связанным со скоростью и устойчивостью определения необходимых величин. Далее, немаловажным условием успешного решения задачи минимизации является поведение функционала невязки. Из представленных в работе примеров было видно, что функционал невязки может иметь множество локальных минимумов и максимумов. Таким образом, разработка стратегии минимизации, позволяющей достичь глобальный минимум функционала, являлась важным этапом исследований. Важной проблемой минимизации является задача вычисления градиента функционала невязки, что позволяет применять градиентные методы для поиска минимума функционала. Этим двум вопросам была посвящена Глава 2 настоящей работы. И последнее, после того как имеются в наличии метод решения прямой задачи, изучены свойства функционала невязки, необходимо проведения исследования основных математических свойств обратной задачи. Глава 3 содержит теоретическое и численное исследование данных свойств обратной задачи. На базе знаний этих свойств могут быть построены алгоритмы решения обратных задач для конкретных постановок.

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

• Прямая задача

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

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

- Не существует ограничений на мощности слоёв: модель среды может содержать как очень толстые, так и очень тонкие слои.

- Идея построения метода решения прямой задачи отличается от известных аналогов.

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

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

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

Функционал невязки

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

- Для минимизации функционала невязки предлагается стратегия минимизации, которая приводит к глобальному минимуму, минуя "ловушки" локальных минимумов.

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

Градиент функционала невязки.

- Метод построения градиента функционала невязки позволяет получить явные аналитические представления для компонент градиента. В работе метод продемонстрирован на примере обратной задачи для СДУ теории упругости для горизонтально-слоистой однородной изотропной среды.

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

Обратная задача

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

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

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

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

• в области интерпретации на реальных данных:

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

• в области математического моделирования:

— создание эффективного алгоритма суммирования по пространственным и временной частотам с целыо получения синтетических сейсмограмм;

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

Выше приведенные пункты определяют направления дальнейших исследований.

Список публикаций по теме диссертации

1. Карчевский A.J1. Поведение функционала невязки для нелинейной коэффициентной гиперболической обратной задачи / Труды международного семинара "Обратные задачи геофизики". Новосибирск: ВЦ СО РАН, 1996, с. 106-109.

2. Karchevsky A.L. Properties of the misfit functional for a nonlinear one-dimensional coefficient hyperbolic inverse problem / Journal of Inverse and Ill-posed Problems, 1997, v. 5, n. 2, p. 139-163.

3. Karchevsky A.L. Finite-Difference Coefficient Inverse Problem and Properties of the Misfit Functional / Journal of Inverse and Ill-posed Problems,

1998, v. 6, n. 5, p. 431-452.

4. Karchevsky A.L., Yakhno V.G. One-dimensional inverse problems for systems of elasticity with a source of explosive type / Journal of Inverse and Ill-posed Problems, 1999, v. 7, n. 4, p. 347-364.

5. Карчевский A. Jl. О поведении функционала невязки для одномерной гиперболической обратной задачи / Сибирский журнал вычислительной математики, 1999, т. 2, № 2, с. 137-160.

6. Karchevsky A.L. Numerical solution of the inverse problem for the system of elasticity for vertically inhomogeneous medium / Bulletin of the Novosibirsk Computing Center, series Mathematical Modeling in Geophysics, n. 5,

1999, p. 63-69.

7. Карчевский А.Л. Численное решение одномерной обратной задачи для системы упругости / Доклады РАН, 2000, т. 375, № 2, с. 235-238.

8. Карчевский А.Л., Фатьянов А.Г. Численное решение обратной задачи для системы упругости с последействием для вертикально неоднородной среды / Сибирский журнал вычислительной математики, 2001, т. 4, № 3, с. 259-269.

9. Karchevsky A.L. Several remarks on numerical solution of the one-dimensional coefficient inverse problem / Journal of Inverse and Ill-posed Problems, 2002, v. 10, n. 4, p. 361-384.

10. Karchevsky A.L. The analytical formulas for the gradient of the residual functional for the coefficient inverse problem for the elasticity system / Journal of Inverse and Ill-posed Problems, 2003, v. 11, n. 6, p. 619-629.

11. Карчевский A.JI. Численное определение функции последействия среды Ц Обратные задачи и информационные технологии, 2003, т. 2, № 1, с. 53-63.

12. Kurpinar Е., Karchevsky A.L. Numerical solution of the inverse problem for the elasticity system for horizontally stratified media Ц Inverse Problems, 2004, v. 20, n. 3, p. 953-976.

13. Karchevsky A.L. Numerical reconstruction of medium parameters of member of thin anisotropic layers Ц Journal Inverse and 111-Posed Problems,

2004, v. 12, n. 5, p. 519-534.

14. Карчевский A.JI. Метод численного решения системы упругости для горизонтально-слоистой анизотропной среды / Геология и Геофизика,

2005, т. 46, Xs 3, с. 339-351.

15. Курпинар Э., Карчевский А.Л. Вычисление градиента при оптимизационном методе решения обратной динамической задачи сейсмики для горизонтально-слоистой среды / Геология и Геофизика, 2005, т. 46, № 4, с. 439-447.

16. Карчевский А.Л. Прямая динамическая задача сейсмики для горизонтально-слоистых сред Ц Сибирские Электронные Математические Известия, 2005, т. 2, с. 23-61 (http://semr.math.nsc.ru/cr2005.html).

¡í

1 è

»202 51

РНБ Русский фонд

2006-4 18688

Подписано к печати 26 сентября 2005г.

Тираж 120 экз. Заказ № 1679. Отпечатано "Документ-Сервис", 630090, Новосибирск, Институтская 4/1, тел. 335-66-00

Содержание диссертации, доктора физико-математических наук, Карчевский, Андрей Леонидович

Введение

0.1 Общая характеристика работы

0.2 Решения поставленных задач на современном этапе

0.2.1 Обратная задача

0.2.2 Прямая задача

0.2.3 Свойства функционала невязки

0.2.4 Градиент функционала невязки

Глава 1 Прямая задача.

1.1 Основные уравнения

1.1.1 Постановка прямой задачи,

1.1.2 Постановка прямой задачи, ф 1.2 Представление решения прямой задачи

1.3 Переход к ДМУР

1.4 Выражение для следа

1.5 Частное решение ДМУР

1.6 Решение матричного уравнения

1.7 Построение матриц Си С.

1.8 Выражение для решения ДМУР

1.9 Полезное упрощение

1.10 Вопросы численной реализации.

1.10.1 Решение Проблемы

• Метод Хичкока

• Априорная информация о решении. о . из физики о . из математики о . из практики

• Численный эксперимент

1.10.2 Решение Проблемы

1.10.3 Решение Проблемы

1.11 Порядок вычислений

1.12 Численное решение прямой задачи

Ф 1.12.1 Формула для

1.12.2 Вычисление матрицанта

1.12.3 Формулы для решения прямой задачи

1.13 Обсуждение и сравнение алгоритмов

Ф 1.14 Примеры

Глава 2 Свойства функционала невязки.

2.1 Постановка обратной задачи

2.1 Свойства решения прямой задачи и функционала невязки.

2.2.1 Модельные среды

• 2.2.2 Влияние параметра а

2.2.3 Влияние параметра v

2.2.4 Поведение функционала невязки

2.2.5 Выводы из численных экспериментов. ф 2.3 Градиент функционала невязки.

2.3.1 Выражение для градиента функционала невязки

2.3.2 Аналитические выражения для V^(0, и,р)

2.3.3 Порядок вычислений

2.3.4 Использование конечных разностей

Глава 3 Обратная задача

3.1 Уточнение постановки обратной задачи

3.2 Изотропная среда

3.3 Орторомбическая среда

3.3.1 Свойства решения прямой задачи и функционала невязки

• Расщепление обратной задачи.

• Ограниченность области спектральных значений

• Овражность функционала невязки.

3.3.1 Численный эксперимент

3.4 Трансверсально-изотропная среда

3.4.1 Свойства решения прямой задачи

3.4.2 Численный эксперимент

3.5 Выводы

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

0.1 Общая характеристика работы

Актуальность. Как правило, использование математического моделирования в практической сейсморазведке ограничено. При этом в широко распространенных обрабатывающих и интерпретационных комплексах программ оно ограничивается акустическими моделями. Это связано с тем, что до настоящего времени основной целью стандартной обработки является получение временных разрезов, которые в некотором приближении могут I соответствовать реакции акустической модели среды на нормальное падение плоской продольной волны. На данном представлении основаны многие интерпретационные подходы, например, псевдоакустический каротаж. Однако, при реализации этого подхода происходит значительная потеря информации о среде, содержащаяся в исходных сейсмограммах. Желание не терять данную информацию, а использовать ее при интерпретации, приводит к необходимости использования способов обработки сейсмических данных, в которых параметры среды определяются по полному волновому полю. В этом случае акустическое приближение становится не адекватным и должно быть заменено моделью упругой среды с возможным учетом анизотропии и поглощения в зависимости о ситуации. Естественно, возникает целесообразность использования не только вертикальной, но и горизонтальных компонент вектора смещений, в частности вовлечения в обработку волн PS (хотя возможность определения упругих характеристик среды по лишь одной из компонент остается интересной и в этом случае).

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

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

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

При решении обратной задачи в частотной области встают две основные проблемы.

Первая проблема — переход к образам интегральных преобразований от реальных сейсмограмм.

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

Для функционала невязки, во-первых, необходимо проведение исследований

• на наличие локальных минимумов и максимумов,

• на выявление характера овражности,

• на выяснение возможности достижения точки глобального минимума,

• на установление характера зависимости функционала невязки от его параметров.

• Какая часть спектральных данных пригодна для численного решения обратной задачи?

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

• создания устойчивого вычислительного алгоритма для решения данной задачи.

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

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

F(t)V5(xhx2,x3-x*3). (1)

Источник вида (1) есть центр давления, который является моделью взрыва.

Для перехода в частотную область к СДУ теории упругости применено преобразование Фурье по пространственным переменным Х\ и х2 и преобразование Лапласа по временнбй переменной t. В случае изотропной среды СДУ теории упругости может быть записана в цилиндрической системе координат, тогда для перехода в частотную область по пространственной переменной г применено преобразование Фурье-Бесселя, а по временнбй переменной t — преобразование Лапласа.

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

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

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

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

1. Прямая задача:

• построено явное аналитическое решение СДУ теории упругости в частотной области для горизонтально-слоистых однородных сред любого!вида анизотропии в любой точке среды, не имеющее ограничений на толщину слоев (модель среды может содержать как очень толстые, так и очень тонкие слои), в основе которого лежит редукция СДУ второго порядка к ДМУР, для которого также построено решение в явной аналитической форме в каждом слое среды;

2. Функционал невязки:

• на основе знания аналитической структуры решения прямой задачи и .численного эксперимента установлен характер поведения функционала невязки в зависимости от своих параметров;

3. Градиент функционала невязки:

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

- объединение в единую постановку прямой задачи для решения СДУ теории упругости и его приращений,

4. Обратная задача: На основе знания аналитической структуры решения прямой задачи и численного исследования зависимости обратной задачи и функционала невязки от своих параметров

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

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

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

Здесь необходимо сделать следующие пояснения. Как уже отмечалось выше, для создания алгоритмов решения обратной динамической задачи сейсмики необходимо решать задачи первой проблемы. Входными данными для обратной задачи являются сейсмограммы. Для того чтобы на практике могли быть использованы предложенные в данной работе методы, необходим алгоритм перехода от сейсмограмм к их образам интегральных преобразований. В работе показано, что частотная область для построения функционала невязки делится на две части: доступную и недоступную. Недоступная область содержит особые точки, которые связаны со спектром дифференциального оператора СДУ теории упругости, доступная область является регулярной частью образа. Данное обстоятельство должно облегчить задачу перехода от сейсмограмм к их образам, поскольку будет необходимо добиться только хорошего совпадения теоретического и практического спектров регулярной части образа. Для решения данной задачи может быть применим метод близких операторов [62].

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

Результаты работы докладывались на следующих научных конференциях:

• II Международная конференция по математическому моделированию, Якутск, 28 июня - 2 июля, 1997.

• III Сибирский конгресс по прикладной и индустриальной математике, Новосибирск, 22-25 июня, 1998.

• Международная конференция "Обратные задачи математической фи-•ф, зики", Новосибирск, 21-25 сентября, 1998.

• VI конференция "Обратные и некорректно поставленные задачи", Москва, 20-21 июня, 2000.

• Международная конференция "Некорректные и обратные задачи фи-^ зики", Новосибирск, 5-9 августа, 2002.

• Международная конференция по Математическим Методам в Геофизике (ММГТ2003), Новосибирск, 8-12 октября, 2003.

• International Conference "Inverse Problems: Modeling and Simulation", Fethiye, Turkey, June 07-12, 2004.

• Международная конференция по вычислительной математике (МКВМ-2004), Новосибирск, 21-25 июня, 2004.

Результаты работы вошли в отчеты по следующим грантам:

• Гранты РФФИ: 05-01-00171, 05-01-00559, 02-01-00818, 02-01-00809, 99-01-00563, 96-01-01937, 93-01-01739.

• Грантам научных школ: НШ-1172.2003.1, 00-15-96183, 96-15-96284.

• Грантам Министерства образования РФ: УР.04.01.026, N-00-1.0-73.

Результаты работы докладывались на научных семинарах:

• • Семинар лаборатории обратных задач математической физики Института математики им. С.Л.Соболева СО РАН (руководитель семинара ^ д.ф.-м.н. Ю.Е.Аниконов).

• Семинар лаборатории численных методов решения обратных задач Института математики им. С.Л.Соболева СО РАН (руководитель семинара д.ф.-м.н! А.Л.Бухгейм).

0.2 Решения поставленных задач на современном этапе

Прежде всего необходимо отметить, что первые работы по многоволновой сейсморазведке появились в России. Большой вклад в эту область был внесен работами Пузырева Н.Н. [139-142] и его учеников (см., например, [125, 126, 325]). В работах Пузырева Н.Н. [143-145] дан обзор становления и развития многоволновой сейсмики в России. Идеи многоволновой сейсмики быстро развивались, создавалась соответствующая техника. За рубежом эти идеи также получили широкое развитие и применение (см., например, [234, 267, 320-322, 327, 383], естественно, что это далеко не полный список работ).

0.2.1 Обратная задача

Систематическая разработка численных методов по решению обратной динамической задачи сейсмики в различных ее постановках началась с работ Алексеева А.С. [4, 5]. Дальнейшее развитие заложенных там идей нашло продолжение в работах Бородаевой Н.М., Антоненко О.Ф., Добринско-го В.И., Чеверды В.А., Авдеева А.В. и др. (см., например, работы [1, 6-9, 15, 16, 32, 45-47, 127, 235-233, 305-307, 315, 318]).

Теоретическому исследованию обратных задач, частным случаем которых является обратная динамическая задача сейсмики, развитию и обоснованию численных методов по их решению были посвящены работы Али-фанова О.М., Артюхина Е.А., Румянцева С.В. [10], Аниконова Ю.Е. [1114], Баева А.В. [21, 22, 242-244], Благовещенского А.С. [25-31], Бухгейма A.JI. [35, 36], Васина В.В. [39-41, 75], Гольдина С.В. [58-65], Денисова A.M. [67, 68, 259], Дмитриева В.И. [72, 73], Искендерова А.Д. [78-80, 264, 265, 274], Кабанихина С.И. [81-100, 275-296] (полный список работ [297]), Карчевско-го A.JI. [101-104, 299-304, 316], Клибанова М.В. (список основных работ

313]), Крейна М.Г. [107, 108], Лаврентьева М.М. [111-115, 317, 318], Ма-датова А.Г., Митрофанова Г.М., Середы В.-А.И. [117-119], Мартакова С.В. [42, 97-99, 120, 348, 255, 256, 323] Романова В.Г. [182, 195], Тихонова А.Н. и его учеников (см. сборник избранных трудов Тихонова А.Н. [213] и приведенную там полную библиографию его работ), Щеглова А.Ю. [223-227], Эпова М.И. в соавторстве с Антоновым Ю.Н., Дашевским Ю.А. Ельцовым И.Н. (см., например, работы [17-20, 66, 74, 124, 228, 260, 261, 391]), Яхно В.Г. [229, 302], их соавторов и многих других. В особенности необходимо отметить работы Романова В.Г. Им был развит метод выделения особенностей и сведения коэффициентной обратной задачи к системе нелинейных уравнений Вольтерра второго рода, что позволило исследовать вопросы корректности широкого круга обратных задач [149-209, 335-354].

Параллельно с российским исследованиями работы по численному решению, теоретическому исследованию обратных задач и обоснованию численных методов по их решению шли и за рубежом: Chavent G., Lailly Р. и их соавторы [23', 245-247, 250, 252, 253] (полный список работ Chavent G. [254]), Clarke T.J. [257], Cook D.A. и Schneider W.A. [258], Hildebrand S.I. и McMechan G.A. [268, 269], Isakov V. [273], Kormendi F., Dietrich M. [314], Rakesh [331-334], Sacks P., Santosa F., Symes W.W. и их соавторы [248, 262, 266, 319, 329, 330, 355-363, 367, 370-374] (полный список работ Santosa F. [364] и Symes W.W. [375]), Song Н., Ma Z., Zhang G. [366], Tarantola А. и его соавторы [376-381], Uhlmann G. (список работ и их файлы [382]), Ursin В. и его соавторы [384-386] (полный список работ Ursin В. [387]).

Особняком от выше перечисленных работ стоят работы Тихонова А.Н. [212], Гласко В.Б. [50-53] и Silvester J., Winerbrenner D., Gylys-Colwell F. [369], которые наследовали вопрос единственности решения поставленных в работах обратных динамических задач в частотной области.

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

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

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

Как правило, в математических постановках обратных задач для определения неизвестных функций использовались комбинированные источники: вертикальный удар в комбинации с центром вращения (Алексеев А.С. [5]), наклонный удар (Благовещенский А.С., [25, 28, 30]; Романов В.Г. [163, 182]; Яхно В.Г. [229]), шнуровые и площадные источники (Яхно В.Г. [229]), направленный взрыв (Карчевский A.JI. и Яхно В.Г. [302]).

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

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

0.2.2 Прямая задача

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

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

Большой вклад в понимание процессов распространения волн в слоистых средах, в развитие методов вычислений для таких сред был внесен такими учеными — как Бреховских JI.M. [33, 34], Молотков Л.А. [122], Петрашень Г.И. [133-138], Ризниченко Ю.В. [148].

Существует масса методов решения СДУ теории упругости, например, метод Томсена-Хаскелла и несколько его модификаций [122, 231] (значительный вклад в развитие и применение на практике матричного метода внесен Молотковым Л.А.), конечно-разностный метод (см., например, монографию Иванова Г.В., Волчкова Ю.М., Вогульского И.О., Анисимова С.А., Кургузова В.Д. [76] и цитируемую там литературу), метод конечных элементов, и прочие. Сочетание различных подходов при решении дифференциальных уравнений теории упругости для вертикально неоднородных анизотропных сред можно найти в работах Мартынова В.Н. и Михайленко Б.Г. [121, 324]. На основе лучевого подхода и неполного разделения переменных значительные результаты для решения прямой задачи для слоистых анизотропных сред получены в работах Петрашеня Г.И. [136, 137], Оболенце-вой И.Р. [128-132]. При решении обратной динамической задачи сейсмики, когда прямая задача решается много раз и пространственные масштабы велики, требуется метод решения, который дает значения требуемых величин быстро. Упрощение модели среды — считается, что среда является горизонтально-слоистой и однородной — помогает построить такой алгоритм.

Одним из первых технологичных алгоритмов послойного пересчета для решения ДУ второго порядка для горизонтально-слоистой однородной среды был алгоритм Тихонова А.Н. и Шахсуварова Д.Н. [211]. Однако он имел некоторые ограничения: при его численной реализации существовали выражения, записанные с участием экспонент, имеющих показатели с положительными действительными частями, что приводило к накоплению ошибок округления при послойном пересчете.

Далее, идея послойного пересчета была реализована в следующем виде. Хорошо известно, что ДУ или СДУ второго порядка могут быть сведены с помощью специальной замены функций к дифференциальному уравнению Риккати (ДУР) или дифференциальному матричному уравнению Риккати (ДМУР). Замечательным является тот факт, что когда коэффициенты ДУР или ДМУР являются постоянными, тогда ДУР или ДМУР имеют решения, которые могут быть записаны в аналитическом виде. По всей видимости, впервые для построения численных алгоритмов, активно применяемых в геофизической практике, этот прием был использован в работе Дмитриева В.И. [70] для ДУ второго порядка для решения прямой задачи электроразведки. Теперь уже стало общепризнано [365], что метод послойного пересчета является наиболее подходящим при численном решении обратной задачи при помощи метода минимизации функционала невязки в случае, когда среда является горизонтально-слоистой и однородной.

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

Пусть имеем следующую прямую задачу: uzz — x2(z)u = 0, z £ [0, оо), (2)

Uz\z=O = uq, и -t 0 (zоо), (3) uz]z=zk = 0, [u]z=Zk = 0. (4)

Функция x(z) является кусочно постоянной функцией. Будем обозначать нщ — значения функции k(z) в интервале [zk~i, zt], точки Zk — точки разрыва среды, к = l,Ni, значение X[Nt+i] будет соответствовать значению функции x{z) в полупространстве [zjyn со). Здесь использовано обозначение [u]Zk = u(zk + 0) — u{zk — 0) для значения скачка функции u(z) в точке Введем функцию s(z) следующим равенством: uz = su. (5)

Подстановка (5) в (2) даст нам ДУР, которому удовлетворяет функция s(z) s'+s2 = >г. (в)

Из условий склейки (3) получаем:

М* = 0- (7)

Таким образом, если мы умеем решать ДУР (6), то значение функции на поверхности z — 0 из первого краевого условия находится немедленно:

I 0 I /о\ u\z=Q = -q , s = s\z=0, (8) s а решение прямой задачи (2)-(4) дается выражением: u(z) = -^ехр < / s{x)dx > . (9)

U J

Известно [298], что ДУР имеет три решения. ДУР с постоянными коэффициентами замечательно тем, что имеет решения, которые могут быть выписаны в аналитическом виде. В частности, эти решения для уравнения (6) в каждом интервале [zk-i,zk], где функция >c{z) постоянна, имеют вид: , [** + + - x[fc]] s(z) = Хщ, S(z) = -Я[к], I где Rеящ < 0, и sk — известное значение функции s(z) в точке Zk- Эти решения не могут совпадать ни в одной точке 2 6 [zk-ъ Zk\.

Метод послойного пересчета заключается в следующем. Положим s(z) = -х[Лгг+:1], г е [zNl, оо). i

Это значение позволит удовлетворить второе краевое условие (3). Далее, используя условие склейки (7), получим

Ni — I |

5 — s\z=zNl-0 — S\z=zNi+o.

Используя sNl как начальное условие для решения ДУР (10), положив z = z^i-1, получим значение s|z=ZNi1+o: дЦ + >сщ]е2х№(Щ-1-Щ) + [sNi Я[Щ] s- [sNi + [sNl X[iV;]] •

Используем снова условия склейки (7), получаем значение

Ni-l — I I

S — S\z=zNil-Q — S|2=ZNj1+o? и вычисляем значение используя решения ДУР (10), положив z = ZNt2.

Таким образом, вычислительный процесс для уравнения (6) с условиями склейки (7) может быть записан коротко в следующем виде: ,-1 „ [** + + [sk X[fc]]

11)

S = ->c[Nl+l].

Нетрудно получить значение u(z) из (9). Пусть z G [zji,zj], тогда интеграл распадается на сумму интегралов

171 1 Л Л s(x)dx = ^ / s(x)dx -f / s(x)dx, (zo = 0).

0 1 zm-1

Очевидно, что s(x)dx = *м [ [si + + V - ^

W UJ J [si + - [Si - X(i]] In ([sj + - [SJ - - + const.

Таким образом нетрудно получить для u(z) из (9) следующее выражение:

- ^ [si + - [si - 6 ' G 1 ri 2x[m]e>f[mi(zm-1"z™) m= jJi [sm + ЛГ[т]]е2хн(^-1"2-) - [s™ - '

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

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

Процесс послойного пересчета (11) имеет также название метода "прогонки", который описан в работе Гельфанда И.М. и Локуциевского О.В. [49] в приложении к разностным методам решения ДУ.

Идеи работы Дмитриева В.И. [70] развивались далее им и его соавторами для других задач электродинамики (см., например, [71]). Для задач теории упругости эта методика была применена в работах Аккуратова Г.В. и Дмитриева В.И. [2, 3] для получения следа решения СДУ теории упругости на поверхности z = 0. Авторы рассмотрели СДУ теории упругости в терминах смещений для горизонтально-слоистой однородной изотропной среды, а затем перешли к СДУ для потенциалов. В данном случае уже СДУ сводилась к ДМУР. Замечательно было то, что полученное ДМУР, у которого матрицы-коэффициенты были постоянны, также имело решение в аналитическом виде.

Далее свое применение для СДУ теории упругости алгоритм послойного пересчета получйл в работах Фатьянова А.Г. и Михайленко Б.Г. [215-217]. В работе [215] рассуждения ведутся для получения следа решения СДУ теории упругости для горизонтально-слоистой однородной изотропной среды с поглощением, в работе [216] — для трансверсально-изотропной среды, когда ось симметрии бесконечного порядка направлена по оси Oz.

Представленный в работе метод решения прямой задачи для СДУ тео

U)2 V2 // рии упругости в частотной области для горизонтально-слоистых однородных анизотропных сред был опубликован в работах [105, 106].

0.2.3 Свойства функционала невязки

По теоретическому исследованию свойств функционала невязки автору известно несколько работ. Это результат Чеверды В.А. [221], работа Каба-нихина С.И. [285] и работы автора [101, 299], в которых были развиты и расширены результаты работы [285].

В работе [221] функционал невязки был записан в области пространственной и временнбй частот в следующем виде:

U)2 Vi п(ш, v)\2dujdv, (12)

Ui 1/1 где n(co,v) — невязка, то есть разность между наблюденными данными и решением прямой задачи, v и и — пространственная и временная круговая частоты соответственно.

В работах [101, 285, 299] функционал невязки был записан в области временнбй переменной в следующем виде: т

N(t)\2dt, (13) о где N(t) — невязка, t — временная переменная.

Результаты работ [285, 299] были распространены на обратные задачи для дифференциально-разностных уравнений в работах Кабанихина С.И. и Баканова Г.В. [286] и автора [300].

Во всех выше приведенных работах было показано, что в предположении существования решения обратной задачи функционал невязки (12) или (13) имеет единственную стационарную точку, которая является точкой глобального минимума функционала невязки и решением обратной задачи. В работах [101, 299, 300] для итерационного метода минимизации функционала невязки типа (13) на основе градиентного метода получена оценка скорости сходимости экспоненциального вида. i /

Тем не менее, для задач практической сейсморазведки данные результаты использовать не удается. На практике области изменения пространственных и временных переменных очень велики, нельзя использовать непрерывI ные области изменения пространственных и временной частот.

В диссертационной работе обратная задача решается в частотной области. Функционал невязки имеет следующий вид:

2, (14) здесь и некоторые конечные числовые множества значений временной и пространственной частот.

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

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

Однако, начиная с работ [102, 103, 301] по численному решению обратной задачи для СДУ теории упругости, было замечено, что функционал невязки имеет различное поведение для различного набора параметров обратной задачи, которые используются для его построения. Это были следующие параметры: значения пространственной и временной частот, число пространственных и временных частот, интервалы из которых они взяты, значение параметра затухания а.

Под параметром затухания а понимается следующее. В работах [102, 103, 301] вместо преобразования Фурье по временной переменной использовалось преобразование Лапласа с параметром преобразования р = —а + гш. Первоначальной причиной, по которой было использовано преобразование Лапласа вместо преобразования Фурье, было наличие параметра а, который являлся дополнительным параметром обратной задачи и находился в наших руках, можно было попытаться, выбирая его различные значения, оказать влияние на поведение функционала невязки.

Дополнительным аргументом к использованию преобразования Лапласа могут служить работы Хачай О.А. с соавторами [218-220], в которых комплексируются методы сейсмо- и электроразведки. В этих работах используется преобразование Лапласа с действительным параметром преобразования р = —а по временной переменной.

Исследование зависимости поведения функционала невязки от своих параметров применительно к обратным задачам для ДУ теории упругости впервые было начато в работе [303] и продолжено в работе [316]. После численных экспериментов предложена стратегия минимизации функционала невязки, которая позволяет при его минимизации достичь точки глобального минимума, не попав в "ловушки" локальных минимумов.

Другой подход для численного решения задачи минимизации функционала, который может иметь локальные минимумы и максимумы, предложен в работах Клибанова М.В. и Тимонова А. [308-312]. Он основан на использовании аппарата карлемановских оценок и на сведении обратной задачи к серии задач минимизации для строго выпуклых функционалов. Метод предложен для ряда обратных задач электроразведки, которые ставятся пока для одного дифференциального уравнения, однако, по всей видимости, не существует ограничений для распространения данного метода на обратные задачи для СДУ.

0.2.4 Градиент функционала невязки

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

В огромном количестве работ, посвященных исследованию и численному решению различных обратных задач, приведены, как часть исследования или решения, алгоритмы вычисления градиента функционала невязки. Общим методом для получения вычислительных формул для градиента служит метод, использующий постановку сопряженной задачи (см., например, [37, 376-381]). Метод является общим и подходит для любых постановок коэффициентных обратных задач. Выражения для градиента функционала невязки получаются в виде интегралов от решений прямой и сопряженной задач. Очевидно, если для выражений, стоящих под интегралами, не удается получить первообразных в аналитическом виде, данные выражения приходится интегрировать численно. Следовательно, при вычислении градиента на точность вычислений будет влиять ошибка метода численного интегрирования. На практике области интегрирования могут быть значительными, и, чтобы получить выражение градиента с приемлемой точностью, придется использовать приемы и методы численного интегрирования, значитель-j но увеличивающие время счета. Необходимо отметить, если используется какой-либо метод сопряженных направлений, точность вычисления градиента сильно влияет на правильность построения сопряженных направлений и в целом на скорость сходимости процесса минимизации [38]. Как правило, в коэффициентных обратных задачах функционал невязки, с одной стороны, имеет овражный характер, а с другой стороны, в окрестности точки глобального минимума является "пологим". Для функционалов с таким поведением точность вычисления градиента имеет принципиальное значение.

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

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

Этот результат является прямым следствием предложенного в этой работе метода решения прямой задачи для СДУ теории упругости для горизонтально-слоистой однородной анизотропной среды. Решение ДМУР дано в виде, который оказался очень удобным для получения аналитических формул для градиента функционала невязки. Впервые идея метода получения аналитических формул для градиента функционала невязки в случае горизонтально-слоистой однородной среды была предложена в работе [303], но не имея общего метода нахождения решения ДМУР, получить аналитические формулы для градиента функционала невязки для коэффициентных обратных задач для более сложных дифференциальных уравнений и систем второго порядка, чем было рассмотрено в [303], не представлялось возможным. Основной результат по получению аналитических формул для градиента функционала невязки в случае СДУ опубликован в работе [304], а затем развит в работе [110].

Заключение Диссертация по теме "Геофизика, геофизические методы поисков полезных ископаемых", Карчевский, Андрей Леонидович

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

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

Выше приведенные пункты определяют направления дальнейших исследований.,

Заключение I

В работе предлагается решать обратную динамическую задачу сейсмики в частотной области методом минимизации функционала невязки. Известно, что скорость решения обратной задачи методом минимизации функционала невязки напрямую зависит от скорости решения прямой задачи, поэтому Глава 1 данной работы посвящена разработке метода решения прямой задачи. Было уделено максимальное внимание вычислительным вопросам, связанным со скоростью и устойчивостью определения необходимых величин. Далее, немаловажным условием успешного решения задачи минимизации является поведение функционала невязки. Из представленных в работе примеров было видно, что функционал невязки может иметь множество локальных минимумов и максимумов. Таким образом, разработка стратегии минимизации, позволяющей достичь глобальный минимум функционала, являлась важным этапом исследований. Важной проблемой минимизации является задача вычисления градиента функционала невязки, что позволяет применять градиентные методы для поиска минимума функционала. Этим двум вопросам была посвящена Глава 2 настоящей работы. И последнее, после того как имеются в наличии метод решения прямой задачи, изучены свойства функционала невязки, необходимо проведения исследования основных математических свойств обратной задачи. Глава 3 содержит теоретическое и численное исследование данных свойств обратной задачи. На базе знаний этих свойств могут быть построены алгоритмы решения обратных задач для конкретных постановок.

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

• Прямая задача

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

Не существует ограничений на мощности слоёв: модель среды может содержать как очень толстые, так и очень тонкие слои.

Идея построения метода решения прямой задачи отличается от известных аналогов.

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

Алгоритм решения прямой задачи может быть легко распараллелен.

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

• Функционал невязки

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

Для минимизации функционала невязки предлагается стратегия минимизации, которая приводит к глобальному минимуму, минуя "ловушки" локальных минимумов.

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

• Градиент функционала невязки.

Метод построения градиента функционала невязки позволяет получить явные аналитические представления для компонент градиента. В работе метод продемонстрирован на примере обратной задачи для СДУ теории упругости для горизонтально-слоистой однородной изотропной среды.

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

• Обратная задача

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

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

Библиография Диссертация по наукам о земле, доктора физико-математических наук, Карчевский, Андрей Леонидович, Новосибирск

1. Модельные представления/Геологил и геофизика, 1991, 10, с. 97-106. [118] Мадатов А.Г., Митрофанов Г.М., Середа В.-А.И. Анпроксимационный подход при динамическом анализе многоканальных сейсмограмм,

2. Оценивание параметров// Геология и геофизика, 1991, 11, с. 117-127. [119] Мадатов А.Г., Митрофанов Г.М., Середа В.-А.И. Апнроксимационный подход при динамическом анализе многоканальных сейсмограмм,

3. Прикладные аспекты/Геология и геофизика, 1992, 4, с. 112-122. [120] Мартаков С В Численное решение одномерной обратной задачи в нестационарном диэлектрическом каротаже/ Труды международного семинара "Обратные задачи геофизики", Новосибирск: ВЦ СО РАН, 1996, с. 126-129. [121] Мартынов В.Н., Михайленко Б.Г. Расчет полных волновых полей в анизотронных средах/Математическое моделирование в геофизике,

4. Clarke T.J. Full reconstruction of a layered elastic medium from P-SV slant-stack data/ Geophysical Journal of Royal Astronomical Society, 1984, v. 78, p. 775-

5. Cook D.A., Schneider W.A. Generalized linear inversion of reflection seismic data. Geophysics, 1983, v. 48, p. 665-

6. Denisov A.M. Solution uniqueness in some inverse problems for heat conduction equation with piece-wise coefficients/ Journal of Galculation Mathematics and Mathematical Physics, 1982, v. 22, p. 858-

7. European Association of Geoscientists and Engineers. [268] Hildebrand S.I., McMechan G.A. 1-D seismic inversion of dual wavefleld data: Part I, Nonuniqueness and stabihty/ Geophysics, 1994, v.59, N 5, p. 782-788.

8. Русский перевод: Камке Э. Справочник по обыкновенным дифференциальным уравнениям. 1976, 576 с. [299] Karchevsky A.L. Properties of the misfit functional for a nonlinear one-dimensional coefficient hyperbolic inverse problem/ Journal of Inverse and Ill-Posed Problems, 1997, V. 5, n. 2, p. 139-163. [300] Karchevsky A.L. Finite-difference coefficient inverse problem and properties of the misfit functional/JowrnaZ of Inverse and Ill-Posed Problems, 1998, v. 6, n. 5, p. 431-452. М.: Наука,

9. European Association of Geoscientists and Engineers. [321] Lou M., Pham D., Lee S. Anisotropic parameters estimation from P- and PS-converted wave ddid.// EAGE 64-th Conference and Exhibition, Florence, Italy, 27-30 May, 2002, Extended Abstract Book, P

10. Русский перевод: Новацкий В. Теория упругости. М.: Мир, 1975, 872 с. [327] Perez М.А., Grechka V., Michelena R.J. Fracture detection in a carbonate reservoir using a variety of seismic methods/ Geophysics, 1999, v. 64, n. 4, p. 1266-1276. [328] Polak E. Computational Methods in Optimization, A Unified Approach. Math, in Science and Engineering, v. 77, Academic Press, New York London, 1

11. Русский перевод: Полак Э. Численные методы оптимизации, amplitude/ Geophysics 2002, v. 67, p. 167-176. единый подход. М.: Мир, 1974, 376 с. [329] Qian, J., and Symes, W.W. Adaptive finite difference method for traveltime and

12. Rakesh, Symes W.W. Uniqueness for an inverse problem for the wave equation/ Communication in PDE, 1988, v. 13, p. 87-

13. Rakesh. Reconstruction for an inverse problem for the wave equation/ Inverse Problems, 1990, V. 6, p. 91-

14. Rakesh. An inverse problem for the wave equation in the half plane/ Inverse Problems, 1993, V. 9, p. 433-

15. Rakesh. An inverse problem for a layered medium with a point source/ Inverse Problems, 2003, v. 19, p. 497-506. [335] Romanov V.G. Problemes inverses pour les systemes hyperboliques/Methodes Mathematiques de Linformatique 7 (etude numerique des grands systemes). 1

16. Iitml [388] Wilkinson J.H. The Algebraic Eigenvalue Problem. Oxford, Clarendon Press, 1

Информация о работе

Карчевский, Андрей Леонидович
доктора физико-математических наук
Новосибирск, 2005
ВАК 25.00.10

Диссертация

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

Автореферат

Похожие работы