Содержание к диссертации
Введение
Глава 1. Численное решение задачи Лэмба для неоднородных неупругих сред Больцмана с экспоненциальными функциями последействия 31
1.1. Постановка задач 31
1.2 Сведение к дифференциальным уравнениям 35
1.3. Конечные интегральные преобразования по пространственной переменной в задачах распространения неупругих волн 37
1.4. Некоторые вопросы сходимости и оценка точности метода 39
Глава 2. Полуаналитический метод решения прямых динамических задач для различных моделей слоистых сред и источников 44
2.1. Конечные интегральные преобразования по пространственной и временной переменным в задаче Лэмба для сред Больцмана с произвольными функциями последействия 44
2.2. Расчет полных волновых полей в неупругих средах 45
2.3. Методы решения краевых задач, полученных после отделения переменных 49
2.4. Полуаналитический метод расчёта волновых полей 55
2.5. Нестационарные. волновые поля в анизотропных средах с поглощением энергии 63
2.6. Постановки задач и некоторые вопросы реализации полуаналитического метода в средах с гравитацией 78
2.7. Волновые поля в слоистых пористых средах 86
2.8. Волновые поля от точечных источников дислокационного типа и источников конечных размеров 97
2.9. Расчёт однократных и монотипных волн без использования коэффициентов отражения и сравнение с лучевым методом 105
2.10. Волновые поля в разномасштабных средах 109
Глава 3. Метод расчета функции Грина в многомерно-неоднородных средах 116
3.1. Постановка задачи 117
3.2. Энергетический метод расчета функции Грина для многомерно-неоднородных моделей сред 119
3.3. Волновые поля в средах с криволинейными границами 126
3.4. Некоторые численные эксперименты и вопросы реализации метода для различных сред 135
Глава 4. Некоторые обратные динамические задачи 147
4.1. Оптимизационный метод одновременного определения скорости и декремента поглощения для акустических сред 147
4.1.1. Постановка задачи 148
4.1.2. Вычисление градиента функционала и построение итерационного процесса 152
4.1.3. Исследование чувствительности функционала и некоторые результаты расчётов 155
4.2. Обратная динамическая задача определения дислокационных параметров 162
4.2.1. Восстановление компонент направленной силы 163
4.2.2. Одновременное определение компонент тензора и входного импульса 166
4.2.3. Некоторые примеры численных расчётов 168
Глава 5. Численное моделирование волновых полей для некоторых моделей неоднородных сред 174
5.1. Технологические вопросы математического моделирования волновых полей 174
5.2. Волновые явления в средах с поглощением энергии 178
5.3. Моделирование волновых процессов в анизотропных средах с поглощением энергии 195
5.4. Волновые поля в сложно построенных средах 207
5.4.1. Моделирование вибросейсмических волновых полей 207
5.4.2. Вертикальное сейсмическое профилирование и межскважинное просвечивание 217
5.4.3. Волновые поля сложно построенных сред Сибири 226
Заключение 232
Литература 234
Акт о внедренииe
- Конечные интегральные преобразования по пространственной переменной в задачах распространения неупругих волн
- Расчет полных волновых полей в неупругих средах
- Энергетический метод расчета функции Грина для многомерно-неоднородных моделей сред
- Вычисление градиента функционала и построение итерационного процесса
Введение к работе
Настоящая работа посвящена разработке методов моделирования волновых полей для широкого класса геофизических сред, написанию комплекса программ и исследованию на этой основе физики распространения сейсмических волн в сложно построенных средах. Совершенствование систем наблюдения и повышение детализации экспериментальной информации приводит в настоящее время к необходимости рассмотрения волновых полей с большей точностью и для всё более возрастающих времён и баз наблюдения. Это, в свою очередь, требует дальнейшего развития методов решения соответствующих задач. Отметим, что методы моделирования с использованием высокопроизводительных компьютеров часто оказываются единственным способом получения информации об исследуемых явлениях, так как возможности теоретического и экспериментального изучения этих явлений ограничены методическими и техническими трудностями.
При изучении законов распространения сейсмических волн фундаментальное значение имеет выбор модели, описывающей основные сейсмические явления. Широкое распространение получила модель упругих сред. Эта физическая идеализация оправдана во многих случаях, но требует усовершенствования при исследовании спектрального состава, затухания колебаний в различных типах волн и т.п. за счёт неидеальной упругости среды. В данной работе поглощение описывается в рамках наиболее общей модели среды Больцмана с упругим последействием [1].
В настоящее время существует большое количество различных методов решения прямых динамических задач сейсмики. Одним из наиболее популярных, является лучевой метод, предложенный в работе [79]. В последствие он развивался во многих работах.
6 Приведём только некоторые из них: [28], [80], [81], [82], [83], [88], [157]. Лучевой метод с учетом только нулевого члена лучевого ряда требует меньше вычислительных затрат чем другие методы. Кроме того, с помощью лучевого метода можно легко учесть вклад тех или иных сейсмических волн в формировании сложного волнового поля. Однако, с развитием высокоточной широкополосной сейсмологической аппаратуры, появились факты регистрации «нелучевых» сейсмических волн, которые не описываются нулевым членом лучевого ряда, и для их вычисления необходимо учесть последующие члены ряда. В работах [141], [143] такие «нелучевые» волны были обнаружены с помощью чсленно-аналитических методов моделирования сейсмических полей.
Традиционно исследование волновых полей в слоистых средах является объектом многочисленных исследований. Это обуславливается многими причинами. Отметим некоторые из них. С точки зрения приложений многие практические ситуации можно достоверно моделировать слоистыми структурами. Слоистая среда, кроме того, может использоваться для тестирования и проверки схем решения прямых задач для сред сложной геометрии. Далее, появление все более высокопроизводительного компьютерного обеспечения позволяет решать обратные задачи оптимизационным методом для реальных баз наблюдений и частот.
В настоящее время существует ряд методов расчета волновых полей в слоистых средах. Укажем на некоторые из них. Еще в 50-х годах Г.И. Петрашень ([21]) построил точное решение для слоисто-однородных моделей сред в замкнутом виде. Полученное в виде многомерных интегралов решение было исследовано асимптотическими методами. Сделано большое число физических выводов о протекании в них физических процессов [22]. После стандартных процедур преобразования физических переменных по времени и пространству посредством преобразования Фурье и Бесселя (в предположении, например, радиальной симметрии) уравнения динамической теории вязко упругости сводятся к системе обыкновенных дифференциальных уравнений, которые решаются для дискретных частот и волновых чисел. Эти уравнения решаются в основном разновидностями матричного метода ([23]-[27]) либо асимптотически ([28]-[30]). Широко распространён также "рефлективный" метод ([31]-[32]) расчёта части волнового поля. Все эти методы имеют свою область применимости. В матричном методе, например, при расчёте волновых полей для высоких частот (тонких слоев) происходит потеря точности. Для преодоления этого приходится рассматривать матрицы более высоких порядков, чем это необходимо по физической постановке задачи [27]. Отметим, что в матричном методе эта проблема не решена принципиально [33]. Асимптотические методы при всей их компактности и элегантности имеют ограниченную лучевыми представлениями область применимости. При этом развитие их на более сложные модели (например, с поглощением энергии) требует значительных усилий. Рефлективный метод позволяет надёжно рассчитывать только часть волнового поля, соответствующую, например, объёмным волнам.
В начале 70-х годов в Вычислительном Центре СО РАН (в последствие институт Вычислительной Математики и Математической Геофизики) стали развиваться численно-аналитические методы решения задач сейсмики, основанные на расщеплении двумерных и трехмерных задач на серию независимых одномерных с помощью интегральных преобразований по горизонтальным координатам и с последующим решением их конечно-разностным методом ([6],[90]). Применение этих методов показало, что для получения приемлемых по точности результатов при численной дискретизации необходимо брать не меньше 20 точек на минимальную длину волны для расчёта волновых полей на расстояния порядка 100 длин волн, что приводило к большим объёмам оперативной памяти. Применение интегрального преобразования Лагерра по времени с последующей дискретизацией пространственных переменных позволяет снизить вычислительные затраты [91-92]. Однако уменьшение точности при расчёте сейсмических волновых полей на большие расстояния приводит к необходимости развития, например, конечно-разностных методов с высоким порядком аппроксимации [158]. Кроме того, применение преобразования Лагерра для сложно построенных сред, например для сред с последействием, приводит к громоздким выкладкам. При этом теряется смысл частоты колебаний, а тем самым наглядность при исследовании физики распространения волн в сложно построенных средах.
В данной работе решение прямых динамических задач осуществляется путём соответствующих аналитических преобразований с последующей дифференциальной факторизацией, скалярной в случае одномерных сред и матричной в случае произвольной геометрии. Это позволяет проводить расчёты без потери точности на расстояния свыше тысячи длин волн для широкого класса геофизических сред. Кроме того, предлагаемые алгоритмы позволяют эффективно решать обратные динамические задачи оптимизационным методом.
Работа состоит из пяти глав, заключения, списка литературы и приложения. В первой главе приводятся постановки прямых динамических задач, рассматриваемых в данной работе. Далее приводится метод решения задачи Лэмба для неоднородных неупругих сред в случае произвольной релаксационной модели. Исходные задачи с использованием экспоненциальности функций последействия сводятся к дифференциальным уравнениям более высокого порядка. Далее, с использованием преобразования Фурье-Бесселя исходная задача сводится к однопараметрическому семейству краевых задач меньшей размерности, которые решаются конечно-разностными методами для различных значений параметра суммирования в формулах обращения преобразования Фурье-Бесселя. Исследуются вопросы сходимости и оценка точности алгоритма.
Во второй главе предложен полуаналитический метод решения прямых динамических задач для различных моделей слоистых сред и источников. Сначала выясняются возможности и условия применения конечных интегральных преобразований Фурье по t для решения задач из первой главы. После отделения пространственной и временной переменных получено двухпараметрическое семейство одномерных краевых задач, которые для дискретного набора параметров к, и лп можно решать численно с использованием современных компьютеров. На выбор численного метода при этом принципиальное значение имеют быстродействие и точность. Алгоритм должен быть устойчивым и мало зависеть от случайных ошибок, которые неизбежны при расчетах волновых полей на большие расстояния и для больших частот. Требуется также правильно передавать решения краевой задачи с кусочно-непрерывными коэффициентами, так как литологически реальные среды имеют блочную структуру и параметры среды, а вместе с ними коэффициенты краевых задач, рвутся на некоторых горизонтах. Кроме того, алгоритм должен быть эффективным в плане затрат компьютерных ресурсов. В силу этого, интегро-интерполяционным методом для сильно контрастных сред строится трехточечная разностная схема со вторым порядком аппроксимации, дополненная известными условиями сопряжения в точках разрыва коэффициентов. Далее для нахождения решения используется метод прогонки.
В случае слоистых сред предлагается алгоритм, позволяющий рассчитывать полные волновые поля для слоисто однородных сред Больцмана с произвольными функциями последействия. Количество слоев, их толщины, параметры среды могут быть произвольными величинами. Вычисления при этом строятся рекуррентно, по однотипным формулам, что допускает естественное распараллеливание на многопроцессорных системах. Заметим, что К. Аки и П. Ричарде ([33]) говорят о нецелесообразности перехода к потенциалам из за возникающих трудностей. В данном методе переход к потенциалам не только не усложняет задачу, но приводит к её значительному упрощению. Полуаналитический метод развит для расчёта любых линейных моделей геофизических сред. В работе рассматриваются нестационарные волновые поля в анизотропных средах с поглощением энергии. В основу теории расчёта волн в анизотропных неупругих средах взяты определяющие соотношения Вольтерра [41], учитывающие влияние упругого последействия. А именно, упругие постоянные си заменяются интегральными операторами q: cvxscttx(t)-cv fh^t-fyWdr, в которых с'у и htj определяют, соответственно, уровень и спектральный состав поглощения энергии. В практике оперируют с
11 коэффициентами поглощения различных типов волн, декрементами затухания и другими аналогичными величинами. Поэтому для моделирования неупругих эффектов разработаны методы определения c'yVi hv по различным физическим параметрам поглощения. В связи с этим вначале рассматриваются спектральные характеристики поглощения для "просто" неупругих сред. В качестве модельных при определении коэффициентов поглощения выбраны одномерные уравнения распространения дги„ д2и„ плоских продольных и поперечных волн: Л—j- = p- dz2 r dt2 'Я2 я2 м—^ = /7-^. В этом случае, часто встречающаяся на практике линейная зависимость коэффициента поглощения от частоты аппроксимировалась с использованием подхода, предложенного в [45]. В этой работе отказываются от частотно независимого поглощения, но с затуханием, обеспечивающем его практическое постоянство в некотором частотном диапазоне. Как показывают многочисленные экспериментальные исследования, в ограниченном диапазоне частот для широкого класса горных пород характерна близкая к линейной зависимость коэффициентов поглощения продольных и поперечных волн от частоты. При этом исключается случай сильно водо-насыщенных пород, в которых эта зависимость близка к квадратичной. Здесь уже более правомерно рассматривать случай пористой среды, заполненной жидким флюидом, что так же рассмотрено в работе. Модель Максвелла не даёт линейной зависимости от частоты [46]. Показано, что это отклонение от линейности тем меньше, чем уже спектр входного сигнала по отношению к доминирующей частоте. Т.е. для входных импульсов с не очень широким спектром модель
Максвелла практически даёт линейность коэффициента от частоты и отсутствие дисперсии, часто наблюдаемую на практике. При введении физических параметров поглощения в анизотропной среде учитывается, что в ней характер распространения волн зависит от угла между направлением распространения и осью Oz. Показано что для замкнутого описания поглощения в рамках линейной теории наследственности требуется знать пять декрементов поглощения: АХР, д1|Я, д^, д||5, аор. Отметим, что теорию поглощения в анизотропной среде можно строить, зная декременты (коэффициенты) поглощения, например, для нескольких направлений. Данный подход, на взгляд автора, представляется более оптимальным, так как позволяет использовать обширную справочную литературу по поглощению в упругих средах. Введение поглощения в анизотропных средах, в первом приближении, можно осуществить следующим образом. Предположим, что декременты поглощения квазипродольных и квазипоперечных волн совпадают с декрементами продольных и поперечных волн: АХР = дпя = Д0/) = Ар, д^ = д||5 = д5.
В этом случае (если нет других данных) для моделирования неупругих анизотропных волн можно использовать достоверные значения декрементов поглощения для упругих сред. Кроме того, в силу наибольшей общности модели с последействием можно добиться заранее заданных частотных зависимостей для разных углов падения. Аналогичным образом в работе определяются неупругие коэффициенты c]j для широко распространённых в геофизике моделей Гуревича, Дерягина и т.д.
Далее во второй главе рассматриваются постановки задач и некоторые вопросы реализации полуаналитического метода в среде с гравитацией. В своей фундаментальной работе [49] Ляв рассчитал статические деформации и малые колебания однородного гравитирующего сжимаемого шара. Выделим работы [50]-[51], в которой охватываются единым методом исследования проблемы общей теории упругих колебаний Земли. В данной работе в рамках слоистой модели Земли в сферической системе координат рассматриваются уравнения упругих колебаний гравитирующего шара [49] с различными (в том числе тензором сейсмического момента) источниками. Решение находится с помощью аналитических преобразований после применения ортогонального векторного разложения [33]. Данный метод применён также для расчёта гравитационных волн цунами в рамках широко известной модели Подъяпольского. Отметим, что здесь удалось одновременно рассчитать как упругие, так и гравитационные волн, что не было известно раньше. Расчёт гравитационных волн требует учитывать то обстоятельство, что их скорость существенно (десятки раз) ниже скорости упругих волн в твёрдой Земле. Это приводило к необходимости измельчать пространственный шаг в конечно-разностных методах, которые ранее использовались для решения этой задачи. Всё это приводило к значительному увеличению объёма вычислительной работы, что делало эту задачу неподъёмной. Полуаналитический метод, в силу его большей эффективности, позволил сделать это. Важное значение в задачах разведочной геофизики имеет численное моделирование распространения упругих волн в слоистых средах с пластами пористых, насыщенных вязкой жидкостью пород. Для описания волновых процессов в таких средах чаще всего используются уравнения теории Био в том или ином виде. В данной работе рассматриваются уравнения Био для расширенного частотного диапазона, приведенные в [59]. Развит метод расчёта волновых полей и проведено моделирование процессов распространения пористых колебаний. Отметим, что слоистая упруго-жидкая модель является частным случаем модели Био [87], число неравных коэффициентов которой в изотропном случае сокращается до трёх [159].
Данный подход применим для рассмотрения источников дислокационного типа и источников конечных размеров. Для описания этих явлений используется понятие тензора сейсмического момента [33]. В работе решена задача моделирования волновых полей для такого источника в анизотропной среде. Прямая динамическая задача в этом случае сведена к четырём задачам типа P-SV и двум для SH волн [70]. Это даёт возможность устойчивого и точного расчёта волновых полей в средах с произвольным количеством слоев и с произвольным отношением длин волн различных типов к мощностям слоев для произвольных тензоров сейсмического момента. Анализ данных, проведённых многими исследователями, показывает, что источник землетрясений часто имеет протяжённые размеры по пространству. Для учёта этого явления в работе сконструирован источник конечных размеров. Показано, что он является пространственным аналогом классического источника типа центра давления. Отметим, что таким образом с помощью развитой методики можно рассчитывать волновые поля для произвольно заданного источника конечных размеров.
В задачах моделирования волновых полей в сложно построенных средах большое, а зачастую и определяющее значение имеют алгоритмы, позволяющие рассчитывать динамику отдельно взятых волн. В настоящее время единственным методом численного анализа волнового поля по частям является асимптотический лучевой метод. Его применение, однако, имеет известные ограничения. В работе впервые предлагается алгоритм расчёта однократных и монотипных волн для слоисто- неоднородных сред с произвольным числом слоев на основе специальных разложений точных решений не имеющий ограничений лучевого метода. В [77] проведено сравнение лучевого метода с полуаналитическим в случае однократных продольных и поперечных волн [78]. Показано, что отличия наблюдаются не только в точке выхода головной волны, как считалось ранее, но и в некоторой области, зависящей от длительности входного импульса. При этом как показывают численные расчёты, продольные волны лучевым методом рассчитываются точнее, чем поперечные.
Неоднородная среда часто устроена так, что её характерные размеры, например, по глубине меняются в широких пределах. Известно, что в этом случае практически невозможно рассматривать такую задачу в рамках детерминированного подхода. Для учёта разномасштабной неоднородности реальной среды в работе использован статистический подход и получено уравнение распространения сферических сейсмических волн в сложно построенных средах. Разработан алгоритм решения уравнения для распространения разномасштабных волн.
Для построения модели разномасштабной неоднородности используется телеграфный случайный процесс [160]. Отличие данного подхода от других исследований заключается в том, что так построенная разномасштабная среда позволяет точно рассматривать среднее поле для мелко- и крупномасштабных неоднородностей. Отметим, также, что в разномасштабной среде наблюдается наличие параметрического резонанса, как и в уравнениях для микронеоднородных сред, полученных в [161].
В третьей главе предлагается метод расчета функции Грина в многомерно-неоднородных средах. Существующие методы расчёта в таких средах в общем можно разделить на
16 аналитические и численные. Все их привести в данной работе представляется невозможным. Отметим только некоторые из них [79]-[92]. В данной главе развивается метод расчёта функции Грина для многомерно-неоднородных моделей сред. Он основан на спектральных разложениях решения по пространственным и временным частотам. Поскольку полного разделения переменных здесь уже не происходит, возникают обыкновенные матричные дифференциальные уравнения. Их решения строятся аналитически. В итоге получен алгоритм, позволяющий моделировать волновые поля для блоковой геометрии сред. Как показали численные эксперименты, после соответствующей аппроксимации, данный алгоритм позволяет также моделировать произвольные границы. Расчёт функции Грина позволяет, в случае необходимости, размножать волновые поля для произвольного количества источников без дополнительных вычислительных затрат, что актуально, например, для задач площадной сейсморазведки и решения обратных задач оптимизационными методами.
Данный метод без принципиальных изменений переносится на все постановки из главы 2. Отметим только, что в матричном виде данный подход в этом случае полностью сохраняет свою структуру. Меняются только в зависимости от рассматриваемой среды соответствующие матрицы. Поскольку сборка решения в физической области осуществляется путём суммирования независимых величин, это даёт возможность эффективной программной реализации метода на многопроцессорных вычислительных комплексах. Это обеспечивает возможность проводить моделирование для любых реальных сред с произвольными геометрией и параметрами [162]-[165].
Энергетический метод расчета функции Грина для многомерно-неоднородных моделей сред с учётом принципа взаимности позволяет учесть наличие слоистой пачки с произвольными величинами мощностей слоев и скоростных параметров [101].
В работе теоретически показано, что число обусловленности соответствующих матриц практически не зависит от свойств среды. Это позволяет говорить о вычислительной независимости метода от структуры и свойств сложно построенной среды. Доказано, что число обусловленности зависит в первой степени от количества удерживаемых членов ряда. Отметим, что при сеточных методах решения аналогичных задач числа обусловленности возникающих при этом матриц величины порядка 0(—) (например, [75]). Это позволяет решать матричные системы стандартными методами без потери точности. В случае произвольного изменения параметров среды по латерали (х) в каждом неоднородном слое соответствующие интегралы вычисляются с помощью процедур интерполяции и сглаживания [90]. В ряде случаев, когда, например, границы рассматриваемой области уходят на бесконечность или когда граничная поверхность имеет геометрически сингулярные точки, например, острые рёбра, можно получить несколько математически корректных решений уравнения, среди которых лишь одно верно описывает исследуемое явление [96]. Дополнительное физическое условие, необходимое для однозначного определения решения, в этом случае известно как условие на ребре [96]. Оно заключается в требовании конечности энергии поля, запасённого в любом конечном объёме в окрестности ребра. В работе доказано, что сингулярные точки для так построенного решения не являются источниками излучения.
В третьем пункте данной главы предлагается метод расчёта волновых полей в средах с криволинейными границами раздела. Такие среды можно рассчитывать методом, описанным выше. Однако для таких границ после применения интегральных преобразований решение удаётся получить в явном виде. Для решения этой задачи вводится мощность слоя. Эта величина вещественна в случае плоской границы и совпадает с обычной мощностью слоя. В случае криволинейной границы расстояние становится комплексным. Использование комплексных величин приводит к тому, что условия ограниченности решения на бесконечности и принцип излучения Зоммерфельда выполняются автоматически. При этом наличие в решении дифракционной составляющей автоматически приводит к выбору нужного решения [166]-[167]. Без принципиальных трудностей данный метод переносится на различные геофизические среды. Учёт, например, анизотропного поглощения производится в соответствие с [34], [42]. Кроме того, можно рассматривать различные распределённые и сосредоточенные источники, например, тензор сейсмического момента [69]. Отметим, что в работе [99] в рассмотрение вводится комплексная фаза, что позволило решить некоторые задачи, к которым лучевой метод в обычной форме неприменим.
Далее в третьей главе рассматриваются некоторые примеры расчёта волновых полей для сложно построенных сред. По выше приведённому алгоритму составлен комплекс программ для проведения численного моделирования волновых полей в средах, возникающих в практических задачах. Расчёты проводились на 64-х разрядных процессорах. В качестве иллюстрации возможностей предлагаемого метода рассматривается моделирование волновых полей для вулканической области. Среда, в том числе и магматическая камера, аппроксимируется с помощью произвольного количества прямоугольных блоков. Так как в итоге находится функция Грина, это даёт возможность моделирования для различных режимов входных сигналов (монохроматическом, свип-сигналов и др.). Приведён пример расчёта волновых полей для очаговой зоны живущего вулкана Эльбрус (Собисевич А.Л. и др., Катастрофические процессы, Москва, 2002). Модель характеризуется пачкой слоев, соответствующей осадочным породам, гранитам и базальтам. Магматическая камера взята в виде квадратного включения. При этом очаговая зона является не полой, а содержит среду с более низкой скоростью.
Рассмотрены дифракционные явления в средах с шероховатыми границами, которые исследуются не только в геофизике, но и других областях современной науки (оптика, акустика, биология и т.д.) [168]. В работе проведено численное исследование волновой картины для сферических волн при изменении основных параметров шероховатости.
На практике требуется проводить моделирование для различных геофизических сред, включая упругие, неупругие, анизотропные и т.д. Данный метод без принципиальных трудностей переносится на все постановки из главы 2. По своей структуре предлагаемый метод состоит в матричном пересчёте, в отличие от скалярного для "чисто" слоистых сред. Иначе говоря, модель среды "сидит" в элементах соответствующих матриц. Таким образом, использование матриц позволяет универсализировать структуру алгоритма для различных моделей сред. В качестве примера в работе приведён расчёт для реальной нефтяной структуры, модель которой взята из промышленного пакета "Tesoral" (Канада).
Для частных геометрий среды использование, например, вместо преобразований Фурье преобразования Бесселя позволяет проводить трёхмерные (3D) расчёты на современной вычислительной технике [102].
Отметим, что разработанные в данной главе методы решения прямых динамических задач в многомерно неоднородных средах позволяют довольно просто рассматривать практические задачи продолжения волновых полей для различных моделей геофизических сред (упругие, неупругие, анизотропные и т.д.). Поскольку в работе решение строится в спектральной области, то для этого достаточно введение соответствующего матричного парциального уравнения, по аналогии со скалярным, рассмотренным Клаербоутом [169]. Это позволяет для произвольных многомерных геофизических сред, в отличие от конечно-разностного метода, без дополнительных аппроксимаций при продолжении волнового поля автоматически учесть, например, только проходящие волны для всех углов падения. Таким образом, в рамках данного подхода естественным образом можно осуществить продолжения волновых полей для любых геофизических сред с произвольной геометрией.
В четвёртой главе рассматриваются некоторые обратные динамические задачи, основанные на методах решения прямых задач приведённых выше. Развитые эффективные методы решения прямых задач позволяют строить решения обратных задач оптимизационным методом. В качестве иллюстрации в данной работе рассмотрены две задачи: определение параметров среды и источника. Рассмотрены обратная коэффициентная задача одновременного определения срорости и декремента поглощения акустических волн и задача определения тензора сейсмического момента для слоистых сред. Задачи определения коэффициентов гиперболических уравнений по некоторой дополнительной информации относятся к некорректным задачам, теория которых была заложена в работах А.Н. Тихонова, В. К. Иванова, М.М. Лаврентьева, А.С. Алексеева, В.Г. Романова и других учёных.
В последнее время в связи с быстрым развитием вычислительной техники одним из наиболее часто используемых методов решения обратных задач является метод минимизации целевого функционала [108-112]. В [113] исследована обратная задача для уравнения акустики в интегральной постановке и получена оценка скорости сходимости в среднем метода наискорейшего спуска.
В работе оптимизационным методом решена задача одновременного определения скорости и декремента поглощения для акустических сред по информации о режиме колебаний на свободной поверхности. Для решения данной задачи использован способ поглощения, описанный в [33]. Как следует из [33], подход Больцмана, если функция последействия получена на основе хорошо установленных механизмов внутреннего трения таким образом, что обеспечивает практически постоянное значение Q на сейсмических частотах, согласуется с логарифмической зависимостью фазовой скорости от частоты. На этой основе получено представление для комплексной скорости в случае модели Максвелла. Искомое решение ищется как точка минимума целевого функционала, характеризующего собой квадратичное отклонение зарегистрированного волнового поля от рассчитанного для текущей модели среды. Получены формулы для градиента целевого функционала в комплексной области и исследована его чувствительность. Это позволило построить итерационный процесс одновременного определения скорости и декремента поглощения в спектральной области. Показано, что за счёт выбора оптимальных диапазонов пространственных частот в зависимости от интервала временных частот, можно повысить чувствительность целевого функционала к низкочастотным компонентам скоростной *' функции и декремента поглощения, как и при восполнении информации за счёт соответствующего выбора пространственных частот в случае отсутствия поглощения [108]-[111]. Для построения точки минимума целевого функционала были использованы комбинации методов сопряжённых градиентов (модификация ф, Флетчера-Ривса) и наискорейшего спуска [117]-[118]. Результаты численных расчётов показали, что наиболее эффективно использовать комбинацию приведённых методов. Приведены результаты численных экспериментов и выяснены основные физические факторы, влияющие на сходимость метода в целом. По зарегистрированным данным на свободной поверхности развит метод определения всех компонент точечного тензора сейсмического момента в предположении, что временная форма сигнала в источнике неизвестна. Решение искомой обратной задачи строится на основе решения прямой, приведённой выше. Доказана устойчивость решения обратной задачи. Теоретические аспекты метода рассматриваются на примере слоистой трансверсально-изотропной неупругой среды. Основные этапы алгоритма рассматриваются на примере
3 реконструкции источника типа точечной направленной силы. Пусть на свободной поверхности z=0 на линии наблюдения в = в0 известен в некоторых пунктах приёма вектор смещения, вызванный действием сосредоточенной силы на некоторой глубине. Обратная задача формулируется следующим образом: по этой информации требуется найти угол действия силы. Решение Ф данной обратной задачи строится оптимизационным методом в ш физической области. Показана положительная определённость гессиана целевого функционала. Это обеспечивает единственность и устойчивость точки глобального минимума.
Далее рассматривается задача определения всех компонент тензора сейсмического момента. Эта задача рассматривалась во многих работах. Отметим только некоторые из них [121-123], в которых среда, как правило, считается однородной. Разработанные в данной работе эффективные алгоритмы построения функции Грина в неоднородных средах позволили снять это ограничение.
В случае если известна форма входного импульса, решение строится полностью аналогично вышеизложенному в физической области. Если входной импульс неизвестен, то решение обратной задачи строится на основе минимизации соответствующего целевого функционала в частотной области. Относительно решения обратной задачи доказана устойчивость процесса нахождения всех компонент тензора. Отметим, что предложенный метод решения обратной задачи позволяет учесть более сложные реальные геологические условия. Для этого достаточно использовать решение соответствующей прямой задачи, полученное выше.
Численное опробование алгоритма решения обратной задачи показывает его хорошую помехоустойчивость к случайным и регулярным помехам (неточное задание параметров среды и источника). Погрешность восстановления компонент тензора за счёт автоматического учёта процедуры накопления при использовании нескольких станций приёма может быть сделана ниже уровня погрешности в данных даже при сильном дефиците исходной информации.
Пятая глава посвящена численному моделированию волновых полей для некоторых моделей неоднородных сред. Здесь приводятся результаты численного моделирования для широкого спектра сред (упругие, неупругие, анизотропные и т.д.) и источников ("центр давления", направленная сила, тензор сейсмического момента).
По алгоритмам, описанным выше, разработан комплекс подпрограмм, позволяющий для сосредоточенных источников широкого класса (центр давления, произвольно ориентированная сила, тензор сейсмического момента) проводить расчёты геофизических полей различной природы в сложно построенных средах на единой программно-алгоритмической базе с учётом поглощения, анизотропии и других усложняющих задачу факторов. В [58] начала осуществляться программная реализация и создание транспортабельного пакета прикладных программ, который может использоваться как самостоятельно для расчёта полных волновых полей при моделировании геофизических процессов в реальных средах, так и в качестве подсистемы решения прямой задачи, встроенной в промышленные обрабатывающие системы. Изначально в ППП в качестве языка программирования был выбран Fortran. Это обеспечивало транспортабельность пакета и возможность в настоящее время без особых трудностей использовать математическое обеспечение, ориентированное на вычисление конкретных специальных функций. Учёт этого, позволяет повысить производительность программ, что имеет большое значение само по себе и особенно возрастает в связи с бурным развитием объектно-ориентированного программирования и его использования для визуализации результата на экране. Моделирование волновых полей проводилось также на современных процессорах (Itanium, Opteron).
На основе разработанных программ исследованы динамические особенности распространения волн в средах с поглощением энергии. Если динамика волн в упругих средах довольно хорошо исследована, то для неупругих сред исследования носят в основном качественный характер [1],[126]-[129].
Экспериментальные исследования приводят, в основном, к независимости декремента поглощения от частоты и практической линейности коэффициента поглощения от частоты [130]-[134]. Эти два факта и взяты в основу при динамическом моделировании сейсмических процессов в неупругих средах.
Приведено количественное сравнение динамических характеристик сейсмических волн для широко распространённых в геофизике моделей вязкого трения, Максвелла, Гуревича, Дерягина для сред с разными декрементами затухания. Экспериментальные данные о поглощении в реальной среде свидетельствуют о большом разбросе д5/д;,. Это отношение может быть как больше, так и меньше единицы [138]. Отметим многочисленные исследования дисперсионных уравнений для неупругих сред в случае различных законов поглощения (например, [134]). Основное отличие данной работы состоит в приведении количественных характеристик отношения доминирующих частот поперечных и продольных волн в зависимости от различных соотношений декрементов поглощения, согласующихся с экспериментальными данными.
Задача отражения и преломления сферических волн и в настоящее время остаётся сложной математической задачей, особенно для тонких слоев. В работе численно исследованы некоторые явления, связанные с распространением неупругих сферических волн. При этом получены новые результаты по динамике волн не только в неупругой, но и упругой средах.
Рассмотрен случай, когда источник колебаний, например, источник типа центра давления, расположен на расстоянии d
Приведены численные снимки (фотографии) волнового поля для упругой и неупругой сред. Исследованы волновые явления в случае двух полупространств при различном расположении источников. В этом случае также возникают интенсивные "нелучевые" волны, обладающие более низкочастотным спектром по сравнению со спектром прямой волны. Проведены численные исследования волновых полей для различных типов источников. Количественно исследованы явления изменения формы колебаний, перераспределения интенсивности между различными типами волн, изменение частотного состава в зависимости от соотношения декрементов затухания.
Рассчитаны интерференционные волновые поля для тонких низкоскоростных слоев, с различным расположением источников. Приведены волновые поля и амплитудный спектр для реальной тонкослоистой структуры. Модель среды представляет собой слоистую пачку из 800 слоев с мощностями порядка 1м. За счёт выбора поглощения удалось добиться побочного максимума на спектре порядка 20гц, наблюдаемого на практике. При расчете волновых полей для таких сред применять сеточные методы нельзя, так как необходимо брать слишком мелкий пространственный шаг. По-видимому, в этом случае единственная возможность это применение выше описанного полуаналитического метода.
Многочисленные экспериментальные наблюдения указывают на наличие поглощения [37] и анизотропии различных типов [38]-[40]. Несомненно, что эти факторы (за исключением структурных) оказывают наиболее существенное влияние на формирование волнового поля. В работе исследованы особенности распространения различных типов волн, а так же проведен анализ анизотропных (классификация М.В. Невского [145]) неупругих эффектов [42]. Обнаружено "нелучевое" явление, что угол раствора "петли" в волновом поле квазипоперечных волн на 8-10 больше угла раствора петли на лучевых индикатрисах соответствующих волн для неупругих анизотропных сред. Исследованы волновые поля для анизотропной упругой и неупругой среды для модели ПИ. Петрашеня ([146]), характеризующейся двумя петлями по осям координат, с различными декрементами поглощения.
Среди систем возбуждения особое место по функциональным возможностям занимают невзрывные источники [147]-[148]. В работе приведены некоторые численные расчёты вибросейсмических полей и проведён анализ полученных результатов для реалистичной сейсмической модели Быстровского полигона (г. Новосибирск). Для этого построен алгоритм и написана программа моделирования виброграмм для реального свип-сигнала. Отметим качественное совпадение реальных и теоретических виброграмм [151].
В экспериментах с мощными вибрационными сейсмоисточниками на больших базах наблюдения впервые обнаружен эффект акустосейсмическои индукции [154]. На расстоянии 20км были зарегистрированы поверхностные волны, возбуждённые акустическим излучением от вибрационного сейсмического источника в атмосфере. В работе данный практический эффект был объяснён численным моделированием в рамках комбинированной модели среды: атмосфера - упругая Земля с распределённым источником из [150].
В последние годы всё большую роль в разведочной нефтяной геофизике приобретает детальное изучение сейсмогеологического разреза по сейсмическим данным. Одними из наиболее результативных методов решения таких задач являются методы вертикального сейсмического профилирования и межскважинного просвечивания. На основе разработанных алгоритмов был создан комплекс программ, позволяющий рассчитывать полные волновые поля на заданной серии вертикальных расстановок приёмников и выдавать моментальные пространственные срезы поля (мгновенные снимки) в произвольные моменты времени. Для задач скважинных исследований нефтегазоносных залежей было проведено моделирование пористых волновых полей для межскважинного просвечивания.
Рассмотрены результаты численного моделирования для сложно построенных сред Юрубчено-Тахомской зоны (правобережье Енисея). В результате численного моделирования удалось добиться хорошего совпадения результатов реального и теоретического ВСП. Разработанные компьютерные технологии моделирования волновых полей для сложно построенных реальных сред впервые позволили выяснить (совместно с ОАО " Енисейгеофизика") причины отсутствия материалов на первичных сейсмограммах.
Каждая глава диссертации имеет независимую нумерацию формул в виде группы из двух чисел: первое число соответствует номеру главы, второе - порядковый номер.
Рисунки, иллюстрирующие рассматриваемые задачи и полученные результаты, располагаются в тексте в местах ссылок на них и имеют соответствующий порядок нумерации. Список цитируемой литературы располагается за заключением, в котором приведены основные результаты диссертации, выносимые на защиту.
Разработанные алгоритмы и программы внедрены в Институте геофизики СО РАН, о также использовались и используются в ОАО "Енисейгеофизика" г. Красноярск, НИРФИ г. Горький, Институт вулканологии г. Петропавловск-Камчатский.
Работа написана по публикациям: [12,16,17,19,20,34], [36,42,53- 58,63,64,69-71,78,94,101,102,108-112,154,162-167]. По материалам, используемым в диссертации, делались доклады на Международных конференциях: "Численные методы интерпретации сейсмических данных" (Суздаль, 1980), "Точные асимптотические и стохастические методы в геофизике" (Санкт- Петербург, 1981), "Модельная оптимизация в исследовательской геофизике"(Западный Берлин, 1991), "Математические методы в геофизике" (Новосибирск, 2003), "Международная конференция по вычислительной математике" (Новосибирск, 2004),"
Вычислительные методы и решение оптимизационных задач" (Новосибирск, 2004), " The International symposium on mathematical modeling of dynamic processes in atmosphere, ocean, and solid Earth" (Novosibirsk, 2004), "International Workshop on Active Monitoring in the зо solid Earth Geophysics (IWAM04)" (Japan, 2004). Полученные результаты обсуждались на семинарах Института вычислительной математики и математической геофизики и института геофизики СО РАН, Института физики земли РАН," С.-Петербургского отделения математического института РАН.
Автор выражает благодарность член-корр. РАН Б.Г. Михайленко за научные консультации и поддержку при выполнении работы, академику РАН А.С. Алексееву за внимание к работе и обсуждение основных результатов.
Конечные интегральные преобразования по пространственной переменной в задачах распространения неупругих волн
В дальнейшем для определения решения используется метод комплексирования, состоящий в отделении переменной г и последующим численным решением соответствующих редуцированных задач, развитый в [6]. Основные этапы метода проиллюстрируем на примере задачи распространения SH волн (1.15)-(1.17).
Для каждого параметра разделения к решение строится конечно-разностным методом. При этом воспользуемся стандартными обозначениями из [8]. Искомая область покрывается сеткой ши, неравномерной по z и равномерной по t.
Так как разностная схема (1.25) пятислойная, то вводится дополнительное условие v„,(z,.,0) = 0, которое получается дифференцированием по t уравнения (1.21) в предположении дополнительной гладкости искомого решения. Для построения разностной аппроксимации краевого условия (1.23) со вторым порядком по z используем стандартный прием повышения порядка аппроксимации с использованием уравнения движения [8]. В итоге получим искомую аппроксимацию, которая здесь не приводится ввиду её громоздкости и очевидности получения. Пусть сетка mh разбивает полосу вдоль оси Oz на М интервалов. Второе краевое условие (1.23) аппроксимируется в этом случае точно v(zw,O = 0.
В этой полосе решается полученная разностная задача для различных значений параметра разделения к. Окончательно решение определяется суммированием в (1.20).
Основные погрешности при вычислении смещения возникают вследствие погрешности разностной схемы и отбрасывания остатка ряда (1.20). Ниже отдельно оценивается каждая из них в предположении однородности среды и равномерности аппроксимации по z.
При практических расчетах для произвольных у условие (1.28) служило ориентиром выбора шагов разностной схемы (1.25). При этом, как показали численные эксперименты, при увеличении / условие устойчивости становится более жестким и для его выполнения требуется мельчить шаг по времени, что приводит к значительному увеличению объема вычислительной работы. Решив одномерные задачи для дискретного набора пространственных частот, искомое смещение находим суммированием ряда (1.20). При этом возникает вопрос о количестве членов ряда, которое нужно удержать при суммировании, чтобы получить требуемую точность.
В дальнейшем будем предполагать, что временная функция в источнике имеет продолжительность от 0 до I и обладает непрерывными производными до к-го порядка включительно. Решение (1.21)-(1.23) записывается в этом случае в виде свертки.
5( Ш)/«о =- )ut).Ju[kx{t)e- dT 4л-// 0J
Из теории интеграла Фурье известно [11], что порядок убывания интегралов Фурье-Бесселя зависит от гладкости fx(t). Пусть fx(t) п раз дифференцируема и ограниченна на промежутке [0,1] и обращается в нуль на его концах вплоть до производных п-1 порядка. Тогда при t l S(0,k,t) C/k", При - оо
При практических расчетах в качестве f(t) брались финитные функции, обладающие достаточным количеством производных для хорошей сходимости. Использовались также функции близкие к финитным, для которых условия сходимости выполнялись сколь угодно точно. Случай z 0 исследован в [12].
Из теоремы 1. следует, в частности, что заглубление линии наблюдения или источника служит дополнительным фактором улучшения сходимости. Отметим, что в дальнейшем это станет физически очевидным при рассмотрении неоднородных волн в полуаналитическом методе. Эти исследования показывают, что при надлежащем выборе шагов аппроксимации разностной схемы и числа гармоник в (1.20) погрешность нахождения искомой функции может быть сделана сколь угодно малой. Результативные погрешности метода оценивались (// = 0) сравнением с [13].
Задача (1.29)-(1.30) решалась квадратурами. При этом относительная погрешность не превышала 2-3% на расстояниях порядка 50я {я- доминирующая длина волны). Данный подход без принципиальных трудностей распространяется на все постановки п. 1.1. Однако для общей релаксационной модели он приводит к многослойным разностным схемам, что существенно затрудняет его реализацию и увеличивает объем вычислительной работы. Есть и принципиальные ограничения. Так, например, широко известная модель Гуревича ([14]) не может быть сведена к дифференциальному уравнению.
Значительно более эффективным оказался метод, основанный на использовании конечных интегральных преобразований по временной и пространственной переменным. Этот подход позволяет не только вычислять волновые поля на расстояниях, превышающих 1000А, НО и решать искомую задачу не только с экспоненциальными, но и с произвольными функциями последействия.
Расчет полных волновых полей в неупругих средах
В случае источника типа горизонтальной силы нет осевой симметрии. В этом случае в волновом поле будут присутствовать все три компоненты вектора смещения. Учитывая, однако, специфику задачи можно подобрать преобразование, которое сводит задачу к одномерной. Для однородных и слоисто однородных сред решение данной задачи в виде контурных интегралов построено в [4]. Следуя работе [18], решение (1.5)-(1.7) будем искать в виде комбинации рядов Фурье-Бесселя с элементарными тригонометрическими функциями.
После отделения пространственной и временной переменных получено двухпараметрическое семейство одномерных краевых задач, которые для дискретного набора параметров к, и хп можно решать численно с использованием современных PC. На выбор численного метода при этом принципиальное значение имеют быстродействие и точность. Алгоритм должен быть устойчивым и мало зависеть от случайных ошибок, которые неизбежны при расчетах волновых полей на большие расстояния и для больших частот. Требуется также правильно передавать решения краевой задачи с кусочно-непрерывными коэффициентами, так как литологически реальные среды имеют блочную структуру и параметры среды, а вместе с ними коэффициенты краевых задач, рвутся на некоторых горизонтах. Кроме того, алгоритм должен быть эффективным в плане затрат компьютерных ресурсов. В силу этого, интегро-интерполяционным методом строится трехточечная разностная схема со вторым порядком аппроксимации, дополненная известными условиями сопряжения в точках разрыва коэффициентов.
Так как для реальных сред уровень поглощения энергии мал, то из (2.28) видно, что устойчивость зависит главным образом от минимальной скорости в среде и частоты. Для вычисления более высокочастотных компонент поля или для расчета сложно построенных сред требуется уменьшать пространственный шаг. При решении задачи (2.22)-(2.23) строится разностная схема со вторым порядком по пространственной переменной [16],[19] с учетом условий сопряжения на границах разрыва параметров. Далее решение краевой задачи для набора соответствующих параметров строится методом прогонки. Получено условие устойчивости прогонки в зависимости от упругих и неупругих параметров.
Аналогично выписываются и решаются разностные уравнения в случае источника типа горизонтальной силы. Отметим преимущества данного подхода по сравнению с алгоритмом из главы 1. Значительно сокращается объем вычислительной работы. Это происходит вследствие того, что здесь для произвольных моделей поглощения решаются однотипные краевые задачи. Метод же сведения к дифференциальным уравнениям приводит уже для однопараметрической релаксационной модели к уравнению третьего порядка. При численной реализации возникает пятислойная разностная схема. Для общей релаксационной модели это приводит к многослойным разностным схемам, что значительно увеличивает объём вычислительной работы. Важным вопросом при численной реализации является выбор системы тестов, контролирующих точность решения задачи. Проверка алгоритма и программы в целом осуществлялась, прежде всего, расчётом волновых полей для экспоненциальных функций последействия. Далее рассчитывались волновые поля SH волн с различными функциями последействия для однородной среды и сравнивались с точным решением. Особое внимание уделялось проверке точности расчётов в случае, когда параметры среды терпят разрыв первого рода на некоторой глубине. В этом случае построено точное (в спектральной области) решение, которое приведено выше.
Ветвь двузначной функции v выбирается из условия совпадения её значения для действительных значений аргумента с арифметическим значением корня. Иначе говоря, выбирается значение, соответствующее главному значению аргумента. Далее решение (2.31) склеивается на границах разрыва исходя из соответствующих прообразов граничных условий. Рассмотрим случай безграничного полупространства в случае источника горизонтальной силы.
Далее искомое волновое поле находится суммированием (2.30). Данная схема нахождения решения проходит в общем случае системы вязко упругости в случае слоисто-однородной среды для различных типов поверхностных и заглублённых источников [20]. Данный алгоритм удобен и эффективен при реализации, когда среда не содержит большого количества слоев. Кроме того, он служил контролем гарантированной точности для задач, приведенных в главе 1.
Исследование волновых полей в слоистых средах является объектом многочисленных исследований. Это обуславливается многими причинами. Отметим некоторые из них. С точки зрения приложений многие практические ситуации можно достоверно моделировать слоистыми структурами. Слоистая среда, кроме того, может использоваться для тестирования и проверки схем решения прямых задач для сред сложной геометрии. Далее, появление все более высокопроизводительного компьютерного обеспечения позволяет решать обратные задачи оптимизационным методом для реальных баз наблюдений и частот.
В настоящее время существует ряд методов расчета волновых полей в слоистых средах. Укажем на некоторые из них. Еще в 50-х годах ПИ. Петрашень ([21]) построил точное решение для слоисто-однородных моделей сред в замкнутом виде. Полученное в виде многомерных интегралов решение было исследовано асимптотическими методами. Сделано большое число физических выводов о протекании в них физических процессов [22]. После стандартных процедур преобразования физических переменных по времени и пространству посредством преобразования Фурье и Бесселя (в предположении радиальной симметрии) уравнения динамической теории вязко упругости сводятся к системе обыкновенных дифференциальных уравнений, которые решаются для дискретных частот и волновых чисел. Эти уравнения решаются в основном разновидностями матричного метода ([23]-[27]) либо асимптотически ([28]-[30]). Широко распространён также "рефлективный метод ([31]-[32]) расчёта части волнового поля. Все эти методы имеют свою область применимости. В матричном методе, например, при расчёте волновых полей для высоких частот (тонких слоев) происходит потеря точности. Для преодоления этого приходится рассматривать матрицы более высоких порядков, чем это необходимо по физической постановке задачи. Отметим, что эта проблема не решена принципиально [33]. Асимптотические методы при всей их компактности и элегантности имеют ограниченную лучевыми представлениями область применимости. При этом развитие их на более сложные модели (например, с поглощением энергии) требует значительных усилий. Рефлективный метод позволяет надёжно рассчитывать только часть волнового поля, соответствующую, например, объёмным волнам. Ниже предлагается алгоритм свободный от данных ограничений. Им можно рассчитывать полные волновые поля для слоисто однородных. сред Больцмана с произвольными функциями последействия. Количество слоев, их толщины, параметры среды могут быть произвольными величинами.
Энергетический метод расчета функции Грина для многомерно-неоднородных моделей сред
Рассмотрим метод расчета функции Грина для неоднородных моделей сред с использованием принципа взаимности. Среда предполагается состоящей из произвольных блоков. При этом её некоторая область может представлять слоистую пачку с произвольным числом слоев. Мощности слоев и скоростные параметры сред могут быть произвольными величинами, как в слоистой пачке, так и в блоковой среде.
В (3.10) t - время расчёта волнового поля в длинах волн, с -параметр, зависящий только от свойств гладкости входного импульса. Из (3.10) следует, что число обусловленности практически не зависит от свойств среды. Знаменатель в (3.10) может теоретически меняться только от единицы до двух. Это позволяет говорить о вычислительной независимости метода от структуры и свойств сложно построенной среды. Перейдём к построению алгоритма решения задачи (3.4)-(3.5).
Организуя этот процесс по числу вертикальных блоков, получаем функцию Грина для произвольной среды. В случае произвольного изменения параметров среды по латерали (х) в каждом неоднородном слое, интегралы из (3.6) вычисляются с помощью процедур интерполяции и сглаживания, разработанных в [90]. После соответствующих преобразований Фурье получаем функцию Грина g в физической области. В практике часто рассматривают среды, содержащие пачки однородных слоев. Ниже рассматривается такой случай. При построении алгоритма в этом случае используется принцип взаимности.
В ряде случаев, когда, например, границы рассматриваемой области уходят на бесконечность или когда граничная поверхность имеет геометрически сингулярные точки, например, острые рёбра, можно получить несколько математически корректных решений уравнения, среди которых лишь одно верно описывает исследуемое явление [96]. Дополнительное физическое условие, необходимое для однозначного определения решения, в этом случае известно как условие на ребре [96]. Оно заключается в требовании конечности энергии поля, запасённого в любом конечном объёме в окрестности ребра. Это равносильно требованию.
Предлагается метод расчёта волновых полей в средах с криволинейными границами раздела. Такие среды можно рассчитывать методом, описанным выше. Однако для таких границ после применения интегральных преобразований решение удаётся получить в явном виде.
Для решения задачи вводится мощность слоя. Эта величина вещественна в случае плоской границы. В случае криволинейной границы расстояние становится комплексным. Использование комплексных величин приводит к тому, что условия ограниченности решения на бесконечности и принцип излучения Зоммерфельда выполняются автоматически.
В задачах моделирования волновых полей, в частности, геофизических определяющее значение имеет выбор тестов, позволяющих контролировать точность используемых алгоритмов.
Точных решений в настоящее время, особенно для многомерных сред существует немного [97]. В тех же немногих случаях, в которых известно строгое решение задачи, оно, как правило, имеет сложный вид, состоящий из рядов, каждый член которых представляется в виде интегралов. Ниже предлагается аналитический метод расчёта волновых полей в средах с произвольными криволинейными границами. Он служил контролем точности для более сложных задач. Данные среды имеют и самостоятельное значение, например, в нефтяной сейсморазведке. При этом аналитическая форма предлагаемого метода отличается относительной простотой. Решение данной задачи строится на основе введения комплексного расстояния. Это приводит к упрощению алгоритма. Кроме того, при введении комплексных величин автоматически учитываются условия Зоммерфельда.
Условия непрерывности теперь должны выполняются на границе разрыва скорости (3.22), которая разбивает полупространство на две области. В силу этого на границах разрыва параметров вводятся известные условия сопряжения.
Искомое решение ищется в виде преобразования Фурье по переменным t их. Пусть ииЬ временные и пространственные частоты, g- соответствующий образ искомого преобразования. В дальнейшем соответствующие значки и индексы для сокращения записи будут опускаться.
Не трудно показать, что в случае плоской границы экстремум (3.23) будет равен толщине слоя. Данный принцип положен в основу нахождения расстояния для сред с нерегулярными границами. Для построения решения в спектральной области задачи (3.1)-(3.3) рассмотрим вначале случай плоской границы.
Физически это означает, что на волновой картине в этом случае будет присутствовать две петли на годографе. При этом если h = h" две петли будут смыкаться в одну. Это подтверждается численными расчётами. Отметим, что из общих соображений следует, что для произвольной нерегулярной границы возможно наличие произвольного количества петель на годографе. Проанализируем выражение (3.25). Для этого рассмотрим случай однородных волн. В этом случае каждая экспонента в (3.25) будет выглядеть следующим образом.
В (3.38) присутствует экспоненциально затухающий множитель, описывающий дифракционную волну. Таким образом, условия ограниченности решения на бесконечности и принцип излучения Зоммерфельда в данном случае будут выполняться автоматически. На взгляд автора это более физично, т.к. наличие в решении дифракционной составляющей автоматически приводит к выбору нужного решения. Волновая картина и в частности количество петель на годографе определяется количеством вещественных корней соответствующего уравнения. Комплексное расстояние описывает возникающие дифракционные явления.
Отметим, что в работе [99] в рассмотрение вводится комплексная фаза. Это позволило авторам [99] решить некоторые задачи, к которым лучевой метод в обычной форме неприменим. Окончательно, используя обратные преобразования Фурье, получим решение (3.1)-(3.3) в физической области. Отметим, что без принципиальных трудностей данный метод переносится на различные геофизические среды.
Вычисление градиента функционала и построение итерационного процесса
В приведенных выше методах n2{z), УФ, / (Z) являются комплексными величинами, и итерационный процесс идёт одновременно по реальной и мнимой частям комплексной скорости.
Результаты численных расчётов показали, что наиболее эффективно использовать комбинацию приведённых методов. Алгоритм был построен таким образом, что после каждых п- шагов по первому методу делался один шаг по методу наискорейшего спуска. Подправляя, таким образом pJt удалось на каждом интервале частот более близко приближаться к точке глобального минимума.
На сходимость процесса существенно влияет выбор пространственных и временных частот, по которым проводилась оптимизация. Отметим, что в [108-111] подробно рассмотрены вопросы, связанные с восполнением информации за счёт соответствующего выбора пространственных частот в случае отсутствия поглощения. Таким образом, на каждом этапе, задавая пределы по к и й, мы можем расширять и интервал спектра восстанавливаемых параметров. Здесь мы не будем подробно останавливаться на этом. Эта технология восполнения отсутствия информации на низких временных частотах, приводящая к неустойчивости определения трендовой составляющей скоростного состава среды, подробно описана в [108-111]. Ясно, при этом, что гладкие функции, имеющие спектры, близкие к финитным, при правильном выборе интервалов частот будут восстанавливаться быстро и точно. По данному алгоритму была написана программа расчёта и проведены методические и практические расчёты.
Параметры искомой среды, как правило, являются кусочно-непрерывными функциями. Для повышения сходимости в этом случае нужно использовать дополнительные регуляризирующие процедуры (например, сглаживание)
Обычный на практике дефицит исходной информации здесь восполняется за счёт привлечения пространственных частот. На взгляд автора этих строк, для более широкого практического использования обратных задач восполнение недостающей информации должно производится не в спектральной, а во временной области. Это не только более наглядно, но и позволяет более корректно использовать априорную информацию. Практически это осуществляется путём рассмотрения дискретных аналогов соответствующих дифференциальных постановок. Решение в этом случае будет строиться с помощью разложений по дискретным ортогональным функциям. В итоге будут получаться аналогичные уравнения, но уже в физической области. Так, например, пространственная частота к будет заменяться на дискретные собственные числа в физической области [120]. Это позволяет восполнять недостающую информацию, без изменения алгоритма, уже в физической области. Это даёт возможность говорить о реальной возможности применения обратных задач при решении конкретных практических задач.
Всё сказанное, с учётом выше развитой теории решения прямых задач, позволяет говорить о реальности создания замкнутого цикла решения сложных практических задач. Отметим, что в [112] решена задача восстановления упругих параметров в случае, когда поглощение в среде считается известным.
В данной главе по зарегистрированным данным на свободной поверхности развит метод определения тензора 9 компонент точечного тензора сейсмического момента в предположении, что временная форма сигнала в источнике неизвестна. Решение искомой обратной задачи строится на основе решения прямой, приведённой выше. Доказана устойчивость решения обратной задачи. Теоретические аспекты метода рассматриваются на примере слоистой трансверсально-изотропной неупругой среды. Т.е. среды, для которой в данной работе построена функция Грина. Алгоритм же решения обратной задачи в данном случае зависит от зарегистрированного волнового поля, которое может быть произвольной природы. Использование данного метода даёт возможность эффективного сейсмомониторинга очаговых зон и районирования сейсмически опасных регионов.
В силу того, что в данном случае v = o в (4.35) будет строгое неравенство. Это означает, что собственные числа матрицы А строго положительны. Таким образом, найденное решение целевого функционала даёт единственную точку глобального минимума.
Это обеспечивает устойчивость процесса нахождения всех компонент тензора. В настоящее время в большинстве работ компоненты тензора сейсмического момента определяются в основном по записям длиннопериодных составляющих волнового поля. Отчасти это объяснялось недостатком соответствующей экспериментальной базы. В настоящее время, в связи с появлением новых технологий, реальные записи можно проводить в широком спектральном диапазоне. Очевидно, что использование такой информации будет способствовать повышению точности решения обратных задач и в, частности, определения компонент тензора. Отметим далее, что предложенный метод решения обратной задачи позволяет учесть более сложные реальные геологические условия (наличие осадочного чехла, складчатое строение земной коры и т.д.). Для этого в (4.40) достаточно учесть решение соответствующей прямой задачи, полученное выше. Как известно [124-125], сейсмический источник наиболее общего типа может быть описан как совокупность дипольных источников, распределённых в пространстве и времени, с некоторым тензорным моментом m,j(x,t)dVdt. Выше рассмотрен метод определения функции Грина для излучающих объектов конечных размеров. Использование вышеизложенного метода даёт возможность эффективного решения такой задачи в слоистых анизотропных средах с поглощением энергии в "ближней" и "дальней" зонах. Это позволяет решать задачу оперативного мониторинга сейсмоопасных зон земной коры и задачу микрорайонирования зон строительства промышленных объектов при помощи имитации землетрясения как природным, так и искусственным мощным источником вибрационного типа.
