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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

• 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−171, 05−01−559, 02−01−818, 02−01−809, 99−01−563, 96−01−1 937, 93−01−1 739.

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

• Грантам Министерства образования РФ: УР.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).

Uzz=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 /о uz=Q = -q, s = sz=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 — sz=zNl-0 — Sz=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 — Sz=zNil-Q — S|2=ZNj1+o? и вычисляем значение используя решения ДУР (10), положив z = ZNt2.

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

11).

N,.

S = ->c[Nl+l].

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

171 1 Л Л s (x)dx = ^ / s (x)dxf / 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.

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

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

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

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

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

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