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

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

В настоящее время существует большое количестно различных методов решения прямых динамических задач сейсмики Одним из наиболее популярных, инляег-<�я лучевой метод, провоженный н работах [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их средах с заданными функциями последействия и долакнея выводы обэффективное&tradeисполь"линия предложенною алгоритма.

В дв}х последующих приложениях приводятся некоторые известные теорет-ческие постулаты и формулы, исполы^емые для построения предлаыемых в дис-сергации алюритмов. В приложении 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раницами. Предложены способы опт ими шиш алюршмон с целыо уменьшения вычислительных трат.