The use of parallel computing for the high-resolution determination of earthquake source parameters

Full Text

ВВЕДЕНИЕ

Исследование особенностей излучения длиннопериодных поверхностных волн мелкофокусным землетрясением, проведенное в работах [Букчин, 2006; Bukchin et al., 2010], показало, что в случае, когда одна из нодальных плоскостей источника субгоризонтальна, небольшие изменения ее угла падения существенно меняют диаграмму поверхностно-волнового излучения такого очага и оценку моментной магнитуды события. Пример такого изменения диаграммы излучения фундаментальной моды волны Лява на периоде 200 с в зависимости от величины угла падения d приведен на рис. 1. Источники представляют собой сдвиги по простиранию с различными значениями угла падения d и фиксированными значениями угла простирания (y = 0°), угла подвижки (l = 0°) и глубины источника (h = 30 км). В верхнем ряду указаны текущие значения угла падения, в среднем ряду приведены соответствующие фокальные механизмы, а в нижнем ряду – соответствующие диаграммы излучения волны Лява. Как видно из рисунка, при изменении угла падения от 90° до 15° диаграмма излучения волны Лява почти не меняется, и лишь на последних десяти градусах, когда одна из нодальных плоскостей становится субгоризонтальной, небольшие изменения фокального механизма существенно меняют диаграмму поверхностно-волнового излучения. Как следствие, применяемые к таким событиям методы оценки углов фокального механизма должны обеспечивать достаточно высокий уровень их разрешения.

 

Рис. 1. Диаграммы излучения фундаментталльной моды волны Лява мелкофокусным источником, представляющим собой сдвиг по простиранию, для различных значений угла падения $\delta$.

 

ОПРЕДЕЛЕНИЕ ФОКАЛЬНОГО МЕХАНИЗМА ЗЕМЛЕТРЯСЕНИЯ, ЕГО ГЛУБИНЫ И СЕЙСМИЧЕСКОГО МОМЕНТА ПО ДЛИННОПЕРИОДНЫМ АМПЛИТУДНЫМ СПЕКТРАМ ПОВЕРХНОСТНЫХ ВОЛН

Наиболее интенсивные смещения земной поверхности, регистрируемые сейсмическими станциями, наблюдаются в поверхностных волнах. Каждая из таких волн является результатом свертки соответствующей функции Грина, определяемой строением среды вдоль пути распространения волны, и тензора плотности сейсмического момента, характеризующего неупругие процессы, происходящие в очаге землетрясения [Backus, Mulcahy, 1976]. При расчете функции Грина мы используем модель среды со слабой горизонтальной неоднородностью [Woodhouse, 1974; Бабич и др., 1976]. Поверхностно-волновая функция Грина для такой модели зависит от параметров среды в окрестности очага и в окрестности точки регистрации, от средней вдоль луча фазовой скорости волны, от ее геометрического расхождения и затухания. При этом амплитудный спектр смещений в поверхностной волне не зависит от ее фазовой скорости. Скорости поверхностных волн в реальной Земле не известны с достаточной точностью. Поэтому, как правило, мы используем для определения параметров источника лишь амплитудные спектры поверхностных волн. Детальное описание метода и примеры применения можно найти в работах [Букчин и др., 1992; Lasserre et al., 2001; Букчин и др., 2015].

Описывая источник в приближении тензора момента, мы рассматриваем мгновенную точечную сдвиговую дислокацию (двойной диполь) на глубине h. Такой источник задается пятью параметрами: его глубиной, фокальным механизмом, определяемым тремя углами (простирания y, падения d и подвижки l) и сейсмическим моментом М0. Четыре первых параметра мы определяем прямым перебором их возможных значений на сетке в параметрическом пространстве, а пятый параметр М0 – минимизируя отличия (невязку e) наблюденных амплитудных спектров от их теоретических значений для каждой текущей комбинации значений остальных параметров. Детальность сетки для углов фокального механизма может достигать одного градуса. Значения параметров, минимизирующие невязку, мы рассматриваем как оценки этих параметров. Для оценки степени разрешения каждого из этих четырех параметров мы строим четыре частные функции невязки: eh (h), ey (y), ed (d) и el (l). Осуществляя перебор возможных значений параметров, мы рассматриваем лишь ту из двух нодальных плоскостей, которая падает круче. Ее угол падения не может быть меньше 45°, и его область значений определяется неравенствами 45° ≤ d ≤ 90°.

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

ПРИМЕНЕНИЕ ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛЕНИЙ ДЛЯ УСКОРЕНИЯ РАСЧЕТА ПАРАМЕТРОВ ОЧАГА ЗЕМЛЕТРЯСЕНИЯ

По описанному выше алгоритму была написана программа с распараллеливанием вычислений (данную программу MomTens_omp можно скачать с сайта ИТПЗ РАН http: //mitp.ru/en/soft). Одним из наиболее популярных средств параллельного программирования для компьютеров с общей памятью в настоящее время является технология OpenMP. Данная технология базируется на традиционных языках программирования, таких как Си, Си++ и Фортран. За основу берется последовательная программа, а для создания ее параллельной версии используется набор директив, функций и переменных окружения. Предполагается, что создаваемая параллельная программа будет переносимой между различными компьютерами с разделяемой памятью, поддерживающими OpenMP. Технология OpenMP нацелена на то, чтобы пользователь имел один вариант программы для параллельного и последовательного выполнения [Антонов, 2009]. Здесь мы рассматриваем возможность применения данной технологии для ускорения процесса определения параметров очага землетрясения по спектрам поверхностных волн. В общем виде алгоритм работы программы представлен на блок-схеме на рис. 2.

 

Рис. 2. Блок-схема работы программы.

 

Как видно из блок-схемы самым внешним циклом является цикл по глубине. Именно по данному циклу проводится распараллеливание. Как описывалось выше, технология OpenMP представляет собой распараллеливание приложений на многопроцессорных системах с общей памятью, в нашем случае – на ядрах одного процессора. Программа начинается с последовательной области – сначала работает один процесс (нить), при входе в параллельную область порождается еще некоторое число процессов, между которыми в дальнейшем распределяются части программы. Итерации цикла по глубине распределяются на максимальное количество ядер и выполняются независимо на каждом из них за исключением критических секций. Под критической секцией подразумевается такая часть программы, которая может исполняться не более, чем одной нитью. Если критическая секция уже выполняется какой-либо нитью, то все другие нити, выполнившие директиву для секции с данным именем, будут заблокированы, пока вошедшая нить не закончит выполнение данной критической секции. Как только работавшая нить выйдет из критической секции, одна из заблокированных на входе нитей войдет в нее. Если на входе в критическую секцию стояло несколько нитей, то случайным образом выбирается одна из них, а остальные заблокированные нити продолжают ожидание [Антонов, 2009]. Мы используем критические секции для работы с глобальными (общими для всех ядер) минимумами невязок. Это необходимо для предотвращения одновременной перезаписи глобальных минимумов несколькими нитями, когда несколько текущих локальных минимумов меньше текущего глобального минимума.

СИЛЬНЕЙШИЙ АФТЕРШОК ЯПОНСКОГО ЗЕМЛЕТРЯСЕНИЯ В ТОХОКУ, 11.03.2011 г., MW = 8.4.

Описанный подход был использован для изучения сильнейшего афтершока землетрясения в Тохоку, произошедшего через полчаса после основного толчка. В результате его записи оказались сильно зашумленными излучением предшествующего главного события. С помощью программы спектрально-временного анализа СВАН мы отфильтровали в полосе периодов от 100 до 200 с и использовали для определения параметров изучаемого события фундаментальные моды Лява и Рэлея, зарегистрированные 13 станциями сетей IRIS, GEOSCOPE и GEOFON. Расположение этих станций приведено на рис. 3а. Шаг сетки для углов фокального механизма был выбран равным 1°. В таблице представлены фокальные механизмы и фазовые невязки для четырех эквивалентных решений, полученных из анализа амплитудных спектров поверхностных волн. Как видно из таблицы, оптимальным является решение 1, приведенное на рис. 3б. Полученные значения углов простирания, падения и подвижки: 33°, 89° и 91°, соответственно. Полученное значение сейсмического момента равно 0.46 × 1022 Нм. Это значение соответствует значению моментной магнитуды Mw = 8.4. Наша оценка глубины наилучшего точечного источника равна 10 км. Частные функ­ции невязки для глубины источника и для углов фокального механизма, характеризующие их разрешение, приведены на рис. 4.

 

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

Номер решения	Угол падения, град	Угол простирания, град	Угол подвижки, град	Нормированная фазовая невязка
1	89	33	91	0.345
2	89	213	91	0.483
3	89	213	–89	0.621
4	89	33	–89	0.717

 

Рис. 3. Распределение точек регистрации поверхностных волн (а) и наилучшее из четырех эквивалентных решений (б), полученных по амплитудным спектрам поверхностных волн. Темные треугольники соответствуют волнам Рэлея, светлые треугольники – волнам Лява.

 

Рис. 4. Частные функции невязки для глубины источника и углов фокального механизма.

 

Обе версии программы (без распараллеливания и с распараллеливанием) были отхронометрированы на использованных записях сильнейшего афтершока Тохоку. Время работы программы без распараллеливания оказалось равным 99 минутам, с распараллеливанием на 2 ядра (каждое ядро – 2 потока) – 41‑й минуте. Как видим, время работы программы уменьшилось больше чем в 2 раза (количество ядер), но меньше чем в 4 раза (количество потоков). Это связано с несколькими факторами. Во‑первых, в программе помимо параллельной части есть еще последовательные, связанные со считыванием данных, критическими секциями, записью и выводом результатов. Во‑вторых, стоит отметить, что потоки не эквивалент­ны ядрам. В отличие от «настоящих» ядер, являющихся полными и независимыми копиями, в случае многопоточности в одном процессоре дублируется лишь часть внутренних узлов, в первую очередь отвечающих за хранение и подготовку данных. Исполнительные же узлы, ответственные за организацию и обработку данных, остаются в единственном числе, и в любой момент времени используются максимум одним из потоков.

СРАВНЕНИЕ РЕЗУЛЬТАТОВ С СМТ-РЕШЕНИЕМ

Сравним полученное решение с решением из Глобального СМТ-каталога (решение ГСМТ). Но вначале рассмотрим проблему выбора спектрального диапазона наблюдений. При перемещении этого диапазона из области периодов, недостаточно длинных для использования приближения точечного мгновенного источника, в область более длинных периодов оценки параметров очага существенно изменяются. Когда же периоды становятся достаточно длинными и указанное приближение – адекватным, оценки параметров принимают свои предельные значения и перестают изменяться с дальнейшим ростом периодов. Соответствующая иллюстрация для сильнейшего афтершока землетрясения в Тохоку приведена на рис. 5а, 5б, 5в. Здесь изображены решения, полученные из одного и того же набора записей поверхностных волн, отфильтрованных в трех разных спектральных диапазонах.

Сравним теперь решение, приведенное на рис. 5в, с ГСМТ-решением, представленном на рис. 5г. Эти решения существенно различаются. Так полученная нами оценка Mw равна 8.4, в то время как ГСМТ-оценка Mw равна 7.9. Заметим, однако, что сравниваемые решения получены по записям поверхностных волн, отфильтрованных в разных спектральных диапазонах. Если же мы сравним ГСМТ-решение (Т > 75 c) с решением, приведенном на рисунке 5а (75 c < Т < 120 c), то увидим, что эти решения весьма близки. Но как было отмечено выше, диапазон периодов 75 c < Т < 120 c является недостаточно длиннопериодным для рассматриваемого события. Это дает нам основание предполагать, что существенные отличия нашего решения от ГСМТ-решения объясняются неадекватным выбором спектрального диапазона наблюдений при расчете последнего.

 

Рис. 5. Результаты инверсии спектров поверхностных волн в трех спектральных диапазонах. Сравнение с ГСМТ-решением.

 

ВЫВОДЫ

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

Авторы благодарят глобальные сейсмические сети IRIS, GEOSCOPE и GEOFON за предоставление сейсмических записей.

About the authors

A. S. Fomochkina

Institute of Earthquake Prediction Theory and Mathematical Geophysics

Author for correspondence.
Email: nastja_f@bk.ru
Russian Federation, 84/32, Profsoyuznaya street, Moscow, 117997

V. G. Bukchin

Institute of Earthquake Prediction Theory and Mathematical Geophysics

Email: nastja_f@bk.ru
Russian Federation, 84/32, Profsoyuznaya street, Moscow, 117997

References

Supplementary files