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

ИСПОЛЬЗОВАНИЕ HPC ДЛЯ РЕШЕНИЯ ЗАДАЧ ОБРАБОТКИ СЕЙСМИЧЕСКИХ ДАННЫХ ПРИ ПОМОЩИ РАЗЛОЖЕНИЯ ПО ВОЛНОВЫМ ПАКЕТАМ

, Андерссон Фредрик

Лундский университет, г. Лунд, Швеция

Современное развитие сейсмических технологии характеризуется постоянным усложнением систем наблюдений и увеличением объемов полевого материала (переход от 2D к 3D съемке, использование трехкомпонентных сейсмоприемников, наблюдения по широкому азимуту и т. д.), так что его обработка неизменно требует использования самых мощных компьютеров.

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

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

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

Для двумерного случая существует пакет Curvelab [4,5], который используется в большинстве приложений в настоящее время. Для трехмерного разложения требуется параллельная реализация на кластерных системах. В пакете Curvelab есть трехмерное разложение, основанное на использовании MPI. Но использовать его для решения практических задач очень затруднительно (см. обсуждение в [6]).

Для трехмерного разложения недавно были предложены новые базисы волновых пакетов [7, 8], для которых базисные функции обладают лучшими свойствами (симметрией и изометричной формой) для обработки и фильтрации данных по сравнению с пакетом Curvelab. Эти новые типы разложения также допускают реализацию в форме алгоритмов на основе быстрых преобразований Фурье на нерегулярных сетках (USFFT, [9,10]). Вычислительная сложность таких алгоритмов имеет такой же порядок, что стандартное быстрое преобразование Фурье (FFT). 

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

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

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

- прямое преобразование Фурье с регулярной на нерегулярную сетку

- обратное преобразование Фурье с нерегулярной на регулярную сетку

Алгоритм USFFT с регулярной на нерегулярную сетку и состоит из следующих частей:

Деление на гауссиан в области изображения Удвоение изображения, дополнение нулями FFT Задание функции на границах используя условие периодичности Интерполяция и свёртка в частотной области

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

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

таких что

где - координаты в трехмерном пространстве. - базисные функции волновых пакетов.

На рисунке 4 показан пример трёхмерного гауссова волнового пакета.

Рисунок 4. Пример трёхмерного гауссова волнового пакета (слева). Срез через центр пакета по оси x2. Для пакета заданы размеры, угол наклона и количество осцилляций вдоль каждой оси.

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

,

где k, \nu, j – параметры гауссова волнового пакета, ,- обозначает спектр Фурье от .

Спектральные функции волновых пакетов подобраны таким образом, чтобы при возведении в квадрат задавать разбиение единицы Фурье области:

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

Ускорение. Наиболее ресурсо-затратными частями алгоритмов являются преобразования Фурье и процессы интерполяции. Для их ускорения использовались высокопроизводительные средства Intel (MKL) и NVidia (CUDA). Интерполяции для CPU и GPU была ускорена при помощи различных стратегий распараллеливания.

Было проведено тестирование прямого и обратного 3D USFFT на процессорах Intel i7 и NVidia Tesla C2050. За время последовательной программы бралось время выполнения (программно) оптимизированного кода на процессоре Intel i7, откомпилированного при помощи стандартного для Visual Studio компилятора cl.

Рассмотрено 2 случая: количество точек нерегулярной сетки равно количеству точек нерегулярной (Nus = N3) и в два раза меньше (Nus = N3/2). (таблицы 1,2).

Таблица 1. Время выполнения прямого и обратного трёхмерного 3D USFFT на CPU (последовательная программа), на CPU с использованием MKL (1 и 4 потока) и на GPU. Nus = N3

N CPU (посл.) MKL (1 поток) MKL (4 потока) GPU пр. обр. пр. обр. пр. обр. пр. обр. 26 15.27 15.47 2.27 2.96 0.99 1.57 0.55 1.22 27 129.35 136.22 20.68 27.92 8.06 13.50 2.7 12.41 28 1070.92 1092.92 205.72 245.52 65.23 118.53 19.7 67.10

Таблица 2. Время выполнения прямого и обратного трёхмерного 3D USFFT на CPU (последовательная программа), на CPU с использованием MKL (1 и 4 потока) и на GPU. Nus = N3/2

N CPU (посл.) MKL (1 поток) MKL (4 потока) GPU пр. обр. пр. обр. пр. обр. пр. обр. 26 7.61 7.65 1.18 1.63 0.51 0.94 0.46 1.02 27 64.93 68.46 10.82 15.47 4.28 7.68 1.8 9.54 28 531.21 546.74 99.32 132.32 39.32 80.22 12,14 50.86

По таблицам видно, что версии с использованием MKL + openMP (4 потока) и с использованием CUDA демонстрируют существенные ускорения. Сравним время выполнения обратного преобразования для Nus=N3 и Nus=N3/2, N=28.

Seq: , MKL4t : , GPU:

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

Полученные ускоренные версии USFFT были использованы в программе разложения сеймических данных по гауссовым волновым пакетам. Было проведено тестирование на процессоре Intel i7 с использованием MKL и OpenMP (4 потока). Для двумерного случая ускорения достигают 5-10 раз. Наиболее существенные ускорения проявляются в трёхмерном случае разложения. В таблицах 3,4 представлено время выполнения прямого и обратного преобразований для входных данных размерностью 1283.

Из таблиц видно, что для данных больших объёмов больше подходит использование GPU. Однако ввиду больших затрат на копирование данных, а также различных инициализаций составных частей видеокарты версия с входными данными малого размера на CPU работает быстрее.

Таблица 4. Время (в секундах) выполнения программы на различных платформах. N=128

Платформа Прямое преобр-е (сек) Обратное преобр-е (сек) Intel i7 (посл.) 16522.54 16822.21 Intel i7 MKL+OpenMP (4 threads) 963.38 1325.12 Tesla C2050 538.43 960.21

Оптимизированные версии программ разложения данных по гауссовым волновым пакетам по сравнению с последовательным кодом демонстрируют существенные ускорения для обеих платформ: CPU при помощи Intel MKL и OpenMP  ~ 17 раз и GPU с использованием NVidia CUDA ~ 30 раз.

Апробация. Программа была внедрена в пакет обработки сейсмических данных Madagascar. Было проведено успешное тестирование алгоритмов на синтетических сейсмических данных.

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

Рис. 10. Сжатие сейсмических данных и подавление шумов: исходные данные со случайным шумом (слева); результат сжатия в 50 раз (CR=.02) и подавления шумов (справа).

Пусть имеется куб сейсмических данных размера N (количество пикселей), с шумом (рис. 10, слева) или пропущенными трассами (рис. 11, слева); на гранях показаны сечения через куб в местах, показанных тонкими линиями. Процедуру подавления помех и интерполяции данных можно реализовать, применив прямое ПВП, сохранение M коэффициентов, абсолютное значение которых больше заданного уровня, обратное ПВП:

.

Здесь соответствует данным с подавленной помехой (рис. 10, справа) и проведенной интерполяцией (рис. 11, справа); С - оператор прямого ПВП (см. (1)). При M < N происходит также сжатие данных. Степень сжатия будем измерять коэффициентом . В примерах на (рис. 10-11) показано сжатие данных в 50 раз (CR=.02) без заметной потери качества изображения целевых волн.

Рис. 11. Интерполяция сейсмических данных: исходные данные с пропущенными трассами (слева); результат сжатия в 50 раз (CR=.02) и интерполяции данных (справа).