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

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

МИХАЙЛОВ Александр Анатольевич

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

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

Автореферат

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

Новосибирск, 2006

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

Научный руководитель: член-корреспондент РАН,

доктор физико-математических наук Михайленко Борис Григорьевич

Официальные оппоненты: доктор физико-математических наук

Имомназаров Холматжон Худайназарович

Ведущая организация: Санкт-Петербургское отделение

Математического института им. В.А. Стеклова РАН.

Защита состоится 19 декабря 2006 года в 15 часов на заседании Диссертационного совета Д 003.061.02 по присуждению ученной степени кандидата физико-математических наук в конференц-зале Института вычислительной математики и математической геофизики СО РАН (630090, Новосибирск, проспект Лаврентьева, 6).

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

Автореферат разослан 19 ноября 2006 года..

| доктор физико-математических наук ; Клем-Мусатов Камилл Давыдович

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

С.Б.Сорокин

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

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

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

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

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

Задачи исследования:

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

2. Построение алгоритмов решения пространственных осесиммет-ричиых задач в смещениях и скоростях / нал ряжен иях для вяз-коупругой слоисто-неоднородной среды.

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

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

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

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

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

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

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

Апробация работы. Результаты, изложенные в диссертационной работе, докладывались и обсуждались на семинаре отдела «Математические задачи геофизике», а также на 5 международных конференциях (1999 -2004 гг.).

Публикации. Основные результаты диссертации опубликованы в 11 печатных работах, в том числе статья в Journal of Computational Acoustics (2001 г.), статья в Journal of Geophysical Prospecting (2003 г.) и статья в Journal of Pure and Applied Geophysics (2003 г.). Список работ помещен в конце автореферата.

Структура и объем работы. Диссертация состоит из введения, трех глав, содержащих 15 параграфов, заключения, двух приложений и списка литературы из 86 наименований. Объем работы—104 машинописных страницы, в работе содержится 26 рисунков.

Благодарности. Хочу выразить благодарность сотрудникам лаборатории «Численное моделирование сейсмических полей» Института вычислительной математики и математической геофизики (ИВМиМГ СО РАН), в частности В.Н. Мартынову и Г.В. Решетовой, за сотрудничество и обмен опытом, накопленным в данной лаборатории, в решении задач математического моделирования.

Особая благодарность моему научному руководителю—заведующему лабораторией, чл.-корр. РАН, д.ф.-м.н. Б.Г. Михайленко, который привлек автора к исследованиям новых методов решения прямых динамических задач сейсмики.

Основное содержание работы

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

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

В п. 1.1 дается математическая постановка задачи для волнового уравнения в полярной системы координат (г,<р). Предлагается рассмотреть описываемый алгоритм декомпозиции областей на примере разбиения исходной области на три участка. На участке [0,Г1] скорость распространения волн задана постоянной. На участке [п, га] скорость задана в виде произвольной двухмерной функции V (г, ¡р). И дня г > гг скорость в среде задана постоянной.

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

При этом на границах пересечения областей вводятся неизвестные функции 01(0 = Д(гь*), <3а(*) = Д(п -Дг,*), <?3(() = Л(г2 + Дг,(), Ф^М = -Кп(г21 *)) которые могут быть определены из условия непрерывности общего решения. После этого для каждого участка среды, используя соответствующий метод (аналитический или численно-аналитический), может быть записано решение с заданными граничными функциями. Затем, используя условие непрерывности, строится система уравнений сшивки найденных решений вида:

Решив данную систему относительно Q^, (¿2, £?3) <?4, можно определить неизвестные граничные функции и в итоге построить общее решение для любой точки исходной области.

В п. 1.3 описывается построение аналитического решения задачи для однородной среды, используя интегральное косинус-Фурье преобразования по пространственной координате <р и преобразования Ханкеля по координате г. В результате этих преобразований на

I

Р(па,<2 г.Оз) =<?1(0. 5„(п - Дг.^ОО =

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

В п. 1.4 излагается численно-аналитический метод решения для неоднородной среды. В основу метода положен алгоритм матричного преобразования системы уравнений, получаемой после косинус-Фурье преобразования по пространственной координате <р и конечно-разностной аппроксимации по координате г. Используя разложение матрицы получаемой системы на собственные вектора и собственные числа, решение задачи на данном участке записывается аналитически по времени в виде интегральной свертки:

L L 1

ЯЮ = Х>мЕ-тг f :(i-r))dr,

i=l j=l J

где Ai — собственные числа, a tjj — элементы матрицы собственных векторов.

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

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

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

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

В п. 1.8 сделаны основные выводы о возможностях предлагаемого алгоритма и его отличительные особенности в сравнении с подобными алгоритмами.

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

В п. 2.1 дается постановка осесимметричной задачи в цилиндрической системе координат для системы уравнений первого порядка в скоростях смещений и напряжениях. Среда задается изотропной вертикально-неоднородной. При этом полагается, что механизм последействия задан в виде интегральных соотношений Больцмана для произвольных функций последействия. В этом случае связь между напряжением и деформацией для произвольной нагрузки о — <r(i) имеет вид:

Здесь Му-упругий модуль, а — функция последействия для этого упругого модуля. В случае изотропной упругой среды ф — О при Мц = /1 для поперечных и Ми = Х + 2/л для продольных волн.

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

оо

Теорема. Пусть две произвольные функции представимы в виде разложения в ряд по функциям Лагерра как

ОО | ОО VI

т = <«)* £ ^-/„(«(ы), фц) = £

Тогда функция = $ /(г)^(г — т)<1т может быть представлена в виде разложения по функциям Лагерра как

где ■

7 1 т

= / (Л*)- = ^ £

О 3=°

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

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

«0 = 1-1;(1-1=1

где количество релаксационных механизмов, а г,| и т„( — времена релаксации деформации и напряжения, соответствующие 2-му механизму.

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

^A(kn) + ^Ejf(kn,m) = P(kn> m-1),

где kn — корни преобразования Бесселя, т — номер гармоники Jla-герра, h — Const—параметр сдвига функций Лагерра.

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

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

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

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

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

циям Бесселя по однородной пространственной координате г вида

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

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

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

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

В п. 2.10 делаются выводы об эффективности применения предлагаемых алгоритмов. Рассматриваются отличительные особенности

. оо

ШЬ;?)} - / {5&::5!}

о

и преимущества в сравнении с другими известными методами решения подобных задач.

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

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

В п. 3.2 описывается теоретический вывод основных формул предлагаемого алгоритма решения поставленной задачи. Алгоритм основывается на сведении исходной интегрально-дифференциальной задачи к системе линейных алгебраических уравнений, используя комплексирование интегральных преобразований и конечно-разностной схемы. Для этого используется интегральное преобразование Лагерра по временной координате и конечное интегральное преобразование Фурье по однородной пространственной координате. В результате исходная задача распадается на N независимых двухмерных задач, где N — количество гармоник в Фурье-

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

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

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

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

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

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

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

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

В приложении 2 показывается принцип построения и вывода основных формул теории вязкоупругости для законов поглощения и дисперсии с заданной функцией последействия.

Заключение

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

На защиту выносятся следующие основные результаты:

1. Совместно с сотрудниками лаборатории «Численного моделирования сейсмических полей> разработаны алгоритмы и лично автором осуществлена программная реализация следующих задач:

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

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

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

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

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

. двум пространственным координатам в комплексе с аналитиче-

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

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

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

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

[1] Конюх Г.В., Михайлов А.А. Об одном алгоритме декомпозиции областей с использованием конечных интегральных преобразований Фурье-Бесселя // Труды ИВМиМГ. Сер. Математическое моделирование в геофизике,—Новосибирск: Изд. ИВМиМГ, 1998.— Вып. 5.— С. 79-90.

[2[ Konyukh G.V., Mikhailenko B.G., Mikhailov А.А. On an algorithm of domain decomposition based on finite Integral Fourier-Bessel transforms // Bulletin of the Novosibirsk Computing Center. Ser, Mathematical Modeling in Geophysics. —Novosibirsk.—1998. —Iss. 4.—P. 103-113.

[3] Конюх Г. В., Михайленко В. Г., Михайлов А. А. Решение прямых динамических задач сейсмшси с помощью интегрального преобразования Лагерра // Об. тр. 8 Всерос. щколы-семияара «Современные проблемы математического моделирования». — Ростов-на-Дону: Изд-ао РГУ, 1999.—С. 107-118.

[4] Konyukh G.V., Mikhailenko B.G., Mikhail о v А.А, Numerical modeling of transient seismic fields in viscoelastic media based on the Laguerre spectral method // Bulletin of the Novosibirsk Computing Center. Ser. Mathematical Modeling in Geophysics. — 2000. — Iss. 6. —P. 31-40.

[5] Конюх Г.В., Михайленко Б.Г., Михайлов А.А, Моделирование нестационарных сейсмических полей в неоднородных вязкоупругих средах на основе преобразования Лагерра // Науки о Земле: современные проблемы сейсмологии.—М.: Вузовская книга, 2001. — С. 25-46.

[6] Конюх Г.В., Михайленко Б.Г., Михайлов А.А. Численное моделирование сейсмических полей в вязкоупругих средах на основе спектрального метода // Математическое моделирование.—2001.—Т. 13, № 2.— С. 61-70.

[7] Konyukh G.V.t Mikhailenko B.G., Mikhailov А.А, Application of the . integral Laguerre transforms for forward seismic modeling // J. of

Computational Acoustics. —2001. —Vol. 9, № 4.—P, 1523-1541.

[8] Mikhailenko B.G., Mikhailov A.A., Reshetova G.V. Seismic wavefield modeling for laterally heterogeneous whole earth models // Proc. of Intern. Conf. on Computational Mathematics. — Novosibirsk, 2002. — P. 632-638.

[9] Mikhailenko B.G., Mikhailov A.A., Reshetova G.V. Numerical viscoelastic modeling by the spectral Laguerre method // J. Geophysical Prospecting. —2003.—Vol. 51.—P. 37-48.

[10] Mikhailenko B.G., Mikhailov A.A., Reshetova G.V. Numerical modeling of transient seismic fields in viscoelastic media based on the Laguerre spectra] method // J. Pure and Applied Geophysics.—2003.—Vol. 160.— P. 1207-1224.

[11] Михайлов А. А. Моделирование сейсмических полей для 2.5D неоднородных вязкоупругих сред // Труды междунар. конференции «Математические методы и геофизике - ММГ 2003».—Новосибирск, 2003.— Ч. 1,—С. 146-152.

МИХАЙЛОВ Александр Анатольевич

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

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

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

Лицензия ИД № 02202 от 30 нюня 2000 г. Подписано в печать 15.11.2006 г.

Формат бумаги 60 у И1/« Объем 1,0 п.л. 0,9 уч.-нзд.я. Тираж 80 экз. Заказ №73 6

ООО «Омега Принт», Новосибирск-90, пр. Лаврентьева, б

ВВЕДЕНИЕ Обзор вычислительных методов при решении задач моделирования сейсмических волновых полей. Актуальность и новизна предлагаемых ашоритмов.

Глава 1. Алгоритм декомпозиции областей с использованием конечных интегральных преобразований Фурье-Бесселя.

§1.1. Постановка задачи для полярной системы координат.

§ 1.2. Теоретические аспекты построения решения.

§ 1.3. Аналитическое решение задачи для однородной среды.

§ 1.4. Численно-аналитическое решение для неоднородной среды.

§ 1.5. 11остроение общей системы для декомпозиции областей и определение решения исходной задачи.

§ 1.6. Некоторые аспекты численных расчетов.

§ 1.7. Результаты численно! о моделирования.

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

В настоящее время существует большое количестно различных методов решения прямых динамических задач сейсмики Одним из наиболее популярных, инляег-<я лучевой метод, провоженный н работах [2], [3], [7] и в последствие развишй и |Ю], |41], (50]. Лучевой метод с учётом нулевого члена лучевою ряда требу о I значительно меньше вычислительных затрат, чем другие методы. Кроме юг о, с помощью лучевою метода можно легко учесть вклад тех или иных сейсмических волн в формировании сложною волновою ноля. Однако, с развитием высокоточной широкополосной сейсмологической аппаратуры, появились факш регистрации "нелучевых" сейсмических волн, коюрые не описываются нулевым членом лучевого ряда, и для их вычисления необходимо учесть последующие члены ряда. В работах [6] и [51| такие "нелучевые" волны были обнаружены с помощью численно-аналитических методов моделирования сейсмических полей.

В начале б()-х годов, с появлением высокопроизводительных ЭВМ, в геофи зике стали использоваться численные методы расчета волновых полей, например, методы конечных разностей и конечных элементов. Возникло новое направление -вычислительная геофизика. Первые результаты применения численных методов в сейсмике показали, что для получения приемлемых но точности результатов при численной дискретизации необходимо 15-20 точек на минимальную длину волны ,1ля расчета волновых полей на расстояние 80-100 длин волн ([34], [35], |37], [53], [54]). Такие требования оказались неприемлемыми, для существовавших в то время ЭВМ, по вычислительным затратам и оперативной памяти даже при решении двухмерных задач сейсмологии и сейсморазведки.

В начале 70-х годов в Вычислительном Центре СО РАН (в последствие ип-сгиту г Вычислительной Математики и Математической Геофизики) стали развиваться численно-аналитические методы решения задач сейсмики, основанные на расщеплении двумерных и трехмерных задач на серию независимых, одномерных с помощью интегральных преобразований по горизонтальным координатам и с после,^ ющим решением их конечно-разностным методом (см. [1], [5], [19], [22], [33]). Такой ио,1\од не требовал большой оперативной памяти, обладал хорошей точностью и был реализован на отечественных ЭВМ типа БЭСМ-б для различных моделей сред: для радиально-неоднородпых - в работе [5], для ани зотропных - в [18], ]С8], ,1ля вязкоупру\их сред - в [30]-[32|. Обобщение вышеуказанного подхода для двумерных и трехмерных моделей сред дано и работах [20], [21], [23], |70]-[72].

Низкая точность и большие вычислительные затраты при расчете сейсмических волновых полей на большие расстояния стандартными численными методами дали толчок к развитию за рубежом конечно-разностных методов с иысоким порядком аппроксимации (см. [36], [42[, [19[, [63]). Кроме тою, получили разви-ше псендоспектральные меюды решения задач сейсмики (см. [60] - [62]). Если в конечно-разностных методах ирои людные аппроксимируются конечно-разностным отношением к дискретных точках фишческою пространства, то в исевдоспек-фальных они определяются на основе разложения и ряды по базисным функциям, которые бесконечное число раз дифференцируемы во всей расчетной обласш. Обычно, в качестве таких базисных функций выбираются триюномсчрические функции, либо полиномы Чебышева [ 18], так как для них существуют эффеюив-ные методы быстрою иреобраювания (БПФ). Важным достоинспюм псевдоспектральных методов является высокая скорость сходимости (теоретически зксионен-ниальная), если решение обладаем достаточной степенью гладкости.

Апа.101 ичной точностью обладаки и спектральные методы, предложенные ранее для решения задач сейсмики в работах [20], [21], [23], [71], [72[. В отличие oi нсевдоспектральных методов в спектральном подходе все вычисления проводятся в спектральном пространстве беï возвращения в каждый момент времени в физическое пространство. Основное время расчета в этом случае идет на вычисления сумм типа свертки, кснорые вычисляются с помощью БПФ или на спецпроцессорах. Сочетание спекiралыюю подхода но юриюнтальным коордишиам и нысокоючпыми разностными схемами с переменным шаюм по вертикальной коордшше и явной разностной схемой но времени оказалось эффективным для решения миошх динамических задач сейсмики (см. [23], [74]).

Все вышеуказанные меюды расчета сейсмических волновых полей основаны на ипюльювании явной разностной схемы вюрого порядка по времени и условно называкнея меюдами расчета во временной области. Дальнейшее развшие ^тих методов связано с попытками использовать более высокий порядок аппроксимации производных но времени. Действительно, спектральный и псевдоспектральный меюды дают очень высокий порядок аппроксимации npocrpaiiciвенных производных, а накопление ошибки при расчею волновых полей на больших расстояниях и временах распространения происходит, преимущественно, за счет второю порядка аппроксимации временных производных. Попытки использован, четверней порядок аппроксимации по времени наылкиваются на большие смраничения на mai дискретизации для тою, чтобы схема была устойчивой. В работе [81| предложено использовать спектральный подход на основе полиномов Чебышева для аппроксимации временною оператора. Меюд обеспечивает высокую ючносп», но требует больших вычислшельных затрат

В работах [13], [55], [56] предложен спектрально-разностный алюритм сведения задачи распространения сейсмических ноли в неоднородных средах к задаче Ко-пш для системы обыкновенных дифференциальных уравнений. Используя меюд матричной декомпозиции, сис юма расщепляется на Дг независимых обыкновенных дифференциальных уравнений, которые решаются анали1ически. Поскольку, декомпошция (диаюнализация) не зависит от пространственною расположения псючника, который задается в правой час1и системы, то эта трудоемкая ироце,;у-ра проводится только один раз для любою положения исючника. Таким образом, отпадает необходимость репинь задачу каждый раз при изменении положения источника. Такой подход является весьма эффективным при моделировании сейсмических волновых нолей в сейсморазведке методом общей ыубинной точки (ОГТ), ко!да необходимо решать задачу для фиксированною строения среды, но при различных (несколько сотен) положениях источника. Как разни 1ие ыкою подхода, в данной диссертации предложен алюритм сшивки аналитических и численно-аиалшическою решения для волновоюу равнения в полярной системы коордшш, резулыаш моделирования, на основе данною алгоритма, были опубликованы в работах [16], [57]

Вторая большая группа методов расчеы сейсмических волновых нолей основана на различных подходах в частотной обла(чи. После применения преобразования Фурье по времени мы получаем уравнения или систему уравнений Гельмюльца. Эти уравнения необходимо решип» для набора временных частот, количество ко-юрых зависит от ширины спектра сшнала в исючнике, затем зги решения просуммировать, чтобы получить импульсную сейсмо1рамму. Наиболее известными меюдами решения в часюшои области является метод Томсона-Хаскелла [17], [84], в дальнейшем развитый в работе [25], а также "рефлективный" метод [16], [77|. В работе [32] предложен иолу аналитический алюритм, основанный на использовании преобразования Фу рье-Бесселя но юри зонтальной координате и введением новых функций, которые сводят решение редуцированной задачи к решению сисчемы уравнений Рикатти. В каждом однородном слое решение преде твлжч-ся аналитически в виде комбинации экспонент только с отрицательным знаком. Просше формулы пересчета волновою поля с одною слоя на друюй выводяк-я на основе 1раничных условий.

Волыное число работ посвящено численным методам решения уравнений ниш Гельмюльца, (см., например, [65]-[67]). Заметим, что численные методы расчет сейсмических волновых полей в частотной области имеют свои особенности. По сравнению с методами расчета во временной области здесь нег 01раничений на ша1 но времени, в то время как во временной области он определяется максимальным перепадом скорос ni распространения сейсмических но.in в выбранной модели среды. Кроме юю, при решении прямых динамических задач но временной области для неу прупгх сред с последействием возникают проблемы, свя данные с наличием в у равнениях интегралов пша свертки. Эти трудности преодолеваю к я путем аппроксимации интегралов конечной суммой и введением дополнительной неременной [38], [79], что приводит к существенному увеличению вычислительных затрат. В то время, как в часки ной области эти проблемы легко преодолеваются введением комплексных j n¡>> i их параметров.

Остановимся на особенностях численною решения прямых динамических задач в частотной области. После применения преобразования Фурье но времени и аппроксимации пространавенных нрои зводных, например, конечно-разностным методом, либо методом конечных элементов задача сводится к системе алгебраических уравнений большой размерности Так как, матрица этой системы зависит от временной частоты, то при решении системы матрицу нужно трансформировав для каждой частот, что предс1авляе1 чрезвычайно трудоемкую нроце,1уру. Чтобы уменьшить вычислительные затраты, используют конечно-разностную аппроксимацию пространственных производных высокою порядка точности. При ном, размерность матрицы уменьшается, но растет количество ненулевых диагоналей в самой матрице В последние юды upoipecc в этой области сними с применением высокоточных компактных разностных схем, с использованием коюрых уменьшается размерность матрицы, но не происходит увеличение ненулевых диа-юиалей (см. [52], [80]).

Кардинальное решение вышеуказанной проблемы связано с применением ин-нчральиою преобразования JIaieppa по времени, вмесю преобразования Фурье. Такой подход предложен в работах [12], [73]. В отличие от преобразования Фурье, применение интегральною преобразования JIaieppa по времени с последющей дискретизацией пространственных переменных позволяет свести исходную задачу к решению алгебраической системы уравнений, в которой параметр разделения т ирису icinyeT только в правой част уравнений и имеет рекуррешиую зависимость. В настоящее время вышеуказанный подход развит л'я упругих изотропных [58], [59] и анизотропных [69] неоднородных сред, а также для вязкоунругих неоднородных сред [21], [75], [76].

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

Главная трудность в осуществлении численных методов в вязкоуируюм моделировании - присутствие интеграла свёртки в уравнении движения. В настоящее время су ществу ют дна основных подхода при численном моделировании распро-(Iранения сейсмическою волновою ноля и вязкоуируюй среде. Первый подход основан на численном моделировании сейсмических полей во временной области. В эюм случае, уравнение движения можем быть записано в дифференциальной форме, вводя дополнительные неременные в исходные формулы, которые позво-ляю1 исключит!» из уравнений интсчральные свертки [39], [43], [15]. Полеченные, мким обра юм, динамические уравнения нязкоупруюсти мо1ут быи» решены, неволь?} я конечно-разностный метод аппроксимации производных (см., например, [44], [85]). Также, используется псевдоспектральный метод, 1де наря,^' с дискретизацией но пространству применяется полиномиальная интерполяция ошоешель-но времени [39], [82]. Псевдоспектральные методы обеспечивают более высокою ючносп» решения, но требуют больших вычислительных зацыь

Второй подход связан с численным моделированием сейсмических полей в вя з-коу пруюй среде в частотной обласнь Эюг спектральный подход основан на применении преобразования Фурье по времени к динамическим уравнениям нязкоупруюсти. Чтобы преодолеть трудность численною нредаавления инкчралов свер1ки, можно преобразовать уравнения в частотную область и моделирован» результирующие уравнения типа уравнения Гельмтльца. Это позволяет пред-с ывип» чаеютио-зависимые нарамефы поыощения в виде комплексных у пру I их параметров [78] Использование конечно-ралностною меюда и меюда конечных 1ементов, для полученных таким обраюм уравнений распространения у пру I ой волны в частотной области, было впервые предложено в работе [05] и позже раз-рабоыны в работах [66], [07]. Полуаналишческий меюд решения в чаеюнюй об-ласш для вязкоупруюй слоистой среды был предложен в работе [32], 1де решение исходной задачи сводится к решению сисгемы уравнений Рикагти. После чет, в каждом у пру юм слое, решение задачи может быть представлено анали1ичееки в виде комбинации убывающих экспоненциальных функций. Общее решение сфо-шея по формулам пересчета волновою поля от нижней к верхней ¡ранице на основе заданных 1раничных условий. В случае произвольно неоднородных сред, после тою, как полученные уравнения в частотной области дискрсчи зированы но прос транс Iвенным координатам, общее решение задачи нредставляе1ся в виде систем уравнений чрезвычайно большою порядка для каждою значения частоты. Численное решение такой задачи 1ребуе1 больших вычислительных зафаь Чю-бы облегчить решение задачи и уменьшить обьем вычислений, в работе [80] предложена новая конечно-разноеIпая схема решения в частотной облает, коюрая о( нонлна на операторах вращения.

Чю касается предлагаемых в диссертации алгоритмов моделирования волновых нолей в вязкоу пру Iих неоднородных средах, то их можно рассматривать как аналог спектрального моделирования, где нмосто временной чагнны ш мы имеем номер т - степень многочленов Латерра. Главное отличие данною меюда от спектрального моделирования заключается в том, что после применения ише-1ральн010 преобразования Лагерра к временным производным и для интегралов свергки по времени с последующей дискр(ми зацией но пространственным координатам, мы получаем систему алгебраических уравнений с матрицей, независимой от параметра т. Так как, только правая часть системы имеет реку рренгну то зависимость от параметра т. Таким обраюм, возможно использовать бьк 1рые методы решения систем линейных алгебраических уравнений с большим числом правых сторон, например, на основе разложений но мето,ту Холецкого В этом случае, численные операции по преобразованию матрицы осуществляю 1ся только только один раз в отличии 01 спектрального подхода, когда необходимо преобразовать матрицу для каждой частотной гармоники и!, что требует больших вычислиIель-ных затрат. Сравнивая интегральные преобразования Фурье и Ла1ерра, можно сказан., чю применение преобразования Лагерра, в этом случае, позволяет существенно сократить вычислительные затраты. Для осуществления такою повода к решению динамических задач вязкоупруюсти потребовалось доказать теорему для преобразовании Ла1срра интегральной свертки. Эта теорема подобна н'ореме преобразования Фурье интегралов свергки. Её использование даёт возможность рассматривать самые общие связи меж,^ тензорами напряжения и деформации, выраженными произвольными релаксационными функциями в интегральных соотношениях Больцмана При этом не требуется введение дополнигельных переменных, которые приводят к появлению дополнительных уравнений и тем самым увеличивают размерность решаемой системы. Предлагаемая методика может быть лег ко обьединена с рядом методов для решения полученной, после преобразования Латерра, дифференциальной системы но пространственным координатам, включая использование конечных разностей или спектральные меюды (Фурье, Фурье-Бесселя или преобразование Лежандра). В предлагаемой диссертации такой подход к решению динамических задач вязкоупругости обсуждается на примере решения пространственных осесиммегричных задач в цилиндрических координатах ,ия вертикально-неоднородных сред и для 2.50 моделей сред в Декартовых координатах с релаксационными функциями для стандартного линейною твердого тела.

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

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

Во в юрой ыаве рассматриваю юя два алюритма моделирования волновых нолей в вязко} пр}I их средах для решения пространственных осесиммегричных задач первого порядка в скоростях смещений и напряжениях и для уравнений второю порядка в смещениях. Среда задаеюя изотропной вер I икал ыю-неоднородной При эгом иолатется, что механизм последействия задан в виде ишегральпых соотношений Больцмана л'я произвольных функций последействия. Главная 01-личительная особенность предлаыемых алгоритмов заключается в иснользова-нии интсчральною преобразования Лаюрра но временной координате Построение алюршмов основывается на комнлексировании шгтральных преобразований и конечно-разностных методов ,Ц1Я сведения исходных инте1рально- дифференциальных систем к системам линейных алюбраических }равнений, решении коюрых мо1}т бьпь найдены наиболее эффективными известными численными меч одами (например, типа Ш— разложения). Приводятся результаты численною моделирования для вязко} нр> I их среде функциями последействия для стандартною линейною твердою тема с несколькими релаксационными механизмами. В ¡аключшельном пункте данной ¡лавы делаются выводы об эффективности применения прельщаемых алюритмов, Рассматриваются отличительные особенности и преим}1цества в сравнении с др}1ими известными методами решения подобных задач.

В 1рс1ьей ыаве рассматривается алгоритм моделирования волновых полей для 2.50 неоднородных вязкоунр}1их сред Построение предла! аемою алюригма основывается на комнлексировании инкчральною преобразования Лаюрра но временной координате и численной конечно-разностной схемы аппроксимации прои з-водных по пространственным координатам ,1дя сведения исходной шгшралыкь дифференциальной системы к системам линейных алюбраических уравнений. Подобное сведение задачи к хорошо обновленной системе алгебраических }равнений с множеством правых частей позволяет использовать быстрые алюритмы решения на основе итерационных меюды шиа сопряженных ¡радистов, сходящиеся к решению задачи всею за несколько итераций. На данном этапе алю-ршм был эффективно распараллелен В частности, была реализована распараллеленная версия меюда сопряженных ¡радиентов. На уровне входных данных, при задании модели среды, это равносильно декомношции исходной обдаст на миожеспю подобласюй, равных количеству процессоров. Это даем вошожность распродо.нчшя намят, как при задании входных параметров модели, ык и при дальнейшей численной реалтации алюршма в подобластях. В $аключительных н>нк1ах ыавы приводятся результаты численною моделирования волновых полги в вя 5ко}и[))1их средах с заданными функциями последействия и долакнея выводы об -эффективное™ исполь«линия предложенною алгоритма.

В дв}х последующих приложениях приводятся некоторые известные теорет-ческие постулаты и формулы, исполы^емые для построения предлаыемых в дис-сергации алюритмов. В приложении 1 описываются некоторые свойпва и приводятся основные формулы для шшчралыюю преобразовании Лаюрра. Расем.и-риваюн'я Iрафики функций Лаюрра в мвисимости от разных параметров и их влияние на спектр коэффициентов разложения но полиномам Лаюрра мданных ф)нкций. В приложении 2 нокаи.1вает(я принцип построения и вывода основных формул к'ории вя?ко>ир>юс1и для шконон иоиющения и дисперсии ( заданной ф> нкцией последействия.

3.5 Основные выводы.

В длиной главе представлен спектральный меюд Лагерра для вязкоупруюю иро-с1ранс1венною моделирования для 2.50 неоднородных сред. Испольювание подхода комплексирования конечною интегральною преобразования Фурье и конечно-разностного метода по пространственным координатам позволяем свести решение исходной задачи к решению независимых систем линейных алгебраических уравнений ,ия каждой гармоники Фурье в отдельности. Такой подход к решению позволяет довольно просю и очень эффективно организовать процесс распараллеливания требуемых вычислений на многопроцессорных вычислительных комплексах. Так как, в этом случае, мы имеем возможность решатI. каждую, полученную таким образом, независимую задачу на отдельно выделенном процессоре. Особенно эффективно применение данного алгоритма вычислений на многопроцессорных вычислительных комплексах кластерного тина, например на подобных юму на котором он и был реализован. Но г ому что проводимые вычисления не требуют частою и большою объема обмена данных между процессорами.

Эффективность использования преобразования Лагерра для задач вязкоунру-госги у же отмечалась в преды,гущей главе. Можно добавить, что такой подход может быть также использован и при решении задач для трёхмерно-неоднородных сред При этом решение получаемой после преобразования Лагерра дифференциальной задачи по пространственным координатам может быть осуществлено, как и в описанной в этой главе способом комплексирования аналитических и конечно-разностным методом, так и с использованием трехмерных пространс Iвенных сеток ,ия разностной аппроксимации [44], [85]. В этом случае, процесс распараллеливания решения производится на этапе нахождения решения системы алгебраических уравнение. Используемый в предлагаемом алгоритме метод сопряженных градиентов позволяет это лет ко сделать Такой подход распараллеливания был опробован на данном ллюргпме при моделировании волновых полей для больших моделей сред. Так как при решение таких задач требуется большой обьем оперативной нами I и. А при данном типе распараллеливания, на каждом процессоре прои зводигея построение вектора решения только для ограниченною пространственною участка среды и следовательно только ,тля определённою участка матрицы решаемой системы. При этом требуется обмен между процессорами найденными значениями вектора решения только в двух соседних пространственных узлах разносIной сетки, так как от параметра Лагерра матрица решаемой системы не зависит, а правая часть уравнений насчитывается независимо для решаемой части системы на данном процессоре. Поэтому спектральный подход решения но времени более экономичен относительно объемов обмена при реал и зации распараллеливания, в отличии от методов использующих численное решения во временной области.

ЗАКЛЮЧЕНИЕ. Основные результаты.

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

На защиту выносятся следующие основные результаты:

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

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

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

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

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

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

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

Исследованы методы повышения ючности иред.Ю/К(Ч1ш»1х алюршмов при расчемч волновых нолей для моделей сред с тонкими слоями и сред с ре$-кокош растными 1раницами. Предложены способы опт ими шиш алюршмон с целыо уменьшения вычислительных трат.

1. Лки К., Ричарде II.

2. Количественная юйсмолошя. Т. 1. Москва, Мир, 1983.

3. Алексеев A.C., Гельчинский 13.51.

4. О лу чевом методе вычисления нолей волн в случае неоднородных сред с криволинейными 1рашщами раздела. Сб. Вопросы динамической теории раепро-( гранения сейсмических волн. Л., 1959, .Vo 3, С. 16 17.

5. Алексеев A.C., Бабич В.М., Гельчинский Б.Я.

6. Лучевой метод вычисления ишенеивности волновых фронтов. Сб. Вопросы динамической теории распространения сейсмических волн. Л., 19G1, .V 5, С. 3-21.

7. Алексеев A.C., Михайленко Б.Г.

8. Решение задач Лэмба для вертикально-неоднородною полупространства. IIш. АН СССР, сер. Фишка Земли, 197G, X« 12, С. 11 25.5| Алексеев A.C., Михайленко Б.Г.

9. Численное моделирование распространения сейсмических волн в радиально-неоднородной модели Земли. Докл. АН СССР, 1977, т. 235, .V 1, С. IG 19.

10. G. Алексеев A.C., Михайленко Б.Г.

11. Нелучевые"-эффекты в юории распространения сейсмических волн. Докл. АН СССР, 1982, т. 2G7, .V 5, С. 1079 1083.7. Бабич В.М., Алексеев A.C.

12. Лучевой метод вычисления ишенсивности волновых фронюв. Изв. АН СССР, Сер. Геофизическая, .V 1, 1958, С. 9 15.

13. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. Москва, Наука, 1987.

14. Горский H.М., Михайленко Б.Г.

15. Решение пространственных задач теории распространения сейсмических волн в скоростях и напряжениях Новосибирск. Препринт ВЦ СО АН СССР, вып 754, 1987, 42 с.10. Грей Э., Мэтбюз Г. Б.

16. Функции Бесим я и их приложения в физике и механике. Москва, II зд иностр. лит., 195311. Коган С.Я.

17. Краткий обзор теории поглощения сейсмических волн. Изв. АН СССР. сер. Фишка Земли. Л* 11, 1966, С. 1 28

18. Конюх Г.В., Михайлеико Б.Г

19. Применение интегральною преобразования Латерра при решении динамических задач сейсмики. Труды ИВМиМГ СО РАН, выи. 5, серия "Математическое моделирование в геофизике", Новосибирск, ИВМиМГ СО РАИ, 1998, С 106 123.

20. Конюх Г.В., Михайлеико Б.Г.

21. Численно-аналитический алюригм ^гя решения прямых динамических задач сейсмики. Труды 6-ой Российской школы-семинара по современным проблемам математического моделирования, Ростовский уииверсиюг, Ростов-на-Дону, 1997, С. 60 70.

22. Конюх Г.В., Михайлеико Б.Г, Михайлов А.А.

23. Численное моделирование сейсмических полей в низкоупру г их средах на основе спектрального метода. Математическое моделирование, Москва, т. 13, .V 2, 2001. С. 01-70

24. Конюх Г. В., Михайлеико Б.Г, Михайлов А. А.

25. Моделирование нестационарных сейсмических полей в неоднородных вязко-упруг их средах на основе преобразования Лагерра. Науки о Земле. Современные проблемы сейсмологии. Москва, Вузовская книга, 2001, С. 25 10

26. Конюх Г. В., Михайлов А. А.

27. Об одном алгоритме декомпозиции областей с использованием конечных ин-нчральных преобразований Фурье-Бесселя. Труды ИВМиМГ СО РАН, вып. 5, серия "Математическое моделирование в геофизике", Новосибирск, ИВМиМГ СО РАН, 1998, С. 91 105.17. Кристенсен Р.

28. Введение в теорию вязкоунругости. Москва, Мир, 197-1.

29. Мартынов В.Н , Михайленко Б.Г.

30. Численное моделирование распространения у пру I их волн в анизотропных неоднородных средах (случай полупространства и сферы). Сб. Математические методы интерпретации юофи шческих наблюдений. Новосибирск, 1979, С. 85 113.19. Михайленко Б.Г.

31. Численное решение задачи Лэмба для неоднородною иолу пространства. Математические проблемы геофизики. Новосибирск, 1973, С. 273 297.20. Михайленко Б.Г.

32. Расчет теоретических сейсмограмм для мноюмерно-неоднородных моделей сред. Сб. Условно-корректные задачи математической физики в интерпретации геофизических наблюдений. Новосибирск, 1978, С. 75 88.21. Михайленко Б.Г.

33. Меюд решения динамических задач сейсмики для двумерно неоднородных моделей сред. Докл. АН СССР, 1979, т. 246, .V 1, С. 47 51.22. Михайленко Б.Г.

34. Сейсмические ноля в сложнопостроенных средах. Новосибирск, Изд. СО РАН, 1988,311 с.21| Михайлов А.А.

35. Моделирование сейсмических полей для 2.5Б неоднородных вязкоу пру г их сред. Труды междун. конференции "Математические методы в геофизике", часть 1. Новосибирск, Изд. ИВМиМГ СО РАН, 2003, С. 116 152.25. Молотков Л.А.

36. Матричный метод и теории распространения волн в слоисгых упругих и жидких средах. Ленинград, Наука, 1984, 270 с.26. Партон В.З., Перлип П.И.

37. Меюды математической теории упругости. Москва, Наука, 1981, 688 с.27| Прудников А.П., Брычков Ю.Л., Млричсв О.И.

38. Интегралы и ряды. Специальные функции. Москва, Наука, 1983.28. Снсмдон И.

39. Преобразования Фурье. Москва, Изд. иностр. лит., 1955.29. Суегин U.K.

40. Классические ортогональные многочлены. Москва, Наука, 1974, 280 с.

41. Фатьянов А.Г., Михайленко Б.Г.

42. Нолу аналитический метод расчета нестационарных волновых нолей л'я слоисю-однородных моделей сред. Сб. Математические методы решения прямых и обратных задач геофи шк. Новосибирск, Изд-во ВЦ СО АН СССР, 1981, С. 92 104.

43. Флгьяиов А.Г., Михайленко Б.Г.

44. Нестационарные сейсмические волновые ноля в неоднородных вя зкоупругих моделях сред Сб Математические проблемы геофизики: модели и численные методы Новосибирск, Изд-во ВЦ СО АН СССР, 1981, С. 82 131.

45. Фа г г.янов А.Г., Михайленко Б.Г.

46. Метод расчета нестационарных волновых полей в неупругих слоисто-неоднородных средах. Докл. АН СССР, 1988, т. 301, .V 4, С. 1021 1027.

47. Alekseev A.S., Mrkhailenko B.C.

48. The solution of dynamic problems of elastic wa\e propagation in inhomogeneoub media by a combination of partial separation of variables and finite-difference method. J. Geophysics, 1980, 48, P. 161 172.

49. Afford R.M., Kelly K.R., Boore D.M.

50. Accuracy of finite-difference modeling of the acoustic wave equation. Geophysics, 1974, 39, P. 834 842.33. Alterman Z., Karal F.G.

51. Propagation of elastic waves layered media by finite-difference methods. Bull. Seism. Soc. Am., 1968, 58, P. 367 398.36 Bay liss A. et.al.

52. A fourth order accurate finite-difference scheme for the computation of elastic waves. Bull., Seis. Soc. Am., 1986, 76, P. 1115 1132.37. Boore D.M.

53. Finite-difference methods for seismic wave propagation in heterogeneous materials. Methods in Computation Physics, 1972, 11, P. 1 37.38. Carcione J.M., Wang Л.Р.

54. A Chebyshev collocation method for the elastodynarnic equation in generalized coordinates. Cornput. Fluid Dynamics J., 1993, 2, P. 269 290.39| Carcione J.M., Kosloff D. and Kosloff R.

55. Wave propagation simulation in a linear viscoelastic medium. Geophys. Л. R. Astr. Soc., 1988, 95, P. 597 -Gil.40. Сепеиу V., Ravindra R

56. Theory of seismic head waves. Toronto Univ. Press, Toronto, 1971, 250 p

57. Cerveny V., Molotkov I.A., Psencik I.

58. Ray method in seismology Prague, Yarlovar. Univ., 1977, 281 p.12. Dablain M.A.

59. The application of high-order differencing to the scale wave equation. Geophysics, 198G, 51, P. 51 -5G.43. Day S.M, Minster J.B.

60. Numerical simulation of attenuated vvavefields using a Pade approximant method. Geophjs. J. II. Astr. Soc., 1984, 78, P. 105 11811. Dong Z., McMechan G.A.

61. Wa\e propagation in three-dimensional spherical sections by the Chebyshev spectral method. Geophys. J. Int., 1999, 136, P. 559 566.19. Iloldberg 0.

62. I Iron F., Mikhailenko B.G.

63. Numerical modeling of nongeometrical effects by the Alekseev-Mikhailenko method. Bull. Soc. Am., 1981, vol. 71, 4, P. 10011 1099.

64. Jo C.I I., Shin C.S., Suh J.II.

65. An optimal 9-point, finite-difference, frequency-space, 2D scalar wa\e extrapolator. Geophysics, 1996, 61, P. 529 537.53. Karal F.C., Keller J.B.

66. Elastic wave propagation in homogeneous and inhomogeneous media. J. Acoust. Soc. Am., 1959, 31, P. 691 705.51. Kelly K.R. et. all.

67. Synthetic seismograms. a finite-difference approach. Geophysics, 1976, 11, P. 2 -27

68. Konyukh G.V., Krivtsov Y.V., Mikhailenko B.G.

69. Numerical-analytical algorithm of seismic wave propagation in inhomogeneous media. Appl. Math. Lett., 1998, 1, vol. 11., P. 23 29.

70. Konyukh G.V., Mikhailenko B.G.

71. Forward seismic modeling based oil combination of finite Fourier transforms with matrix decomposition method. Bulletin of the Novosibirsk Computing Center, series Mathematical Modeling in Geophysics, Novosibirsk, 1998, 4, P. 93 102.

72. Konyukh G.V., Mikhailenko B.G., Mikhailov A. A.

73. On an algorithm of domain decomposition based on finite Fourier-Bessel transforms. Bulletin of the Novosibirsk Computing Center, series Mathematical Modeling in Geophysics, Novosibirsk, 1998, 4, P. 103 113.

74. Konyukh G.V., Mikhailenko B.G., Mikhailov A. A.1.tegral Laguerre transform as applied to forward seismic modeling. Bulletin of the Novosibirsk Computing Center, series Mathematical Modeling in Geophysics, Novosibirsk, 1999, 5, P. 71 91.

75. Konyukh G.V., Mikhailenko B.G., Mikhailov A.A.

76. Application of the integral Laguerre transforms for forward seismic modeling Journal of Computational Acoustics, 2001, vol. 9, .V 4, P. 1523 1511.1. GO. Koslofr D., Baysal E.

77. Forward modeling by a Fourier method. Geophysics, 1982, 47, P. 1402 1412.61. Kosloff D , et. al.

78. Solution of equations of dynamic elasity by a Chebyshev spectral method. Geophysics, 1990, 55, P. 731 718.

79. Kosloff D., Reshef M., Loevventhal D.

80. Clastic wave calculations by the Fourier method. Bull. Seis Soc. Am., 1984, 74, P. 875 891.03. Levander A.R.

81. Fourth order velocity-stress finite-difference scheme. Proc. 57"' SFG Annual Meeting New Orleans, 1987, P. 234 245.

82. Gl. Liu H.P, Anderson D.L., Kanamori H.

83. Velocity dispersion due to anelasticity; implications for seismology and mantle composition. Geophys. J. R. Astr. Soc., 47, 197G, P 41 58.63. Lysmor B , Drake N.

84. A finite-element method for seismology. In Bolt B.A. lids., Methods in computational physics, Seismology: Surface waves and Earth oscillations, Academic Press Inr , 1972, 11, P 181 216.66. Marfurt K.J.

85. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave (»(illations. Geophysics, 1981, 19, P. 533 549.67. Marfurt K.J., Shin C.S.

86. The future of iterative modeling in geophysical exploration In Eisner E. Ed., Handbook of geophysical exploration: Seismic exploration, supercomputers in seismic exploration, Pergainon Press, 1989, 21, P. 203 228

87. Martynov V N., Mikhailenko B.G.

88. Numerical modeling of propagation of elastic waves in anisotropic inhomogeneous media for the half-space and the sphere. Geophys. J. II. Astr. Soe., 1981, 76, P 53-63

89. Martynov V.N., Mikhailenko B.G.

90. Two algorithms for calculation of theoretical seismograms for anisotropic media Bulletin of ICNNG, Novosibirsk, 1999, 5, P. 105 115.

91. Mikhailenko B.G., Korneev V.I.

92. Calculation of synthetic seismograms for complex subsurface geometries by a combination of finite integral Fourier transforms and finite-difference techniques. J Geophys., 1981, 51, P. 195 206.71. Mikhailenko B.G.

93. Synthetic seismograms for complex 3D geometries using an analytical-numerical algorithm. Geophys. J. R. Astr. Soc , 1984, 79, vol. 3, P. 963-986.72. Mikhailenko B.G.

94. Numerical experiments in seismic investigations. J. Geophys., 1985, 58, P. 101 -121.73. Mikhailenko B.G.

95. Spectral Laguerre method for the approximate solution of time dependent problems. Appl. Math. Lett., 1999, 12, P. 105 110.74. Mikhailenko B.G.

96. Seismic modeling by the spectral-finite difference method. Physics of the Earth and Planetary Interiors, 2000, 119, P. 133 147.

97. Mikhailenko B.G , Mikhailov A.A , Reshetova G.V.

98. Numerical modeling of transient seismic fields in viscoelastic media based on the Laguerre spectral method. Pure apll. geophys., 2003, 1G0, P. 1207-1224.170J Mikhailenko B.G., Mikhailov A.A., Reshetova G.V.

99. Numerical viscoelastic modeling by the spectral Laguerre method. Geophysical Prospecting, 2003, 51, P. 37-18.77. MullerG.

100. The reflectivity method; a tutorial Geophysics, 1985, 58, P. 153 174.78. MullerG.

101. Rheological properties and \elocity dispersion of a medium with power-law dependence of Q on frequency. Geophysics, 1983, 54, P. 20 29.

102. Robertsson J., Blanch J., Syines W.

103. Viscoelastic finite-difference modeling. Geophysics, 1999, 61, P. 1444 115G80. Stekl I., Pratt R.G.

104. Accurate viscoelastic modeling by frequency-domain finite-difference using rotated operators. Geophysics, 1998, G3, P. 1779 1791.

105. Tal-E/er H., Kosloff D , Koren Z.

106. An accurate scheme for seismic forward modeling. Geophys. Prosp., 1987, 35, P. 479 190.

107. Tal-Ezer H., Carcione J.M, Kosloff D.

108. An accurate and efficient scheme for wave propagation in linear viscoelastic media. Geophysics, 1990, 55, P. 13G6 1379.83. Tarantola A.

109. Theoretical background for the inversion of seismic waveforms, including elasticity and attenuation. Pure and Applied Geophysics, 1988, 128, P 3G5 -399.84. Thomson W.T.

110. Transmission of clastic wa\es through a stratified solid. J. Appl. Pliys., 1950, 21, P. 89 9385. Xu T., McMechan G.A.

111. Efficient 3-D viscoelastic modeling with application to near-surface land seismic data. Geophysics, 1998, 63, P. 601 612.86. Zener C M.

112. Elasticity and Anelasticity of Metals. The University of Chicago Pres, 1918.