Тема 10: ОПТИМАЛЬНЫЕ ЛИНЕЙНЫЕ ФИЛЬТРЫ.

Как много дел считалось невозможными, пока они не были осуществлены.

Гай Плиний Секунд (философ).

Специалисты в науке подобны старателям. Стоит одному найти крупинку золота, как другие выроют в этом месте котлован. А тема оптимальности, это вообще золотое Эльдорадо, можно копать в любом направлении.

Володя Старцев (Геофизик).

Содержание: 10.1. Модели случайных процессов и шумов. Белый шум. Фильтрация белого шума. 10.2. Критерии построения оптимальных фильтров. Среднее квадратическое отклонение. Амплитудное отношение сигнал/шум. Энергетическое отношение сигнал/шум. 10.3. Фильтр Колмогорова-Винера. Условие оптимальности фильтра. Система линейных уравнений фильтра. Частотная характеристика фильтра. Эффективность фильтра. 10.4. Оптимальные фильтры сжатия сигналов. Частотная характеристика. Условие оптимальности. 10.5. Фильтр прогнозирования. 10.6. Фильтр обнаружения сигналов. Частотная характеристика. Система линейных уравнений. Эффективность фильтра. Согласованный фильтр. Обратный фильтр. 10.7. Энергетический фильтр. Критерий оптимальности. Расчет векторов операторов фильтров. Литература.

10.1. Модели случайных процессов и шумов /л12/.

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

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

Белый шум q(t) можно представлять как случайную по времени (аргументу) последовательность дельта - импульсов d(t_i) со случайными амплитудными значениями a_i:

q(t) = S_i a_id(t-t_i), (10.1.1)

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

R_q(t) = c²d(t), W_q(w) = c²,

где с² = Ns_a², N - число импульсов на интервале Т реализации случайного процесса, s_a² -дисперсия амплитуд импульсов.

Фильтрация белого шума. Если на входе фильтра с импульсным откликом h(t) действует белый шум q(t), то сигнал на выходе фильтра:

g(t) = h(t) _* q(t) = h(t) _* a_id(t-t_i) = a_i h(t-t_i), (10.1.2)

т.е. выходной сигнал будет представлять собой последовательность сигналов импульсной реакции фильтра h(t) с амплитудой a_i, при этом автокорреляционная функция и спектр мощности выходного потока также становятся подобными ФАК и спектру мощности импульсной реакции фильтра и в первом приближении определяются выражениями:

R_g(t) N d_a² R_h(t) = c² R_h(t), (10.1.3)

W_g(w) N d_a² |H(w)|² = c² |H(w)|². (10.1.4)

Этот результат известен как теорема Кэмпбелла.

10.2. Критерии построения оптимальных фильтров.

В практике обработки геофизических данных используются три основных критерия построения оптимальных фильтров: минимум среднего квадратического отклонения профильтрованного сигнала от его действительного или заданного значения, максимум отношения сигнал/шум и максимум энергетического отношения сигнал/шум на выходе фильтра. При анализе и синтезе фильтров используется аддитивная модель входного сигнала: x(k) = s(k)+q(k), где s(k) - полезная составляющая сигнала, q(k) - составляющая помех. Синтез оптимальных фильтров производится с максимальным использованием известной априорной информации как о сигналах, которые необходимо выделить, так и о шумах и помехах. Как правило, используется информация о природе полезного сигнала и шума, об их спектральном составе, о корреляционных и взаимно корреляционных характеристиках. Наличие определенных особенностей (различий) в характеристиках сигнала и шума позволяет реализовать фильтр вообще и оптимальный фильтр в частности. Если такие особенности отсутствуют, постановка задачи становится некорректной.

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

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

y(k) = h(n) _* x(k-n) (10.2.1)

отличается от s(k) на величины e(k) = y(k)-s(k), которые являются абсолютными значениями погрешности воспроизведения полезного сигнала по координатам k. Качество фильтра оценивается средним значением квадрата величины e(k):

. (10.2.2)

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

. (10.2.2')

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

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

r_а = y_экс/s,

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

, (10.2.3)

где y² - средний квадрат амплитуды сигнала в пределах его формы.

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

r = Е_sy/E_qh, (10.2.4)

где Е_sу и E_qh - энергия соответственно сигнала и шума на выходе фильтра.

10.3. Фильтр Колмогорова-Винера.

Условие оптимальности фильтра. Подставим (10.2.1) в выражение (10.2.2') и получим среднее квадратическое отклонение e² формы выходного сигнала y(k) = h(n)_*x(k-n) от оптимальной формы z(k) по всем точкам массива данных (в частном случае z(k) = s(k)):

. (10.3.1)

Минимум функции (10.3.1) определяет значения коэффициентов h(n) оптимального фильтра. Для нахождения их значений продифференцируем выражение (10.3.1) по коэффициентам фильтра и приравняем полученные уравнения нулю. В итоге получаем:

,

где - корреляционная функция входного сигнала, - взаимная корреляционная функция сигналов z(k) и x(k). Отсюда:

h_nR(m-n) = B(m), n = m = 0,1,2, ... , M. (10.3.2)

В краткой форме записи:

h(n)_*R(m-n) = B(m). (10.3.3)

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

. . . . . . . . . . . . .

h_oR(M) + h₁R(M-1)+ h₂R(M-2)+ ....... + h_MR(0) = 0.

При расчете коэффициентов фильтра значение s(0) обычно принимается равным 1.

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

10.5. Фильтры прогнозирования.

Если в правой части уравнения (10.3.3) желаемым сигналом задать входной сигнал со сдвигом на величину kDt, то при этом B(m) = R(m+k) и уравнение принимает вид:

h(n) _* R(m-n) = R(m+k). (10.5.1)

При k > 0 фильтр называется фильтром прогнозирования и вычисляет будущие значения сигнала по его предшествующим значениям. При k < 0 фильтр является фильтром запаздывания. Реализация фильтра заключается в решении соответствующих систем линейных уравнений для каждого заданного значения k. Фильтр может использоваться для интерполяции геофизических полей, в том числе в наперед заданные точки, а также для восстановления утраченных данных.

10.6. Фильтр обнаружения сигналов.

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

Частотная характеристика. Для расчета фильтра требуется задать известную форму полезного сигнала s(k) ó S(w) и функцию автокорреляции или спектр мощности помех R_q(m) ó W_q(w). Полный входной сигнал принимается по аддитивной модели: x(t) = s(t)+q(t). На выходе проектируемого фильтра h(n) ó H(w) для составляющих выходного сигнала имеем:

y(t) = H(w) S(w) exp(jwt) dw, (10.6.1)

s² = |H(w)|² W_q(w) dw, (10.6.2)

где s - средняя квадратическая амплитуда выходной помехи. Значения (10.6.1, 10.6.2) используются для задания критерия максимума отношения сигнал/шум (10.2.3) для произвольной точки t_i:

r = [y(t_i)]²/d². (10.6.3)

Исследование функции (10.6.3) на максимум показывает, что он достигается при частотной характеристике фильтра:

H(w) = exp(-jwt_i) S*(w) / W_q(w), (10.6.4)

Без потери общности можно принять t_i=0:

H(w) = S*(w)/W_q(w) = |S(w)|exp(jj_s(w)) / W_q(w). (10.6.5)

При переходе во временную (координатную) область:

H(w)W_q(w) = S*(w) ó h(n) _* R_q(n-m) = s(-m). (10.6.6)

Система линейных уравнений для расчета фильтра:

h_oR_q(0)+ h₁R_q(1)+ h₂R_q(2)+ h₃R_q(3)+ ...+ h_MR_q(M) = S(-M),

h_oR_q(1)+ h₁R_q(0)+ h₂R_q(1)+ h₃R_q(2)+ ...+ h_MR_q(M-1)= S(-M+1),

h_oR_q(2)+ h₁R_q(1)+ h₂R_q(0)+ h₃R_q(1)+ ...+ h_MR_q(M-2)= S(-M+2),

. . . . . . . . . . . . .

h_oR_q(M)+ h₁R_q(M-1)+ h₂R_q(M-2)+ ..... + h_MR_q(0) = S(0).

Эффективность фильтра. Из выражения (10.6.5) можно видеть, что фильтр имеет максимальный коэффициент передачи на частотах доминирования сигнала и минимальный коэффициент передачи на частотах доминирования помех. Кроме того, фазовая характеристика фильтра j(w) = -j_s(w) обеспечивает синфазность всех частотных составляющих выходного сигнала и соответственно максимальную его амплитуду в заданный момент времени t_i = 0:

y(0)ó S(w) H(w) = = .

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

Согласованный фильтр. При помехах типа белого шума W_q(w) = s² и H(w) = S*(w)/s². Множитель 1/s² не влияет на отношение сигнал/помеха и может быть опущен. Частотная характеристика фильтра определяется только спектром сигнала, при этом:

h(n) = s(-n). (10.6.7)

Фильтр получил название согласованного (по частотной характеристике со спектром сигнала). Он мало эффективен при коротком импульсном или длинном гармоническом сигнале.

Обратный фильтр. Допустим, что помеха имеет такой же частотный состав, что и полезный сигнал, т.е.:

W_q = s²|S(w)|².

Выделение полезного сигнала в таких условиях весьма сомнительно. Тем не менее, определим оптимальный фильтр:

H(w) = S*(w) / [s²|S(w)|²] = 1 / [s²S(w)]. (10.6.8)

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

W_q = a²|S(w)|²+b².

Подставляя это выражение в (10.6.5), с точностью до множителя получаем:

H(w) = S*(w) / [|S(w)|²+g²], (10.6.9)

где g = b/a - отношение дисперсий шума и сигнала. Фильтр стремится к согласованному при больших g, и к обратному (фильтру сжатия) при малых.

10.7. Энергетический фильтр.

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

Критерий оптимальности. Энергия сигнала на выходе фильтра:

E_sh = S_k s_k² = S_k (S_n h_n s_k-n)² = S_k h_k S_n h_n R_s(k-n), (10.7.1)

где R_s- функция автокорреляции сигнала. В векторной форме:

E_sh = . (10.7.2)

Аналогично, выражение для энергии помех на выходе:

E_qh = S_k h_k S_n h_n R_q(k-n) = , (10.7.3)

где R_q - функция автокорреляции помех. При некоррелированной помехе E_qh = s².

Подставим (10.7.2, 10.7.3) в выражение (10.2.4):

r= / . (10.7.4)

Расчет векторов операторов фильтров. Для определения значений вектора продифференцируем r по и приравняем производную к нулю:

.

. (10.7.5)

В системе уравнений (10.7.5) неизвестны собственные значения r матрицы и значения коэффициентов h_n, при этом система имеет N+1 ненулевых решений относительно значений r и соответствующих этим значениям векторов . Для определения коэффициентов фильтра приравнивается к нулю и решается относительно r определитель матрицы , после чего максимальное значение r_max подставляется в (10.7.5) и система уравнений решается относительно коэффициентов h_i вектора . При фильтрации сигнала вектор обеспечивает выделение первой по мощности главной компоненты сигнала, т.е. составляющей сигнала, которая имеет наибольшую энергию и отношение сигнал/шум. В сложных геофизических полях такая компонента, как правило, соответствует региональному фону.

В принципе, расчет может быть продолжен и для других значений r<r_max и определены значения коэффициентов векторов , и т.д., с использованием которых могут выделяться вторая и прочие компоненты сигнала. Наиболее эффективно такой метод используется для разделения сигналов (полей) при некоррелированных помехах. В этом случае корреляционная матрица помех является единичной (единицы по диагонали, остальное - нули) и уравнение (10.7.5) имеет вид:

. (10.7.6)

В развернутой форме:

h_o(R_s(0)-r)+ h₁R_s(1)+ h₂R_s(2)+ h₃R_s(3)+ ...+ h_MR_s(M) = 0,

h_oR_s(1)+ h1(Rs(0)-r)+ h₂R_s(1)+ h₃R_s(2)+ ...+ h_MR_s(M-1) = 0,

h_oR_s(2)+ h₁R_s(1)+ h₂ (Rs(0)-r)+ h₃R_s(1)+ ...+ h_MR_s(M-2) = 0,

. . . . . . . . . . . . .

h_oR_s(M)+ h₁R_s(M-1)+ h₂R_s(M-2)+ ..... + h_M (Rs(0)-r) = 0.

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

В заключение отметим, что расчеты оптимальных фильтров могут производиться с использованием алгоритма Левинсона.