Стандартный кросскорреляционный алгоритм

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

 

Рисунок 2 - Стандартный кросскорреляционный алгоритм вычисления вектора скорости по изображениям частиц в потоке

 

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

В стандартном алгоритме предполагается, что все частиц внутри области совершают одинаковое перемещение. В реальности это не совсем так, в области может существовать ненулевой градиент скорости, приводящий к ошибке определения смещения. Влияние эффекта градиента скорости будет рассмотрено ниже. Таким образом, интерпретация пространственной корреляционной функции в терминах смещения трассеров возможна только для несжимаемого потока[1]PTV (Particle Tracking Velocimetry). с однородным распределением «идеальных» частиц – трассеров. Если эти условия не выполняются, положение максимума корреляционной функции определяется не только движением потока, но и распределением частиц внутри него. В этом случае обработка может быть проведена с использованием другого метода, например

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

где f и g – интенсивности первого и второго кадров изображений, представляет собой свертку. Метод вычисления свертки зависит от конкретной реализации алгоритма обработки. Обычно используют либо непосредственное вычисление сумм (DCC – Direct Cross-Correlation), либо корреляционную теорему Винера – Хинчина:

где J и J^-1 - прямое и обратное преобразование Фурье соответственно (FFT-CC – Fast Fourier Transform Cross-Correlation). При вычислении преобразования Фурье используется алгоритм быстрого преобразования Фурье (БПФ), который дает преимущества во времени обработки перед прямым расчетом корреляционной функции.

Преобразование Фурье подразумевает периодическую трансляцию расчетной области на бесконечность таким образом, что на границах возникает разрыв данных. Для уменьшения краевых эффектов на границах областей, а также для уменьшения клада в корреляционную функцию частиц, находящихся на границе области, для которых нет парных частиц на втором кадре изображения (эффект «потери пары»), может применятся техника весовых функций (рисунок 3). Метод дополнения нулями или свертки с функцией «прямоугольного окна» обязательно должен применятся с вычитанием средней интенсивности по области, дополненной нулями, чтобы убрать искажения корреляционной функции, связанные с резким изменением значений интенсивности на границе окна. При этом также может применятся нормировка корреляционной функции для приведения области ее значений к интервалу (-1;1). Корреляционная функция, записанная в терминах преобразования Фурье, будет иметь следующий вид:

где  W=W(x,y)- весовая функция, сигма-квадрат- дисперсия значений интенсивности внутри расчетной области.

 

Рисунок 3 - Периодическое продолжение расчетной области при вычислении корреляции через преобразование Фурье: краевой эффект связан с частицами на продолжении области;

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

б- дополнение нулями (заштрихованная область) для устранения эффекта периодичности, краевой эффект пропадает

 

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

По определенному наиболее вероятному смещению частиц d и интервалу времени между вспышками лазера дельта-т рассчитывается наиболее вероятная скорость частиц внутри расчетной области:

---