Д.А. Маловичко
Пермский университет
Псевдоспектральные методы моделирования сейсмического волнового поля
С начала 80-х гг. для численного решения уравнений динамической теории упругости активно используются псевдоспектральные методы [2]. По своей сути данные методы близки к разностным методам – в области решения вводится некоторая сетка, значения выражений в уравнении принимаются постоянными в пределах шага сетки. Различие между методами заключается в способе аппроксимации пространственных производных. В разностных методах они заменяются алгебраическими комбинациями значений выражений в узлах сетки, в псевдоспектральных - вычисляются с помощью умножения на волновое число в области Фурье-преобразований. Например, вычисление выражения ¶U(x)/ ¶x в псевдоспектральных методах осуществляется по следующей схеме:
   ¶U(x)/ ¶x= FT-1 ( ikx FT{ U(x) } ),
где FT и FT-1 – соответственно прямое и обратное преобразования Фурье, kx – волновое число в направлении x, i – мнимая единица. Использование  алгоритма  быстрого  преобразования Фурье (БПФ) делает псевдоспектральные методы сопоставимыми по скорости вычислений с разностными методами. Ограничения вычислительных мощностей современных компьютеров позволяют использовать псевдоспектральные методы только для решения одно- и двумерных задач. 
Основным недостатком псевдоспектральных методов является невозможность явного задания в них граничных условий. Применение разложения в ряд Фурье для вычисления пространственных производных предполагает автоматическое выполнение условий периодичности на границах сетки. Поэтому, с помощью описанной выше схемы вычисления пространственных производных некорректно решать уравнения теории упругости в цилиндрических координатах (z, r, q). На оси r=0 становится невозможно задать источник, т.к. вычисление производных по координате r подразумевает выполнение условия ¶U/ ¶r½r=0 = ¶U/ ¶r½r=R , где R – граница сетки в направлении r. Учитывая это, приведенная выше схема вычислений производных позволяет решать уравнения динамической теории упругости только в плоской постановке. 
Эффективный способ устранения условий периодичности приведен в [3]. Используемая в этой работе методика решения уравнений динамической теории упругости для полуплоскости позволяет приближенно описать граничное условие свободной от напряжений поверхности. Это достигается применением разложения в ряд по полиномам Чебышева вместо разложения в ряд Фурье при вычислении производных по вертикальным координатам. Использование аппроксимации полиномами Чебышева приводит к неявному заданию непериодических граничных условий на горизонтальных границах сетки.
Решения уравнений динамической теории упругости на плоскости не в полной мере характеризуют процесс распространения сейсмических волн в реальных средах, поскольку не описывают расхождения энергии в азимутальных направлениях. Данный недостаток устраняется, если уравнения упругости решать в осесимметричной постановке, а при вычислении производных по координате r вместо разложения в ряд Фурье использовать разложение в ряд по полиномам Чебышева. Получаемый таким образом псевдоспектральный метод решения уравнений динамической теории упругости для случая осевой симметрии программно реализован автором. Работоспособность метода демонстрируют результаты тестовых расчетов. 
В однородном упругом полупространстве (скорости продольных волн 4.6 км/с, поперечных – 2.7 км/с, плотность среды 2 кг/м3) был задан вертикальный импульсный источник на глубине 10 км. Длительность синусоидального сигнала в источнике составляла 2 с. На рис. 1 изображена дивергентная часть сейсмического волнового поля (колебания связанные с изменением объема) в различные моменты времени после излучения сигнала. Метками P и PP обозначены прямая и отраженная продольные волны, R – волна Рэлея и, наконец, Ph – головная волна и просачивающаяся мода, образующаяся в результате падения на свободную поверхность поперечной волны SV. Сейсмограмма вертикальных компонент скоростей смещений для приемников, расположенных на поверхности с шагом 10 км, приведена на рис.2. 
Полученная сейсмограмма правильно описывает изложенные в [1] основные закономерности формирования сейсмического волнового поля для тестовой модели среды.
Библиографический список
1. Аки К., Ричардс П. Количественная сейсмология: теория и методы. Т.1. М.: Мир, 1983. 520 с.
2. Kosloff D., Baysall E. Forward modeling by a Fourier method // Geophysics. 1982. 47.  P. 1402-1412.
3. Kosloff D., Kessler D., Filho A. Solution of the equations of dynamic elasticity by a Chebychev spectral method // Geophysics. 1990. 55. P. 734-748. 

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