CC BY
40
Вычислительные технологии
Том 15, № 4, 2010
Корреляционная коррекция в методе слежения за частицами в потоках*
Е.К. Ахметбеков, Д.М. Маркович, М.П. Токарев Институт теплофизики СО РАН, Новосибирск, Россия e-mail: yerbol@itp .nsc.ru, mtokarev@itp .nsc.ru
Исследован алгоритм построения полей скорости в методе слежения за частицами с применением коррекции, основанной на корреляционной обработке изображений частиц. На основе искусственно сгенерированных изображений частиц изучено влияние диаметра и концентрации частиц, а также величины их смещения на погрешность при расчете поля скорости. Алгоритм коррекции апробирован при обработке экспериментально полученных данных. Исследовано влияние алгоритма коррекции на пространственное разрешение метода PTV. Достоинством алгоритма является то, что с его использованием коррекция может осуществляться после применения как стандартного PTV, так и гибридных подходов PIV/PTV. В работе предложен метод отсева ошибочных векторов скорости в нерегулярном векторном поле, основанный на фильтре скользящего среднего.
Ключевые слова: PTV, PIV, корреляция, пространственное разрешение, им-пактная струя.
Введение
В современных методах цифровой траееерной визуализации применяются две основные группы подходов для обработки изображений частиц, полученных с короткой экспозицией в два момента времени: кросскорреляционный анализ с регулярным разбиением полного изображения на элементарные измерительные области (Particle Image Veloeimetry, PIV) и методы слежения за каждой отдельно взятой частицей (Particle Tracking Veloeimetry, PTV).
Корреляционный подход, реализуемый в методе PIV, представляет собой эффективный и удобный способ расчета вектора смещения. Высокая субпиксельная точность определения максимума корреляционной функции и соответственно смещения частиц, а также наличие надежных способов выделения полезного сигнала (высокое отношение сигнал/шум) являются важными достоинствами метода PIV [1]. Однако применение данного метода накладывает определенные ограничения на условия проведения эксперимента, к чему относится, например, необходимость обеспечения некоторого минимального количества траееерных частиц в элементарной измерительной области. Это требование приводит к невозможности радикального повышения пространственного разрешения метода даже с использованием частичного перекрытия элементарных расчетных областей (наиболее часто применяемое перекрытие составляет 50 % и позволяет
* Работа проводилась в рамках реализации ФЦП "Научные и научно-педагогические кадры инновационной России" на 2009-2013 годы.
увеличить количество векторов скорости в 4 раза, однако не приводит к повышению точности при определении ряда значимых характеристик турбулентных потоков). Для достижения большего пространственного разрешения в методе PIV необходимо уменьшать размеры исследуемой области, что не позволяет измерять мгновенные поля скорости в потоке в целом.
Еще одно ограничение традиционного PIV — это невозможность эффективно обрабатывать изображения с существенно неравномерным распределением частиц. Значительные усилия в этом направлении были предприняты авторами работы [2], где предложено использовать информацию об отдельных частицах для выбора положения точек измерения скорости и оптимального значения размера окна в зависимости от локальной концентрации частиц.
Одним из преимуществ метода PTV является большое пространственное разрешение полей скорости, при этом локальность полученного вектора скорости имеет максимально возможное значение. Кроме того, нет необходимости в дополнительных мерах при обработке изображений с неравномерным распределением частиц. Метод PTV включает следующие стадии: определение положения центров частиц, нахождение пар частиц и расчет векторов смещения, отсев ошибочных векторов. Основные проблемы при применении метода PTV связаны с распознаванием неярких и перекрытых образов частиц, а также с недостаточной точностью определения положений центров частиц. Первая проблема ведет к уменьшению количества рассчитанных векторов скорости и соответственно к ухудшению пространственного разрешения. Высокая точность определения положения центра частицы и, в конечном счете, смещения частицы является необходимым условием для расчета и дальнейшего анализа статистических моментов скорости высших порядков в турбулентных потоках.
Изложенный в работе [3] метод корреляции с гауссовой маской PMC (Particle Mask Correlation) в сравнении с применяемым обычно методом простой бинаризации по порогу интенсивности позволяет значительно улучшить способность метода PTV к идентификации затененных частиц и частиц с перекрытыми образами,
В /тухкадровой реализации метода PTV, после определения положения центров частиц на первом и втором кадрах для последующего расчета вектора смещения необходимо каждой частице на первом кадре найти пару во втором кадре, т, е, образ той же самой частицы, переместившейся за заданный промежуток времени в другую точку пространства. Простейшая реализация метода поиска пар заключается в нахождении ближайшей частицы на втором кадре, которой присваивается метка пары текущей частицы на первом кадре. Этот стандартный двух кадровый метод PTV [4] имеет следующее ограничение: смещение частицы на двух последовательных кадрах не должно превышать половины среднего расстояния между частицами на одном кадре изображения, Фактически это является ограничением на пространственное разрешение или на динамический диапазон измеряемых скоростей в потоке, поскольку на плотность изображений частиц в кадре накладывается ограничение сверху, связанное с максимальным смещением частицы. Для другой модификации метода PTV [5, 6] необходимо использование четырех кадров изображений частиц, полученных с одинаковыми временными интервалами между ними. Треки по четырем изображениям идентифицируются в соответствии с критерием расходимости. Трек с минимальной расходимостью считается верным для текущей частицы. Метод легко реализуем, но не применим при существенных значениях третьей, перпендикулярной лазерному ножу, компоненты скорости, Основная идея еще одной модификации метода PTV (Enhanced PTV—EPTV [7])
заключается в использовании информации о размерах трассеров для поиска пар частиц, Преимуществом метода является то, что он теоретически может быть применим к потокам любого типа, в том числе с достаточно большими градиентами скорости, К его недостаткам относится сложность выполнения некоторых условий эксперимента, а именно, обеспечения однородности интенсивности излучения по толщине лазерного ножа. Кроме того, достаточно широкое распространение получили комбинированные PIV/PTV-подходы, Метод Super Resolution PIV (SE PIV) [8-12] основывается на корреляционной обработке изображений частиц на первом этапе и уточнении вектора смещения каждой отдельной частицы на втором этапе,
В данной работе описан комбинированный подход, в котором применены корреляционные методы и методы слежения за частицей. Это позволяет избежать большинства недостатков каждого из методов, сохраняя их достоинства — высокое пространственное разрешение и точность. Основная идея подхода заключается в использовании коррекции, основанной на корреляции изображений индивидуальной частицы в обоих кадрах после нахождения пары и построения вектора скорости. Такой метод коррекции называется Individual Particle Correlation (IPC) [13], Применение этого метода позволяет уменьшить размер корреляционного окна до размеров изображения отдельной частицы. Стандартный PIV-метод не позволяет успешно определять смещения частиц по расчетным областям малого размера, так как нельзя гарантировать наличие и достаточную концентрацию частиц в данной области. Предварительно рассчитанное нерегулярное поле PTV дает возможность использовать изображения отдельных найденных частиц. Подход IPC позволяет увеличить точность расчета смещения частиц в стандартном PTV-методе,
1. Описание алгоритма
Процесс обработки изображений частиц методом PTV с коррекцией состоит из четырех стадий: 1 — определение положения центров частиц, 2 — нахождение пар частиц, 3 — корреляционная коррекция смещения частиц, 4 — отсев ошибочных векторов. Положения центров частиц определяются при помощи корреляционного поиска Particle Mask Correlation, Метод PMC позволяет распознать частицу, основываясь на форме ее образа (форме распределения интенсивности), которая может быть аппроксимирована двумерным распределением Гаусса, при этом абсолютное значение максимума интенсивности не имеет значения. Нормализованный коэффициент корреляции области изображения вблизи текущего пикселя с предварительно построенной гауссовой маской рассчитывается для каждого пикселя. Таким образом, получается поле корреляций, для которого проводится пороговая бинаризация, В полученных доменах ведется поиск максимума, центр его и соответствует положению частицы. Частицы с частично перекрытыми образами на изображении различаются автоматически,
1.1. Нахождение пар частиц
В данной работе для поиска пар частиц использовался релаксационный метод. Этот метод, описанный в [13], основан на оценке вероятности каждой из возможных пар. Для каждой частицы на первом кадре ведется поиск пары на втором кадре в пределах окружности с радиусом, равным максимально возможному перемещению Tm, Каждой из возможных пар на втором кадре присваивается одна и та же начальная вероятность, кроме того рассматривается отдельная вероятность для состояния поте-
ри пары, которая также равна начальной вероятности остальных вариантов. Сумма вероятностей всех исходов равна единице. Вероятность каждого исхода итерационно переечитываетея. При оценке вероятности текущей пары рассматриваемой частицы у соседних частиц ведется поиск пар, совершающих близкое по величине смещение, и в случае нахождения их вероятности дают вклад в вероятность текущей пары. Этот метод позволяет определять пары частиц без дополнительных предварительных расчетов и может быть применен для обработки потоков с высокой концентрацией частиц (вплоть до 0,05 частиц/пике2). Однако описанный метод требует значительных вычислительных ресурсов. Предварительно обработанные методом PIV данные могут в существенной мере уменьшить количество частиц — кандидатов в пары и увеличить производительность. Ниже предлагаемый подход будет называться гибридным релаксационным, Оригинальный метод имеет три параметра. Поиск кандидатов в пары для каждой частицы в первом кадре ведется среди частиц во втором кадре в пределах окружности с радиусом Tm и с центром, совпадающим с положением частицы в первом кадре. Второй параметр Tn определяет радиус окрестности текущей частицы, в которой соседние частицы показывают схожее движение с небольшими отклонениями. Величина такого отклонения является третьим параметром Tq.
Предварительно рассчитанные методом PIV поля скорости на регулярной сетке могут быть использованы несколькими способами. Одним из таких способов является поиск ближайшей частицы к положению, предсказанному вектором смещения, который, в свою очередь, получен методом PIV, Этот метод называется SR PIV (SuperResolution PIV), Такой подход является расширением метода PIV и обладает теми же недостатками, что и PIV, — сложностью в обработке при наличии больших градиентов скорости и невозможностью разрешения мелкомасштабных структур, В нашей работе использовался релаксационный алгоритм слежения за частицей. Основной задачей было уменьшить число кандидатов, отбросив очевидно неверные пары. Таким образом, гибридный метод имеет три дополнительных параметра, которые используются вместо параметра Tm. Новые параметры необходимы для определения области поиска пары на втором изображении. Основной параметр — это отношение радиуса области поиска к длине предсказанного вектора, который был интерполирован в местоположение центра текущей частицы из PIV-данных. Оставшиеся два параметра — это минимальное и максимальное абсолютное значение этого радиуса,
1.2. Коррекция методом корреляции изображения частицы
Важными характеристиками любого метода определения положения частиц являются коэффициент распознавания частиц и субпиксельная точность определения их положения, Метод корреляции с гауссовой маской обладает хорошей способностью к распознаванию частиц.
Для достижения высокой субпиксельной точности при определении положения центров частиц, в том числе и методом РМС, чаще всего используется трехточечная и двумерная гауссова интерполяция. Однако результаты тестирования на искусственных изображениях показали, что точность определения положения частиц, полученного таким образом, недостаточна и данные методы не позволяют избежать ошибки, связанной с эффектом группировки рассчитываемых величин (положения центра частицы и вектора смещения) у целых значений пикселей (peak-locking). Это обусловлено систематической ошибкой, связанной с дискретизацией изображения.
Рис. 1. Схема работы метода коррекции
Дня уменьшения эффекта peak-locking в данной работе используется алгоритм коррекции векторов смещения методом IPC, Коррекция применяется поело нахождения вектора смещения и до процедуры отсева ошибочных векторов. Цель алгоритма коррекции — увеличение точности результирующего вектора смещения. Таким образом, совместное применение метода корреляции с гауссовой маской и алгоритма коррекции позволяет достичь высокой субникселыюй точности и значительной величины коэффициента восстановления.
Дня применения коррекции но изолированным изображениям одной и той же частицы па обоих кадрах рассчитывается поло кросекорреляции (рис, 1), Положение максимума па корреляционной плоскости определяет значение вектора коррекции.
Перед применением корреляционной коррекции малые области па изображениях первого и второго кадров восстанавливаются путем интерполяции в новых системах координат. Началом координат в них являются дробные положения центров частиц па исходном изображении. Интерполяция проводится дня того, чтобы сместить положение корреляционного максимума в последующей процедуре корреляционной коррекции ближе к началу координат, где субникеелытя оценка положения максимума имеет меньшую ошибку. Иными словами, чем меньше значение вектора коррекции, тем меньше величина систематической ошибки. Дня интерполяции изображений использовалась двумерная кубическая В-силайн схема, которая но данным тестирования |1| имеет наилучшее соотношение точность/время расчета применительно к изображениям частиц. Коррекцию IPC можно применять итерационно несколько раз, в каждом случае используя уточненное па предыдущем шаге значение вектора смещения.
Дня определения положения максимума па корреляционной плоскости с субник-селыюй точностью применяется трехточечпая гауссова интерполяция, имеющая такие же ограничения субникеелыюй точности, как и при определении субникеелыюго положения центра частицы, однако подобная процедура проводится один раз дня каждого вектора смещения, тогда как без применения коррекции она проводится дважды. Кроме того, с каждой итерацией величина ошибки должна уменьшаться (в продело стремиться к пуню), поскольку если значение вектора коррекции близко к целому (в пашем случае — к пуню), то ошибки методов интерполяции значительно снижаются.
1.3. Отсев ошибочных векторов
Автоматическая обработка изображений неизбежно приводит к появлению ошибочных векторов смещения. Такие векторы могут появиться из-за низкого качества изображения частиц, наличия больших градиентов скорости, существенных значений компоненты скорости, перпендикулярной плоскости лазерного ножа, и других эффектов. Поиск и удаление таких векторов до последующей обработки является важным шагом, поскольку их влияние на итоговые данные существенно, особенно на старшие статистические моменты в турбулентных течениях [15], и связанные с этим погрешности полностью устранить весьма сложно [15, 16]
1.3.1. Алгоритм скользящего среднего для поля скорости на регулярной сетке
Этот метод основан на предположении о локальной гладкости векторного поля. Значение каждого вектора сравнивается со значением соседних векторов. Вектор лг(х,у) рассматривается как ошибочный, если выполняется соотношение ¡V(х,у) — у(х,у)| < к, где у(х,у) — средняя локальная скорость. Пороговое значение к рассчитывается следующим образом: к = а тахх,у ¡V(х, у) — у(х,у)|, где 0 < а < 1 [17].
1.3.2. Алгоритм скользящего среднего для поля скорости на нерегулярной сетке
Сложность применения алгоритма скользящего среднего к нерегулярным данным состоит в том, что в этом случае невозможно подобрать оптимальный размер области осреднения, которая, безусловно, будет содержать соседние векторы. Если рассматривать достаточно большую область, заведомо включающую соседние векторы, метод будет работать неточно, так как среднее значение скорости может сильно отличаться от локального значения. Для решения этой проблемы было предложено рассматривать большую область с учетом весовой функции по Гауссу. Таким образом, оригинальные выражения были преобразованы к следующему виду:
где N(г0) — множество векторов из области осреднения. Описанный метод может быть применен несколько раз итерационно. Каждая итерация удаляет векторы, наиболее отличные от среднего вектора, рассчитанного с весовой функцией.
Использование двух дополнительных параметров позволяет избежать другой проблемы, которая проявляется, когда поле скорости уже не содержит ошибочных векторов. Особенность этого метода в том, что он находит ошибочные векторы всегда, помечая наиболее отличимый от среднего вектор скорости как неверный. Дополнительные параметры ограничивают нижнюю границу отличия, достаточного для определения вектора как неверного. Первый параметр — абсолютная нижняя граница, второй — нижняя граница, рассчитываемая относительно средней длины вектора.
пеМ (го)
пеМ (го)
к = а тах V — v(r0)|
гаем (го)
2. Тестирование алгоритма на синтетических изображениях
Дня тестирования алгоритма РТУ с коррекцией 1РС были созданы синтетические изображения. Сравнивались результаты обработки методами РТУ, стандартным ЕЕТ-СС Р1У и итерационным Р1У (С\¥Б В-БрНпе Р1У) (рис, 2, а). При этом учитывались данные |16, 181,
Отдельные тесты использовались дня оценки различных свойств алгоритма коррекции: определение центра частицы, нахождение пар частиц, систематическая и случайная ошибки. Размер синтетических изображений — 1024 х 1024 пике, глубина оттенков серого составляла 16 бит. Тепловой шум цифровой матрицы не моделировался. Компонента скорости, перпендикулярная плоскости изображения, была равна нулю. Во всех тестах, кроме теста па нахождение пар частиц (рис, 2,6), количество частиц па изображении составляло 20000 (0,019 частиц/пике2).
Тестирование способности алгоритма к определению центров частиц проводилось на искусственных изображениях с частицами гауссовой формы и диаметрами в диапазоне от 1 до 7 пике с шагом 0,1 пике. Диаметр ¿р образа частицы на изображении определялся шириной гауссова распределения по уровню интенсивности I = 10 ехр-2, оде 10 — интенсивность в центре распределения, согласно |16|, Равномерное смещение частиц составляло й = 0.5 пике, так как при этом значении алгоритм имеет максимальную погрешность.
На рис, 2, а представлена зависимость случайной ошибки определения смещения от диаметра частиц ¿р. Результаты, полученные по стандартным Р1У- и РТУ-алгоритмам, имеют минимумы при йр = 2.5 и 1,9 пике соответственно, что хорошо согласуется с данными [16, 18, 19], Кривая РТУ имеет второй минимум. Согласно [18], с увеличением ¿р наблюдаются и другие минимумы. Применение алгоритма коррекции позволяет увеличить диапазон диаметров частиц, при котором уровень шума остается относительно низким.
График, соответствующий РТУ с коррекцией, повторяет форму графика дня адаптивного Р1У-а.нгоритма. В диапазоне 1-1.5 пике оба РТУ-метода показывают неустойчивое поведение, что может быть связано с частотной фильтрацией в ходе определения
Рис. 2. Случайная ошибка в зависимости от диаметра образа частицы (а) и количество верно определенных векторов для различных концентраций частиц на синтетическом изображении (б)
положения частиц методом РМС. В ходе тестирования способности метода к нахождению пар рассчитывалось количество полученных пар Njmeas для различных концентраций частиц Nj (см, рис, 2, б). Однородное смещение частиц в исследуемых изображениях составляло d = 0.5 пике. При малых концентрациях число найденных пар совпадает с числом сгенерированных частиц. Начиная с 4000 частиц на изображении (концентрация 0,0076 частиц/пикс2) появляются неопределенные пары. При концентрации в 40000 (0,038 частиц/пикс2) определяется 70% частиц. Таким образом, с увеличением концентрации частиц эффективность нахождения нар снижается.
Оценка точности алгоритмов проводилась при помощи синтетических изображений с равномерным горизонтальным смещением в продолах от 0 до 2 пике и шагом 0,05 пике. На рис, 3 представлены средние значения систематической и случайной ошибок дня РТУ- и PIV-методов. Графики ошибок имеют периодическую форму, зависящую от величины смещения частиц |20|, Ошибка уменьшается дня РТУ с коррекцией, при этом оставаясь выше уровня адаптивного PIV, Это может быть объяснено небольшим размером корреляционного окна, используемого в РТУ с коррекцией, но сравнению с размером окна в случае PIV, Кроме того, па рис, 3 приведены зависимости дня РТУ с коррекцией с размером корреляционных окоп 3 и 5 пике.
Следующий тест был проведен дня подтверждения того, что пространственное разрешение метода РТУ после применения корреляционной коррекции не ухудшается. Искусственные изображения были сгенерированы с одномерным синусоидальным смещением и различными волновыми числами дня вертикальной компоненты скорости вдоль оси x v = 3sin(2nwx/5), оде ш = F/А = 5/А = 0.005 — 1 (шаг 0,05) — безразмерное волновое число, нормированное на размер области F, Оптимальный размер окна для алгоритма коррекции F = 5 пике. Передаточная функция T(ш) = v(w)/v0(w) была рассчитана во всех точках поля скорости для текущего волнового числа ш [21]:
N
8 (ш)2
N
J^(vQi - T(ш^0г)2; T(ш) = 1
i= 1
\
N
8(ш)2/х; vQi, T(ш) < 1
i=1
На рис, 4, а представлены несколько передаточных функций дня текущего теста: прямоугольный фильтр скользящего среднего (top-hat moving average filter — THMA),
0.012
0.006
S
s с
cd
-0.006
-0.012
0.01
8 к в
0.001
0.0001 2 0
d, пике
Рис. 3. Средние значения систематической (а) и случайной (б) ошибок как функции от смещения частиц
(О ю
Рис. 4. Пространственное разрешение но сравнению с прямоугольным фильтром скользящих) ереднмх) (ТНМА) (а) и полная ошибка (б) методов Р1У и РТУ для одномерншх) еинусоидаль-но!'о смещения в зависимости от наложенной пространственной частоты
и графики для стандартного РТУ и РТУ с коррекцией при размере окна 3, 5, 7, 9, 11, 13, 15 пике. Результаты ТНМА приведены дня сравнения, так как они схожи с данными кроеекорреляциошюго анализа PIV, Передаточные функции дня РТУ и РТУ с окном коррекции в 3 и 5 пике различается незначительно. Кривые с другими размерами окон показывают меньшее пространственное разрешение, так как имеют меньшую частотную границу отсечения. Из рис. 3, б" следует, что среднеквадратичное значение ошибки дня окна размером в 5 пике ниже но сравнению с окном в 3 пике. Таким образом, оптимальное значение размера корреляционного окна было выбрано равным 5 пике. Для ш < 0.5 на рис. 4, a, РТУ с коррекцией показывает несколько лучшие результаты но сравнению со стандартным РТУ. Рисунок 4, б подтверждает этот вывод. Таким образом, алгоритм корреляционной коррекции с небольшим размером корреляционного окна практически не уменьшает локальность метода РТУ в изображениях со стандартной дня PIV концентрацией частиц.
3. Эксперимент по исследованию турбулентной осесимметричной импактной струи
3.1. Параметры проведения эксперимента
В ходе эксперимента определялись характеристики затопленной осесимметричной импактной струи. Полученные изображения частиц использовались дня тестирования метода IPC. Течение было организовано в прямоугольной емкости из органического стекла высотой 40 см, шириной 20 см и длиной 20 см (рис. 5, а). Во время эксперимента температура воды поддерживалась постоянной 32 ± 0.5 °С. Режим турбулентного струйного течения характеризовался числом Re = 19 000, которое определялось по средней расходной скорости течения U0 = 1.028 м/с и диаметру профилированного сопла ^nozzle = 15 мм. Измерения проводились с помощью системы ПОЛИС под управлением программного обеспечения ActnalFlow |22|, Расстояние между соплом и импактной поверхностью составляло три диаметра сопла. Двойной импульсный Xd:YAG-,na3ep NewWave Solo III с энергией в импульсе, равной 50 мДж, использовался дня освещения
х, мм
Рис. 5. Схема эксперимента (а) и типичное мгновенное ноле скорости для экспериментов № 1 и 2 с размером исследуемой области соответственно 35 х 25 мм и 6.5 х 5 мм (б); жирные линии — горизонтальное и вертикальное сечения, вдоль которых проводилось наблюдение
Импактная поверхность
. м/с
потока. Лазерный нож формировался при помощи цилиндрической линзы. Толщина лазерного ножа в области исследования составляла 0.8 мм. Для регистрации изображений частиц применялась цифровая ПЗС-камера Видеоскан 205Д с разрешением 1360 х 1024 пике. Для регулировки положения камеры использовалось прецизионное координатное устройство. Были проведены два эксперимента с разным размером исследуемой области: №1 — 35 х 25 мм и .№ 2 — 6.5 х 5 мм с высоким пространственным разрешением. В экспериментах использовались полиамидные частицы размером 20 и 5 мкм соответственно. Плотность частиц равнялась 1.05 кг/м3, В обоих случаях верхняя граница области измерения совпадала с имнактпой поверхностью. Для получения достоверной статистической информации в каждом эксперименте зарегистрировано и обработано 3000 пар изображений.
3.2. Параметры обработки данных
Экспериментальные данные были обработаны гибридным релаксационным методом РТУ как с коррекцией, так и без коррекции IPC. Результаты РТУ-метода сравнивались с результатами итерационного PIV-алгоритма с нереразбиением сетки между итерациями, Перед поиском смещения частиц проводилась предварительная обработка изобра-
жений, 11мпак I нам поверхность в верхней части изображения маскировалась черным цветом. Далее осредненное по ансамблю изображение частиц вычиталось из каждого исходного изображения. Ниже представлены параметры, использованные в первом эксперименте с большим размером измерительной области. Число частиц на изображении составляло около 20000, что соответствовало плотности N = 0.02 частиц/пикс2.
Гибридная схема релаксационного метода была применена со следующими параметрами: Тп = 15 пике, Т = 3 пике, три итерации. Радиус предсказанной области находился в интервале [5; 8] пике и определялся относительным коэффициентом Я = 0.5 (т.е. радиус области поиска равнялся половине длины предсказанного вектора). Алгоритм коррекции применялся с использованием корреляционного окна размером 5 пике с четырьмя итерациями. Последующий отсев ошибочных векторов проводился при следующих значениях параметров: радиус осреднения Яаугё = 20 пике, ширина гауссовой весовой функции = 10 пике, а = 0.99, две итерации, абсолютная и относитель-
ная нижние границы отсева равны соответственно 1 пике и 10%, Статистические характеристики рассчитывались на базе оригинальных РТУ-данных без использования промежуточной интерполяции на регулярную сетку.
Конечный размер корреляционного окна для алгоритма Р1У был выбран Я =16 пике, расстояние между узлами сетки равнялось 8 пике (50%-е перекрытие областей). Начальный размер элементарной области Я = 32 пике. Две итерации адаптивного Р1У-ал-горптма с непрерывным смещением и деформацией областей проводились с начальным размером элементарной области. После переразбиения сетки выполнялись еще две итерации с конечным размером областей. Между итерациями проводился отсев ошибочных векторов в соответствии с процедурой, описанной в [23], Все идентифицированные ошибочные векторы заменялись на интерполированные. После заключительной итерации замена векторов не проводилась. При расчете кросскорреляции для увеличения эффективного пространственного разрешения элементарные области изображения рассматривались с учетом весовой функции — двумерного распределения Гаусса с полушириной 0.5Я
двумерный кубический В-сплайн,
Второй эксперимент с зуммированием исследуемой области и с числом частиц на изображении, равным 7500 (плотность частиц N2 = 0.007 чаетиц/пикс2), был обработан со следующими параметрами. Для гибридной релаксационной схемы РТУ минимальный и максимальный радиусы предсказанной области равнялись соответственно 2 и 8 пике. Поиск ошибочных векторов проводился с параметрами: ширина гауссовой весовой функции = 15 пике, а = 0.99, одна итерация, абсолютная нижняя грани-
ца отсева равна 0 пике, относительная граница составляла 20%, Остальные параметры аналогичны используемым в первом эксперименте за исключением того, что размеры окон были в два раза больше, т. е. Я = 64 и Я = 32 пике для начального и конечного размера расчетной области соответственно,
3.3. Результаты эксперимента
На рисунке 5, б" приведены примеры мгновенных распределений скорости для обоих экспериментов, Профили компонент средней скорости и компонент кинетической энергии турбулентности (КЭТ), приведенные на рис, 6-8, построены вдоль сечений, показанных жирными линиями на рис, 5, б.
На рис, 6, а представлена средняя скорость в горизонтальном сечении (у = 5 мм). Ось симметрии струи совпадает с осью абсцисс. Максимальное значение скорости па выходе сопла равнялось Vy = 1.028 м/с. Распределения средней продольной скорости, полученные тремя методами (Р1У, РТУ, РТУ с коррекцией), с хорошей точностью совпадают друг с другом. Исключением является область в диапазоне т/О от 0,2 до 0,45 в слое смешения, где средняя скорость, полученная методом Р1У, меньше скорости, рассчитанной методом РТУ, примерно па 4%, Это, очевидно, связано с влиянием пространственного осреднения в методе Р1У, которое в областях с большими градиентами скорости может давать дополнительную ошибку.
На рис, 6, б" приведены распределения средней компоненты скорости Ух в вертикальном сечении с большей детализацией вблизи имнактпой поверхности (эксперимент У2 2), Статистика дня РТУ-даппых была рассчитана па сетке с элементарной ячейкой размером в 16 пике без перекрытия. Шаг сетки для Р1У равен 16 пике (Я = 32 пике с 50%-м перекрытием). Разница между кривыми Р1У и РТУ на рис, 6, б для у' = 0.2 мм составляла 3%, что, как и в эксперименте .У 1, может быть связано с недооценкой значений скорости, полученных методом Р1У, за счет высоких градиентов скорости в данной области.
Распределения компонент кинетической энергии турбулентности у^2 и у'у2 в горизон-
у=5
различия между этими величинами, рассчитанными различными методами, наблюдаются в центре сдвигового слоя, где уровень турбулентных пульсаций достигает максимума, Видно, что значения КЭТ, полученные с помощью метода РТУ, находятся выше данных, полученных методом Р1У, Причиной такого различия может быть как ошибка используемого алгоритма РТУ/Р1У, так и эффект пространственного разрешения, который определяется общим размером измерительной области, размером элементарной расчетной ячейки дня Р1У или характерным расстоянием между частицами на изображении дня РТУ |24|,
Погрешности алгоритма в первую очередь влияют на измеренный уровень турбулентности в эксперименте .У 2 с большим пространственным разрешением (уменьшеп-
а б
Рис. 6. Средняя скорость в горизонтальном сечении (у = 5 мм) для осевой компоненты скорости, рассчитанная различными методами, эксперимент № 1 (а) и средняя скорость в вертикальном сечении (ж = 3.4 мм) для компоненты скорости вдоль импактной поверхности, эксперимент № 2 с большим пространственным разрешением (б)
х/П хШ
Рис. 7. Радиальная (а) и осевая (б) компоненты КЭТ в горизонтальном сечении (у = 5 мм), рассчитанные различными методами; эксперимент № 1
пая область измерения). Дня эксперимента № 1 с меньшим пространственным разрешением (большая область измерения) ошибка алгоритма и эффекты пространственного разрешения оказывают примерно равнозначное влияние на измеренные значения ну.ньсадионных характеристик. Измеренные значения турбулентных характеристик (статистические моменты высших порядков, диссипация КЭТ и др.) при невысоком пространственном разрешении, как правило, существенно недооценены в методах цифровой трассерной визуализации ввиду того, что сетка с конечным значением измерительной ячейки играет роль фильтра низкой частоты с передаточной функцией, ослабляющей вклад высоких частот. Таким образом, вклад высокочастотной части турбулентного спектра при расчете характеристик оказывается недостаточно учтенным. Повышение пространственного разрешения (уменьшение размера элементарной расчетной ячейки) вплоть до величин нескольких колмогоровских масштабов |25| может существенно уменьшить эту погрешность, особенно дня статистических моментов первого и второго порядка. С другой стороны, алгоритмическая ошибка приводит, как правило, к завышению оцененных статистических моментов турбулентных пульсаций скорости.
На рис. 8 приведены распределения КЭТ дня вертикального сечения, рассчитанные но мгновенным нолям скорости, полученным в экспериментах № 1 и 2 с использованием методов РТУ, РТУ с коррекцией и адаптивного Р1У. Видно, что результаты РТУ дня эксперимента № 1 дают значение компоненты КЭТ < у'х2 >, близкое к полученному методом Р1У в эксперименте У2 2 с большим пространственным разрешением. Таким образом, можно сделать вывод о том, что РТУ-алгоритмы позволяют учесть вклад мелкомасштабных флуктуаций скорости даже дня эксперимента с низким пространственным разрешением в отличие от адаптивного Р1У, который недооценивает уровень флуктуаций в тех же условиях. Анализ зависимостей дня КЭТ вдали от имнактной поверхности (у' > 5 мм), где значения ее компонент так же, как и осредненной скорости, уменьшаются, асимптотически приближаясь к пуню, показал, что данные, полученные методом РТУ с коррекцией, близки к результатам адаптивного РГУ-метода, обеспечивающего максимальную точность в зоне, где отсутствуют значительные градиенты скорости, Данные, полученные но стандартному РТУ-алгоритму, демонстрируют больший уровень ошибки, особенно дня поперечной компоненты КЭТ (рис, 8,6),
4 б
у', мм
4 6
у', мм
Рис. 8. Распределения радиальной (а) и осевой (б) компонент КЭТ в вертикальном сечении вблизи импактпой поверхности (ж = 3.4 мм), рассчитанные различными методами; PIV, PTV, РТУ IPC эксперимент № 1, PIV (увел.) № 2
0.02 "Jj 0.015 0.01 0.005
0.005
0.025
0.015
Таким образом, описанный и реализованный алгоритм коррекции IPC позволяет обеспечить уровень погрешности в методе РТУ, сопоставимый с погрешностью адаптивного PIV, при возможности достижения высокого пространственного разрешения без уменьшения размеров полной исследуемой области потока в продолах одного эксперимента.
Заключение
В статье исследован метод РТУ с корреляционной коррекцией IPC. Алгоритм коррекции протестирован на синтетических изображениях, изучено влияние диаметра образов частиц, их концентрации, а также вони чины смещения на точность рассчитанного ноля скорости. Определены оптимальные параметры дня изображений частиц и параметры алгоритма обработки. Показано, что применение метода IPC позволяет использовать в эксперименте более широкий диапазон размеров частиц при том же уровне ошибки определения их смещения, что и без применения коррекции. Результаты тестов на синтетических изображениях с равномерно смещенными частицами показали уменьшение полной ошибки в три раза.
Для типичного PIV-изображения с плотностью частиц около 0.02 частиц/пикс2 определен оптимальный размер корреляционного окна, равный 5 пике, практически не снижающий пространственное разрешение метода РТУ.
Проведено экспериментальное исследование гидродинамики осесимметричпой импактпой струи в двух режимах с различным пространственным разрешением. Турбулентные статистические характеристики течения вблизи импактпой поверхности были получены при помощи стандартного РТУ, РТУ с коррекцией IPC и адаптивного PIV методов.
Результаты эксперимента продемонстрировали преимущества исследуемого подхода. Важным качеством метода РТУ с корреляционной коррекцией IPC является небольшая величина ошибки, определяемой эффектом группировки значений смещений частиц у целых величин (peak-locking), сопоставимая с ошибкой в адаптивном PIV, При этом сохраняются преимущества, связанные с более высоким пространственным разрешением и низким уровнем шумов.
Предложен метод отсева ошибочных векторов на основе фильтра скользящего среднего, который может быть использован при обработке полей скорости на нерегулярной сетке с неравномерной концентрацией частиц на изображении. Адаптивный характер данного алгоритма позволяет сохранить высокое пространственное разрешение метода PTV,
Список литературы
[1] Токарев М.П., Маркович Д.М., Бильский А.В. Адаптивные алгоритмы обработки изображений частиц для расчета мгновенных полей скорости // Вычисл. технологии. 2007. Т. 12, № 3. С. 109-131.
[2] Theunissen R., Scarano F., Riethmuller M.L. An adaptive sampling and windowing interrogation method in PIV // Meas. Sci. Technol. 2007. Vol. 18. P. 275-287.
[3] Etoh R., Takehara K. The particle mask correlation method // Proc. 8th Intern. Svmp. on Flow Visualization. Italy, Sorrento, 1998. P. 283.
[4] Uemura Т., Yamamoto F., Ohmi K. A high speed algorithm of image analysis for real time measurement of two-dimensional velocity distribution // Flow Visualization / Eds. B. Khalighi, M. Braun, C. Freitas. N.Y.: ASME, 1989. P. 129-134.
[5] Kobayashi Т., Saga Т., Segawa S. Multipoint velocity measurement for unsteady flow field by digital image processing // Flow Visualization. 1989. Vol. 5. P. 197-202.
[6] Hassan Y.A., Canaan R.E. Full-field bubbly flow velocity measurements using a multiframe particle tracking technique // Exp. Fluids. 1991. Vol. 12. P. 49-60
[71 Mikheev A.V., Zubtsov V.M. Enhanced particle-tracking velocimetrv (EPTV) with a combined two-component pair-matching algorithm // Meas. Sci. Technol. 2008. Vol. 19. P. 1-16.
[8] Cowen E., Monismith S. A hybrid digital particle tracking velocimetrv technique // Exp. Fluids. 1996. Vol. 22. P. 199-201.
[9] Ishikawa XL. mljral Y., Wada A. et al. A novel algorithm for particle tracking velocimetrv using the velocity gradient tensor // Ibid. 2000. Vol. 29. P. 19-31.
[10] Takehara K., Adrian R., Etoh G., Christensen K. A Kalman tracker for superresolution PIV // Ibid. 2000. Vol. 29. P. 34-41.
[11] Wernet M. New insights into particle image velocimetrv data using fuzzv-logic-based correlation/particle tracking processing // Ibid. 2000. Vol. 30. P. 434-447.
[12] Stitou A., Riethmuller M.L. Extention of PIV to super resolution using PTV // Meas. Sci. Technol. 2001. Vol. 12. P. 1398-1403.
[13] Theunissen R., Stitou A. and Riethmuller M.L. A novel approach to improve the accuracy of PTV methods // Proc. 12th Intern. Svmp. on Appl. of Laser Techniques to Fluid Mechanics. Lisbon, Portugal, 2004.
[14] Влек S.J., Lee S.J. A new two-frame particle tracking algorithm using match probability // Exp. Fluids. 1996. Vol. 22. P. 23-32.
[15] Heinz O.M., Ilyushin B.B., Markovich D.M. Application of a PDF's method for the statistical processing of experimental data // Intern. J. Heat Fluid Flow. 2004. Vol. 25. P. 864-874.
[16] Rappel XL. Willert С.Е., Wereley S.T., Kompenhans J. Particle Image Velocimetrv. A Practical Guide. Springer, 2007.
[17] Host-Masden A., McCluskey D.R. On the Accuracy and Reliability of PIV measurements // Proc. of the Seventh Intern. Svmp. on Appl. of Laser Techniques to Fluid Mechanics. Lisbon, Portugal, 1994.
[18] Estrada-Perez C.E., Dominguez Ontiveros E.E., Hassan Y.A., Ortiz-Villaf-uerte J. Near-wall study of turbulent channel flows with an improved PTV algorithm // Proc. 6th Intern. Svmp. on Particle Image Velocimetrv Pasadena. California, USA, 2005.
[19] Westerweel J. Theoretical analisvs of the measurement precision in particle image velocimetrv // Exp. Fluids. 2000. Vol.'29 (Suppl.). P. S3 SI2.
[20] Astarita Т., Cardone G. Analysis of interpolation schemes for image deformation methods in PIV // Ibid. 2005. Vol. 38. P. 233-243.
[21] Astarita T. Analysis of interpolation schemes for image deformation methods in PIV: effect of noise on the accuracy and spatial resolution // Ibid. 2006. Vol. 40. P. 977-987.
[22] Ахметбеков E.K., Бильский А.В., Ложкин Ю.А. и др. Система управления экспериментом и обработки данных, полученных методами цифровой трассерной визуализации (ActualFlow) // Вычисл. методы и программирование. 2006. Т. 7. С. 79-85.
[23] Westerweel J., Scarano F. Universal outlier detection // Exp. Fluids. 2005. Vol. 39. P. 1096-1100.
[24] Cowen E.A., Monismith S.G. A hybrid digital particle tracking velocimetrv technique // Ibid. 1997. Vol. 22. P. 199-211.
[25] Saarenrinne P., pllrto XL. Eloranta H. Experiences of turbulence measurement with PIV // Meas. Sci. Technol. 2001. Vol. 12. P. 1904-1910.
Поступила в редакцию 30 сентября 2008 г., с доработки — 3 февраля 2010 г.
