Декомпозиция и распараллеливание вычислительноёмких функций обработки данных микросейсмического мониторинга Текст научной статьи по специальности «Математика»

АПВПМ-2019

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

К, С, Алеынбаев1

1 НИИ прикладной информатики и математической геофизики БФУ имени И.Канта, 236001, Калининград

УДК 519.683:550.8

Б01: 10.24411/9999-016А-2019-10004

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

Ключевые слова: микросейсмический мониторинг, кластеризация, тетраэдризация, воксель, параллельные вычисления, распознавание.

Введение

Программное обеспечение обработки данных микросейсмического мониторинга можно условно разделить на две части. Первая представляет собой обработку исходных оцифрованных сигналов от поверхностных датчиков. Обычно посредством решения обратной кинематической задачи вычисляются локации сейсмических событий и преобразуются в представление набора (облака) трёх- или четырёх-мерных параметризованных точек, соответствующих событиям (далее будем называть их точками-событиями). Среди параметров — время проявления сигнала сейсмической эмиссии, энергия, скорость распространения сейсмических волн в среде и др. Во второй части решаются задачи выделения полезного сигнала посредством пространственно-временной фильтрации, выявления характерных особенностей каждой точки-события и их совокупностей для интерпретации происходящих в околоскважинной области процессов, создаётся трёхмерная динамическая визуализация полученных результатов. Одна из главных целей — выявление или подчёркивание особенностей процессов развития трещиноватости. В коллективе автора статьи созданы средства визуального подчёркивания тех или иных особенностей точек-событий для экспертной оценки картины данных пользователем [1]. Разрабатываются экспериментальные средства выявления уплотнений и образований в облаке точек-событий с использованием средств геометрического моделирования.

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

Есть алгоритмы, объём вычислений которых растёт приблизительно линейно от объёма обрабатываемых данных, и вычисления локализованы для анализируемых точек, например, в задаче вычисления плотности точек (чем больше в окрестности точки соседей, тем выше плотность). Очевидно, такие алгоритмы легко распараллеливаются. Задача сильно усложняется, когда в алгоритме рассматриваемая точка может быть

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 16-2915090).

!ЯВ.\ 978-5-901548-42-4

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

Распараллеливание, даже если оно возможно, не решает кардинально задачи ускорения алгоритмов с временной сложностью больше 0(N2 log N). Одним из подходов может быть изменение алгоритма таким образом, чтобы разделить обрабатываемые данные на фрагменты с ограниченным количеством элементов в каждом и организовать их раздельную обработку, чтобы время счета алгоритма для любого фрагмента не выросло до неприемлемых значений. При таком способе разделения обработки данных необходимо решить задачу последующего объединения локальных результатов, что может представлять непростую проблему.

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

1 Модификация алгоритма кластеризации

1.1 Построение дерева для точек-событий

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

На практике могут использоваться бинарные деревья (например, BSP-деревья [2]), или, в зависимости от размерности пространства, октодеревья и 16-деревья. В наших реализациях выбран второй подход, поскольку теория бинарных деревьев развита для оптимизации последующего поиска, а в нашем случае поиск отсутствует, основные временные затраты при построении дерева и обратной операции, важнее оптимизация для сложных операций, ориентированных по координатам.

Используется рекурсивный способ построения дочерних узлов. Теоретически временная сложность такого алгоритма в вырожденных случаях, когда большинство точек сконцентрированы в небольшой части рассматриваемого пространства, может достигать 0(N2). Но облако точек-событий микросейсмических данных содержит в целом равномерно распределённые точки, «сгущения» представляют умеренно плотные образования. Построение дерева вглубь ограничиваются минимально допустимым размером ребра области, отвечающей вершине, поэтому количество V вершин дерева и глубина дерева Н ограничены. Поскольку известна оценка временной сложности рекурсивного построения октодеревьев, равная 0(NH), то практическая временная сложность может быть оценена между 0(N) и 0(N log N).

1.2 Декомпозиция алгоритма кластеризации

Кластеризация точек-событий в трёх- или четырёх-мерном пространстве необходима для некоторых задач пространственно-временной фильтрации и для выявления «сгущений» событий перед геометрическим моделированием выявляемых образований. Задача кластеризации ставится при изначально неизвестном количестве кластеров, критерии принадлежности — допустимое максимальное расстояние D между соседними точками кластера, расстояние может быть эвклидовым, или же неэвклидовым с построением неизотропной метрики на основе тензоров сейсмического момента [1]). В качестве основы выбран алгоритм кластеризация точек на регулярной сети [3], легко модифицирующийся для случаев трех- и четырех- измерений и для различных видов используемой метрики. Алгоритм позволяет выделять связные компоненты типа сгущений или лент, которые наиболее пригодны для распознавания разломов или зон трещиноватости. Однопроцессорная реализация без фрагментации данных имеет временную сложность от 0(N2 log N) до 0(N3) [4], где N — количество точек-событий в облаке.

Модификация алгоритма заключается в следующих этапах:

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

превышать заданное значение Z или очередной шаг разбиения не приведёт к заданному предельно

малому размеру области.

Рис. 1: А. Укрупнённый фрагмент тестового множества, попадающий на несколвко соседних листьев. Б. Пример порожденных вершин октодерева для теста. Имеются 2 ограничивающие точки (на нижней плоскости ближе к наблюдателю и вдали от наблюдателя) для формирования требуемых размеров октодерева

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

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

3. Проводится иерархическая склейка элементов разбиения вместе с полученными в них кластерами.

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

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

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

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

Есть 2 типа шагов: 1) поиск в дереве всех вершин, соседствующих с пустыми областями, однотипными по форме и размеру и укрупнение дерева путём операций слияния всех пар соседних областей. На рис. 3 приводится пример процесса такого укрупнения некоторой области дерева до приведения к соседству только вершин. 2) поиск и операции слияния всех пар однотипных по размеру вершин с одновременной склейкой соседствующих по границам кластеров. Результатом слияния пар элементов дерева является прямоугольный параллелепипед, необязательно куб. Шаги типов 1) и 2) попеременно чередуются, пока не останется един-

Рис. 3: Укрупнение элементов дерева в случае соседства вершин и пустых узлов

ствеппая вершина дерева, которая, собственно, и содержит все кластеры, включая склеенные. Внутренние циклы шагов как типа 1), так и типа 2) могут распараллеливаться.

Будем обозначать множество точек большой латинской буквой, например Р, а точку в этом множестве как Рг, кластер, к которому принадлежит точка и номер этого кластера, как к(Рг). При построении дерева возникают множества точек, содержащиеся в прямоугольной области пространства, соответствующее вершине при построении и модификации дерева. Обозначим как Р множество точек, соответствующее первой из вершин, участвующих в склейке, и как Я — множество точек, соответствующее второй вершине.

Алгоритм склейки кластеров соседних вершин шага типа 2) осуществляется следующим образом.

• Обеспечивается уникальность номеров кластеров в объединяемом множестве кластеров путём увеличения номеров кластеров в Я на значение количества кластеров в Р.

• Выбирается в отдельный массив С вершины из Д, которые отстоят от границы с областью Р не более чем на Б (рис. 4 А).

• Об л ас ти Я и Р сливаются в одну область в с объединением списочных структур.

• Для каждой точки Рг области Р производится просмотр отобранных граничных точек из С. Если встретится точка такая, что расстояние от Рг до С? меньше Д список кластера к((3]) подсоединяется к списку кластера к(Рг).

Легко видеть, что после просмотра всех точек из Р все кластеры в, расстояние между какими либо точками которых меньше И, будут объединены (рис. 4 Б, В). На рис. 4 А и Б приведены иллюстрации результатов склейки для теста большого размера, па рис. 5 В — для реального облака.

Сложность по времени алгоритма необходимо оценивать по-отделыюсти для трёх этапов.

Выше мы оцепили сложность построения дерева для облака точек-событий микросейсмических данных между 0(М) и 0(М N).

Кластеризация точек одной вершины немодифицированным алгоритмом оценивается от 0(М2 1og М) до 0(М3) [4], где М — число точек в области одной вершины, а время последовательной кластеризации

Рис. 4: А. Склейка кластеров соседствующих вершин: зелёным цветом окрашено множество приграничных точек, которое просматривается для точек области левой вершины. Б. Кластеры, прошедшие кластеризацию раздельно для вершин (цветом выделены результаты локальной кластеризации). В. Кластеры исходного облака, полученные в результате нескольких шагов иерархической склейки

Рис. 5: А. Большой тест (облако в 770 ООО точек, 426 вершин) после фрагментарной кластеризации. Б. Облако после иерархической склейки. В. Реальное облако (452 тысячи точек), кластеризованное модифицированным алгоритмом.

точек всех вершин примерно линейно зависит от числа вершин. Поскольку M ограничена константой Z, то сложность этапа кластеризации можно оценить в 0(N ), тел и M всегда меньше Z. Однако это условие достигается не всегда, так как оно может конфликтовать с условием минимального размера ребра порождаемых областей, соответствующих вершинам. Практическая оценка показывает временную сложность этапа кластеризации точек всех вершин примерно в 0(N log N). Алгоритм этапа кластеризации организуется в один цикл просмотра всех вершин дерева и легко распараллеливается.

Сложность по времени алгоритма этапа иерархической склейки точно оценить непросто. При операции склейки имеет место цикл по всем точкам множества Р, с перебором для каждой точки множества G и проверкой условия близости. Время может быть 0(1) в случае пустого множества G и О(Мд) в случае непустого G, где g — количество точек в G. Очевидно, в среднем g << M. Есть ещё внешние циклы с перебором вершин, что для прохода для верхнего уровня иерархии с числом вершин V даёт верхнюю оценку 0(MgV) = O(Ng), поскольку значение MV близко N. Примерно такую же оценку даёт цикл склейки каждого уровня иерархии дерева. Экспериментальные оценки показали, что сложность по времени не выше 0(N log N). Заметим, что этап иерархической склейки распараллеливается не так эффективно, как этап кластеризации, так как по мере объединения вершин дерева объём точек в их областях увеличивается, а количество вершин уменьшается. В целом можно оценить сложность по времени всего модифицированного алгоритма от 0(N) до 0(N log N).

Были проведены сравнительные эксперименты. Использовался компьютер с процессором Intel Xeon W3565 (создаётся до 8 потоков), 3.2 Ггц, память 12 Гбт. Время выполнения немодифицированной программы для 770 ООО точек тестового множества составило 12635 сек. (3.5 часа). Модифицированная программа с распараллеливанием выполнилась за 20 сек. (выигрыш в 632 раза). При этом кластеризация выполнялась 9 секунд, склейка — 11 секунд. Модифицированная программа без распараллеливания выполнилась за 60 сек. (выигрыш в 211 раз). При этом кластеризация выполнялась 44 секунды, склейка — 16 секунд. Таким

Рис. 6: Двумерная иллюстрация декомпозиции функции тетраэдризации

образом, если немодифицированный алгоритм может выполняться за приемлемое время (несколько десятков минут) для множества не более 50-100 тысяч точек, то модифицированный — для десятков миллионов точек.

2 Декомпозиция и распараллеливание воксельного алгоритма

В геоинформатике известен подход дополнения друг другом векторных и растровых представлений двухмерных объектов. В трехмерных задачах микросейсмики аналогичный приём по сочетанию векторных и воксельных представлений пространства точек-событий показывает свою результативность [1]. Покажем подход создания воксельного алгоритма на примере задачи моделирования объёмных тел.

В векторном случае построение моделей тел обычно опирается на тетраэдризацию Делоне с последующим исключением тетраэдров, у которых радиус описанной сферы превышает заданный параметр U, для моделирования невыпуклых тел с внутренними полостями. Тетраэдризация имеет те же проблемы, что и алгоритм кластеризации: трудности по декомпозиции данных для распараллеливания, квадратичный или кубический рост вычислительной сложности от количества точек. Используемый алгоритм тетраэдризации имеет сложность 0(N2 log N). Переход от векторного представления облака точек-событий к вексельному и обратно позволяет решить эти проблемы.

При вокселизации исходных данных создаётся воксельное пространство, где воксели могут иметь определенные свойства: «событие» (то есть содержат исходную точку-событие), «фон», «тело» и др. В алгоритме по воксельному пространству перемещается маска кубической формы, с ребром равным или больше 2U. Шаг движения маски равен U. Маска движется построчно и послойно. Тетраэдризация проводится только в пределах текущего положения маски для вокселов типа «событие». Тетраэдры, удовлетворяющие выбранному выше критерию, заполняются вокселями со значением «тело». Результаты заполнения тетраэдров для очередного положения маски накладываются на перекрывающиеся области ранее вычисленных положений маски с приоритетом вокселей со значением «тело» или «событие». При этом размер маски и шаг её перемещения в результате наложения позволяют корректно определять свойства вокселей. На рис. 6 приведена

Рис. 7: Слияние с наложением тетраэдизованных зон вексельного пространства

двумерная иллюстрация четырёх шагов алгоритма: по два шага для двух соседних строк движения маски, на рис. 7 — слияние тетраэдризованных участков.

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

Заключение

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

Список литературы

[1] Ерохин Г.Н., Алсынбаев К.С., Брыксин В.М., Савеленко В.В., Строков В.П., Козлов А.В., Козлов М.В. Алгоритмы интерпретации и визуализации результатов обработки данных постоянно действующего мониторинга месторождений углеводородов / / Труды международной научной конференции «Марчуков-ские научные чтения — 2017», издательство Института вычислительной математики и математической геофизики Сибирского отделения РАН (Новосибирск), Новосибирск, 2017, с.289-295.

[2] Computational Geometry: Algorithms and Applications. — Springer Science and Business Media, 2008. — P. 259. ISBN 978-3-540-77973-5.

[3] Кластеризация точек на регулярной сети, http://habraliabr.ru/post/138185/

[4] Алсынбаев К.С. Алгоритмы определения тел объёмных объектов в трехмерном нерегулярном облаке точек // Вестник БФУ им. И.Канта, 2015 г., № 10. С. 159-165. ISSN 2223-2095, е- ISSN 2310-3698.

Алсынбаев Камил Салихович — к.т,.н., зав. лаб. НИИ прикладной информатики и математической геофизики Балт,ийского федерального университ,ет,а имени И.Кант,а;

e-mail: [email protected]. Дата поступления — 30 апреля 2019 г.