CC BY
207
Вычислительные технологии
Том 12, № 3, 2007
АДАПТИВНЫЕ АЛГОРИТМЫ ОБРАБОТКИ ИЗОБРАЖЕНИЙ ЧАСТИЦ ДЛЯ РАСЧЕТА МГНОВЕННЫХ ПОЛЕЙ СКОРОСТИ
М.П. Токарев, Д. М. Маркович, А. В. Бильский Институт теплофизики СО РАН, Новосибирск, Россия e-mail: [email protected]
State-of-the-art of research in the field of numerical algorithms for processing images obtained by Particle Image Velocimetry (PIV) development is presented. Implementations of the procedures, allowing to receive 2D2C velocity fields, are described. On the basis of the generated synthetic images of particles the analysis of accuracy and spatial resolution of the mentioned procedures is performed.
Введение
Методы визуализации для изучения потоков были известны задолго до появления электронных вычислительных машин. Первые наблюдения за течением жидкости в водоемах при помощи естественных природных трассеров были описаны еще Леонардо да Винчи. Людвиг Прандтль (1875-1953) использовал взвесь из частиц слюды на поверхности воды для анализа обтекания цилиндров, призм и профилей крыла в экспериментальном канале.
Постепенно от качественных наблюдений произошел переход к измерению количественных характеристик течений, и к 60-м годам XX столетия сформировалось широкое направление в диагностике, известное как "стробоскопическая визуализация". Принцип стробоскопической трассерной визуализации заключается в измерении смещения трассеров в заданном сечении потока жидкости или газа за известный интервал времени. Областью измерения служит плоскость, освещаемая световым ножом. Результатом измерения является мгновенное поле скорости в измерительной плоскости. Одной из наиболее показательных отечественных работ того времени является статья [1], в которой приведены результаты измерения полей скорости и турбулентных пульсаций в пограничном слое на пластине, полученные путем ручной обработки трассерных картин, зарегистрированных на фотографической пленке. Однако ручная обработка была трудоемка, занимала длительное время и не позволяла получать достаточного объема данных для расчета статистических характеристик. Данная ситуация являлась типичной для всех исследовательских групп в мире, занимавшихся количественной визуализацией потоков.
Появление термина PIV (Particle Image Velocimetry) — международное название метода цифровой трассерной визуализации, связывают с работой [2], в которой метод PIV был выделен как частный случай метода лазерной спеклометрии LSV (Laser Speckle Velocimetry) [3]), базирующегося на оптическом преобразовании Фурье яркостных картин.
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.
Следует отметить, что методы LSV развивались изначально применительно к анализу деформаций при нагружении твердотельных образцов. Метод PIV, вследствие более низкой концентрации трассеров, при которой различимы отдельные частицы, более подходит для анализа структуры турбулентного течения, так как влияние второй фазы на поток в этом случае незначительно.
Развитие цифровой и компьютерной техники и, как следствие, применение цифровых методов регистрации изображений и обработки данных для стробоскопической трассерной визуализации еще более укрепило позиции PIV среди оптических методов исследования потоков, сократив время регистрации и обработки на порядки. Применение корреляционного анализа для определения смещения частиц [4], а также изобретение цифровой кросскорреляционной камеры для записи пар изображений [5] обеспечили еще больший прогресс и создали условия для дальнейшего совершенствования алгоритмов обработки.
Метод цифровой трассерной визуализации, относящийся к классу бесконтактных оптических методов, позволяет регистрировать мгновенные поля скоростей в плоскости изме-рения1. Одним важнейших его преимуществ является отсутствие возмущающего влияния на поток. К достоинствам метода можно также отнести широкий динамический диапазон измеряемых скоростей (порядка 500:1), что позволяет использовать его для исследования сложных турбулентных течений. Еще одно преимущество PIV состоит в возможности набирать и обрабатывать на обычном персональном компьютере значительный объем экспериментальных данных для расчета статистических характеристик течения.
К ограничениям PIV можно отнести неидеальность трассирующих частиц (размер, плотность), приводящую к тому, что трассеры не всегда точно следуют за потоком. Кроме того, диаметр используемых частиц ограничивает размер элементарной области измерения снизу, а использование более мелких частиц сопряжено с влиянием броуновского движения на их смещение, а также накладывает ограничения на мощность источника излучения и чувствительность регистрирующей аппаратуры. Конечный объем элементарной измерительной области, для которого рассчитывается вектор скорости, приводит к пространственному усреднению по размеру ячейки подобно тому, как это происходит при численном моделировании2.
Область применения метода PIV довольно обширна. Она включает в себя фундаментальные научные исследования, направленные на изучение динамики и масштабов вихревых структур в потоках жидкости и газа, получение дифференциальных и статистических характеристик, а также оценку достоверности математического моделирования и коррекцию численных моделей для турбулентных потоков [6, 7]. Помимо этого, метод PIV широко применяется и в прикладных исследованиях. Здесь можно выделить задачи оптимизации обтекания летательных аппаратов и судов в авиастроении, кораблестроении, конструкций промышленных агрегатов в энергетике и нефтегазовой промышленности, изучение процессов в двигателях внутреннего сгорания. Постепенно увеличивается число приложений PIV для задач исследования микропотоков, таких как управление процессами синтеза веществ в микрофлюидных аналитических системах (^TAS), физическое моделирование работы искусственных сосудов и клапанов в медицине и др.
На протяжении последних 20 лет развития и применения метода обозначилось несколь-
ХВ настоящее время начинают активно развиваться модификации метода, позволяющие измерять поля скорости в объеме (голографический PIV, томографический PIV).
2Данные проблемы в существенной степени преодолеваются применением современных модификаций алгоритмов обработки, в том числе описанных в данной работе адаптивных методов, а также комбинацией методов PIV и PTV (Particle Tracking Velocimetry).
ко направлений, связанных с оптимизацией и совершенствованием методов обработки данных с целью улучшения качества получаемых данных. Прогресс в развитии компьютерной техники позволяет разрабатывать и совершенствовать методы обработки с сохранением приемлемого для исследователей времени расчета. Это увеличение динамического диапазона измеряемых скоростей [8], улучшение пространственного разрешения метода [9], применение кросскорреляционного анализа для течений с существенным градиентом скорости, таких как пограничный слой и слой смешения [10]. К этому можно добавить общую проблему выявления источников систематической и случайной погрешностей расчетных процедур и минимизации этих погрешностей [11]. Одним из таких источников можно назвать проблему группировки рассчитанных смещений частиц на изображении около целочисленных значений (peak-locking), что приводит к существенным искажениям реальной формы распределений компонент скорости по полю [12]. Существует проблема детектирования и отсева ошибочных векторов, появление которых связано с влиянием шума и ошибками в определении положения максимума корреляционной функции [13], неизбежными при автоматизированной обработке больших массивов данных. Все указанные проблемы тесно связаны между собой и не могут рассматриваться отдельно друг от друга. Целью данной работы является представление адаптивного похода к кросскорреляционной обработке изображений частиц, который позволяет устранить или заметно уменьшить часть указанных выше проблем.
В работе освещаются современные подходы к реализации адаптивных методов и связанные с ними ограничения. В части 1 и 2 обсуждаются базовые принципы корреляционного анализа PIV-изображений и ограничения стандартного алгоритма. Часть 3 посвящена обсуждению адаптивного алгоритма, его возможностям в решении проблем, присущим стандартному алгоритму. В части 4 обсуждаются ограничения адаптивного алгоритма.
1. Стандартный кросскорреляционный алгоритм
На начальном этапе развития методов цифровой обработки PIV-изображений для определения смещений частиц широко использовался базовый кросскорреляционный алгоритм [8], который в дальнейшем был назван стандартным. Стандартный алгоритм состоит из следующих основных операций (рис. 1):
— разбиение пары изображений на элементарные расчетные области равного размера;
— расчет кросскорреляционной функции;
— нахождение максимума на корреляционной функции;
— подпиксельная интерполяция максимума корреляционной функции.
Приведем подробное описание каждой из операций. На вход алгоритма подается два изображения образов частиц, полученных в два последовательных момента времени с известной задержкой между ними. Необходимо получить мгновенное поле скорости, соответствующее мгновенной скорости потока в измерительной области. Для этого вся измерительная область разбивается на элементарные ячейки (расчетные области размером M х N) таким образом, чтобы в каждую расчетную область попало несколько частиц. Для стандартного алгоритма размер и количество расчетных областей в течение вычисления всего поля скорости остается постоянным. Смещение частиц между первым и вторым кадрами изображения может быть найдено путем вычисления функции пространственной корреляции.
Рис. 1. Стандартный кросскорреляционный алгоритм вычисления вектора скорости по изображениям частиц в потоке
В стандартном алгоритме предполагается, что все частицы внутри области совершают одинаковое перемещение. В реальности это не совсем так, в области может существовать ненулевой градиент скорости, приводящий к ошибке определения смещения. Влияние эффекта градиента скорости будет рассмотрено ниже. Таким образом, интерпретация пространственной корреляционной функции в терминах смещения трассеров возможна только
3 (( 11
для несжимаемого потока3 с однородным распределением идеальных частиц — трассеров. Если эти условия не выполняются, положение максимума корреляционной функции определяется не только движением потока, но и распределением частиц внутри него. В этом случае обработка может быть проведена c использованием другого метода, например PTV [14].
Если принять во внимание тот факт, что в цифровой обработке изображений сигналы являются дискретными, функция пространственной корреляции
M N
Щт,п) = f 0 д = Y^ f (i,j+ m,j + n) i=0 j=0
где f и д — интенсивности первого и второго кадров изображений, представляет собой свертку. Метод вычисления свертки зависит от конкретной реализации алгоритма обработки. Обычно используют либо непосредственное вычисление сумм (DCC — Direct Cross-Correlation), либо корреляционную теорему Винера—Хинчина:
R(m,n) = Ö-1[Ö[f ]Ö* [д]],
где ö и ö-1 — прямое и обратное преобразования Фурье соответственно (FFT-CC — Fast Fourier Transform Cross-Correlation). При вычислении преобразования Фурье используется алгоритм быстрого преобразования Фурье (БФП), который дает преимущества по времени обработки перед прямым расчетом кросскорреляционной функции.
Преобразование Фурье подразумевает периодическую трансляцию расчетной области на бесконечность таким образом, что на границах возникает разрыв данных. Для уменьшения краевых эффектов на границах областей, а также уменьшения вклада в корреляционную функцию частиц, находящихся на границе области, для которых нет парных
3Подходы для PIV-диагностики высокоскоростных сжимаемых газовых потоков развиваются в ряде аэродинамических центров и университетов Европы и США. Их рассмотрение выходит за рамки данной статьи.
частиц на втором кадре изображения (эффект "потери пары"), может применяться техника весовых функций (рис. 2). Метод дополнения нулями или свертки с функцией "прямоугольного окна" обязательно должен применяться с вычитанием средней интенсивности по области, дополненной нулями, чтобы убрать искажения корреляционной функции, связанные с резким изменением значений интенсивности на границе окна. При этом также может применяться нормировка корреляционной функции для приведения области ее значений к интервалу (—1; 1). Корреляционная функция, записанная в терминах преобразования Фурье, будет иметь следующий вид:
-(Wf )]S*[Wg -(Wg)]] R(m,n) =-. -,
л/о^ТмШ '
где W = W(x,y) — весовая функция, a2 — дисперсия значений интенсивности внутри расчетной области.
Далее по схеме (см. рис. 1) после вычисления корреляционной функции производится поиск максимального корреляционного пика, отвечающего наиболее вероятному смещению в области. Поскольку изображение представляется в виде конечного числа элементов, найденное положение максимума будет иметь целочисленные координаты в пикселях. Для того чтобы определить смещение частиц с подпиксельной точностью, используют интерполяцию корреляционного максимума кривой (по одной координате) или поверхностью некоторой формы и затем определяют положение максимума численно.
По определенному наиболее вероятному смещению частиц d и интервалу времени между вспышками лазера At рассчитывается наиболее вероятная скорость частиц внутри расчетной области: v = MAt, где M — масштабный коэффициент увеличения системы регистрации. t
2. Ограничения стандартного кросскорреляционного алгоритма
Описанный выше стандартный метод обладает рядом недостатков. Первая проблема, связанная с ограничением динамического диапазона измеряемых скоростей [8], практически
СТ
о"' 1 о" 1
Р
L <
J 1 1
о
о
Рис. 2. Периодическое продолжение расчетной области при вычислении корреляции через преобразование Фурье: краевой эффект связан с частицами на продолжении области; а — смещенный прямоугольник — одно из положений области при вычислении кросскорреляции; б — дополнение нулями (заштрихованная область) для устранения эффекта периодичности, краевой эффект пропадает
а
полностью решена на сегодняшний день. Под динамическим диапазоном здесь понимается отношение максимально возможной определяемой скорости к минимально возможной в пределах расчета одного мгновенного поля скорости. Минимально возможная измеряемая скорость для данной конфигурации эксперимента (размер измерительной области Ь, интервал времени между двумя последовательными вспышками Д£) определяется погрешностью кросскорреляцонного алгоритма определения смещения, порядка 0.1 пикс. Для минимальной реализации кросскорреляционного алгоритма, которую мы будем называть стандартным кросскорреляционным алгоритмом, максимальное значение скорости определяется размером области изображения, по которому ведется расчет [15]. Размер области не может быть увеличен без потери в пространственном разрешении. Решение указанной проблемы было найдено в смещении расчетной области на одном из кадров на предварительно вычисленный вектор скорости, полученный данным же методом, но с большим размером расчетной области. Благодаря такой технике стало возможным вычислять большие по величине смещения без потери в пространственном разрешении (это отражено на рис. 6, о котором речь пойдет дальше).
Второе ограничение связано с точностью расчета полей скорости по изображениям. Под точностью понимается суммарный уровень систематической в и случайной погрешности а выходных данных:
1 N _ в = - = и - ио;
г=1
а
\
1
N
^^(иг - и)2; 5 =
г=1
\
1
N
N
]>>о - иг)2
г=1
где и0 — истинное значение скорости, иг — значение скорости для г-го измерения, и — среднее значение скорости по N измерениям. Полная ошибка 5 связана с систематической и случайной погрешностями (среднеквадратичным отклонением) соотношением 52 = в2 + а2.
Уровень погрешности определяется качеством входных данных, выбором численных методов, применяющихся в процессе обработки данных, а также выбором свободных параметров процедур обработки. Перечислим вначале особенности изображений частиц, которые влияют на погрешность стандартного алгоритма, —
• качество системы регистрации изображений и освещения области измерения:
— дисторсия и аберрация оптического тракта системы регистрации,
— шум цифровой камеры,
— неравномерное распределение интенсивности на изображении,
— отражение, преломление лучей от поверхностей экспериментальных установок;
• качество засева измерительного объема частицами и характеристики потока:
— неоднородная и недостаточная концентрация частиц по всему изображению,
—1еоптимальный размер частиц на изображении (оптимальный диаметр 1.5... 2.5 пикс, здесь и далее в тексте диаметр определяется по уровню е-2 интенсивности изображения),
— уровень градиента скорости в измерительной области потока выше определенного значения,
— существенная трехмерность потока и недостаточная толщина лазерного ножа (приводит к уменьшению количества коррелирующих пар частиц на первом и втором кадрах изображения за счет выхода из области освещения).
В работе [5] анализируются погрешности стандартного метода. Расчет корреляционной функции по измерительной области конечного размера приводит к систематической недооценке определяемого смещения [11]. Так, вклад в главный пик на корреляционной
плоскости вносят частицы, присутствующие в данной расчетной области и на первом, и на втором кадрах изображения. Чем больше величина смещения частиц, тем меньше площадь подобласти, присутствующей на обоих кадрах.
Устранить систематическую ошибку при свертке областей равного размера можно путем поэлементного деления корреляционной функции на результат свертки весовых функций Ш(х,у) = (1 — |х|/Е)(1 — |у|/Е) (для функции "прямоугольного окна"), Е Е"
координаты в области корреляционной функции; Е — размер
2 ' 2
где x, y £ области.
Такой метод увеличивает значение корреляционной функции по краям области в 2 раза, поэтому в случае, если отношение высоты главного пика ко второму по высоте пику на корреляционной плоскости меньше двух, главный пик может определиться неверно. Вектор скорости, соответствующий неверному главному пику, называется ошибочным вектором (международные термины — outlier, false vector). Даже небольшое количество ошибочных векторов способно внести существенную погрешность при расчете статистических данных [16].
Другой метод устранения неявной свертки с весовой функцией заключается в выборе области поиска большего размера, например, в 2 раза больше, тогда свертка областей будет иметь равномерный участок в диапазоне смещений
" F F' F F"
- 2~ 2~
Применение адаптивного метода с повторным расчетом корреляционной функции, который будет обсуждаться в следующей части работы, также является эффективным способом устранения систематической ошибки, связанной с неявной сверткой с весовой функцией. Метод позволяет сместить главный пик ближе к центру корреляционной плоскости, где влияние весового множителя незначительно.
Систематическая и случайная погрешности численных процедур определения смещения частиц в данной работе оценивалась по синтетическим PIV-изображениям. Для этой цели была создана 21 пара искусственных PIV-изображений частиц со следующими параметрами: размер изображения 512 х 512 пикс2, разрядность 8 бит, однородное горизонтальное смещение частиц в диапазоне от 0 до 2 пикс с шагом 0.1 пикс, гауссова форма распределения интенсивности образов частиц. Диаметр частиц на изображении ~ 3 пикс, плотность 0.075 частиц на пикс2. Средняя интенсивность частиц составляла 75 (из 256) со стандартным отклонением интенсивности 15. Шум на изображения не добавлялся для получения явной формы зависимости ошибки от величины смещения частиц.
На рис. 3 представлены зависимости систематической и случайной погрешности от величины используемого смещения. Видно, что коррекция систематической погрешности позволяет добиться уровня ошибки, близкой к нулю. Краевые эффекты, связанные с периодичностью сигнала, при вычислении корреляционной функции через преобразование Фурье дают осциллирующее поведение систематической ошибки (рис. 3). Случайная погрешность прямого метода расчета корреляционной функции несколько меньше, чем для случая с методом расчета через корреляционную теорему.
Многие авторы указывают на неоднородное распределение систематической ошибки в зависимости от величины смещения частиц. При этом величины компонент смещения группируются около целочисленных значений (рис. 4, а). Этот феномен назван peak-locking, он проявляется на стадии определения подпиксельной части смещения. Влияние
эффекта peak-locking на среднее значение скорости и значение пульсаций скорости для турбулентных течений было исследовано в работах [17, 18]. Результаты показали, что данный эффект приводит к переоценке статистических моментов, начиная с моментов второго порядка.
В работе [19] были определены основные источники, связанные с этим феноменом:
1) геометрия сенсора и разрешение. Диаметр частиц на изображении меньше 2 пикс приводит к потере информации при дискретизации изображения и вследствие этого — к увеличению уровня систематической ошибки;
2) обрезание изображений частиц на границе расчетной области;
3) алгоритм подпиксельной интерполяции.
Первый источник ошибки может быть уменьшен, но не устранен полностью, так как дискретизация изображения — это неотъемлемый атрибут регистрации изображения в цифровых методах. Минимизация ошибки достигается выбором коэффициента увеличения оптической системы, при котором диаметр образов частиц становится больше 2 пикс.
Влияние обрезания изображений частиц на границе области можно уменьшить, при-
0.075 0.05
^0.025
•Г
I о
п.
-0.025 -0.05 -0.075
■ —DCC
' ftr Тз
t А А А—А-^ .АЛЛА,,
-- A D ♦ FI . FI --.- СС с коррек Т-СС Т-СС с кор -.- дией рекцией -.- 'оооо0
0.04
0.03
Л, 0.02
0.01
0.5
1
dfmncc]
1.5
Рис. 3. Систематическая (а) и случайная погрешности (б) стандартного алгоритма в зависимости от величины дробного смещения частиц в диапазоне от 0 до 2 пикс для различных способов расчета корреляционной функции с коррекцией систематической ошибки и без нее
Нет интерполяции ■в— Центрирование Параболическая Кривая Гаусса
-4 -3 VX [пикс]
Рис. 4. Плотность распределения осевой компоненты скорости для свободной затопленной струи (тест A, PIV Challenge 2003, Westerweel [23]): сравниваются результаты стандартного и адаптивных PIV алгоритмов CWS, CWD (а); среднеквадратичная ошибка определения подпиксельного смещениия зависит от величины дробного смешения [24] (б)
а
менив адаптивный алгоритм с дробным смещением расчетных областей [20]. При этом количество частиц, попавших на границу области, уменьшается за счет восстановления изображения в промежутках между узлами. Последний метод требует применения эффективных схем интерполяции изображения [21, 22].
Наиболее исследованные типы подпиксельной интерполяции применительно к методу Р1У — трехточечная интерполяция параболой и кривой в виде распределения Гаусса [15]. Минимальную погрешность имеет интерполяция кривой Гаусса (рис. 4, б). Это связано с тем, что образы частиц на изображении (в случае, когда эффекты, связанные с дифракцией, малы) могут быть аппроксимированы поверхностью Гаусса с наилучшей точностью [15], а свертка двух двумерных гауссиан представляет собой поверхность той же формы. Поэтому локализацию положения главного пика на корреляционной плоскости в окрестности его максимума наиболее естественно проводить трехточечной подгонкой кривой Гаусса:
Аи = 1п Я(г — 1,?) — 1п Я(г + 1,?)
Ау
21п Я(г — 1.) — 41п Я(г.) + 21п Я(г + 1, ?):
1п Я(г — 1,?) — 1п Я(г + 1, ?) 21п Я(г — 1,?) — 41п Я(г.) + 21п Я(г + 1, ?):
где (и, у) = (г,?) + (Аи, Ау) — соответственно полное значение смещения, его целочисленное значение и дробное значение.
Трехточечная подгонка кривой эффективно работает для оценки смещения по изображениям с диапазоном диаметров частиц 1.5... 2.5 пикс [11, 15]. Меньшие значения диаметров образов частиц приводят к уменьшению ширины главного корреляционного пика, при этом значения, прилегающие к максимуму справа и слева, становятся меньше уровня шума. Увеличение среднего диаметра частиц приводит к тому, что различия между тремя точками, участвующими в интерполяции, становятся незначительными и определяются случайным шумом. Для больших диаметров частиц (более 3 пикс) эффективно подходит метод центрирования пика [15].
Еще один используемый метод оценки подпиксельного смещения — это градиентный метод [25, 26]. Подпиксельное смещение таким методом определяется путем приближенного решения уравнения для некоторой области размером М х N относительно дробных значений смещения: f (х,у) = д(х + г + Аи,у + ? + Ау). При этом интенсивность второго кадра разлагается в ряд Тейлора в точке, связанной с целочисленным смещением: f (х, у) = д(х + г, у + ?) + дх(х + г, у + ?)Аи + ду(х + г, у + ?)Ау, где дх,ду — пространственные производные, вычисленные, например, по схеме экстраполяции Ричардсона. Подпик-сельные значения смещения (Аи, Ау) определяются из решения системы линейных уравнений по методу наименьших квадратов. Метод градиентов дает меньшую систематическую погрешность: в < 0.03 пикс [25, 26], по сравнению с трехточечной подгонкой, где в < 0.1 пикс при использовании его в стандартном алгоритме. С увеличением размера частиц до О = 10 пикс значение погрешности увеличивается незначительно, как и для метода центрирования пика.
Последний метод подпиксельной интерполяции, который будет рассмотрен здесь, — это подгонка корреляционной функции путем минимизации разницы е(и,у) между основной
R и вспомогательной R' корреляционными функциями:
e(u,v) = ^^(R - R'(u,v))2,
R'(u,v) = f 0 f',
f (x,y) = f (x + u,y + v)>
R(u,v) = f 0 g.
Здесь построение функции R' требует процедуры интерполяции изображения для восстановления значений интенсивности в точках с дробными значениями координат. Такой метод позволяет добиться уменьшения погрешности до уровня ß < 0.01 пикс.
Пространственное разрешение является одним из важных параметров в методе PIV — это связано с пространственным усреднением смещения частиц внутри отдельной измерительной области. В работе [18] было показано, что уровень пространственного разрешения — как размера измерительной области, выбираемого во время проведения эксперимента, так и размера элементарной расчетной области во время обработки изображений частиц — существенно влияет на величину измеренных статистических моментов пульсаций скорости и уровень диссипации турбулентной кинетической энергии. Особенно существенное влияние пространственное разрешение оказывает на характеристики, связанные с дифференцированием исходных полей скорости [27].
По теореме Найквиста, максимально разрешаемый масштаб в потоке при измерениях методом PIV ограничен снизу удвоенным размером элементарной расчетной области. С другой стороны, ограничение на размер расчетной области сверху возникает вследствие уменьшения соотношения сигнал/шум для областей с низкой концентрацией частиц. Экспериментально было установлено [28], что для надежного определения сигнального пика на корреляционной плоскости относительно составляющей шума количество образов частиц внутри расчетной области должно быть не меньше 10... 12 с характерным размером частиц в 2 пикс. Таким образом, характерный минимальный размер расчетной области для стандартного алгоритма составляет 8 пикс. Одним из методов увеличения пространственного разрешения в обработке PIV-изображений является перекрытие расчетных областей. В результате перекрытия уменьшается фактическое расстояние между соседними значениями, однако, как будет показано ниже, перекрытие без применения специальных весовых функций (окон) [9] приводит к спектру сигнала, несколько отличающемуся от случая с аналогичным разрешением, но без применения техники перекрытия.
Передаточная функция T(и) = u(u)/u0(u), где u0 — истинное значение скорости, и = F/Л и Л — длина волны, соответствующая характерному изменению скорости в потоке, для стандартного кросскорреляционного алгоритма c весовой функцией "прямоугольного окна" совпадает с передаточной характеристикой КИХ-фильтра скользящего среднего (THMA — top hat moving average), которая в одномерном случае выглядит как T(и) = sin (пи)/пи [22, 27, 29].
Передаточная функция стандартного алгоритма с произвольной весовой функцией W(x,y) совпадает с передаточной характеристикой соответствующего фильтра Tw(и) (рис. 5), таким образом, вычисленное значение скорости стандартным алгоритмом для различных гармоник равно v = Tw (u)v0.
На рис. 5 видно, что передаточные характеристики для весовых функций THMA (32) и THMA (16) имеют область отрицательных значений на частотах ш=1,5 и 3, что может привести к изменению знака соответствующих гармоник измеряемой скорости. Для устранения этих явлений в корреляционной PIV-обработке часто применяются фильтры
0.6
х £
0.4
/ /> \\ // \ ТНМА(32) к
!, // - ТНМА(7б) 1 -- gauss — - LFC | 1 | \
-16
-8 0
х[пикс]
а
16
ТНМА(32)
- ТНМА(7б)
- gauss
Рис. 5. Весовые функции в горизонтальном сечении (а) и соответствующие им передаточные функции для стандартного алгоритма (б) (теоретические кривые и экспериментальные данные)
с положительным значением передаточной функции на всем интервале частот, например, гауссово окно и ЬЕС-фильтр [9, 29]. Экспериментальные значения передаточной функции получены на наборе синтетических изображений частиц с теми же характеристиками изображения, что и для тестирования точности алгоритма (см. выше), с синусоидальным смещением у0(и) = $\п(2пшх/Г), для Е = 32 пикс, гармоник и = 0.05... 1 с шагом 0.05. Значения передаточной функции вычислялись по формуле
1 N
= (V0i - TMvoi)2; T(u) = 1 i= 1
-
N
*M2/£t,0i, T(u) < 1.
i=1
Пространственный спектр стандартного алгоритма представляет собой квадрат указанной выше передаточной функции. Один из способов расширения диапазона пространственного разрешения — это применение весовых функций для модуляции передаточной функции кросскорреляционного анализа [9]. Таким образом, можно ослабить ограничение на пространственное разрешение, связанное с размером расчетной области Am > 2F. В этом случае минимально возможный измеряемый масштаб определяется расстоянием между соседними узлами поля скорости Am > 2Д либо расстоянием между трассерами Am > 2т в зависимости от их величины.
В заключение стоит отметить, что ошибка в PIV-измерениях имеет спектр белого шума, т.е. равномерно распределена для всех пространственных частот [27]. Дальнейшее повышение уровня пространственного разрешения на стадии численной обработки данных возможно за счет локального (адаптивного) увеличения разрешения. Мгновенное распределение концентрации частиц в потоке неоднородно (они не ведут себя как пассивная примесь, а за счет своей неидеальности, разброса по размерам и других свойств по-разному реагируют на локальные пульсации скорости). Таким образом, некоторые области потока характеризуются более высокой концентрацией частиц и, следовательно, большим соотношением сигнал/шум. Примеры использования концентрации частиц как критерия для построения сетки с неоднородным пространственным разрешением на основе рекурсивного алгоритма двоичного разбиения пространства (BSP — binary space partition) для детального разрешения областей потока с большим градиентом скорости приведены в работах [30, 31].
Также для более точного разрешения областей с градиентом скорости возможно использование метода трассировки отдельных частиц — PTV (Particle Tracking Velocimetry) в комбинации с PIV-алгоритмом [32]. Следует отметить, что в своей традиционной реализации методы PTV в настоящее время применяются все реже. В рамках тестирования методов обработки изображений для цифровой трассерной визуализации на международном тест-симпозиуме PIV Challenge 2005 [33] было продемонстрировано, что несмотря на заявленное высокое пространственное разрешение все реализации PTV-алгоритмов, учавство-вавшие в тестировании, имеют значительный уровень шума на высоких пространственных частотах. Следовательно, метод не дает фактического улучшения пространственного разрешения. В настоящее время известны и другие подходы для обработки трассерных картин с высокой разрешающей способностью. Одним из них является Optical Flow Method (OF), позволяющий добиться пространственного разрешения в 1 пикс [34]. Обсуждение достоинств и ограничений этого метода выходит за рамки данной работы.
3. Адаптивный кросскорреляционный алгоритм
Как было указано выше, для расширения динамического диапазона рассчитываемых скоростей и повышения точности вычисления скорости применяются адаптивные подходы обработки PIV-изображений.
Сущность адаптивной схемы обработки состоит в изменении поведения алгоритма в зависимости от особенностей входных данных. После оценки параметров входных данных оценивается и предлагается наиболее оптимальный путь обработки. Адаптивный метод, описываемый здесь, является итерационным. На каждой итерации метода оценивается уточняющая поправка v^ к полю скорости vk. С использованием текущего приближения vk = vk-1 + vCk, k £ N, v0 = 0 компенсируется смещение трассеров на исходных изображениях и вычисляется остаточное смещение, уточняющее текущую аппроксимацию. Таким образом, можно добиться сходимости результата обработки и увеличить точность получаемых данных.
Наиболее важный элемент итерационного алгоритма — компенсация смещения частиц на изображении. Она достигается выбором расчетных областей на первом и втором кадрах изображения, смещенных друг относительно друга, следующим образом (рис. 6, а) (для итерации k):
т.е. симметрично смещенными относительно центра области (x0,y0), где (x,y) — координаты элементов изображения внутри расчетной области.
Выбор симметричного смещения областей не является обязательным, однако в работе [35] было показано преимущество симметричного смещения областей (CDI — central difference interrogation) как метода, имеющего второй порядок точности по времени, по сравнению с несимметричным смещением (FDI — forward difference interrogation), имеющим первый порядок точности. Таким образом, данная техника предпочтительна для случая больших временных задержек At между регистрацией первого и второго кадров.
g
/ V
0-го порядка
1-го порядка
2-го порядка
Рис. 6. Прямое и симметричное смещение областей (а), компенсация смещения частиц (б): слева — нулевой порядок точности (без учета локального градиента); в центре — первый порядок точности (линейный учет локального градиента); справа — второй (и выше) порядок точности
Смещение расчетных областей позволяет существенно уменьшить эффект "потери пары" за счет лучшего совпадения трассерных картин и тем самым повысить отношение сигнал/шум. Кроме того, высокая степень "совпадения" частиц на изображениях дает возможность (при достаточной концентрации частиц) уменьшить размеры конечной элементарной области, что позволяет повысить разрешающую способность метода без ущерба для качества получаемых данных.
Компенсация смещения позволяет построить поправку V (корректирующее поле) к текущему приближению векторного поля путем вычисления кросскорреляционной функции между областями f и д. В общем случае скорость в любой точке \(х,у) внутри расчетной области можно разложить в ряд Тейлора относительно центра:
v(x,y) = v(xo ,yo) +
dv
dx
d v
(x - x0) + —
(xo,yo) dy
(y - y0) +
(xo,yo)
dv2
dx2
(x - xo)2+
(xo,yo )
+
5 v2
dxdy
(x - xo) (y - y0) +
5 v2
(xo,yo)
dy2
(y - Уо)2
(xo,yo)
+ ... + o(x - xo,y - yo)3
Выбирая подходящий порядок малости поправки смещения, можно получить различные модификации адаптивных алгоритмов. Поправка нулевого порядка соответствует одновременному смещению всех элементов области на один и тот же вектор (рис. 6, б). Известно два метода, соответствующих поправке смещения нулевого порядка: с целым значением вектора смещения [8] (DWS — Iterative algorithm with discrete window shifting) и дробным значением вектора смещения [20] (CWS — Iterative algorithm with continuous window shifting). Метод CWS требует интерполяции интенсивности изображения в точке с дробными значениями координат.
Следующий — первый порядок малости — позволяет учитывать градиент скорости в потоке и соответствует алгоритму CWD (Iterative algorithm with continuous window shifting and deformation) [24]. Первый и второй порядки малости требуют оценки первой и второй
производной скорости в центре расчетной области. Исследование процедуры с использованием поправки второго порядка можно найти в работе [12].
Типичная последовательность действий для итерации k в адаптивном алгоритме АШБ выглядит следующим образом.
1. Получить корректирующее поле скорости V для сетки из областей размером Г алгоритмом с компенсацией смещения частиц, описанным выше, с полем предсказания Vй-1 (если это первая итерация, то поле предсказания нулевое — V0 = 0).
А. Для каждой расчетной области:
I) вычислить компоненты градиента скорости Vv(x0,у0) по центрально-разностной
схеме;
II) для каждого элемента расчетных областей f и g вычислить соответствующие им положения на изображениях 11,12 и интенсивность, используя v(x, у). Интенсивность изображения в точке с дробными координатами находится путем интерполяции по соседним элементам;
III) стандартным кросскорреляционным алгоритмом для областей f и д определить корректирующее значение вектора скорости ^(х,у).
Б. Если это последняя итерация, то vk = Vй-1 + v/k — результат обработки.
2. Провести отсев ошибочных векторов текущего приближения одним из методов фильтрации. Провести интерполяцию текущего приближения vk. Выполнить фильтрацию (сглаживание) векторного поля vk.
3. Провести детализацию сетки, в ином случае перейти к пункту 1. Уменьшить размер расчетной области в 2 раза Гк = Гк-1/2 и интерполировать поле скорости vk на детализированную сетку билинейной интерполяцией. Перейти к пункту 1.
Здесь следует отметить, что интерполяция поля скорости не проводится на последней итерации. Таким образом, по окончании процедуры получаются не интерполированные и не усредненные данные.
Данная процедура выполняется k раз. Как правило, сходимость адаптивных методов достаточно быстрая ^ = 3 ... 5). При этом можно использовать итерации без детализации сетки. Размер области на этих итерациях не изменяется. Таких шагов без детализации сетки может быть несколько. Они служат для улучшения сходимости результатов, но при этом требуют дополнительных вычислительных затрат. В идеале концом обработки можно считать момент полного совпадения первого и второго кадров пары изображений. В действительности, полного совпадения добиться не удается в связи:
• с локальным градиентом скорости (необходима подходящая компенсация градиентов);
• с наличием третьей компоненты скорости (частицы выходят из плоскости лазерного ножа);
• с различием в освещенности изображений.
Описанный адаптивный метод расчета мгновенных полей скорости по Р1У-изобра-жениям был реализован в качестве алгоритма обработки данных в программном обеспечении ActualFlow измерительной системы "ПОЛИС", разработанной в Институте теплофизики СО РАН, основанной на оптических методах измерения Р1У, РТУ, ЫГ, РЫТ [36].
4. Ограничения адаптивных алгоритмов
Адаптивный алгоритм позволяет преодолеть некоторые ограничения стандартного алгоритма. Однако ему также присущи определенные ограничения, наиболее важные из которых рассматриваются ниже.
Увеличение точности в адаптивном алгоритме достигается наилучшим совпадением положений частиц за счет относительного смещения каждого элемента расчетной области на не целое значение методом интерполяции исходного изображения. В работе было исследовано влияние типа интерполяционной схемы на погрешность расчета скорости. На рис. 7, а и б представлены уровни систематической в и случайной а погрешностей в зависимости от типа интерполяционной схемы в сравнении с результатами для стандартного алгоритма без коррекции систематической погрешности, связанной с неявной сверткой с весовой функцией. Видно, что по уровню систематической ошибки интерполяционные схемы Б-эрНпе, э1пс, билинейная интерполяция с коррекцией [37] имеют преимущество перед обычной билинейной и бикубической интерполяциями, погрешность которых больше погрешности стандартного алгоритма. Если рассматривать случайную ошибку, можно заметить, что все схемы по уровню погрешности находятся в последовательности, приведенной выше. Однако существенное отличие в том, что адаптивный алгоритм с применением любой интерполяционной схемы по величине случайной ошибки выигрывает на порядок по сравнению со стандартным алгоритмом.
0.05
-0.05
0.01
0.001
0.0001
0.1
0.01
•— Bilinear Sin corr.
Sine (5x5) ■♦— B-spline (3rd order)
0.001
0.0001
ааа
■ ■
0.5 1
d [пике]
0.5 1
d [пике]
0.01
0.001
♦ FFT-CC a DCC □ Bicubic (4x4) о Bilinear (2x2) a Sine (5x5)
♦ Bilin. Sin corr. (;
♦ B-spline (5x5)
x2)
0.5 1
d [пике]
10 t[c]
1.5
100
Рис. 7. Систематическая погрешность (а), случайная погрешность (б), полная ошибка (в) адаптивного алгоритма в зависимости от величины дробного смещения частиц в диапазоне от 0 до 2 пикс для различных интерполяционных схем в сравнении со стандартным алгоритмом; время расчета поля скорости и полная ошибка в зависимости от интерполяционной схемы (г)
Полная ошибка 8 представлена на рис. 7, в. Из приведенных зависимостей можно сделать вывод, что уменьшение ошибки достигается главным образом за счет уменьшения систематической погрешности. Вторым параметром, который зависит от типа используемой схемы, является время работы алгоритма. Чем больше исходных точек изображения используется при интерполяции и чем "сложнее" его интерполятор, тем медленнее работает алгоритм. На рис. 7, г представлено время расчета поля скорости из 79x63 узлов по изображению размером 1280x1024 пикс2, полученному цифровой кросскорреляцион-ной камерой "Видеоскан" 205/Д-2001 с размером расчетной области 32 пикс, перекрытием 50 % и двумя итерациями для итерационных методов.
Применение адаптивных алгоритмов существенно уменьшает ошибку, связанную с эффектом peak-locking (см. рис. 4, а). Была обнаружена взаимосвязь интенсивности peak-locking c систематической ошибкой, как на рис. 7, а. Аналогично, использование деформации расчетной области (поправки первого порядка) дает небольшое уменьшение эффекта.
Пространственное разрешение адаптивного кросскорреляционного алгоритма исследовалось по экспериментально полученным передаточным функциям. В общем случае теоретическое выражение для скорости итерационного алгоритма можно записать следующим образом: vk = vk—1 + TW (w)(v0 - vk-1) [22, 29, 38]. Второе слагаемое в правой части выражения означает компенсацию смещения, т. е. на k-й итерации кросскорреляционным алгоритмом вычисляется остаточное смещение. Таким образом, передаточная функция для итерации k простейшего итерационного алгоритма выглядит следующим образом:
Ta = 0, k k
k k i k i ^ Tk = 1 - (1 - Tw)k.
rpk rpk — 1 i rp ¡л rpk —1\ A v W )
T A = TA + TW(1 - Ta )
При k =1 она равна передаточной характеристике стандартного алгоритма: T\ (w) = Tw(w). Итерационный алгоритм может быть дополнен цифровой фильтраций поля предсказания fp, корректирующего поля скорости fc и компенсации смещения fsc:
vk = fp(w)vk—1 + fc(w)Tw(w)(vo - fsc(w)vk—1).
При этом можно моделировать передаточные функции фильтров для оптимизации передаточной характеристики адаптивного алгоритма. Легко получить выражение для передаточной функции такого алгоритма:
k-1
Tk = fcTwY, (fp - fscfcTw )г при k > 1.
г=0
С использованием указанной выше теоретической зависимости было изучено поведение адаптивного алгоритма по критериям пространственного разрешения и сходимости. В табл. 1 показано влияние фильтрации на параметры итерационного метода, где 0 означает, что фильтр практически не влияет на поведение алгоритма. Видно, что наиболее выраженным влияние фильтрации оказывается для фильтра ТЩ, который позволяет повысить пространственное разрешение, а также для фильтра fsc, который ухудшает сходимость метода. Ухудшение сходимости сглаживающего fsc выражается в переоценке скорости ук > у0, (т.е. Тд > 1), что также наблюдалось на мгновенных полях скорости, рассчитанных по реальным Р1У-изображениям при сглаживании поля предсказания на этапе коррекции смещения.
Таблица 1. Влияние фильтрации гауссовым фильтром на пространственное разрешение и сходимость адаптивного алгоритма
Тип фильтра ¡V ¡с ¡вс
Сглаживающий разр. 0 разр. - разр. — разр. 0
сход. + сход. 0 сход. 0 сход. --
Выделяющий разр. ++ разр. 0 разр. + разр. +
сход. + сход. 0 сход. — сход. —
На рис. 8, а представлены передаточные характеристики адаптивного алгоритма СШБ после третьей итерации для двух вариантов обработки с весовой гауссовой функцией и без нее в сравнении с результатом после первой итерации. Характеристики синтетических изображений приведены в части 2 статьи при тестировании пространственного разрешения. Видно, что для результата без оконной функции отличие передаточных характеристик больше, чем при применении окна. В обоих случаях адаптивный алгоритм дает более точное восстановление гармоник. Идеальная передаточная характеристика стремится к единице.
Дополнительно внимание было уделено вопросу сходимости СШБ алгоритма. В работе [38] показано, что трех-четырех итераций достаточно для достижения сходимости результата. На рис. 8, б показана зависимость случайной погрешности расчета от числа шагов итерационного алгоритма. Для данного теста концентрация частиц была снижена до уровня 0.018 частиц на 1 пикс2 и смоделирован шум, связанный с выходом частиц из плоскости изображения путем удаления 15% частиц со второго кадра изображения. Также на изображения был наложен случайный шум с гауссовым распределением на уровне 15 разрядов с дисперсией 5 разрядов, имитирующий тепловой шум матрицы цифровой камеры. По графику видно, что уровни погрешности для всех случаев на первой итерации равны — это связано с одинаковой передаточной характеристикой для всех методов на первой итерации. Далее, с увеличением номера итерации, погрешность БШБ частично увеличивается, а СШБ и СШБ — уменьшается. Для последних двух методов сходимость
Рис. 8. Передаточная характеристика итерационного алгоритма в сравнении со стандартным кросскорреляционным алгоритмом и КИХ-фильтром скользящего среднего (а); зависимость случайной погрешности адаптивных алгоритмов DWS, CWS, CWD от номера итерации для гармоники и = 0.6 (б)
наблюдается уже на четвертой итерации. Преимущество метода CWD видно по асимптотической величине случайной погрешности: 0.69 пикс вместо 0.94 и 0.99 пикс для CWS и DWS соответственно.
Еще одним ограничением стандартных кросскорреляционных алгоритмов, упоминавшимся во введении, является неспособность измерять поле скорости с существенным градиентом скорости внутри измерительной области. В работе [5] максимальный измеряемый градиент скорости для расчета стандартным алгоритмом составил 0.5 пикс/пикс (15 с-1). В монографии [15] приведена зависимость случайной погрешности от градиента скорости для аргумента в диапазоне 0.. .0.15 пикс. Заявленный авторами работы [10] максимальный градиент скорости, измеренный алгоритмом с деформацией области, составил 1.1 пикс/пикс (55.3 с-1) с уровнем случайной погрешности а = 0.25 пикс.
Тест на способность измерения в области с градиентом скорости проведен на синтетическом изображении частиц в турбулентном потоке в канале с исходным полем скорости, полученным DNS-расчетом (тест В, PIV Challenge 2003, Okamoto [23]). Величина градиента в области пограничного слоя составляет дvx/dy = 0.125 пикс/пикс для размера расчетной области F =16 пикс. На рис. 9, а показаны среднее значение скорости и турбулентные пульсации скорости, а также спектр пульсаций скорости (рис. 9, б) в горизонтальном сечении y = 32 пикс вблизи стенки, рассчитанные адаптивным и стандартным алгоритмами. Для сравнения на рис. 9, б дан спектр пульсаций исходного DNS-поля скорости. В работе [23] показано, что степень совпадения спектральных пульсаций скорости с исходными DNS-данными вблизи стенки коррелирует с отклонением рассчитанных пульсаций скорости от истинного значения.
На рис. 9, а пульсации скорости при расчете стандартным алгоритмом выше, чем при применении адаптивного алгоритма вследствие большей случайной погрешности вблизи стенки. На рис. 9, б результат алгоритма CWD находится ниже, что свидетельствует о меньшей погрешности восстановления поля скорости в пограничном слое. В работе [40] показано, что при превышении градиента скорости внутри расчетной области
Ух [пикс]
а б
Рис. 9. Профили x-компоненты средней скорости и пульсаций скорости в сечении x = 304 пикс (а), спектр пульсаций скорости для компоненты вдоль потока в сечении y = 32 пикс (б): сравниваются результаты для FFT-CC- и CWD-алгоритмов для тестового случая B, PIV Challenge 2003, Okamoto (турбулентный поток в канале) [23]
400 800 1200 400 800 1200 400 800 1200
х [пике] х [пике] х [пике]
а б в
Рис. 10. Пример расчета мгновенного поля скорости. Объект — вихрь за транспортным самолетом (тест B, PIV Challenge 2001 Kahler [39]): а — PIV-изображение частиц (два кадра объединены на одном изображении); результат стандартного (б) и адаптивного (в) кросскорреляционного алгоритмов (градацией серого цвета обозначен модуль скорости)
0.025 пикс/пикс (у < 100 пикс для теста на рис. 9) определяющую роль в погрешности начинает играть случайная ошибка и выбор типа интерполяции практически не влияет на ее уровень. Зависимость ошибки адаптивного алгоритма с деформацией области от величины шума на изображении исследована в работе [22]. Там же показано уменьшение влияния типа интерполяционной схемы на общую ошибку при увеличении случайного шума ПЗС-матрицы и шума, связанного с выходом частиц из плоскости лазерного ножа.
В заключение приведем пример расчета мгновенного поля скорости, представляющего собой сильный вихрь, образующийся за транспортным самолетом в режиме посадки [39] (рис. 10). Тест обработан стандартным и адаптивным (ОШВ) кросскорреляционны-ми алгоритмами, так что адаптивный алгоритм на последней итерации имеет один набор параметров расчета со стандартным алгоритмом. На рисунке видно, что поле (в) имеет больший динамический диапазон определения скорости (разрешение скорости в центре вихря) по сравнению с (б), а также меньшее количество ошибочных векторов.
Заключение
В работе приведен краткий обзор истории развития численных методов обработки изображений для расчета мгновенных полей скорости 2D PIV-эксперимента. Описан стандартный кросскорреляционный алгоритм расчета мгновенного поля скорости. Показаны основные ограничения стандартного алгоритма, связанные с динамическим диапазоном, точностью, эффектом peak-locking и пространственным разрешением.
Рассмотрены основные источники погрешности метода и представлены зависимости систематической и случайной погрешностей от величины смещения частиц в диапазоне от 0 до 2 пикс на синтетических изображениях частиц с однородным смещением без добавления шума.
Проанализирована составляющая систематической ошибки стандартного алгоритма, возникающая в результате неявного произведения корреляционных коэффициентов на весовую функцию при свертке двух расчетных областей равного размера. Приведены три источника эффекта peak-locking, связанного с искажением исходной формы плотности вероятности для компонент скорости, проанализированы основные подходы по его устра-
нению.
Пространственное разрешение кросскорреляционного алгоритма рассмотрено в терминах цифровой фильтрации данных и передаточных функций.
Продемонстрировано действие стандартного кросскорреляционного PIV-алгоритма как двумерного фильтра с конечной импульсной характеристикой, сходной с фильтром скользящего среднего по расчетной области [22, 27, 29]. Получена модельная передаточная характеристика стандартного алгоритма на наборе искусственных изображений частиц с вертикальным смещением в виде гармонической функции с различными частотами w = F/A = 0.05 ... 1.
В работе описана реализация адаптивного кросскорреляционного алгоритма, который позволяет улучшить качество получаемых данных по сравнению со стандартным алгоритмом, а именно: динамический диапазон измеряемых скоростей, точность и пространственное разрешение. В описании адаптивные алгоритмы классифицированы по порядку малости к поправке скорости на стадии компенсации смещения частиц [24].
Описаны ограничения адаптивного алгоритма по точности, уровню пространственного разрешения и сходимости, а также диапазону измеряемого градиента скорости. Исследовано поведение систематической и случайной ошибок в зависимости от однородного смещения частиц и типа интерполяционной схемы, применяемой при деформации изображения. Показано, что наиболее эффективной из рассмотренных интерполяционных схем по точности является B-spline третьего порядка, а эффективной с точки зрения времени и точности — билинейная интерполяция с коррекцией [21, 37]. Адаптивный алгоритм с деформацией области позволяет заметно уменьшить нежелательный эффект peak-locking, причем интерполяционная схема, дающая наименьшую систематическую ошибку алгоритма, одновременно дает наименьшую выраженность и по эффекту peak-locking.
Пространственное разрешение адаптивного алгоритма исследовано в терминах передаточной функции. Получена теоретическая передаточная характеристика для адаптивного алгоритма с учетом фильтрации. На модельной передаточной характеристике показано, что итерационный алгоритм без уточнения сетки может улучшать эффективное пространственное разрешение. Определено, что достаточно четырех итераций для сходимости реализованных методов CWS и CWD. Рассмотрена проблема расчета полей скорости для потоков со значительным градиентом скорости.
Дальнейшие исследования 2D PIV-алгоритмов будут направлены в сторону адаптивного пространственного разрешения в зависимости от концентрации частиц в локальной области потока, в том числе с совмещением подходов PIV и PTV. Необходимо также проведение исследований по расчету скорости на границе изображения и маскированных областей, где точность адаптивного алгоритма в настоящее время не превышает точности стандартного метода.
Список литературы
[1] хабахпашева е.м., Перепелица б.в. Поля скоростей и турбулентных пульсаций при малых добавках к воде высокомолекулярных веществ // ИФЖ. 1968. T. 14, № 4. C. 598.
[2] Adrian r.j. Scattering particle characteristics and their effect on pulsed laser measurements of fluid flow: speckle velocimetry vs. particle image velocimetry // Appl. Opt. 1984. Vol. 23. P. 1690-1691.
[3] Meynart R. Convective flow field measurement by speckle velocimetry // Rev. Phys. Appl. 1982. Vol. 17. P. 301-330.
[4] Adrian R.J., Yao C.S. Development of pulsed laser velocimetry (PLV) for measurement of fluid flow // Proc. of the 8th Biennial Symp. on Turbulence, Rolla, Missouri, Sept., 1984. P. 170-186.
[5] Willert C.E., Gharib M. Digital particle image velocimetry // Exp. Fluids. 1991. Vol. 10. P. 181-193.
[6] Алексеенко С.В., Бильский А.В., Маркович Д.М. Применение метода цифровой трас-серной визуализации для анализа турбулентных потоков с периодической составляющей // Приборы и техника эксперимента. 2004. №. 5. C. 145-153.
[7] Бильский А.В. Гидродинамическая структура осесимметричной импактной струи: Дис. ... канд. физ.-мат. наук. Новосибирск, 2006. 184 с.
[8] Scarano F., Riethmuller M.L. Iterative multigrid approach in PIV image processing with discrete offset // Exp. Fluids. 1999. Vol. 26. P. 513-523.
[9] Noguiera J., Lecuona A., Rodrigues P., Alfaro J.A., Acosta A. Limits on the resolution of correlation PIV iterative methods. Practical implementation and design of weighting functions // Exp. Fluids. 2005. Vol. 39. P. 314-321.
[10] Huang H.T., Feilder H.F., Wang J.J. Limitation and improvement of PIV, part II. Particle image distortion, a novel technique // Exp. Fluids. 1993. Vol. 15, P. 263-273.
[11] Westerweel J. Fundamentals of digital particle image velocimetry // Meas. Sci. Technol. 1997. Vol. 8. P. 1379-1392.
[12] Chen J., Katz J. Elimination of peak-locking error in PIV analysis using the correlation mapping method // Meas. Sci. Technol. 2005. Vol. 16. P. 1605-1618.
[13] Westerweel J., Scarano F. Universal outlier detection // Exp. Fluids. 2005. Vol. 39. P. 10961100.
[14] Baek S.J., Lee S.J. A new two-frame particle tracking algorithm using match probability // Exp. Fluids. 1996. Vol. 22. P. 23-32.
[15] Raffel M., Willert C., Kompenhans J. Particle Image Velocimetry. A Practical Guide. Springer, 1998.
[16] Heinz O., Ilyushin B., Markovich D. Application of a PDF based method for the experimental statistical processing of experimental data // Int. J. Heat and Fluid Flow. 2004. Vol. 25. P. 864-874.
[17] Christensen K.T. The influence of peak-locking errors on turbulence statistics computed from PIV ensembles // Exp. Fluids. 2004. Vol. 36. P. 484-497.
[18] Alekseenko S.V., Bilsky A.V., Dulin V.M., Markovich D.M., Tokarev M.P. Stereo PIV measurements of turbulence characteristics in axisymmetric jet // Submitted to Exp. Fluids. 2007.
[19] Noguiera J., Lecuona A., Rodrigues P.A. Identification of a new source of peak locking, analysis and its removal in conventional and super-resolution techniques // Exp. Fluids. 2001. Vol. 30. P. 309-316.
[20] Scarano F., Riethmüller M.L. Advances in iterative multigrid PIV image processing // Exp. Fluids. 2000. Suppl. P. S51-S60.
[21] Astarita T., Cardone G. Analysis of interpolation schemes for image deformation methods in PIV // Exp. Fluids. 2005. Vol. 38. P. 233-243.
[22] Astarita T. Analysis of interpolation schemes for image deformation methods in PIV: effect of noise on the accuracy and spatial resolution // Exp. Fluids. 2006. Vol. 40. P. 977-987.
[23] Stanislas M., Okamoto K., KAhler C.J. Main results of the second international PIV Challenge // Exp. Fluids. 2005. Vol. 39. P. 170-191.
[24] Scarano F. Iterative image deformation methods in PIV. Review article // Meas. Sci. Technol. 2002. Vol. 13. P. R1-R19.
[25] Sügii Y., Nishio S., Oküno T., Okamoto K. A highly accurate iterative PIV technique using gradient method // Meas. Sci. Technol. 2000. Vol. 11. P. 1666-1673.
[26] Bing P., Hüi-min X., Bo-qin X., Fü-long D. Performance of sub-pixel registration algorithms in digital image correlation // Meas. Sci. Technol. 2006. Vol. 17. P. 1615-1621.
[27] FoüCAüt J.M., Carlier J., Stanislas M. PIV optimization for the study of turbulent flow using spectral analysis // Meas. Sci. Technol. 2004. Vol. 15. P. 1046-1058.
[28] Kean R.D., Adrian R.J. Optimization of particle image velocimeters, Part 1: double pulsed system // Meas. Sci. Technol. 1990. Vol. 1. P. 1202-1215.
[29] Nogüeira J., Lecuona A., Rodriguez P.A. Local field correction PIV: on the increase of accuracy of digital PIV systems // Exp. Fluids. 1999. Vol. 27. P. 107-116.
[30] Susset A., Most J.M., Honore D. A novel architecture for a super-resolution PIV algorithm developed for the improvement of the resolution of large velocity gradient measurements // Exp. Fluids. 2006. Vol. 40. P. 70-79.
[31] Theunissen R., Scarano F., Riethmüller M.L. An adaptive sampling and windowing interrogation method in PIV // Meas. Sci. Technol. 2007. Vol. 18. P. 275-287.
[32] Stitou A., Riethmüller M.L. Extention of PIV to super resolution using PTV // Meas. Sci. Technol. 2001. Vol. 12. P. 1398-1403.
[33] Stanislas M., Okamoto K., KAhler C.J., Westerweel J. Main results of the third international PIV challenge // To appear in Exp. Fluids. 2007.
[34] Quenot G.M., Pakleza J., Kowalewsky T.A. Particle image velocimetry with optical flow // Exp. Fluids. 1998. Vol. 25. P. 177-189.
[35] Wereley S.T., Meinhart C.D. Second-order accurate particle image velocimetry // Exp. Fliuds. 2001. Vol. 31. P. 258-268.
[36] АхмЕТБЕков Е.К., Бильский А.В., Ложкин Ю.А., Маркович Д.М., Токарев М.П., Тюрюшкин А.Н. Система управления экспериментом и обработки данных, полученных методами цифровой трассерной визуализации (ActualFlow) // Вычисл. методы и программирование. 2006. Т. 7. C. 79-85.
[37] Piirto m., Eloranta h., Saarenrinne p., Karvinen r. A comparative study of five different PIV interrogation algorithms // Exp. Fluids. 2005. Vol. 39. P. 571-588.
[38] Scarano f. On the stability of iterative PIV image interrogation methods // Proc. 12th Int. Symp. on Application of Laser Techniques to Fluid Mechanics Lisbon, Portugal, 12-15 July, 2004.
[39] Stanislas m., Okamoto k., KAhler c. Main results of the First International PIV Challenge // Meas. Sci. Technol. 2003. Vol. 14. P. R63-R89.
[40] Kim b.j., Sung h.j. A further assessment of interpolation schemes for window deformation in PIV // Exp. Fluids. 2006. Vol. 41. P. 499-511.
Поступила в редакцию 2 марта 2007 г.
