Обработка сигналов и изображений в пространственно частотной фильтрации

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

Вейвлетный базис пространства L²(R), R(-, ), целесообразно конструировать из финитных функций, которые должны стремиться к нулю на бесконечности. Чем быстрее эти функции стремятся к нулю, тем удобнее использовать их в качестве базиса преобразования при анализе реальных сигналов. Допустим, что такой функцией является psi - функция ??t?, равная нулю за пределами некоторого конечного интервала и имеющая нулевое среднее значение по интервалу задания. Последнее необходимо для задания локализации спектра вейвлета в частотной области (???????0 при ?=0 и ?>??. На основе этой функции сконструируем базис в пространстве L²(R) с помощью масштабных преобразований независимой переменной.

Функция изменения частотной независимой переменной в спектральном представлении сигналов отображается во временном представлении растяжением/сжатием сигнала. Для вейвлетного базиса это можно выполнить функцией типа ?(t) => ?(a^mt), a = const, m = 0, 1, … , M, т.е. путем линейной операции растяжения/сжатия, обеспечивающей самоподобие функции на разных масштабах представления. Однако локальность функции ?(t) на временной оси требует дополнительной независимой переменной последовательных сдвигов функции ?(t) вдоль оси, типа ?(t) => ?(t+k), для перекрытия всей числовой оси пространства R(-, ). C учетом обеих условий одновременно структура базисной функции может быть принята следующей:

?(t) => ?(a^mt+k). (1.2)

Для упрощения дальнейших выкладок значения переменных m и k примем целочисленными. При приведении функции (1.2) к единичной норме, получаем:

?_mk(t) = a^m/2 ?(a^mt+k). (1.3)

Если для семейства функций ?_mk(t) выполняется условие ортогональности:

?_nk(t), ?_lm(t) =?_nk(t)·?*_lm(t) dt =?_nl·?_km, (1.4)

то семейство ?_mk(t) можно использовать в качестве ортонормированного базиса пространства L²(R). Произвольную функцию этого пространства можно разложить в ряд по базису ?_mk(t):



s(t) =S_mk ?_mk(t), (1.5)

где коэффициенты S_mk - проекции сигнала на новый ортогональный базис функций, как и в преобразовании Фурье, определяются скалярным произведением

S_mk = s(t), ?_mk(t) =s(t)??_mk(t) dt, (1.6)

при этом ряд равномерно сходиться:

||s(t) -S_mk ?_mk(t),|| = 0.

При выполнении этих условий базисная функция преобразования ?(t) называется ортогональным вейвлетом.

Простейшим примером ортогональной системы функций такого типа являются функции Хаара. Базисная функция Хаара определяется соотношением

?(t) = ( 1.7)

Легко проверить, что при а = 2, m = 0, 1, 2, ..., k = 0, 1,2, … две любые функции, полученные с помощью этого базисного вейвлета путем масштабных преобразований и переносов, имеют единичную норму и ортогональны. На рис. 1.5 приведены примеры функций для первых трех значений m и b при различных их комбинациях, где ортогональность функций видна наглядно.

Функции Хаара

Вейвлетный спектр, в отличие от преобразования Фурье, является двумерным и определяет двумерную поверхность в пространстве переменных m и k. При графическом представлении параметр растяжения/сжатия спектра m откладывается по оси абсцисс, параметр локализации k по оси ординат - оси независимой переменной сигнала. Математику процесса вейвлетного разложения сигнала в упрощенной форме рассмотрим на примере разложения сигнала s(t) вейвлетом Хаара с тремя последовательными по масштабу m вейвлетными функциями с параметром а=2, при этом сам сигнал s(t) образуем суммированием этих же вейвлетных функций с одинаковой амплитудой с разным сдвигом от нуля, как это показано на рис.

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

Для начального значения масштабного коэффициента сжатия m определяется функция вейвлета (?1(t) на рис. 1.6), и вычисляется скалярное произведение сигнала с вейвлетом ?1(t), s(t+k) с аргументом по сдвигу k. Для наглядности результаты вычисления скалярных произведений на рис. 1.6 построены по центрам вейвлетных функций (т.е. по аргументу k от нуля со сдвигом на половину длины вейвлетной функции). Как и следовало ожидать, максимальные значения скалярного произведения отмечаются там, где в сигнале s_k локализована эта же вейвлетная функция.

После построения первой масштабной строки разложения, меняется масштаб вейвлетной функции (?2 на рис. 1.6) и выполняется вычисление второй масштабной строки спектра, и т.д.

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

На рис. 1.7 приведен пример графического отображения вейвлетной поверхности реального физического процесса /4/. Вид поверхности определяет изменения во времени спектральных компонент различного масштаба и называется частотно-временным спектром. Поверхность изображается на рисунках, как правило, в виде изолиний или условными цветами. Для расширения диапазона масштабов может применяться логарифмическая шкала.

Пример вейвлетного преобразования

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

· Вейвлет-функции (t), как psi-функции времени с нулевым значением интеграла и частотным фурье-образом (щ). Этой функцией, которую обычно и называют вейвлетом, выделяются локальные особенности сигнала. В качестве вейвлетов обычно выбираются функции, хорошо локализованные и во временной, и в частотной области. Пример временного и частотного образа функцийй приведен на рис. 1.4.

· Масштабирующей функции ?(t), как временной скейлинг-функции phi с единичным значением интеграла, которой выполняется грубое приближение (аппроксимация) сигнала.

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

Непрерывное вейвлет-преобразование (CWT- Continious Wavelet Transform). Допустим, что мы имеем функции s(t) с конечной энергией в пространстве L²(R), определенные по всей действительной оси R(-, ). Для финитных сигналов с конечной энергией средние значения сигналов должны стремиться к нулю на ±.

Непрерывным вейвлет-преобразованием (или вейвлетным образом) функции s(t) L²(R) называют функцию двух переменных:

С(a,b) = s(t), ?(a,b,t) = s(t)??(а,b,t) dt, a, b R, a ? 0. (2.1)

где вейвлеты ?(a,b,t) ?_ab(t) - масштабированные и сдвинутые копии порождающего вейвлета ?(t) L²(R), совокупность которых создает базис пространства L²(R).

Порождающими функциями могут быть самые различные функции с компактным носителем - ограниченные по времени и местоположению на временной оси, и имеющие спектральный образ, локализованный на частотной оси. Базис пространства L²(R) целесообразно конструировать из одной порождающей функции, норма которой должна быть равна 1. Для перекрытия функцией вейвлета всей временной оси пространства используется операция сдвига (смещения по временной оси): ?(b,t) = ?(t-b), где значение b для CWT является величиной непрерывной. Для перекрытия всего частотного диапазона пространства L²(R) используется операция временного масштабирования вейвлета с непрерывным изменением независимой переменной: ?(a,t) = |а|^-1/2?(t/а). На рис. 1.4. видно, что если временной образ вейвлета будет расширяться (изменением значения параметра 'а'), то его "средняя частота" будет понижаться, а частотный образ (частотная локализация) перемещаться на более низкие частоты. Таким образом, путем сдвига по независимой переменной (t-b) вейвлет имеет возможность перемещаться по всей числовой оси произвольного сигнала, а путем изменения масштабной переменной 'а' (в фиксированной точке (t-b) оси) "просматривать" частотный спектр сигнала по определенному интервалу окрестностей этой точки.

С использованием этих операций вейвлетный базис функционального пространства образуется путем масштабных преобразований и сдвигов порождающего вейвлета??(t):

?(a,b,t) = |а|^-1/2?[(t-b)/а], a, b R, a ? 0, ?(t) L²(R). (2.2)

Нетрудно убедиться, что нормы вейвлетов ?(a,b,t) равны норме ?(t), что обеспечивает нормировочный множитель |а|^-1/2. При нормировке к 1 порождающего вейвлета ?(t) все семейство вейвлетов также будет нормированным. Если при этом выполняется требование ортогональности функций, то функции ?(a,b,t) образуют ортонормированный базис пространства L²(R).

Понятие масштаба ВП имеет аналогию с масштабом географических карт. Большие значения масштаба соответствуют глобальному представлению сигнала, а низкие значения масштаба позволяют различить детали. В терминах частоты низкие частоты соответствуют глобальной информации о сигнале, а высокие частоты - детальной информации и особенностям, которые имеют малую протяженность, т.е. масштаб вейвлета, как единица шкалы частотно-временного представления сигналов, обратен частоте. Масштабирование расширяет или сжимает сигнал. Большие значения масштабов соответствуют расширениям сигнала, а малые значения - сжатым версиям. В определении вейвлета коэффициент масштаба а стоит в знаменателе. Соответственно, а > 1 расширяет сигнал, а < 1 сжимает его.

Процедура преобразования стартует с масштаба а=1 и продолжается при увеличивающихся значениях а, т.e. анализ начинается с высоких частот и проводится в сторону низких частот. Первое значение 'а' соответствует наиболее сжатому вейвлету. При увеличении значения 'а' вейвлет расширяется. Вейвлет помещается в начало сигнала (t=0), перемножается с сигналом, интегрируется на интервале своего задания и нормализуется на 1/. Результат вычисления С(a,b) помещается в точку (a=1, b=0) масштабно-временного спектра преобразования. Сдвиг b может рассматриваться как время с момента t=0, при этом координатная ось b повторяет временную ось сигнала. Для полного включения в обработку всех точек входного сигнала требуется задание начальных и конечных условий преобразования (определенных значений входного сигнала при t<0 и t>t_max на полуширину окна вейвлета). При одностороннем задании вейвлетов результат относится, как правило, к временному положению средней точки окна вейвлета.

Затем вейвлет масштаба а=1 сдвигается вправо на значение b и процедура повторяется. Получаем значение, соответствующее t=b в строке а=1 на частотно-временном плане. Процедура повторяется до тех пор, пока вейвлет не достигнет конца сигнала. Таким образом получаем строку точек на масштабно-временном плане для масштаба а=1.

Для вычисления следующей масштабной строки значение а увеличивается на некоторое значение. При CWT в аналитической форме ?b0 и ?a0. При выполнении преобразования в компьютере выполняется увеличение обоих параметров с определенным шагом. Тем самым осуществляется дискретизация масштабно-временной плоскости.

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

На рис. приведен пример модельного сигнала и спектра его непрерывного вейвлет-преобразования.

При непрерывных значениях параметров 'а' и 'b' сигналу, определенному на R, соответствует вейвлетный спектр R Ч R. Отсюда следует, что вейвлетный спектр НПВ имеет огромную избыточность.

Обратное преобразование. Так как форма базисных функций ?(a,b,t) зафиксирована, то вся информация о сигнале в переносится на значения функции С(a,b). Точность обратного интегрального вейвлет-преобразования зависит от выбора базисного вейвлета и способа построения базиса, т.е. от значений базисных параметров a, b. Строго теоретически вейвлет может считаться базисной функцией L²(R) только в случае его ортонормированности. Для практических целей непрерывного преобразования часто бывает вполне достаточна устойчивость и "приблизительность" ортогональности системы разложения функций. Под устойчивостью понимается достаточно точная реконструкция произвольных сигналов. Для ортонормированных вейвлетов обратное вейвлет-преобразование записывается с помощью того же базиса, что и прямое:

s(t) = (1/C_?) (1/a²) С(a,b) ?(a,b,t) da db. (2.3)

где C_? - нормализующий коэффициент:

C_? = (|?(?)|² /?) d? < . (2.4)

Условие конечности C_? ограничивает класс функций, которые можно использовать в качестве вейвлетов. В частности, при щ=0, значение ?(?) должно быть равно нулю. Это обеспечивает условие компактности фурье-образа вейвлета с локализацией вокруг некоторой частоты ?_o - средней частоты вейвлетной функции. Следовательно, функция ?(t) должна иметь нулевое среднее значение по области его определения:

?(t) dt =0.

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

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

s(t) = (1/C_?)(1/a²) С(a,b) ?^#(a,b,t) da db, (2.3')

где индексом ?^#(a,b,t) обозначен ортогональный "двойник" базиса ?(a,b,t), о котором будет изложено ниже.

Таким образом, непрерывное вейвлет-преобразование представляет собой разложение сигнала по всем возможным сдвигам и сжатиям/растяжениям некоторой локализованной финитной вейвлет-функции. При этом переменная 'a' определяет масштаб вейвлета и эквивалентна частоте в преобразованиях Фурье, а переменная 'b' - сдвиг вейвлета по сигналу от начальной точки в области его определения, шкала которого повторяет временную шкалу анализируемого сигнала. Вейвлетный анализ является частотно-пространственным анализом сигналов.

В качестве примера рассмотрим вейвлет-преобразование чистого гармонического сигнала s(t), приведенного на рис. 2.2. На этом же рисунке ниже приведены вейвлеты ?_a(t) симметричного типа разных масштабов.

Скалярное произведение (2.1) "просмотра" сигнала вейвлетом определенного масштаба 'a' может быть записано в следующей форме:

C_a(b)= s(t), ?_a(t+b) =s(t)??_a(t+b) dt. (2.5)



Но выражение (2.5) эквивалентно взаимной корреляционной функции R_a(b) сигналов s(t) и ?_а(t). Если сигнал s(t) представляет собой гармонику, а второй сигнал симметричен, задан на компактном носителе и имеет нулевое среднее значение, то, как известно, форма взаимной корреляционной функции таких сигналов также является центрированным гармоническим сигналом. В частотной области скалярное произведение двух функций отображается произведением Фурье-образов этих функций, которые приведены на рисунке в правом столбце спектров. Масштабы спектров ?_a(?) и R_a(?) для наглядности сопоставления нормированы к спектру s(t). Максимальная амплитуда гармоники R_а(b) будет наблюдаться при совпадении средней частоты локализации вейвлета ?_а(t) определенного масштаба 'а' в частотной области с частотой сигнала s(t), что и можно видеть на рис. 2.2 для функции R_a(b) при масштабе вейвлета a=20. Результирующий вейвлетный спектр непрерывного вейвлет-преобразования гармоники приведен на левом нижнем графике и показывает точное положение на временной оси 'b' максимумов и минимумов гармонического сигнала.

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

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



a = а_о^-m, b = k·а_о^-m, a_o > 1, m, k I,

где I - пространство целых чисел {-, }, m - параметр масштаба, k - параметр сдвига. Базис пространства L²(R) в дискретном представлении:

?_mk(t) = |а_о|^m/2?(а_о^mt-k), m,k I, ?(t) L²(R). (2.6)

Вейвлет-коэффициенты прямого преобразования:

C_mk =s(t)??_mk(t) dt. (2.7)

Значение 'a' может быть произвольным, но обычно принимается равным 2, при этом преобразование называется диадным вейвлет-преобразованием. Для диадного преобразования разработан быстрый алгоритм вычислений, аналогичный быстрому преобразованию Фурье, что предопределило его широкое использование при анализе массивов цифровых данных.

Обратное дискретное преобразование для непрерывных сигналов при нормированном ортогональном вейвлетном базисе пространства:

s(t) = C_mk??_mk(t). (2.8)

Число использованных вейвлетов по масштабному коэффициенту m задает уровень декомпозиции сигнала, при этом за нулевой уровень (m = 0) обычно принимается уровень максимального временного разрешения сигнала, т.е. сам сигнал, а последующие уровни (m < 0) образуют ниспадающее вейвлет-дерево. В программном обеспечении вычислений для исключения использования отрицательной нумерации по m знак 'минус' обычно переносится непосредственно в (2.6), т.е. используется следующее представление базисных функций:

?_mk(t) = |а_о|^-m/2?(а_о^-mt-k), m,k I, ?(t) L²(R). (2.6')

Устойчивость дискретного базиса определяется следующим образом.

Функция ?(t) L²(R) называется R-функцией, если базис на ее основе по (2.6) является базисом Рисса (Riesz). Для базиса Рисса существуют значения А и В, 0 < A ? B < , для которых выполняется соотношение

A||C_mk||² ? || C_mk??_mk(t)||² ? B||C_mk||²,

если энергия ряда C_mk конечна. При этом для любой R-функции существует базис ?^#_mk(t), который ортогонален базису ?_mk(t). Его называют ортогональным "двойником" базиса ?_mk(t), таким, что

?_mk(t), ?^#_nl(t) = ?_mn·?_kl.

Если A = B = 1 и а_о = 2, то семейство базисных функций {?_mk(t)} является ортонормированным базисом и возможно полное восстановление исходного сигнала, при этом ?_mk(t) ? ?^#_mk(t) и для реконструкции сигналов используется формула (2.8). Если ?(t) не ортогональный вейвлет, но имеет "двойника", то на базе "двойника" вычисляется семейство ?^#_mk(t), которое и используется при обратном преобразовании вместо ?_mk(t), при этом точное восстановление исходного сигнала не гарантировано, но оно будет близко к нему в среднеквадратическом смысле.

Как и для непрерывного вейвлет-преобразования, обратное дискретное преобразование (2.8) не может выполнить восстановление нецентрированных сигналов в силу нулевого первого момента вейвлетных функций и, соответственно, центрирования значения вейвлет-коэффициентов C_mk при прямом вейвлет-преобразовании. Поэтому при обработке числовых массивов данных дискретные вейвлеты используются, как правило, в паре со связанными с ними дискретными скейлинг-функциями. Скейлин-функции имеют с вейвлетами общую область задания и определенное соотношение между значениями, но первый момент скейлин-функций по области определения равен 1. Если вейвлеты рассматривать, как аналоги полосовых фильтров сигнала, в основном, высокочастотных при выделении локальных особенностей в сигнале, то скейлин-функции вейвлетов представляет собой аналоги низкочастотных фильтров, которыми из сигнала выделяются в отдельный массив составляющие, не прошедшие вейвлетную фильтрацию. Так, например, порождающая скейлинг-функция вейвлета Хаара задается следующим выражением:

При обозначении скейлинг-функций индексом ?_mk(t) аналитика скейлин-функций повторяет выражения (2.6-2.7) и образует дополнительный базис пространства L²(R). Сумма вейвлет-коэффициентов и скейлинг-коэффициентов разложения сигналов соответственно дает возможность выполнять точную реконструкцию сигналов, при этом вместо (2.8) используется следующее выражение обратного вейвлет-преобразования:

s(t) =Сa_k ?_k(t) + Сd_mk??_mk(t), (2.9)

где Ca_k - скейлин-коэффициенты, которые обычно называют коэффициентами аппроксимации сигнала, Cd_mk - вейвлет-коэффициенты или коэффициенты детализации.

Частотно-временная локализация вейвлет-анализа. Реальные сигналы, как правило, конечны и принадлежат пространству L²(R). Частотный спектр сигналов обратно пропорционален их длительности. Соответственно, достаточно точный низкочастотный анализ сигнала должен производиться на больших интервалах его задания, а высокочастотный - на малых.

Если считать, что каждый вейвлет имеет определенную "ширину" своего временного окна, которому соответствует определенная "средняя" частота спектрального образа вейвлета, обратная его масштабному коэффициенту а, то семейства масштабных коэффициентов вейвлет-преобразования можно считать аналогичными семействам частотных спектров оконного преобразования Фурье, но с одним принципиальным отличием. Масштабные коэффициенты изменяют "ширину" вейвлетов и, соответственно, "среднюю" частоту их фурье-образов, а, следовательно, каждой частоте соответствует своя длительность временного окна анализа, и наоборот. Так малые значения параметра а, характеризующие быстрые составляющие в сигналах, соответствуют высоким частотам, а большие значения - низким частотам. За счёт изменения масштаба вейвлеты способны выявлять различия на разных частотах, а за счёт сдвига (параметр b) проанализировать свойства сигнала в разных точках на всём исследуемом временном интервале. Многоразмерное временное окно вейвлет-преобразования адаптировано для оптимального выявления и низкочастотных, и высокочастотных характеристики сигналов. При этом, на высоких частотах лучше разрешение по времени, а на низких - по частоте. Для высокочастотной компоненты сигнала мы можем точнее указать ее временную позицию, а для низкочастотной - ее значение частоты.

Изменение частотно-временного окна вейвлета определяет угол влияния значений функции в произвольных точках t_i на значения коэффициентов С(а,b). И наоборот, угол влияния из точки С(a_i,b_i) на ось t определяет интервал значений функции, которые принимают участие в вычислении данного коэффициента С(a_i,b_i) - область достоверности. Схематически это показано на рис. 2.3.

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

Достоинства и недостатки вейвлетных преобразований.

· Вейвлетные преобразования обладают всеми достоинствами преобразований Фурье.

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

· Вейвлетные базисы, в отличие от преобразования Фурье, имеют много разнообразных базовых функций, свойства которых ориентированы на решение различных задач. Базисные вейвлеты могут реализоваться функциями различной гладкости.

· Недостатком вейвлетных преобразований является их относительная сложность.

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

Непрерывные вейвлеты дают несколько более наглядное представление результатов анализа в виде поверхностей вейвлет-коэффициентов по непрерывным переменным. На рис. 2.4 анализируемый сигнал состоит из двух модулированных гауссианов. Преобразование вейвлетом Морлета четко показывает их пространственную и частотную локализацию.

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

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

Пространственная фильтрация

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

Существует четкое соответствие между временными / частотными областями и пространственными / волновыми областями. Рассмотрим простой «треугольный» фильтр, при значениях дискретизации 1, 2, 3, 2, 1. Во временной области он (и его амплитудные спектры) выглядят следующим образом:

Временные дискреты (через каждые 4 мсек) преобразуются в функцию «типа синхронизации» в частотной области. Теперь рассмотрим ту же функцию в пространстве:

Кроме осей и того, что в пространстве мы рассматриваем и положительные, и отрицательные волновые числа, эти два рисунка идентичны.

Чтобы проверить и временные и пространственные частоты одновременно, в обработке используется такой тип анализа, как анализ FK (F = временная частота, K = пространственная частота или волновое число), для разработки некоторых видов пространственных фильтров (обычно, это FK-фильтры).

Преобразование блока трасс из времени и пространства в «F» и «K» довольно (или почти!) просто математически, но занимает время. Каждая трасса преобразуется из времени в частоту через преобразование Фурье, а затем каждый частотный компонент (сложный термин) по всем трассам преобразуется через преобразование Фурье, чтобы получить спектр FK.



Этот рисунок показывает пару секунд (хотите верьте, хотите нет!) синтетической сейсмограммы, которая содержит всего одну частоту, падающую слева направо при постоянном угле наклона. Если мы посмотрим вниз по одной трассе на этом рисунке, мы увидим только одну частоту. Если вы посмотрите поперек трасс (например, вдоль линии времени), мы увидим порядок пиков и падений, которые опять представляют всего одну (пространственную) частоту. А это спектр FK сейсмограммы, приведенной выше. Вертикальная шкала от 0 до 125 Гц (наши данные записаны при шаге дискретизации 4 мсек). Горизонтальная шкала от -20 до +20 циклов на километр (наши трассы расположены через 25 м). Цвета представляют амплитуду, причем все данные вышеприведенной записи сжаты в один маленький крошечный «пичок», показанный желтым. Т.к. этот рисунок пытается показать трехмерную поверхность в двухмерном измерении, вы должны представить, что желтый пик торчит из экрана в вашем направлении! Если бы нам нужно было удалить пик из спектра FK (просто вырезав его при редакции), то когда мы преобразовывали бы обратно во время и пространство, у нас были бы мертвые трассы! Заметьте, что мы обычно не смотрим на фазовую часть спектра FK (хотя она, конечно, рассчитывается). Мы обычно предполагаем, что любая выполняемая редакция FK, оставит (очень сложные) фазы неизмененными. Еще одна синтетическая запись, на этот раз она выглядит более запутанной. На деле, она просто состоит из двух сложенных различных временных частот (две косинусные волны), падающие под двумя различными углами наклона. Как мы могли надеяться, спектр FK показывает и частоты и углы наклона как отдельные точки на спектре FK



Если мы «хирургически» удалим участок внутри красного круга в вышеприведенном преобразовании, а затем преобразуем обратно во временную / пространственную области, мы полностью удалим косинусоидальную волну с положительным углом наклона (вниз вправо). Все, что нам останется, это частотное проявление (временное) чуть повыше с отрицательным углом наклона (вверх вправо).

Помните наше правило угла наклона об увеличении времени с увеличением расстояния означает что на этой «записи» я предположил, что расстояние увеличивается в правую сторону.

Еще одна синтетическая сейсмограмма, содержащая несколько больше смысла! У нас имеется всего одно проявление с ограниченным диапазоном частот и одним углом наклона. (Я построил ее, просто добавив пик во все трассы при одном времени, фильтруя (сворачивая) его с единичным широкополосным фильтром, а затем применил статику, чтобы добавить угол наклона.)

Анализ FK также показывает угол наклона! Помните наше уравнение:

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

Для данного расстояния между трассами, мы можем нарисовать постоянный угол наклона на нашем анализе FK (это иногда для вас делается). Вот набор углов наклона от -20 до +20 миллисекунд на трассу при расстоянии между трассами 25 метров.

Любые углы наклона более 4 мсек (один временной дискрет) на трассу в волновом числе будут подвергаться действию аляйсинга.

Этот рисунок пытается показать кое-что из потенциального аляйсинга (я включил этот рисунок потому, что он красиво выглядит).

Если вы проследите за самой нижней линей справа (угол наклона +20 мсек / трассу) от исходной точки (центр), вы заметите, что она «отражается» или искажается (под действием аляйсинга) дважды в показанном диапазоне частот (от 0 до 125 Гц).

Еще один последний синтетический пример - три падающих проявления и их анализ FK (внизу). Самое крутое проявление (с самым широким диапазоном частот) довольно хорошо подвергается эффекту аляйсинга.



Рисунок ниже показывает результаты обратного преобразования после редакции двух углов наклона, близких к нулю (два проявления, которые пересекаются на линии времени).

Заметьте, что преобразование на концах работает неправильно - быстрое преобразование Фурье предполагает, что данные повторяются в пространстве, а мы усекли преобразование до всего 64 трасс.

Теперь мы можем вычесть трассы (внизу) из разреза (справа) и (более или менее) удалить большой угол наклона, подвергшийся аляйсингу.

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

А теперь давайте обсудим некоторые формы фильтрации углов налона.



Линейная фильтрация Тау-P

Когда мы преобразуем данные в область FK, наклонные проявления (во времени / пространстве) становятся наклонными проявлениями (в частоте / волновом числе). Преобразование Tau-P это еще один способ взглянуть на наши данные, который преобразует углы наклона (эффективно) в отдельные точки. Тогда мы можем редактировать данные в этой «наклонной» области и преобразовывать ее обратно во временную область

.

Вот проблема того типа, что обычно попадаются в тестах на IQ и им подобных. Даны общие цены (в любой предпочтительной валюте) в различных рядах и колонках, можете ли вы восстановить пропущенные? (Я не буду оскорблять ваш интеллект, давая готовый ответ!)Это просто графическое представление 4 одновременных уравнений для трех неизвестных. Уравнения все внутренне последовательны, поэтому есть только один возможный ответ. Если «11» внизу рисунка было изменено на что-то еще (например, на «12»), уравнения были бы внутренне непоследовательны, и нам нужно было бы перейти к решению наименьших квадратов, чтобы получить соответствующий ответ.

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

Как обычно, для демонстрации - синтетическая сейсмограмма.

Эта «полевая сейсмограмма» (один ПВ) имеет три наклонных проявления, расширяющиеся от ближних выносов (слева) по направлению к бoльшим выносам (справа). Трассы расположены через 25 м, а линии времени через 500 мсек.

0.08

«Горизонтальная скорость», показанная выше, - это просто скорость, с которой «проявление» движется по регистрирующей расстановке, расстояние между трассами разделить на угол наклона. Интервальное время пробега волны (термин, часто используемый в преобразованиях Tau-P) - это просто обратная величина этой скорости (обычно измеряется в секундах на километр).



Это линейное преобразование Tau-P вышеприведенной сейсмограммы. Цвета фона соответствуют углам наклона трех разноцветных проявлений, показанных выше.

Время этих «петель» - это время при нулевом выносе (грубо говоря, три трассы влево от записи, показанной выше - ближний вынос трассы составляет порядка 75 м).

Я использовал почти в три раза больше трасс в преобразовании, чем их присутствует в исходной записи. Это хорошая идея, т.к. она предоставляет больше статистики для обратного преобразования.Теперь у нас есть рисунок, который показывает время (вертикально) относительно угла наклона (горизонтально). Каждое из наших проявлений свернулось (более или менее) в маленькое «пятно» на преобразовании Tau-P, и они легко выделяются в этой области.



Чтобы сохранить редакцию в этой области простой (на практике мы можем «вырезать» любые участки, какие нам захочется), я просто преобразую обратно во временную область, используя только отрицательные углы наклона в вышеприведенном преобразовании (я сместил все вправо от сине-зеленой области - 0-вой угол наклона).

Я не стал раскрашивать этот рисунок, но должно быть очевидно, что я (более или менее) сдвинул одно проявление (C), которое имело положительный угол наклона («вниз» и вправо).

Несколько остатков от этого проявления остались на конечных трассах. Этот «краевой эффект» часто встречается во всех типах пространственной фильтрации - у фильтра нет больше данных, чтобы продолжить расчеты.

Чтобы минимизировать некоторые ошибки, связанные с фильтрацией Tau-P, мы зачастую преобразуем обратно то, что хотим удалить (например, два проявления с отрицательными углами наклона, показанные выше), а затем удаляем их из исходной записи. Любое «сглаживание», присущее преобразованиям, тогда менее заметно.

Мы посмотрим на несколько реальных примеров фильтрации Tau-P (и некоторых других методов в области Tau-P) через пару страниц. Мы также вернемся к другой форме фильтрации Tau-P, когда будем говорить об удалении кратных отражений.

А теперь, тем не менее, мы посмотрим еще один пример пространственной фильтрации - деконволюцию FX!

Деконволюция FX

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

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

Мы начнем с формального определения деконволюции FX:-

Если мы преобразуем сейсмоданные из времени (t) и расстояния (x) в частоту (f) и расстояние (x), временной срез преобразуется в частотный срез, каждый дискрет в преобразованных данных будет комплексным числом с реальными и мнимыми компонентами.

Проявления с похожими углами наклона появятся как синусоидальный комплексный сигнал (cos wt + i sin wt) вдоль данного частотного среза, а, следовательно, они предсказуемы.

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

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

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

Этот рисунок показывает истинный компонент преобразования Фурье каждой трассы в вышеприведенном «разрезе».

Горизонтальная шкала теперь - это частота, а вертикальная - «пространство» (исходные номера трасс).

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

Опять-таки, четко виден тренд.

(Эти два рисунка следовало бы нарисовать как комплексные числа на одном рисунке. Тогда потребовался бы четырехмерный рисунок, что, к сожалению, немного превосходит мои художественные способности!)

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

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

В данном случае, однако, каждый термин уравнения является комплексным числом, а решением является комплексный фильтр!

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

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

Мы посмотрим на эти результаты (вместе с другими типами пространственной фильтрации)

Примеры пространственной фильтрации (1)

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

Я применил к данным некоторые исходные поправки (включая статику), а панель справа показывает часть данных, выделенных (очень слабо) красным.

Выносы в вышеприведенном разрезе приращены на 25 м, давая нам пространственный диапазон (500/25) 20 циклов на километр.

Это анализ FK вышеприведенного ПВ, временные частоты снизу вверх, а пространственные - поперек сверху. Шкала амплитуды в дБ, от белого через голубой до желтого и красного.

Вы можете видеть (временную) полосу, ограниченную природой этих данных (от примерно 10 до 90 Гц), некоторые сильные углы наклона (возможно, первые срывы), и даже некоторый пространственный аляйсинг (желтые окончания линии упирающиеся в букву «W» «Волнового числа»).



Тот же анализ FK после сильных «углов наклона» или «веерного» фильтра.

Мы удалили все проявления с углами наклона более 2,6 мсек / трассу (горизонтальная скорость 25/0,0026 = 9615 м/с).

А вот так вышеприведенный фильтр выглядит в трех измерениях - фильтр треугольной формы (эти фильтры иногда называют «веерными фильтрами») с острыми краями.

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

Обратное преобразование отфильтрованного спектра.



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

Один известный способ снижения эффектов тяжелого фильтра FK - это «обратное смешение» некоторых исходных (не отфильтрованных) данных с конечным результатом.



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

Мы определяем это, как всегда, путем тестирования!

Другая копия нашего анализа FK, на этот раз мы подняли каждую точку в спектре до степени «1,3».

Эта степенная экспонента FK усиливает высокие значения в спектре (постоянные углы наклона ) за счет более слабых значений.

Например, 21,3 равняется примерно 2,46, в то время как 41,3 равняется примерно 6,06, увеличивая разницу от 2:1 до почти 2,46:1.

Еще раз результат обратного преобразования.

Легкое когерентное обогащение лишь слегка усилило когерентные данные (включая поверхностную волну).

Опять-таки, довольно эффективный процесс, который дает некоторое улучшение.



Литература

1. Сагадеев Г.И., Седельников Ю.Е. Пространственно-частотная фильтрация широкополосных сигналов в системах с фазированными антенными решетками. «Инфокоммуникационные технологии» Том 6, № 1, 2008, с. 99 - 104.

2. Федосов В. П. Пространственно-временная обработка сигналов в исследованиях и разработках кафедры теоретических основ радиотехники ТРТУ

- Таганрог: ТРТИ, 1991. - С. 161 - 162.

3. Баскаков С.И. Радиотехнические цепи и сигналы? Учебник для вузов. - М.? Высшая школа, 1988.

4. Даджион Д., Мерсеро Р. Цифровая обработка многомерных сигналов. - М.: Мир, 1988. - 488 с.

5. Дмитриев В.И. Прикладная теория информации: Учебник для вузов. - М.: Высшая школа, 1989.

6. Игнатов В.А. Теория информации и передачи сигналов. - М.: Советское радио, 1979.

7. Купер Дж., Макгиллем А. Вероятностные методы анализа сигналов и систем. - М.: Мир, 1989.

8. Лосев А.К. Линейные радиотехнические цепи: Учебник для вузов. - М.: Высшая школа, 1971.

9. Оппенгейм А.В., Шафер Р.В. Цифровая обработка сигналов. - М.: Связь, 1979. - 416 с.

10. Сергиенко А.Б. Цифровая обработка сигналов. / Учебник для вузов. - СПб.: Питер, 203. - 608 с.

11. Колесник В.Д., Полтырев Г.Ш. Курс теории информации. - М.: Наука, 1982. - 416 с.

12. Даджион Д., Мерсеро Р. Цифровая обработка многомерных сигналов. - М.: Мир, 1988. - 488 с.

Левкович-Маслюк Л, Переберин А. Введение в вейвлет-анализ: Учебный курс. - Москва, ГрафиКон'99, 1999.

14. Алексеев К.А. Очерк "Вокруг CWT". http://support.sibsiu.ru/MATLAB_RU/wavelet/book3/ index.asp.htm.

15. Переберин А.В. О систематизации вейвлет-преобразований. - Вычислительные методы и программирование, 2002, т. 2, с. 15-40.

16. Новиков Л.В. Основы вейвлет-анализа сигналов: Учебное пособие. - СПб, ИАнП РАН, 1999, 152 с.

17. Polikar R. Введение в вейвлет-преобразование. Пер. Грибунина В.Г. - СПб, АВТЭКС. - http://www.autex.spb.ru.

18. prodav.narod.ru - Персональный сайт профессора Давыдова А.В.

Размещено на Allbest.ru