Локализация особенностей газодинамических полей при помощи комплексных ортогональных вейвлет-разложений
( Localization of singularities of gas dynamics fields with the aid of complex orthogonal wavelet transform
Preprint, Inst. Appl. Math., the Russian Academy of Science)

1.    ВВЕДЕНИЕ

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

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

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

Известен ряд прямых подходов к  задаче о выделении сингулярностей см. [7]. Другим из возможных подходов к этой задаче  является переход к вейвлет-преобразованию векторного или скалярного поля с последующим анализом уже этого объекта. Этот подход успешно используется в задачах локализации разного рода особенностей, в частности, краев при анализе изображений – см., например, [1]. Возможность использования фазы непрерывного комплекснозначного вейвлет-преобразования в качестве параметра, позволяющего локализовать особенности сигнала, также известна (см. [1],[2]). Основным достоинством такого подхода является его многомасштабность, позволяющая классифицировать сингулярности по их интенсивности и итеративно уточнять их положение. Однако, непрерывное вейвлет-преобразование, особенно для двух- и трехмерных полей, требует огромного объема вычислений. Поэтому в данной работе для анализа газодинамических полей предлагается использование дискретного  комплекснозначного  ортогонального вейвлет-преобразования (ОВП), сочетающего такой чувствительный индикатор нарушения гладкости, как фаза, и быстрый (линейный по объему входных данных) алгоритм вычислений.  Конечно, для обобщенного решения газодинамической задачи, заданной в виде системы законов сохранения, мы не имеем априорного описания всех возможных сингулярностей решения, а численный метод может вносить в решение свои возмущения. Предлагаемый в данной работе подход фиксирует наличие сингулярностей, но  оставляет открытым вопрос об их характере.

Разные методы выделения сингулярностей обладают различной чувствительностью и могут быть адаптированы для определения интересующих нас сингулярностей   выделенного типа. В отличие от  средств классического гармонического анализа (см. [5]), применяющихся в вычислительной математике, здесь необходимым этапом анализа является вычислительный эксперимент. Например, при использовании вейвлетных разложений приходится сталкиваться с необходимостью выбора оптимального материнского вейвлета, по которому проводится разложение.

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

2.    Комплексное ортогональное вейвлет-преобразование

Одним из наиболее широко используемых вариантов вейвлет-преобразования является  преобразование, использующее классические ортогональные вейвлеты Добеши с компактными носителями [8]. Известно, что размер носителя и гладкость этих вейвлетов взаимосвязаны и обычно используются построенные в [8] семейства «наиболее гладких» или «наиболее симметричных» вейвлетов. Однако конструкция многомасштабной аппроксимации, используемая для построения вейвлетов Добеши, позволяет получить и другие семейства ортогональных вейвлетов, обладающих различными дополнительными свойствами. В работе [9] таким способом были построены новые семейства ортогональных вейвлетов с компактными носителями, среди которых есть и комплекснозначные. Их значение для анализа данных состоит в том, что появляется возможность использовать такой дополнительный информативный параметр, как аргумент (фаза) вычисленных коэффициентов.

Эти вейвлеты обладают и обычными достоинствами ортогональных вейвлетов с компактным носителем. Коэффициенты многомасштабных разложений вычисляются при помощи двух сопряженных квадратурных фильтров (Conjugate Quadrature Filters), один из которых низкочастотный, а другой высокочастотный, а разложение осуществляется при помощи быстрого алгоритма Малла [1]. Точное восстановление данных может быть выполнено с помощью этих же фильтров, но для задачи выделения особенностей это несущественно.

На рисунке 1 приведены частотные характеристики использованных в наших экспериментах фильтров из семейства, описанного в [9] и обозначаемого DCOLA (Daubechies Complex Orthogonal Least-Asymmetric), а именно фильтров DCOLA8, содержащих по 8 коэффициентов.

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

,

                                                                  (1)

,

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

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

, .                                 (2)

Графики этих функций на рис. 1 (слева) показывают, что фильтры  h и g , примененные по формулам (1), действительно разделяют сигнал на сглаженную (низкочастотную) часть и высокочастотные «детали».

Двумерный вариант этого алгоритма на каждом шаге расщепляет массив, полученный на предыдущем шаге на 4 матрицы-компоненты. При этом шаг одномерного алгоритма применяется сначала к строкам, а потом к столбцам.

3.    Описание экспериментов

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

К этой функции был применен один шаг комплексного ОВП с вейвлетом DCOLA8. На рис. 3 показаны фаза и амплитуда одной из трех высокочастотных составляющих преобразованного массива. Видно, что ненулевые значения фазы сконцентрированы в точности вдоль отрезка, по которому проходит разрыв градиента (узкая темная полоса в нижней части возникла из-за периодизации при вычислении ОВП, и не учитывается при анализе). Положение этих точек позволяет сразу же идентифицировать положение разрыва. В изображении амплитуды этот разрыв тоже четко обозначен светлой линией максимальных значений.

На рис. 4 и 5 показаны аналогичные результаты для функции, имеющей разрыв первого рода вдоль той же линии. Отметим, что в этом случае скачок локализуется при помощи распределения фазы не только высокочастотной, но и низкочастотной составляющей.

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

Перейдем к результатам обработки массива данных, полученных авторами благодаря любезности А.Е. Луцкого. Рассматривается нестационарная  задача о дифракции сильной ударной волны, набегающей на зону прогретого газа, расположенную в прямоугольнике (x,y)[1,5][0,0.5].  Расчет проводился по методу Годунова, т.е. без выделения ударных волн (см. [6]).

На рис. 6  изображено распределение значений изучаемой функции – плотности газа в области течения в некоторый фиксированный момент времени. В начальный момент плотность газа в темной зоне составляла 0.03, а в светлой 1. На рисунке визуально отслеживается  сложная структура разрывов, включая тройную точку в набегающей волне.

Эти линии четко проявляются в распределении фазы комплексного ОВП. На рис. 7 показаны распределение фазы и модуля одной из высокочастотных составляющих.

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

Очевидно и наличие структур, более тонких, чем изолированный разрыв функции или ее градиента. Физическая интерпретация этих структур – одна из задач дальнейших исследований. Однако ясно, что они тесно связаны с интересными и неочевидными особенностями течения. Например, имеется довольно точно локализованная граница области «хаотического» изменения фазы на рис. 7 и, если верхняя и задняя граница области имеют очевидную физическую интерпретацию, то физическая причина наличия «переднего фронта» этой зоны пока до конца не ясна.

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

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

В данной задаче наиболее интенсивные разрывы (особенности течения) локализуются по распределению фазы низкочастотной части преобразования. При применении к распределению фазы низкочастотной части вейвлетного преобразования детектора краев Канни (Canny edge detector, [7])  получен набор кривых, отвечающих особенностям течения (рис. 8). Эти кривые соответствуют расположению физических особенностей течения. Например, вертикальные белые линии на рис. 8 отвечают положению набегающей ударной волны, горизонтальная линия соответствует контактному разрыву (скачку плотности в начальных данных). Четко локализуются две косые ударные волны (скачки уплотнения в системе координат, связанной с набегающей ударной волной)  и четверная точка в контактном разрыве.

На рисунке 9  локализованы сингулярности в распределении давления. На линии невозмущенного контактного разрыва при x[~4.3,5] давление меняется гладко. Если учесть законы сохранения, то из рисунков 8-9 нетрудно усмотреть наличие слабых разрывов. Структуры аналогичного вида появляются и при анализе компонент вектора скорости, однако здесь имеется еще большее число вопросов, требующих дополнительного анализа.  Интересной задачей является изучение возможности алгоритмического различения сильных и слабых ударных волн с помощью подбора адекватного материнского вейвлета.

4.    Заключение

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

Авторы пользуются приятной возможностью поблагодарить А.Е. Луцкого не только за предоставленные расчетные материалы, но и за плодотворные обсуждения. Работа была поддержана Программой 3.2 ОМН РАН.

5.    ЛИТЕРАТУРА

[1] Stéphane Mallat. A wavelet tour of signal processing, Academic Press, 1997.

[2] Matthias Holschneider. Wavelets: an analysis tool, Calderon Press, Oxford, 1995.

[3] С.К. Годунов, А.В. Забродин и др. Численное решение многомерных задач газовой динамики. Наука 1976

[4] Теоретические основы и конструирование численных методов решения задач математической физики. Под ред. К.И. Бабенко. М., «Наука», 1979.

[5] К.И. Бабенко. Основы численного анализа. М., «Наука», 1984.

[6] V.L.Bychkov, A.I.Klimov, A.E.Lutsky. Numerical Simulation of Shock Wave Propagation in the Gas with Temperature Discontinuities with Regard of Vibration-Translation Energy Transfer. //The 2nd Workshop on Magnetoplasma Aerodynamics in Aerospace Applications - April 5-7 2000 -Moscow - IVTAN - P.289-296.

[7] John Canny. "A Computational Approach to Edge Detection", IEEE Transactions on Pattern Analysis and Machine Intelligence, 1986,Vol. PAMI-8, No. 6, pp. 679-698.

[8] И. Добеши, «Десять лекций по вейвлетам», М., «Регулярная и хаотическая динамика», 1999.

[9] Carl Taswell, “Constraint-Selected and Search-Optimized Families of Daubechies Wavelet Filters Computable by Spectral Factorization”, Journal of Computational and Applied Mathematics (Special Issue on Numerical Analysis in the 20th Century), 121:179-195, June, 2000.