Применение суперкомпьютерных технологий для решения задачи выделения космических объектов в цифровом изображении Текст научной статьи по специальности «Компьютерные и информационные науки»

| ISSN 2221-1179 Вестник Концерна ПВО «Алмаз - Антей» | №2, 2015

| Информатика |

УДК 004.383.2::528.854.2/.4

© А. П. Коновальчик, 2015

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

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

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

Введение

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

повышением разрешающей способности радиолокационных и оптических средств воздушно-космической обороны, необходимой частоты обновления информации, количества гипотетических целей;

увеличением сложности законов эволюции целевой обстановки.

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

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

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

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

Исследования алгоритма обработки изображения проводились на разработанном в ОАО «Концерн ПВО «Алмаз - Антей» реконфигурируемом проблемно-ориентированном вычислительном комплексе (РПВС-К) «Орфей-К» (рис. 1) [3], который является типичным представителем суперЭВМ с ускорительными платами на базе программируе-

Рис. 1. Реконфигурируемый проблемноориентированный вычислительный комплекс «Орфей-К»

82

| ISSN 2221-1179 Вестник Концерна ПВО «Алмаз - Антей» | №2, 2015

/------\

Рис. 2. Структура суперЭВМ «Орфей-К»

мых логических интегральных схем (ПЛИС) (рис. 2).

Структура суперЭВМ «Орфей-К» является гибридной и включает (см. рис. 2):

а) универсальную компоненту, реализованную как обычный вычислительный кластер из 16 специализированных универсальных вычислительных реконфигурируемых вычислительных блоков (СУРВБ) на основе двух процессоров Intel Xeon, соединённых высокоскоростной сетью Infiniband;

б) специализированную компоненту, состоящую из 32 реконфигурируемых вычислительных блоков (РВБ), попарно подключенных к СУРВБ по шине PCI Express; каждый РВБ включает 8 ПЛИС Xilinx XC6VSX475T.

Основные характеристики суперЭВМ «Орфей-К»:

число СУРВБ.............................16

число процессорных ядер в решающем поле (Intel Xeon

E5620 2,4 ГГц 4 ядра)..................128

число ПЛИС на РВБ........................8

число ПЛИС в решающем поле.............256

скорость чтения/записи в системе хранения данных, Мбайт/с..............................~ 250

пропускная способность сети Infiniband, Гбайт/c.4

суммарный объём оперативной памяти в решающем

поле, Гбайт............................192

емкость системы хранения данных, Тбайт..12

пиковая производительность универсального процес-

сора:

float (1ядро), Гфлопс..................19,2

double (1ядро), Гфлопс..................9,6

float (РПВС-К в целом), Тфлопс.........2,46

double (РПВС-К в целом), Тфлопс........1,23

Постановка задачи

Известный алгоритм обработки изображения для выделения космических объектов основан на замене задачи выделения протяжённых следов задачей обнаружения фрагментов фиксированной длины с последующим их группированием [3]. Фрагменты следов и звёзды обнаруживаются в скользящем по изображению прямоугольном пиксельном окне с классификацией по признаку «неподвижный объект/ фрагмент следа объекта». Накопление сигналов вдоль гипотетических фрагментов следов проводится в пределах окна. Максимальное число наиболее правдоподобных объектов каждого типа ограничивается некоторой априори задаваемой константой. Число группируемых фрагментов при этом значительно меньше числа группируемых отметок при традиционном подходе, основанном на пороговой обработке сигналов в пикселях.

Решение задачи

На вход алгоритма подаётся прямоугольное мозаичное изображение (16 бит серого) в виде матрицы размером NX (по горизонтали) на NY

83

| Информатика |

| ISSN 2221-1179 Вестник Концерна ПВО «Алмаз - Антей» | №2, 2015

| Информатика |

(по вертикали) целых чисел (16 бит без знака): ||j i=0...NX- 1, j=0...NY- 1, где NX, NY - размер входного изображения (фиксированные).

В пространственном прямоугольном пиксельном окне размером WLх WL (WL - нечётное целое, фиксированное для каждой отдельно взятой реализации алгоритма) вокруг выбранного пикселя изображения выполняется оценка амплитуды сигнала, уровня фона, вычисление решающей статистики (максимума логарифма функции правдоподобия) и выбор оптимальной гипотезы, приводящей к максимуму функционал (решающая статистика) [3]:

L j(k, a b) =

WL-1 WL-1

P

z

WL-1

2

+p, J+p

WL-1

2

[a ■ K,

+ , О)

где к = 0...K - гипотезы;

K = 2-WL;

WL = 3, 5, 7, 9, 11;

ND - максимальное число обнаруживаемых фрагментов (фиксированное);

а - амплитуда сигнала;

b - амплитуда фона;

Fkq - матрицы опорных отметок обнаруживаемых космических объектов (шаблонов) размером WLxWL для K = 2WL гипотез, к=0... K-1 (рассчитываются заранее для каждой отдельно взятой реализации алгоритма).

Оптимальная гипотеза к* для ij-пикселя изображения Y^ при расчете выражения (1) выбирается по следующим критериям:

к * = arg max {{L*0 }u \l\ \ ak > o}},

где L*k - максимальное значение функционала Ljk, a, b), полученное нахождением максимизирующих правую часть выражения (1) значений a = ак* и b = Ьк*.

Значения а** и Ьк* для каждого ij-пикселя изображения Yiq вычисляются по формулам:

а* = Ска • (SYFk • S - SY • SFk ), k > 0, a* = 0;

b* = C* • (57 - a* • ), Ik > 0, *

5 ’

где C

_____1______.

SfLS^SFLSF^ ’

ИХ-1

tfL-1

^ ' = Z Z<^, >2;

J^L-1 J^L-1

p=---------g=---------

2 2

WL-1

WL-1

S = £ Z1 = (WL )2;

WL-1 WL-1

p=-------q=--------

2 2

WL-1

2

WL-1

2

SF = ^ Z FL-

WL-1 WL-1

p=—-— q=

2

WL-1

2

WL-1

SYF‘ = ^ Z Yp.,+p ■ FL) •

WL -1 WL -1

p=-------q=-------

2 2

WL-1

2

WL-1

2

sr =2 ZY

WL-1

-q=-

ck _ 1

4 _ S'

WL-1

2

2+P, 1 + P ,

P

2

Выражение для L*k с учетом значений ак* и Ьк* может быть представлено в виде:

L -

-Ч ~

а при к

2 • а* • SYFk - а* • а* • SYFk + b*k • b* • S, = 0

L

SY • SY S

Из всего набора гипотез, получивших положительную оценку амплитуды сигнала ак, и гипотезы отсутствия полезного сигнала (к = 0) по вышеприведенной формуле выбирается гипотеза с максимальной величиной Lk.

В результате такого попиксельного выбора оптимального значения формируются две матрицы промежуточных значений размером (NX-WL + 1)х(ЖК - WL + 1). Первая из них содержит номера выбранных гипотез к*, вторая - величины функционала Lk * для каждого пикселя, соответствующие выбранным к*:

к

arg(4,!к*,a*,b*)) ;

к *

Л

к„(к а\ Ъ')

где n, m = 0...(NX - WL + 1).

Ввиду того, что не всегда возможно определить порог предварительно, не проведя об-

84

| ISSN 2221-1179 Вестник Концерна ПВО «Алмаз - Антей» | №2, 2015

/------\

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

На основе массива величин Л формируется неупорядоченный список. Полученный список сортируется по убыванию. В качестве порогового значения TD выбирается (ND+^-й элемент упорядоченного списка (ND-элемент при индексации от 0).

Отсев подпороговых значений Для каждого элемента матрицы Л производится сравнение:

Л, > т.

Если это условие не выполняется, соответствующий элемент матрицы k зануляется.

Операция гарантирует, что число обнаруженных отметок не превысит ND.

Результат пороговой обработки Результат обработки формируется из прошедшей пороговую обработку матрицы k дополнением до матрицы размером NX на NY нулевыми столбцами и строками по (WL-1)/2 с каждой из четырех сторон.

Минимизация объёма выходных данных

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

Реализация на кристалле адаптивной пороговой обработки позволяет отказаться от вывода с кристалла ПЛИС матрицы Л, что даже в грубой реализации вывода массива k (то есть в полном его объёме) преуменьшает объём выходных данных по сравнению с входом.

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

если в процессе сортировки запоминать сопутствующие величинам статистик координаты. Фрагментация кадра

На первом этапе обработки вычисления, выполняемые в прямоугольных пиксельных окнах для двух соседних пикселей изображения, взаимно независимы. Пересекаются лишь входные массивы данных, и фрагментация неизбежно приведёт к дублированию входных данных на границах окон. Объём дублируемых данных приблизительно (WL-1)N элементов для квадратного фрагмента размером N*N, что соотносится с площадью фрагмента как (WL-1)/N.

Второй этап обработки связывает результаты первого этапа процедурой сортировки. Однако в силу физических особенностей решаемой задачи в определённых пределах фрагментация допускается. За грубую оценку минимального размера фрагмента можно выбрать 32^32. При допущении концентрации ложных отметок в 1/16 от площади фрагмента такой выбор позволит всегда обнаруживать след «яркого» объекта целиком. В целях же обнаружения «тусклых» объектов размер фрагмента стоит увеличить хотя бы до 64x64. Одновременно это позволит обнаруживать более одного объекта во фрагменте кадра.

Таким образом, имеет смысл производить фрагментирование кадра дважды: на вход кристалла ПЛИС подавать фрагменты большего размера (в целях сокращения накладных расходов дублирования данных), а пороговую обработку производить над фрагментами меньшего размера.

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

Анализ вычислительных особенностей алгоритма

Упрощенная блок-схема реализации алгоритма на суперЭВМ может быть разбита на четыре

85

| Информатика |

| ISSN 2221-1179 Вестник Концерна ПВО «Алмаз - Антей» | №2, 2015

| Информатика |

этапа:

загрузка изображения; осуществление вычислений, не зависящих от изображения;

обработка изображения; адаптивная обработка.

Проведя анализ каждого из этапов, можно сделать вывод, что алгоритм удовлетворяет основным требованиям, предъявляемым для реализации на ПЛИС, а именно:

наличие в алгоритме параллельных компактных подалгоритмов, в которых сосредоточена основная вычислительная сложность алгоритма;

невысокая ёмкостная сложность подалгоритмов;

возможность распараллеливания реализуемого фрагмента алгоритма и достаточный уровень потенциального параллелизма;

отсутствие или незначительное число ветвлений в алгоритме;

высокая вычислительная интенсивность реализуемого подалгоритма.

Оценка числа операций Основная трудоёмкость приходится на блок обработки изображения. Общее число операций

- (NX-WL)*(NY-WL)*K.

Оценим количество операций на 1 байт данных (входных и выходных) для всех шаблонов.

Оценка объёмов памяти

Входные данные составляют ~25% (250 тыс. пикселей) от суммарного объёма входных и выходных данных.

Объем памяти под шаблоны K* WL2 слов. Для 100 шаблонов 11*11 потребуется 12 КБайт при однобайтовом элементе, 24 КБайт - при двухбайтовом.

При наличии в ПЛИС XC6 VSX475T

— 4 Мбайт блочной памяти возможно использование до 0,5 Мбайт для обрабатываемого изображения, т. е. 1 ПЛИС может обрабатывать фрагменты изображения до 512*512 пикселей, а для обработки изображения 4008*2672 пикселей потребуется от 64 ПЛИС (т. е. 8 РВБ) [2]. Оценка производительности

При задействовании всей суперЭВМ «Орфей-К» и пропорциональном уменьшении

размеров фрагмента (в 4 раза) скорость обработки может достигнуть 340-600 фрагментов на ПЛИС.

Устоявшаяся скорость обмена данными с 1 ПЛИС (при работе со всеми 16 ПЛИС) составляет ~ 80 Мбайт/с в каждую сторону. Поскольку выходные данные в 7 раз превышают входные, именно это является ограничивающим фактором.

С одной ПЛИС можно считывать до —6*106 данных для пикселей изображений и загружать данные для —42*106 пикселей, т. е. обрабатывать —23 фрагмента изображений в секунду (если считать только загрузку, то 160). При уменьшении фрагмента и увеличении числа задействованных ПЛИС пропорционально возрастает достижимая скорость обмена и обработки фрагментов соответственно. На «Орфей-К» в полной конфигурации (256 ПЛИС) соответственно возможна обработка со скоростью в 4 раза большей (92 и 640 фрагментов/с на ПЛИС соответственно).

Пропускная способность системы хранения данных позволяет извлекать из неё в секунду не более 10-15 изображений размером 20-32 Мбайт и соответственно сохранять данные обработки не более 2-3 изображений.

Пропускная способность сети Infiniband 32 Гбит/с позволяет загружать из сервера загрузки, управления и мониторинга (СЗУМ) данные не более чем 128-200 изображений в секунду и сохранять результаты обработки не более чем 18-30 изображений в секунду.

При сокращении объемов возвращаемой информации скорость обработки может быть увеличена вплоть до значений, лимитируемых производительностью ПЛИС: порядка 100 изображений в секунду для 100 шаблонов 11*11 или 15*15.

Декомпозиция задачи и оценка возможности реализации на ПЛИС

Параллельность выполнения задачи естественным образом достигается разбиением исходного изображения на фрагменты. При обработке на суперкомпьютере изображение делится между ПЛИС, а внутри ПЛИС - между логическими вычислительными устройствами (ВУ), работающими параллельно. При проведении

86

| ISSN 2221-1179 Вестник Концерна ПВО «Алмаз - Антей» | №2, 2015

/------\

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

Внутри ПЛИС массив входного изображения (ввиду большого объёма) должен храниться в блоке памяти (BRAM), который позволяет хранить минимально 210 точек изображения (например, квадрат размером 32*32).

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

Так как протокол обмена данными с BRAM последовательный, то для ВУ предполагается естественным следующий порядок работы: для каждой точки последовательно выполнять сличение квадрата изображения с центром в данной точке с шаблоном из списка; далее в одном блоке ВУ выполняются вычисления, требующие многократного чтения из памяти (это суммы SY и SYF), а во втором - вычисления, не требующие обращений к памяти, т. е. досчитывать функционал L(k, i, j). Такая организация вычислений позволяет блокам работать параллельно (с разными шаблонами).

Время обработки одной пары «пикселшаблон» в первом блоке определяется временем вычисления суммы SYF (более простая сумма SY вычисляется попутно). При считывании из памяти по одному элементу массивов изображения и шаблонов сумма вычисляется за WL2 тактов. Время обработки пары «точка - шаблон» во втором блоке можно, при распараллеливании процесса, свести к 9-10 тактам (минимальной длине периода работы первого блока). Тогда время обработки в ВУ фактически будет определяться временем работы первого блока.

Объем ВУ не менее двух 16-Кбит BRAM для хранения ||1Р|| (поскольку к одному фрагменту изображения может обращаться одновременно до четырёх ВУ, а элементы блочной памяти имеют 2 порта; иначе можно хранить в одном ВУ фрагмент, охватывающий все необходимые для анализа точки; в данном случае, его размер 42*42); один блок памяти для хра-

нения результатов (k, Л); несколько (не более 15, зависит от реализации) DSP-блоков для выполнения умножений и делений.

Дальнейшие вычисления в ПЛИС сводятся к отбору результатов с наибольшими значениями. Для каждого фрагмента размером 32*32 достаточно выбрать около 10 наилучших результатов: это выполняется на фоне остальной обработки посредством вставки очередного значения в 10-элементный список и требует нескольких тактов.

Окончательно списки, выработанные различными ВУ, должны быть объединены в общий список из примерно 50 элементов. При такой реализации в одной ПЛИС Virtex6-475SXT можно разместить не менее 128 ВУ. При рабочей частоте 100 МГц время обработки фрагмента размером 32*32 на одном ВУ и при времени обработки всеми (128) ВУ изображения размером 217 точек оценивается как ~ 0,06 с.

Таким образом, обработка изображения размером 225 точек за 0,1 с достижима на комплексе из 256 ПЛИС.

Структурная схема элементарного канала обработки в ПЛИС

Элементарный канал обработки выполняет подсчёт функционала (1) для небольшого фрагмента изображения (ориентировочно 32*32 точки).

Структурная схема канала обработки в ПЛИС приведена на рис. 3.

Реализация алгоритма на архитектуре «Орфей-К» состоит из следующих этапов:

рассылка частей изображения по вычислительным узлам;

«нарезка» изображения на фрагменты; загрузка фрагментов в ПЛИС; обработка фрагмента в ПЛИС; выгрузка результатов обработки из ПЛИС;

сборка результатов, преобразование форматов и сортировка;

сборка частей результата.

Реализация алгоритма на архитектуре «Орфей-К» позволила обработать изображение размером 10 мегапикселей за 35 мс, что более чем в 3900 раз быстрее, чем на универсальном процессоре.

87

| Информатика |

| ISSN 2221-1179 Вестник Концерна ПВО «Алмаз - Антей» | №2, 2015

| Информатика |

Рис. 3. Структурная схема канала обработки в ПЛИС

На рис. 4 представлены исходное (а) и обработанное (б) изображения с выделенным треком (следом) реального космического объекта. Изображения получены с телескопа VT-53 с апертурой 12 см (ПАО «МАК «Вымпел»), размер изображения 4008*2672 пикселей.

Таким образом, на суперЭВМ «Орфей-К» целесообразно решать задачи с высокой вычислительной сложностью, которые за необходимое время не могут быть решены на автоматизированном рабочем месте на базе современной ПЭВМ. В качестве ориентира может быть принят уровень в 27*1015 операций,

а б

Рис. 4. Исходное (а) и обработанное (б) изображения с выделенным треком (следом)

реального космического объекта

88

| ISSN 2221-1179 Вестник Концерна ПВО «Алмаз - Антей» | №2, 2015

/------\

что соответствует более 16 ч счётной работы на современной ПЭВМ.

Выводы

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

2. Функциональный параллелизм на верхнем структурном уровне рассматриваемого алгоритма обеспечивает совмещение во времени вычислений в ПЛИС с вычислениями в универсальных процессорах и с обменами данными, а наличие в подалгоритмах функционального параллелизма - возможность организации макроконвейерной обработки в ПЛИС. Это позволяет реализовать алгоритм автоматического обнаружения в 10-мегапиксельном цифровом изображении перемещающихся за время экспозиции малоконтрастных космических объектов с неизвестными орбитами за 35 мс, что приблизительно в 3900 раз быстрее,

чем на универсальном процессоре.

3. Специализированный вычислительный комплекс «Орфей-К» на базе ПЛИС обеспечивает реализацию рассмотренного алгоритма в режиме реального времени и может быть использован в качестве базового вычислителя для перспективных оптических средств мониторинга космического пространства Минобороны России, Роскосмоса и др.

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

1. Промежуточный научно-технический отчет. НИЭР «КИМС-Вымпел-2014». М.: ОАО «МАК «Вымпел», 2014. 435 с. Инв. 2657.

2. Научно-технический отчет. ОКР «Орфей-К». М.: ОАО «Концерн ПВО «Алмаз - Антей», 2011. 889 с. Инв. 91-214-К.

3. Колесса А. Е, Репин В. Г. Робастный адаптивный алгоритм выделения отметок от целей в цифровом изображении // Космические ин-формационно-управляющие системы. 2009. Вып. 3. С. 203-210.

Поступила 27.07.15

Коновальчик Артем Павлович - советник генерального директора ОАО «Концерн ПВО «Алмаз - Антей», г. Москва.

Область научных интересов: суперкомпьютерные технологии, цифровая обработка сигналов.

89

| Информатика |