Оптимальный фильтр Калмана—Бьюси. Дискретное обобщенное H-оптимальное управление и фильтрация в линейных непрерывных объектах Бирюков Руслан Сергеевич

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2011 Управление, вычислительная техника и информатика № 3(16)

В.И. Смагин, С.В. Смагин

ФИЛЬТРАЦИЯ В ЛИНЕЙНЫХ ДИСКРЕТНЫХ НЕСТАЦИОНАРНЫХ СИСТЕМАХ С НЕИЗВЕСТНЫМИ ВОЗМУЩЕНИЯМИ

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

Ключевые слова: линейные дискретные нестационарные системы, фильтр Калмана, неизвестные возмущения.

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

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

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

1. Постановка задачи

Рассматривается дискретная система, которая описывается следующими разностными уравнениями:

x(k +1) = A(k)x(k) + f + q(k), x(0) = x0 , (1)

где x(k) e Rn - вектор состояния; A(k) - nxn-матрица; f - неизвестный постоянный вектор; q(k) - белая гауссовская случайная последовательность с характери-

M {q(k)} = 0 , M{q(k)qT (j)} = Q(k)bk ] . (2)

Канал наблюдений имеет вид

y(k) = S (k) x(k) + v(k), (3)

y(k) e R1 - вектор измерений; S(k) - матрица размерности l x n ; v(k) - белая гаус-

совская случайная последовательность ошибок измерений, с характеристиками:

М{у(к)} = 0 , М{д(к)ут (])} = 0, М{у(к)ут (у)} = V(к)81} ; (4)

для матриц (£(к), А(к)) выполняются условия наблюдаемости. Вектор х0 является случайным и не зависит от от процессов д(к) и у(к), при этом

М{х(0)} = хо, М {(х(0) - Хо)(х(0) - Хо)т } = Ро.

Для системы (1) и канала наблюдений (3) требуется синтезировать фильтр, вычисляющий оценку вектора состояния, не использующий оценки неизвестной постоянной составляющей возмущений.

2. Синтез фильтра

Преобразуем дискретную систему (1). Исключаем постоянную составляющую возмущений / из описания объекта посредством вычитания из уравнения (1) такого же уравнения, но со сдвигом на один такт:

х(к) = А (к -1) х(к -1) + / + q(k -1). (5)

В результате получаем следующее уравнение:

х(к +1) = (А (к) + Еп) х(к) - А (к -1) х(к -1) + q(k) - q(к -1). (6)

Расширим пространство состояний системы путем добавления к уравнению (6) тождества х(к) = х(к). Обозначим

X (к) = (хГ„) «к) = (q"k)- ■>). ()

Систему (1) представим в векторно-матричной форме

X(к +1) = А (к)X (к) + д (к), X (0) = Х0, (8)

где А (к) - 2п х 2п -матрица имеет следующую блочную структуру:

К) = (А<кЕ+ Еп -А<0 - 0 ^. (

Случайный вектор X0 = (х^ х-1)т имеет следующие характеристики:

М{X(0)} = X0, М {(X0 -X0)^0 -X0)т} = Р0, (10)

где X0 = (х0т х-)т. Отметим, что здесь дополнительно вводится п-мерный вектор х-1, который является независимым от д(к) и у(к), а характеристики (10) могут быть получены по априорной информации об объекте (1).

Отметим, что в рассмотренной модели (8) процесс д (к) не является белой гауссовской последовательностью, процессы д (к) и д (к -1) будут коррелированны:

Q(k), если у = к,

Q(k -1), если ] = к -1, (11)

0, если 0 < ] < к -1,

М{д (к) д т (у)} =

где Q{k) = |"«к> + «к-1) 0), Q ¡¡). (12)

Представим канал наблюдений для расширенной системы (8) в виде

у(к) = 5 (к) X (к) + v(k), (13)

где 5 (к) = (5(к) 0), v(k) - случайная последовательность ошибок измерений с характеристиками (4).

В качестве уравнения для вычисления оценки вектора состояния расширенной системы выберем уравнение, по своей структуре совпадающее с фильтром Кал-мана:

Х(к +1) = Л(к)Х(к) + К (к)(y(k +1) - 5 (к +1) Л (к)Х(к)), Х(0) = Х0. (14)

Учитывая (8) и (14), получим следующее уравнение для ошибки е(к) = Х(к) - X(к):

е(к + 1) = (Л(к) - К (к) 5 (к +1) Л(к))е(к) + К (к)м(к + 1) + (К (к) 5 (к + 1) - Е2п)д (к). (15)

В силу (11) и (15), матрица Р (к) = М{е(к)ет (к)} определится из следующего разностного уравнения:

Р(к +1) = (Л(к) - К (к)5(к +1)Л (к))Р(к)(Л(к) - К (к)5(к +1)Л (к))т +

+(К (к) 5 (к +1) - Е2п Йк)(К (к)5 (к +1) - Е2п)т + К (к) V (к +1) Кт (к) +

+(Л(к) - К (к)5 (к +1) Л(к))(К (к -1)5 (к) - Е2п) х

х0(к -1)(К (к) 5 (к +1) - Е2 п)т + (К (к)5 (к +1) - Е2п) х

х0(к - 1)(К (к -1) 5 (к) - Е2п)т (Л(к) - К (к)5 (к + 1)Л(к))т, Р(0) = Р0. (16)

Оптимизируемый критерий зададим в виде

3 (к +1) = ЪР (к +1). (17)

Оптимальные коэффициенты передачи фильтра К(к) определяются из условия

Учитывая (17) и правую часть уравнения (16), применяя правила матричного дифференцирования следа от матрицы , получим из условия (18) уравнение для определения матрицы К(к):

Л (к) Р (к) Л(к)т 5 (к + 1)т + К (к) 5 (к +1) Л (к) Р (к) Л(к)т 5 (к + 1)т +

К (к) 5 (к + 1)^(к)5 (к)т - &(к) 5 (к + 1)т - К (к) 5 (к + 1)0(к -1) х х5 (к)т К (к - 1)т Л(к)т 5 (к + 1)т + К (к) 5 (к + 1)0(к -1) Л(к)т 5 (к + 1)т -

К (к) 5 (к +1) Л(к) К (к -1) 5 (к)0(к -1)5 (к + 1)т +

К (к)5 (к +1) Л(к)0(к -1) 5 (к + 1)т + 0(к -1)5 (к)т К (к - 1)т х хЛ(к)т 5(к + 1)т - 0(к -1)Л(к)т 5(к + 1)т -Л(к)0(к -1)5(к + 1)т +

Л(к) К (к -1) 5 (к)0(к -1) 5 (к + 1)т + К (к V (к +1) = 0. (19)

Решение последнего уравнения относительно К(к) дает следующий результат:

К (к) = Р(к)5(к + 1)т (5(к +1)Р(к)5(к + 1)т + V(к +1))-1, (20)

где P(к) = Л (к)P(к)Л(к)т + Q(k - 1)(E2n - S(к)т K(к - 1)т)Л(к)т +

A(k)(Eln - K(к -1)5(к))Q(k -1) + Q(k). (21)

Отметим, что для вычисления коэффициентов передачи (20), в силу (21), необходимо задать начальные значения коэффициентов K(-1).

Подставив в уравнение (16) выражение для оптимального коэффициента передачи (20), получим уравнение

P(к +1) = (E2n - K (к)S (к +1))P(к), P(0) = Р0. (22)

Основной результат сформулируем в виде теоремы, учитывая симметричность и блочное представление матриц P(к) и P(к):

P(к) = f p (к) (к) 1, P(k) = f p1(к) p2T (к) 1, (23)

IР 2 (к) p з(к)) У Р2(к) Рз(к))

блочные структуры матриц Л(к), Q(k), Q(k), S(к) и представление матрицы K(к) в виде

k (к >=(K%). <24)

Теорема. Пусть процесс с неизвестным постоянным возмущением определяется уравнениями (1) и канал наблюдений имеет вид (3). Тогда оптимальный алгоритм фильтрации определится следующими разностными уравнениями: x(k +1) = (A (k) + En) x(k) - Л (к -1) x(k -1) + K1 (к)(y(k +1) -

S (к +1)[(Л (к) + En) x(k) - Л (к -1) x(k -1)] (25)

с начальными условиями

x(0) = x0, x(1) = M{x(1)} = x . (26)

Матрица Kx (k) в (25) определяется по формуле

K (к) = рх (к) S (к + 1)т (S (к +1) Р (к) S (к + 1)т + V (к +1))-1, (27)

где матрица р (к) в^гчисляется из системы уравнений

Р(к) = (Л (к) + En)р1(к)(Л (к) + En)т - Л(к -1) Р2(к)(Л (к) + En)т -

-(Л(к) + En) рТ (к) Л(к - 1)т + Л (к -1) р3(к) Л (к - 1)т + Q(к -1) S (к)т Kj(k - 1)т х х(Л (к) + En)т - Q(k -1) S (к)т K2 (к - 1)т Лт (к -1) +

+(Л(к) + En) Kj (k -1) S (k) Q(k -1) - Л(k -1) K2 (k -1) S (k) х xQ (k -1) - (Л(к) + En)Q(k -1) - Q(k -1)(Л(к) + En)T + Q(k) + Q(k -1),

Р2 (к) = #(k)(Л(к) + En)T - p2 (к)Л(к - 1)T +

K1 (k - 1)S(k)Q(k -1) - Q(k -1), ^3 (k) = p1 (k),

Р1 (k +1) = (En - K (k)S(k +1))p (к), Р1 (0) = Р10,

Р2 (k +1) = - K 2 (k) S (k +1) P (k) + p 2 (к), Р2 (0) = Р20,

Р3 (k +1) = -K2 (k)S(k +1) p2 (к) + Р3 (к), Р3 (0) = p^ ,

K2 (k) = p 2 (k)S (k + 1)T (S (k +1) p (k) S (к + 1)T + V (к +1))-1. (28)

В (28) начальные условия р10, р2 0, р3 0, являются соответствующими блоками матрицы Р0. Отметим, что для выполнения расчетов в (28) необходимо задать начальные условия для КД-1) и К2(-1).

Замечание. Управляемый объект

х(к +1) = А(к)х(к) + В(к)и(к) + / + д(к), х(0) = х0, (29)

при исключении неизвестного постоянного возмущения / объекта, необходимо преобразовать к виду, который будет отличаться от (8) одним слагаемым:

X (к +1) = А (к) X (к) + В(к)(и(к) - и(к -1) + д (к), X (0) = Х0, (30)

где матрица А(к) приведена в формуле (9), д (к) имеет характеристики (11), (12). В (30) матрица В (к) имеет вид

В(к) =(В0к)). (31)

Тогда уравнения фильтра будут следующими:

Х(к +1) = (А (к) + Еп) Х(к) - А (к -1) Х(к -1) + В(к)(и(к) - и (к -1)) + К1 (к)(у(к +1) --Б(к +1)[(А(к) + Еп)Х(к) - А(к -1)Х(к -1) + В(к)(и(к) -и(к -1))], (32)

с начальными условиями (26), а матрица К1(к) определяется в соответствии с (27) и (28).

3. Результаты вычислительного эксперимента

Рассмотрим применение алгоритма фильтрации для модели второго порядка вида (1) и канала наблюдений (3) со следующими значениями параметров:

(0 1 А _ (0,01 0 А ТЛ

() = ^0,05 0,925 + 0,Ыи(0,01к)) ’ ® = [ 0 0,02} ’ = , ’

х = (1 1); Х0 =(иР0 =(100 100} (

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

у(к +1) = 5Х(к +1) + у(к +1) = £А(к)Х(к) + Б/ + 5д(к) + у(к +1). (34)

Рекуррентные уравнения оценивания неизвестного вектора / имеют вид /(к +1) = /(к) + К/ (к)(у(к +1) - £А(к)Х(к) - Б/(к)), Д0) = /0,

Кг (к) = Рг (к) Бт (БРГ (к) Бт + ^т + V)-1,

Р/ (к +1) = (Е2 - К/ (к)Б)Р/ (к), Р/ (0) = Р/0, (35)

где М{/} = /0, М{(/ - /0)(/ - /0)т } = Р/0. (36)

Оценка вектора состояния для объекта с неизвестным постоянным входом задается уравнением:

х(к +1) = Л(к)х(к) + /(к) + Кх (к)(у(к +1) - БЛ(к)х(к) - Б/(к)) , (37)

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

(01 Р =Г1,0 01 ,0 ] , Л { 0 1,0 ].

Применение расширенного фильтра Калмана для данного примера (в этом случае уравнение (1) расширяется путем добавления уравнения /(к+1) = /(к)) приводит к необходимости построения фильтра Калмана для дискретной системы со следующими матрицами динамики, канала наблюдений и интенсивностей аддитивных возмущений:

(Л (к) К21 ((б 01

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

п > т и I > т. (40)

В неизвестное возмущение определяется в виде / = Ой, где й - неизвестный т-мерный вектор, О - п х т -известная матрица. В рассмотренном примере О = Е2, п = 2, т = 2, I = 1, а это означает, что условия (40) не выполняются.

Применение алгоритма фильтрации исследовалось также для неизвестного переменного возмущения с тремя возможными значениями компонент вектора /:

1, если 0 < к < 9,

/1(к) = /2(к) = < -1, если 9 < к < 25,

1, если 25 < к < 50.

На рис. 1 приведены реализации процессов и их оценок для трех сравниваемых фильтров. Отметим, что при реализации алгоритма фильтрации (25), начальные значения К1(-1) и К2(-1) задавались нулевые.

Рис. 1. Реализации процессов и оценок (1 - реализация х(к); 2 - оценка, построенная по алгоритму (25); 3 - оценка, построенная по двухэтапному алгоритму; 4 - оценка для расширенного фильтра Калмана)

На рис. 2 приведены ошибки оценивания компонент вектора состояния.

Рис. 2. Графики ошибок фильтрации (1 - ошибка для оценки, построенной по алгоритму (25); 2 - ошибка для оценки, построенной по двухэтапному алгоритму; 3 - ошибка для расширенного фильтра Калмана)

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

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

Средние значения среднеквадратических ошибок для компонент вектора состояния

Алгоритм (25) Двухэтапный алгоритм Расширенный фильтр Калмана

е1>Ср = 0,0912 е1,ср = 0,3128 Єі,ср = 0,4103

Є2,ср = 0,0945 е2,ср = 0,2917 е2,ср = 0,4296

Заключение

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

ЛИТЕРАТУРА

1. Astrom K., EykhoffP. System identification. A survey // Automatica. 1971. V. 7. P. 123-162.

2. FriedlandB. Treatment of bias in recursive filtering // IEEE Trans. on Automat. Contr. 1969. V. AC-14. P. 359-367.

3. Chen J., Patton R. J. Optimal filtering and robust fault diagnosis of stochastic systems with unknown disturbances // IEE Proc. Control Theory Appl. 1996. V. 143. P. 31-36.

4. Darouach M., Zasadzinski M. Unbiased minimum variance estimation for systems with unknown exogenous inputs // Automatica. 1997. V. 33. P. 717-719.

5. Darouach M., Zasadzinski M., Xu S. J. Full-order observers for linear systems with unknown inputs // IEEE Trans. on Automat. Contr. 1999. V. AC-39. P. 606.

6. Gillijns S., Moor B. Unbiased minimum-variance input and state estimation for linear discrete-time systems // Automatica. 2007. V. 43. P. 111-116.

7. Hou M., Patton R. Optimal filtering for systems with unknown inputs // IEEE Trans. on Automat. Contr. 1998. V. AC-43. P. 445-449.

8. Hsieh C.-S. A unified solution to unbiased minimum-variance estimation for systems with unknown inputs // Proc.17th World Congress The International Federation of Automatic Control. Seoul. Korea. July 6 - 11, 2008. P. 14502-14509.

9. Hsieh C.-S. Robust two-stage Kalman filters for systems with unknown inputs // IEEE Trans. on Automat. Contr. 2000. V. AC-45. P. 2374-2378.

10. Hsieh C.-S. Extension of the optimal unbiased minimum-variance filter for systems with unknown inputs // Proc. 15th IEEE International Workshop on Nonlinear Dynamics of Electronic Systems. Tokushima. Japan. 2007. P. 217-220.

11. Hsieh C.-S. Robust parameterized minimum variance filtering for uncertain systems with unknown inputs // Proc. American control conference. New York. 2007. P. 5118-5123.

12. Kalman R.E., Busy R. A new results in linear filtering and prediction theory // Trans. ASME J. Basic Engr. 1961. V. 83. P. 95-108.

13. Браммер К., ЗиффлингГ. Фильтр Калмана - Бьюси. М.: Наука, 1972. 200 с.

14. Пугачев В.С., Синицин И.Н. Стохастические дифференциальные уравнения М.: Наука, 1990. 630 с.

15. Смагин С.В. Фильтрация в линейных дискретных системах с неизвестными возмущениями // Автометрия. 2009. Т. 45. № 6. C. 29-37.

16. Амосов А.А., Колпаков В.В. Скалярно-матричное дифференцирование и его применение к конструктивным задачам теории связи // Проблемы передачи информации. 1972. № 1. С. 3-15.

Смагин Валерий Иванович

Смагин Сергей Валерьевич

Томский государственный университет

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА 2011 Управление, вычислительная техника и информатика № 3(16) УДК 517.511 В.И. Смагин, С.В. Смагин ФИЛЬТРАЦИЯ В ЛИНЕЙНЫХ ДИСКРЕТНЫХ НЕСТАЦИОНАРНЫХ СИСТЕМАХ С НЕИЗВЕСТНЫМИ ВОЗМУЩЕНИЯМИ Рассматривается алгоритм синтеза оптимального фильтра, определяющего оценку вектора состояния дискретной линейной нестационарной динамической системы с аддитивными возмущениями, содержащими неизвестную постоянную составляющую. Приводятся результаты вычислительного эксперимента. Ключевые слова: линейные дискретные нестационарные системы, фильтр Калмана, неизвестные возмущения. В работах многих авторов большое внимание уделяется разработке алгоритмов калмановской фильтрации для класса систем с неизвестными аддитивными возмущениями и параметрами, которые могут использоваться в качестве моделей реальных физических систем, моделей объектов с неизвестными сбоями. Известные методы вычисления оценок вектора состояния базируются на алгоритмах, использующих оценки неизвестного возмущения . В работах рассматриваются алгоритмы расширения пространства состояний (к основной модели объекта добавляется модель ненаблюдаемого возмущения) и алгоритм двухэтапной фильтрации, уменьшающий вычислительные затраты за счет декомпозиции задачи. В работах изучены алгоритмы рекуррентной оптимальной фильтрации, использующие оценки неизвестного возмущения, имеющие достаточно жесткие условия их разрешимости. В настоящей работе для дискретного нестационарного объекта с неизвестной постоянной составляющей возмущений предлагается метод оптимальной фильтрации, не использующий оценки неизвестного возмущения. Метод базируется на преобразовании модели и сведении к задаче линейной калмановской фильтрации . В настоящей статье обобщаются результаты на случай решения задачи для нестационарного дискретного объекта. 1. Постановка задачи Рассматривается дискретная система, которая описывается следующими разностными уравнениями: x(k + 1) = A(k) x(k) + f + q (k), x(0) = x0 , (1) где x(k) ∈ R n – вектор состояния; A(k) – n×n-матрица; f – неизвестный постоянный вектор; q(k) – белая гауссовская случайная последовательность с характеристиками M {q (k)} = 0 , M{q(k)q Τ (j)} = Q(k)δk , j . (2) Канал наблюдений имеет вид y (k) = S (k) x(k) + v(k) , (3) y (k) ∈ R l – вектор измерений; S(k) – матрица размерности l × n ; v(k) – белая гаус- В.И. Смагин, С.В. Смагин 44 совская случайная последовательность ошибок измерений, с характеристиками: M{v(k)} = 0 , M{q (k)v Τ (j)} = 0 , M{v(k)v Τ (j)} = V (k)δi , j ; (4) для матриц (S(k), A(k)) выполняются условия наблюдаемости. Вектор x0 является случайным и не зависит от от процессов q(k) и v(k), при этом M{x(0)} = x0 , M {(x(0) − x0)(x(0) − x0)Τ } = P0 . Для системы (1) и канала наблюдений (3) требуется синтезировать фильтр, вычисляющий оценку вектора состояния, не использующий оценки неизвестной постоянной составляющей возмущений. 2. Синтез фильтра Преобразуем дискретную систему (1). Исключаем постоянную составляющую возмущений f из описания объекта посредством вычитания из уравнения (1) такого же уравнения, но со сдвигом на один такт: x(k) = A(k − 1) x(k − 1) + f + q(k − 1) . (5) В результате получаем следующее уравнение: x(k + 1) = (A(k) + En) x(k) − A(k − 1) x(k − 1) + q (k) − q(k − 1) . (6) Расширим пространство состояний системы путем добавления к уравнению (6) тождества x(k) = x(k) . Обозначим x(k) ⎞ ⎛ q(k) − q(k − 1) ⎞ . X (k) = ⎛⎜ ⎟ ⎟ , q (k) = ⎜ 0 ⎝ ⎠ ⎝ x(k − 1) ⎠ Систему (1) представим в векторно-матричной форме X (k + 1) = A(k) X (k) + q (k), X (0) = X 0 , (7) (8) где А(k) – 2n × 2n -матрица имеет следующую блочную структуру: ⎛ A(k) + En A(k) = ⎜ En ⎝ − A(k − 1) ⎞ ⎟. 0 ⎠ (9) Случайный вектор X 0 = (x0Τ x−Τ1)Τ имеет следующие характеристики: M{ X (0)} = X 0 , M {(X 0 − X 0)(X 0 − X 0)Τ } = P0 , (x0Τ (10) x−Τ1)Τ где X 0 = . Отметим, что здесь дополнительно вводится n-мерный вектор x−1 , который является независимым от q(k) и v(k) , а характеристики (10) могут быть получены по априорной информации об объекте (1). Отметим, что в рассмотренной модели (8) процесс q (k) не является белой гауссовской последовательностью, процессы q (k) и q (k − 1) будут коррелированны: если j = k, ⎧ Q (k), ⎪ M{q (k)q (j)} = ⎨Q (k − 1), если j = k − 1, ⎪ 0, если 0 ≤ j < k − 1, ⎩ (11) Q(k) + Q(k − 1) 0 ⎞ ⎛ −Q(k − 1) 0 ⎞ . Q(k) = ⎛⎜ ⎟ , Q (k − 1) = ⎜ 0 0 0 0 ⎟⎠ ⎝ ⎠ ⎝ (12) Τ где Фильтрация в линейных дискретных нестационарных системах 45 Представим канал наблюдений для расширенной системы (8) в виде y (k) = S (k) X (k) + v(k) , (13) где S (k) = (S (k) 0) , v(k) − случайная последовательность ошибок измерений с характеристиками (4). В качестве уравнения для вычисления оценки вектора состояния расширенной системы выберем уравнение, по своей структуре совпадающее с фильтром Калмана: Xˆ (k + 1) = A(k) Xˆ (k) + K (k)(y (k + 1) − S (k + 1) A(k) Xˆ (k)) , Xˆ (0) = X . (14) 0 Учитывая (8) и (14), получим следующее уравнение для ошибки e(k) = Xˆ (k) − X (k) : e(k + 1) = (A(k) − K (k) S (k + 1) A(k))e(k) + K (k)v(k + 1) + (K (k) S (k + 1) − E2 n)q (k) . (15) В силу (11) и (15), матрица P (k) = M{e(k)eΤ (k)} определится из следующего разностного уравнения: P (k + 1) = (A(k) − K (k) S (k + 1) A(k)) P (k)(A(k) − K (k) S (k + 1) A(k))Τ + +(K (k) S (k + 1) − E2 n)Q (k)(K (k) S (k + 1) − E2 n)Τ + K (k)V (k + 1) K Τ (k) + +(A(k) − K (k) S (k + 1) A(k))(K (k − 1) S (k) − E2 n) × ×Q (k − 1)(K (k) S (k + 1) − E2 n)Τ + (K (k) S (k + 1) − E2 n) × ×Q (k − 1)(K (k − 1) S (k) − E2 n)Τ (A(k) − K (k) S (k + 1) A(k))Τ , P (0) = P0 . (16) Оптимизируемый критерий зададим в виде J (k + 1) = trP (k + 1) . (17) Оптимальные коэффициенты передачи фильтра K(k) определяются из условия dJ (k + 1) =0. (18) dK (k) Учитывая (17) и правую часть уравнения (16), применяя правила матричного дифференцирования следа от матрицы , получим из условия (18) уравнение для определения матрицы K(k): − A(k) P (k) A(k)Τ S (k + 1)Τ + K (k) S (k + 1) A(k) P (k) A(k)Τ S (k + 1)Τ + + K (k) S (k + 1)Q (k) S (k)Τ − Q (k) S (k + 1)Τ − K (k) S (k + 1)Q (k − 1) × ×S (k)Τ K (k − 1)Τ A(k)Τ S (k + 1)Τ + K (k) S (k + 1)Q (k − 1) A(k)Τ S (k + 1)Τ − − K (k) S (k + 1) A(k) K (k − 1) S (k)Q (k − 1) S (k + 1)Τ + + K (k) S (k + 1) A(k)Q (k − 1) S (k + 1)Τ + Q (k − 1) S (k)Τ K (k − 1)Τ × × A(k)Τ S (k + 1)Τ − Q (k − 1) A(k)Τ S (k + 1)Τ − A(k)Q (k − 1) S (k + 1)Τ + + A(k) K (k − 1) S (k)Q (k − 1) S (k + 1)Τ + K (k)V (k + 1) = 0 . (19) Решение последнего уравнения относительно K(k) дает следующий результат: K (k) = P (k) S (k + 1)Τ (S (k + 1) P (k) S (k + 1)Τ + V (k + 1)) −1 , (20) 46 В.И. Смагин, С.В. Смагин где P (k) = A(k) P (k) A(k)Τ + Q (k − 1)(E2 n − S (k)Τ K (k − 1)Τ) A(k)Τ + + A(k)(E2 n − K (k − 1) S (k))Q (k − 1) + Q (k) . (21) Отметим, что для вычисления коэффициентов передачи (20), в силу (21), необходимо задать начальные значения коэффициентов K(−1). Подставив в уравнение (16) выражение для оптимального коэффициента передачи (20), получим уравнение P (k + 1) = (E2 n − K (k) S (k + 1)) P(k) , P (0) = P0 . (22) Основной результат сформулируем в виде теоремы, учитывая симметричность и блочное представление матриц P (k) и P (k) : ⎛ p (k) P(k) = ⎜ 1 ⎝ p2 (k) ⎛ p (k) p2Τ (k) ⎞ , P (k) = ⎜ 1 p3 (k) ⎟⎠ ⎝ p2 (k) p2Τ (k) ⎞ , p3 (k) ⎟⎠ (23) блочные структуры матриц A(k), Q(k), Q (k), S (k) и представление матрицы K (k) в виде ⎛ K (k) ⎞ K (k) = ⎜ 1 ⎟ . (24) ⎝ K 2 (k) ⎠ Теорема. Пусть процесс с неизвестным постоянным возмущением определяется уравнениями (1) и канал наблюдений имеет вид (3). Тогда оптимальный алгоритм фильтрации определится следующими разностными уравнениями: xˆ (k + 1) = (A(k) + En) xˆ (k) − A(k − 1) xˆ (k − 1) + K1 (k)(y (k + 1) − − S (k + 1)[(A(k) + En) xˆ (k) − A(k − 1) xˆ (k − 1)] (25) с начальными условиями xˆ(0) = x0 , xˆ(1) = M{x(1)} = x1 . Матрица K1 (k) в (25) определяется по формуле (26) K1 (k) = p1 (k) S (k + 1)Τ (S (k + 1) p1 (k) S (k + 1)Τ + V (k + 1)) −1 , где матрица p1 (k) вычисляется из системы уравнений (27) p1 (k) = (A(k) + En) p1 (k)(A(k) + En)Τ − A(k − 1) p2 (k)(A(k) + En)Τ − −(A(k) + En) p2Τ (k) A(k − 1)Τ + A(k − 1) p3 (k) A(k − 1)Τ + Q(k − 1) S (k)Τ K1 (k − 1)Τ × ×(A(k) + En)Τ − Q(k − 1) S (k)Τ K 2 (k − 1)Τ AΤ (k − 1) + +(A(k) + En) K1 (k − 1) S (k)Q(k − 1) − A(k − 1) K 2 (k − 1) S (k) × ×Q(k − 1) − (A(k) + En)Q(k − 1) − Q(k − 1)(A(k) + En)Τ + Q(k) + Q(k − 1) , p2 (k) = p1 (k)(A(k) + En)Τ − p2Τ (k) A(k − 1)Τ + + K1 (k − 1) S (k)Q(k − 1) − Q(k − 1) , p3 (k) = p1 (k) , p1 (k + 1) = (En − K1 (k) S (k + 1)) p1 (k) , p1 (0) = p1,0 , p2 (k + 1) = − K 2 (k) S (k + 1) p1 (k) + p2 (k) , p2 (0) = p2,0 , p3 (k + 1) = − K 2 (k) S (k + 1) p2Τ (k) + p3 (k) , p3 (0) = p3,0 , K 2 (k) = p2 (k) S (k + 1)Τ (S (k + 1) p1 (k) S (k + 1)Τ + V (k + 1)) −1 . (28) Фильтрация в линейных дискретных нестационарных системах 47 В (28) начальные условия p1,0 , p2,0 , p3,0 , являются соответствующими блоками матрицы P0 . Отметим, что для выполнения расчетов в (28) необходимо задать начальные условия для K1 (−1) и K 2 (−1) . Замечание. Управляемый объект x(k + 1) = A(k) x(k) + B(k)u (k) + f + q(k), x(0) = x0 , (29) при исключении неизвестного постоянного возмущения f объекта, необходимо преобразовать к виду, который будет отличаться от (8) одним слагаемым: X (k + 1) = A(k) X (k) + B (k)(u (k) − u (k − 1) + q (k), X (0) = X 0 , (30) где матрица A(k) приведена в формуле (9), q (k) имеет характеристики (11), (12). В (30) матрица B (k) имеет вид B (k) ⎞ B (k) = ⎛⎜ ⎟. ⎝ 0 ⎠ Тогда уравнения фильтра будут следующими: (31) xˆ (k + 1) = (A(k) + En) xˆ (k) − A(k − 1) xˆ (k − 1) + B(k)(u (k) − u (k − 1)) + K1 (k)(y (k + 1) − − S (k + 1)[(A(k) + En) xˆ (k) − A(k − 1) xˆ (k − 1) + B(k)(u (k) − u (k − 1))] , (32) с начальными условиями (26), а матрица K1 (k) определяется в соответствии с (27) и (28). 3. Результаты вычислительного эксперимента Рассмотрим применение алгоритма фильтрации для модели второго порядка вида (1) и канала наблюдений (3) со следующими значениями параметров: 0 1 0 ⎞ ⎞ ; Q = ⎛ 0, 01 ; V = 0,9 ; A(k) = ⎛⎜ ⎟ ⎜ 0 0, 02 ⎟⎠ ⎝ ⎝ 0, 05 0,925 + 0,1sin(0, 01k) ⎠ 1, 0 1, 0 0 ⎞ S = (1 1) ; x0 = ⎛⎜ ⎞⎟ ; P0 = ⎛⎜ (33) ⎟. ⎝ 1,5 ⎠ ⎝ 0 1, 0 ⎠ Вычисление оценок вектора x(k) можно выполнить, используя двухэтапный алгоритм фильтрации . Модель измерений в этом случае с учетом (1) представляется в виде y (k + 1) = Sx(k + 1) + v(k + 1) = SA(k) x(k) + Sf + Sq(k) + v(k + 1) . (34) Рекуррентные уравнения оценивания неизвестного вектора f имеют вид fˆ (k + 1) = fˆ (k) + K (k)(y (k + 1) − SA(k) xˆ (k) − Sfˆ (k)) , fˆ (0) = f , 0 f Τ Τ Τ −1 K f (k) = Pf (k) S (SPf (k) S + SQS + V) , где Pf (k + 1) = (E2 − K f (k) S) Pf (k), Pf (0) = Pf0 , (35) M{ f } = f 0 , M{(f − f 0)(f − f 0)Τ } = Pf0 . (36) В.И. Смагин, С.В. Смагин 48 Оценка вектора состояния для объекта с неизвестным постоянным входом задается уравнением: xˆ (k + 1) = A(k) xˆ (k) + fˆ (k) + K (k)(y (k + 1) − SA(k) xˆ (k) − Sfˆ (k)) , (37) x где матрица K x (k) определяет коэффициенты передачи фильтра Калмана. При моделировании используем 0 1, 0 0 ⎞ f 0 = ⎛⎜ ⎞⎟ , Pf0 = ⎛⎜ (38) ⎟. ⎝0⎠ ⎝ 0 1, 0 ⎠ Применение расширенного фильтра Калмана для данного примера (в этом случае уравнение (1) расширяется путем добавления уравнения f(k+1) = f(k)) приводит к необходимости построения фильтра Калмана для дискретной системы со следующими матрицами динамики, канала наблюдений и интенсивностей аддитивных возмущений: Q 0⎞ ⎛ A(k) E2 ⎞ , (S 0) , ⎛⎜ (39) ⎟. ⎜ 0 E2 ⎟⎠ ⎝ 0 0⎠ ⎝ Использование в данном примере методов, описанных в работах , невозможно в силу невыполнения условий существования оптимальных оценок неизвестного входного вектора : n≥m и l≥m. (40) В неизвестное возмущение определяется в виде f = Gd , где d – неизвестный m-мерный вектор, G – n × m -известная матрица. В рассмотренном примере G = E2 , n = 2 , m = 2, l = 1 , а это означает, что условия (40) не выполняются. Применение алгоритма фильтрации исследовалось также для неизвестного переменного возмущения с тремя возможными значениями компонент вектора f: ⎧ 1, если 0 ≤ k ≤ 9, ⎪ f1 (k) = f 2 (k) = ⎨ −1, если 9 < k < 25, ⎪ 1, если 25 ≤ k ≤ 50. ⎩ На рис. 1 приведены реализации процессов и их оценок для трех сравниваемых фильтров. Отметим, что при реализации алгоритма фильтрации (25), начальные значения K1 (−1) и K 2 (−1) задавались нулевые. x1(k) x1(k) x2(k) x2(k) 2 10 0 –10 0 3 4 20 30 40 k –10 0 4 1 0 1 10 3 10 2 10 20 30 40 k Рис. 1. Реализации процессов и оценок (1 – реализация x(k); 2 – оценка, построенная по алгоритму (25); 3 – оценка, построенная по двухэтапному алгоритму; 4 – оценка для расширенного фильтра Калмана) Фильтрация в линейных дискретных нестационарных системах 49 На рис. 2 приведены ошибки оценивания компонент вектора состояния. e1(k) 4 2 e2(k) 4 3 1 0 –2 –4 –6 0 2 2 3 1 0 2 –2 10 20 30 40 k –4 0 10 20 30 40 k Рис. 2. Графики ошибок фильтрации (1 – ошибка для оценки, построенной по алгоритму (25); 2 – ошибка для оценки, построенной по двухэтапному алгоритму; 3 – ошибка для расширенного фильтра Калмана) Как видно из рисунков для рассмотренного примера, качество оценок, полученных с помощью фильтра (25), лучше, чем для двухэтапного алгоритма фильтрации и расширенного фильтра Калмана, использующих оценки неизвестного возмущения. Отметим также, что для алгоритма фильтрации (25) не нужно задавать априорную информацию о характеристиках распределения начальных значений f 0 и Pf0 . Ниже, в таблице, приведены средние значения среднеквадратических ошибок оценивания для трех рассматриваемых методов, рассчитанных по 50 реализациям. Как видно из таблицы, предложенный метод фильтрации (25) обеспечивает среднюю ошибку в 3 – 4 раза меньшую, чем другие методы. Средние значения среднеквадратических ошибок для компонент вектора состояния Алгоритм (25) e1,ср = 0,0912 Двухэтапный алгоритм e1,ср = 0,3128 Расширенный фильтр Калмана e1,ср = 0,4103 e2,ср = 0,0945 e2,ср = 0,2917 e2,ср = 0,4296 Заключение Разработан алгоритм синтеза дискретного оптимального нестационарного фильтра для объекта, возмущения которого содержат неизвестную постоянную составляющую. Алгоритм построен на основе расширения пространства состояния и исключения из модели неизвестной составляющей. В отличие от классического фильтра Калмана, предложенный фильтр использует рекуррентные оценки, построенные на двух предыдущих тактах. Как показали результаты вычислительного эксперимента, алгоритм может быть применен для кусочно-постоянной неизвестной аддитивной составляющей возмущений. ЛИТЕРАТУРА 1. Astrom K., Eykhoff P. System identification. A survey // Automatica. 1971. V. 7. P. 123−162. 2. Friedland B. Treatment of bias in recursive filtering // IEEE Trans. on Automat. Contr. 1969. V. AC-14. P. 359−367. 3. Chen J., Patton R. J. Optimal filtering and robust fault diagnosis of stochastic systems with unknown disturbances // IEE Proc. Control Theory Appl. 1996. V. 143. P. 31–36. 50 В.И. Смагин, С.В. Смагин 4. Darouach M., Zasadzinski M. Unbiased minimum variance estimation for systems with unknown exogenous inputs // Automatica. 1997. V. 33. P. 717–719. 5. Darouach M., Zasadzinski M., Xu S. J. Full-order observers for linear systems with unknown inputs // IEEE Trans. on Automat. Contr. 1999. V. AC-39. P. 606. 6. Gillijns S., Moor B. Unbiased minimum-variance input and state estimation for linear discrete-time systems // Automatica. 2007. V. 43. P. 111–116. 7. Hou M., Patton R. Optimal filtering for systems with unknown inputs // IEEE Trans. on Automat. Contr. 1998. V. AC-43. P. 445–449. 8. Hsieh C.-S. A unified solution to unbiased minimum-variance estimation for systems with unknown inputs // Proc.17th World Congress The International Federation of Automatic Control. Seoul. Korea. July 6 – 11, 2008. P. 14502–14509. 9. Hsieh C.-S. Robust two-stage Kalman filters for systems with unknown inputs // IEEE Trans. on Automat. Contr. 2000. V. AC-45. P. 2374–2378. 10. Hsieh C.-S. Extension of the optimal unbiased minimum-variance filter for systems with unknown inputs // Proc. 15th IEEE International Workshop on Nonlinear Dynamics of Electronic Systems. Tokushima. Japan. 2007. P. 217–220. 11. Hsieh C.-S. Robust parameterized minimum variance filtering for uncertain systems with unknown inputs // Proc. American control conference. New York. 2007. P. 5118–5123. 12. Kalman R.E., Busy R. A new results in linear filtering and prediction theory // Trans. ASME J. Basic Engr. 1961. V. 83. P. 95–108. 13. Браммер К., Зиффлинг Г. Фильтр Калмана – Бьюси. М.: Наука, 1972. 200 с. 14. Пугачев В.С., Синицин И.Н. Стохастические дифференциальные уравнения М.: Наука, 1990. 630 с. 15. Смагин С.В. Фильтрация в линейных дискретных системах с неизвестными возмущениями // Автометрия. 2009. Т. 45. № 6. C. 29−37. 16. Амосов А.А., Колпаков В.В. Скалярно-матричное дифференцирование и его применение к конструктивным задачам теории связи // Проблемы передачи информации. 1972. № 1. С. 3−15. Смагин Валерий Иванович Смагин Сергей Валерьевич Томский государственный университет E-mail: [email protected]; [email protected] Поступила в редакцию 6 декабря 2010 г.

Нахождение оптимального фильтра Винера основывалось на использовании интегрального уравнения Винера - Хопфа, при решении которого стационарные случайные процессы рассматривались в частотной области. В 1960 г. Р. Калман и Р. Бьюси рассмотрели проблему линейной фильтрации во временной области и, используя концепцию «пространства состояний», предложили новый эффективный метод синтеза оптимальных систем по критерию минимума математического ожидания квадрата случайной ошибки, применимый как для стационарных, так и для нестационарных марковских случайных процессов. Так как в основе используемой Калманом и Бьюси концепции «пространства состояний» лежит предположение о том, что случайный процесс является марковским, то их подход к синтезу оптимальных линейных систем иногда называют марковской теорией оптимальной линейной фильтрации.

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

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

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

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

Для вывода уравнения одномерного оптимального фильтра Калмана - Бьюси существенным является то, что случайный процесс должен быть сначала представлен дифференциальным уравнением первого порядка следующего вида:

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

Корреляционные функции нестационарных случайных процессов имеют вид

где - непрерывные, непрерывно дифференцируемые функции времени, причем

В частном случае для стационарных случайных процессов их корреляционные функции

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

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

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

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

Прежде чем определить функции находят некоторую функцию времени равную математическому ожиданию квадрата случайной ошибки (дисперсии ошибки):

она определяется как решение следующего дифференциального уравнения Риккати:

Для решения (9.143) нужно знать начальное значение при Обычно поэтому

После нахождения функции определяют функцию по формуле

и функцию по формуле

Наиболее сложным этапом синтеза оптимальных фильтров методом Калмана - Бьюси является решение уравнения Риккати (9.143). В общем случае оно требует применения ЭВМ.

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

Учитывая (9.146), уравнение оптимального фильтра Калмана-Бьюси иногда записывают в следующем виде:

Дифференциальному уравнению (9.140) соответствует структурная схема оптимального фильтра, показанная на рис. 9.19, а; дифференциальному уравнению (9.147) соответствует структурная схема, показанная на рис. 9.19, б. Таким образом, оптимальный фильтр Калмана-Бьюси можно рассматривать как некоторую динамическую систему с обратной связью, имеющую структурную схему, приведенную либо на рис. 9.19, а, либо на рис. 9.19, б. Естественно, что обе эти структурные схемы эквивалентны.

Для нестационарных случайных процессов функции зависят от времени и оптимальный фильтр Калмана-Бьюси получается нестационарным.

Для стационарных случайных процессов функции а также в установившемся режиме функции не зависят от времени, поэтому оптимальный фильтр Калмана-Бьюси в этом случае является стационарным, определяемым дифференциальным уравнением с постоянными коэффициентами

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

Естественно, что для стационарных процессов результаты, полученные методом Калмана-Бьюси и методом Винера, совпадают. Уравнение (9.148), полученное во временной области, эквивалентно оптимальному фильтру Винера, определяемому в частотной области уравнением (9.125).

Остановимся кратко на очень существенном для фильтров Калмана-Бьюси вопросе о возможности представления случайного процесса в виде дифференциального уравнения (9,135).

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

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

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

Учитывая (9.149), можно записать

откуда частотная передаточная функция формирующего фильтра

Подставляя в последнее выражение получаем выражение для передаточной функции формирующего фильтра

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

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

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

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

где - случайные величины с известными статистическими характеристиками.

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

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

Пример 9.8. На входе линейной следящей системы действует стационарный случайный процесс спектральная плотность которого

и случайная помеха типа «белый шум», имеющая спектральную плотность

Числовые значения коэффициентов

Определить методом Калмана-Бьюси оптимальную передаточную функцию системы, обеспечивающую минимум средней квадратической ошибки.

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

В соответствии с (9.149) представляем выражение для спектральной плотности в виде произведения комплексно-сопряженных сомножителей

и находим

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

3. Находим передаточную функцию формирующего фильтра:

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

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

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

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

Случайны процессы и шумы

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

Белый шум

Белый шум является стационарным случайным процессом , у которого автокорреляционная функция описывается дельта - функцией Дирака и, соответственно, спектральная плотность мощности не зависит от частоты и имеет постоянное значение W q (f) = σ 2 {\displaystyle W_{q}(f)=\sigma ^{2}\,\!} , равное дисперсии значений q (t) {\displaystyle q(t)\,\!} . Другими словами, все спектральные составляющие белого шума имеют одинаковую мощность (как белый цвет содержит все цвета видимого спектра). По существу, это идеализированный случайный процесс с бесконечной энергией. Но в случае постоянства спектральной плотности мощности случайного процесса в конечном диапазоне частот введение такой идеализации позволяет разрабатывать достаточно легко реализуемые оптимальные методы фильтрации. Многие помехи в радиотехнике, в технике связи и в других отраслях, в том числе в информатике, рассматривают как белый шум, если эффективная ширина спектра сигналов B s {\displaystyle B_{s}\,\!} много меньше эффективной ширины спектра шумов B q {\displaystyle B_{q}\,\!} , а спектральная плотность мощности шумов слабо изменяется в интервале спектра сигнала. Понятие "белый шум" определяет только спектральную характеристику случайного процесса и под это понятие подпадают любые случайные процессы, имеющие равномерный энергетический спектр и различные законы распределения.

Если частотный диапазон спектра, на котором рассматриваются сигналы и помехи, равен 0 − B {\displaystyle 0-B\,\!} , то спектральная плотность шума задается в виде:

W q (f) = σ 2 , 0 ≤ f ≤ B ; W q (f) = 0 , f > B (1.1) {\displaystyle W_q(f) = \sigma^2, 0 \le f \le B; W_q(f)=0, f>B \qquad \color{Maroon} (1.1) \,\!}

при этом корреляционная функция шума определяется выражением:

R q (τ) = σ 2 B ⋅ sin ⁡ (2 π B τ) 2 π B τ (1.2) {\displaystyle R_q(\tau) = \sigma^2 B \cdot \frac{\sin(2 \pi B \tau)}{2 \pi B \tau} \qquad \color{Maroon} (1.2) \,\!}

Эффективный интервал корреляции:

T k = 2 ∫ 0 ∞ | R q (τ) | d τ R q (0) (1.3) {\displaystyle T_k = 2 \frac{\int\limits_{0}^{\infty} |R_q (\tau)| d\tau}{R_q(0)} \qquad \color{Maroon} (1.3) \,\!}

Реальный интервал корреляции целесообразно определять по ширине главного максимума функции R q (τ) {\displaystyle R_{q}(\tau)\,\!} (значения при первых пересечениях нулевой линии), в котором сосредоточена основная часть энергии шумов, при этом T k = 1 B {\displaystyle T_{k}={\tfrac {1}{B}}\,\!} и B T k = 1 {\displaystyle BT_{k}=1\,\!} .

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

Рис. 1.1. Функции корреляции белого шума в частотном интервале 0-В.

Модель белого шума q (t) {\displaystyle q(t)\,\!} можно формировать как случайную по времени (аргументу) последовательность дельта - импульсов δ (t i) {\displaystyle \delta (t_{i})\,\!} со случайными амплитудными значениями :

q (t) = ∑ i a i δ (t − t i) (1.4) {\displaystyle q(t) = \sum_{i}^{} a_i \delta (t-t_i) \qquad \color{Maroon} (1.4) \,\!}

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

W q (ω) = c 2 = N σ a 2 {\displaystyle W_{q}(\omega)=c^{2}=N\sigma _{a}^{2}\,\!} где N {\displaystyle N\,\!} - число импульсов на интервале T {\displaystyle T\,\!} реализации случайного процесса, σ a 2 {\displaystyle \sigma _{a}^{2}\,\!} - дисперсия амплитуд импульсов.

Фильтрация белого шума

Спектральное описание белого шума оказывается удобным при учете влияния на него амплитудно-частотных характеристик различных устройств. Если на входе фильтра с импульсным откликом действует белый шум q (t) {\displaystyle q(t)\,\!} , то сигнал на выходе фильтра:

g (t) = h (t) × q (t) = h (t) × ∑ i a i δ (t − t i) = ∑ i a i h (t − t i) (1.5) {\displaystyle g(t) = h(t) \times q(t) = h(t) \times \sum_{i}^{} a_i \delta (t-t_i) = \sum_{i}^{} a_i h (t-t_i) \qquad \color{Maroon} (1.5) \,\!}

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

R g (τ) ≈ N σ a 2 R h (τ) = c 2 R h (τ) (1.6) {\displaystyle R_g(\tau) \approx N \sigma_a^2 R_h(\tau) = c^2 R_h(\tau) \qquad \color{Maroon} (1.6) \,\!} W g (ω) ≈ N σ a 2 | H (ω) | 2 = c 2 | H (ω) | 2 (1.7) {\displaystyle W_g(\omega) \approx N \sigma_a^2 |H(\omega)|^2 = c^2 |H(\omega)|^2 \qquad \color{Maroon} (1.7) \,\!}

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

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

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

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

Среднее квадратическое отклонение

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

y (k) = h (n) × x (k − n) (2.1) {\displaystyle y(k) = h(n) \times x(k-n) \qquad \color{Maroon} (2.1) \,\!}

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

ξ 2 ¯ = [ y (k) − s (k) ] 2 ¯ (2.2) {\displaystyle \overline{\xi^2} = \overline{^2} \qquad \color{Maroon} (2.2) \,\!}

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

ξ 2 ¯ = [ y (k) − z (k) ] 2 ¯ (2.2 ′) {\displaystyle \overline{\xi^2} = \overline{^2} \qquad \color{Maroon} (2.2") \,\!}

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

Амплитудное отношение сигнал/шум

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

ρ a = y e k s σ {\displaystyle \rho _{a}={\frac {y_{eks}}{\sigma }}\,\!} где y e k s {\displaystyle y_{eks}\,\!} - экстремальное (максимальное или минимальное) значение амплитуды сигнала, s {\displaystyle s\,\!} - средний квадратический уровень амплитудных значений помех.

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

ρ 2 ¯ = y 2 σ 2 ¯ (2.3) {\displaystyle \overline{\rho^2} = \overline{\frac{y^2}{\sigma^2}} \qquad \color{Maroon} (2.3) \,\!} где y 2 {\displaystyle y^{2}\,\!} - средний квадрат амплитуды сигнала в пределах его формы.

Энергетическое отношение сигнал/шум

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

ρ 2 = E y E q (2.4) {\displaystyle \rho^2 = \frac{E_y}{E_q} \qquad \color{Maroon} (2.4) \,\!} где E y {\displaystyle E_{y}\,\!} и E h {\displaystyle E_{h}\,\!} - энергия соответственно сигнала и шума на выходе фильтра.

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

Условие оптимальности фильтра

Фильтр Колмогорова-Винера является оптимальным фильтром формирования из входного сигнала выходного сигнала при известной форме полезного сигнала , который содержится во входном сигнале в сумме с шумами. В качестве критерия его оптимизации используется среднее квадратическое отклонение сигнала y (t) {\displaystyle y(t)\,\!} на выходе фильтра от заданной формы сигнала z (t) {\displaystyle z(t)\,\!} . Подставим уравнение свертки (2.1) {\displaystyle \color {Maroon}(2.1)\,\!} в раскрытой форме весового суммирования в выражение (2.2 ′) {\displaystyle \color {Maroon}(2.2")\,\!} и получим отклонение ξ 2 {\displaystyle \xi ^{2}\,\!} выходного сигнала y (k) = h (n) × x (k − n) {\displaystyle y(k)=h(n)\times x(k-n)\,\!} от заданной формы выходного сигнала z (k) {\displaystyle z(k)\,\!} по всем точкам массива данных:

ξ 2 = [ ∑ n h (n) x (k − n) − z (k) ] 2 ¯ (3.1) {\displaystyle \xi^2 = \overline{[ \sum_{n}^{} h(n)x(k-n)-z(k) ]^2} \qquad \color{Maroon} (3.1) \,\!}

В частном случае воспроизведения формы полезного сигнала в качестве функции z (k) {\displaystyle z(k)\,\!} принимается функция s (k) {\displaystyle s(k)\,\!} . Минимум выражения определяет значения коэффициентов оптимального фильтра. Для нахождения их значений продифференцируем выражение (3.1) {\displaystyle \color {Maroon}(3.1)\,\!} по коэффициентам фильтра и приравняем полученные уравнения нулю. В итоге получаем:

d (ξ 2) d h (n) = − z k x k − m ¯ + ∑ n h n x k − m x k − n ¯ {\displaystyle {\frac {d(\xi ^{2})}{dh(n)}}={\overline {-z_{k}x_{k-m}}}+\sum _{n}^{}h_{n}{\overline {x_{k-m}x_{k-n}}}\,\!} где x k − m x k − n ¯ = R (m − n) {\displaystyle {\overline {x_{k-m}x_{k-n}}}=R(m-n)\,\!} - корреляционная функция входного сигнала, z k x k − m ¯ = B (m) {\displaystyle {\overline {z_{k}x_{k-m}}}=B(m)\,\!} - взаимная корреляционная функция сигналов z (k) {\displaystyle z(k)\,\!} и x (k) {\displaystyle x(k)\,\!} . ∑ n h (n) R (m − n) = B (m) , n = m = 0 , 1 , 2 , . . . , M (3.2) {\displaystyle \sum_{n}^{} h(n)R(m-n) = B(m), n = m = 0,1,2, ... , M \qquad \color{Maroon} (3.2) \,\!}

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

h (n) × R (m − n) = B (m) (3.3) {\displaystyle h(n) \times R(m-n) = B(m) \qquad \color{Maroon} (3.3) \,\!}

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

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

Выражение (3.2) {\displaystyle \color {Maroon}(3.2)\,\!} может быть записано в виде системы линейных уравнений - однострочных равенств левой и правой части для фиксированных значений координаты m коэффициентов фильтра. При расчете симметричных фильтров, не сдвигающих фазы выходного сигнала, фильтр может вычисляться только одной половиной своих значений:

m = 0: h 0 R (0) + h 1 R (1) + h 2 R (2) + h 3 R (3) + . . . + h M R (M) = B (0) (3.3 ′) {\displaystyle m=0: h_0 R(0) + h_1 R(1) + h_2 R(2) + h_3 R(3) + ... + h_M R(M) = B(0) \qquad \color{Maroon} (3.3") \,\!}

M = 1: h 0 R (1) + h 1 R (0) + h 2 R (1) + h 3 R (2) + . . . + h M R (M − 1) = B (1) {\displaystyle m=1:h_{0}R(1)+h_{1}R(0)+h_{2}R(1)+h_{3}R(2)+...+h_{M}R(M-1)=B(1)\,\!}

M = 2: h 0 R (2) + h 1 R (1) + h 2 R (0) + h 3 R (1) + . . . + h M R (M − 2) = B (2) {\displaystyle m=2:h_{0}R(2)+h_{1}R(1)+h_{2}R(0)+h_{3}R(1)+...+h_{M}R(M-2)=B(2)\,\!}

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . {\displaystyle ........................................................................\,\!}

m = M: h 0 R (M) + h 1 R (M − 1) + h 2 R (M − 2) + . . . + h M R (0) = B (M) {\displaystyle m=M:h_{0}R(M)+h_{1}R(M-1)+h_{2}R(M-2)+...+h_{M}R(0)=B(M)\,\!}

Решение данной системы уравнений относительно значений дает искомый оператор фильтра.

При расчете каузальных (односторонних) фильтров выходной сигнал z (t) {\displaystyle z(t)\,\!} следует задавать со сдвигом вправо по оси координат таким образом, чтобы значимая часть функции взаимной корреляции B (m) {\displaystyle B(m)\,\!} полностью располагалась в правой части системы уравнений (3.3 ′) {\displaystyle \color {Maroon}(3.3")\,\!} .

Отметим, что R (m) = R s (m) + R q (m) {\displaystyle R(m)=R_{s}(m)+R_{q}(m)\,\!}

где - функция автокорреляции сигнала, - функция автокорреляции шума,

а B (m) = B z s (m) + B z q (m) {\displaystyle B(m)=B_{zs}(m)+B_{zq}(m)\,\!}

где B z s {\displaystyle B_{zs}\,\!} - функция взаимной корреляции сигналов z (k) {\displaystyle z(k)\,\!} и s (k) {\displaystyle s(k)\,\!} , B z q {\displaystyle B_{zq}\,\!} - функция взаимной корреляции сигнала z (k) {\displaystyle z(k)\,\!} и помех q (k) {\displaystyle q(k)\,\!} .

Подставляя данные выражения в , получаем:

h (n) × [ R s (m − n) + R q (m − n) ] = B z s (m) + B z q (m) (3.4) {\displaystyle h(n) \times = B_{zs}(m) + B_{zq}(m) \qquad \color{Maroon} (3.4) \,\!}

Частотная характеристика фильтра

Частотная характеристика фильтра находится преобразованием Фурье левой и правой части уравнения (3.4) {\displaystyle \color {Maroon}(3.4)\,\!} :

H (w) [ W s (w) + W q (w) ] = W z s (w) + W z q (w) {\displaystyle H(w)=W_{zs}(w)+W_{zq}(w)\,\!} H (w) = W z s (w) + W z q (w) W s (w) + W q (w) (3.5) {\displaystyle H(w) = \frac{W_{zs}(w)+W_{zq}(w)}{W_s(w) + W_q(w)} \qquad \color{Maroon} (3.5) \,\!} где W s (ω) ⟺ R s (m) {\displaystyle W_{s}(\omega)\iff R_{s}(m)\,\!} и W q (ω) ⟺ R q (m) {\displaystyle W_{q}(\omega)\iff R_{q}(m)\,\!} - энергетические спектры (плотности мощности) сигнала и помех, W z s (ω) ⟺ B z s (m) {\displaystyle W_{zs}(\omega)\iff B_{zs}(m)\,\!} - взаимный энергетический спектр входного и выходного сигналов, W z q (ω) ⟺ B z q (m) {\displaystyle W_{zq}(\omega)\iff B_{zq}(m)\,\!} - взаимный энергетический спектр выходного сигнала и помех.

В геофизической практике обычно имеет место статистическая независимость полезного сигнала, а, следовательно, и сигнала z (k) {\displaystyle z(k)\,\!} , от шумов, при этом B z q = 0 {\displaystyle B_{zq}=0\,\!} и фильтр называют оптимальным по сглаживанию шумов при заданной форме выходного сигнала:

H (ω) = W z s (ω) W s (ω) + W q (ω) (3.6) {\displaystyle H(\omega) = \frac{W_{zs}(\omega)}{W_s(\omega)+W_q(\omega)} \qquad \color{Maroon} (3.6) \,\!}

Фильтр (3.6) {\displaystyle \color {Maroon}(3.6)\,\!} оптимален в том смысле, что максимизирует отношение мощности сигнала к мощности шума по всему интервалу сигнала, но не в каждой индивидуальной точке.

Выражения (3.5 − 3.6) {\displaystyle \color {Maroon}(3.5-3.6)\,\!} , как правило, всегда имеют решение. Однако это не означает возможность формирования фильтром любой заданной формы выходного сигнала. Из чисто практических соображений можно сразу предполагать, что если спектр заданного сигнала z (t) {\displaystyle z(t)\,\!} больше значимой части спектра полезного сигнала s (t) {\displaystyle s(t)\,\!} , то оператор фильтра попытается сформировать требуемые высокие частоты заданного сигнала из незначимых частот спектра полезного сигнала, что может потребовать огромных коэффициентов усиления на этих частотах, под действие которых попадут и частотные составляющие шумов. Результат такой операции непредсказуем. Эти практические соображения можно распространить и на все частотные соотношения входного и выходного сигналов линейных фильтров: значимые гармоники спектров выходных сигналов должны формироваться из значимых гармоник спектров входных сигналов.

Если заданная форма сигнала z (k) {\displaystyle z(k)\,\!} совпадает с формой полезного сигнала s (k) {\displaystyle s(k)\,\!} , то B (m) = B s s = R s {\displaystyle B(m)=B_{ss}=R_{s}\,\!} и фильтр называют фильтром воспроизведения полезного сигнала:

H (ω) = W s (ω) W s (ω) + W q (ω) (3.7) {\displaystyle H(\omega) = \frac{W_s(\omega)}{W_s(\omega)+W_q(\omega)} \qquad \color{Maroon} (3.7) \,\!}

Выражения (3.6 − 3.7) {\displaystyle \color {Maroon}(3.6-3.7)\,\!} достаточно наглядно демонстрируют физический смысл формирования передаточной функции фильтра. При воспроизведении сигнала частотная функция взаимной корреляции входного сигнала с выходным (плотность взаимной мощности) повторяет частотную функцию автокорреляции (плотность мощности сигнала). Плотность мощности статистических шумов распределена по частотному диапазону равномерно, в отличие от плотности мощности сигнала W s {\displaystyle W_{s}\,\!} , которая, в зависимости от формы сигнала, может занимать любые частотные интервалы спектрального диапазона. Частотная передаточная функция фильтра воспроизведения сигнала формируется отношением W s (ω) W s (ω) + W q (ω) {\displaystyle {\frac {W_{s}(\omega)}{W_{s}(\omega)+W_{q}(\omega)}}\,\!} . На частотах, где сосредоточена основная энергия сигнала, имеет место W s (ω) ≫ W q (ω) {\displaystyle W_{s}(\omega)\gg W_{q}(\omega)\,\!} и (как минимум, больше 0.5). Там, где значение становится меньше W q {\displaystyle W_{q}\,\!} , коэффициент передачи фильтра становится меньше 0.5, и в пределе на всех частотах, где полностью отсутствуют частотные составляющие сигнала.

Аналогичный процесс имеет место и при формировании произвольного сигнала z (t) {\displaystyle z(t)\,\!} на выходе фильтра, только в этом случае из частот входного сигнала устанавливаются на выделение и усиление частоты взаимной мощности входного и выходного сигнала, необходимые для формирования сигнала z (t) {\displaystyle z(t)\,\!} , причем коэффициент на этих частотах может быть много больше 1, а подавляться могут не только шумы, но и частоты основного сигнала, если их нет в сигнале z (t) {\displaystyle z(t)\,\!} .

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

Задание мощности шумов

Следует внимательно относиться к заданию функции шумов . При полной неопределенности шума необходимо, как минимум, выполнять оценку его дисперсии и распространять на весь частотный диапазон с соответствующей нормировкой на его величину ( 2 ∫ 0 Ω W q (ω) d ω = σ 2 {\displaystyle 2\int \limits _{0}^{\Omega }W_{q}(\omega)d\omega =\sigma ^{2}\,\!} ), т.е. считать его белым шумом. При известной функции полезного сигнала в зарегистрированной реализации значение дисперсии шумов в первом приближении может быть оценено по разности дисперсий реализации и функции полезного сигнала. Можно выполнить и выделение шумов из входного сигнала в отдельный шумовой массив, например, вейвлетным преобразованием. Однако использовать выделенный шум непосредственно для вычисления функции W q (ω) {\displaystyle W_{q}(\omega)\,\!} допустимо только по достаточно представительной реализации при условии стационарности и эргодичности шума. В противном случае функция Wq(w) будет отображать только распределение шумов в зарегистрированной реализации сигнала, а соответственно фильтр будет оптимален только к этой реализации, что не гарантирует его оптимальности к любой другой реализации. Но для обработки единичной зарегистрированной реализации сигнала такой метод не только вполне допустим, но и может существенно повысить точность формирования выходного сигнала.

Эффективность фильтра

Из выражений (3.5 − 3.7) {\displaystyle \color {Maroon}(3.5-3.7)\,\!} следует, что с позиции минимального искажения полезного сигнала при максимальном подавлении шумов фильтр Колмогорова-Винера эффективен в тем большей степени, чем больше отношение сигнал/шум на входе фильтра. В пределе, при W q (ω) ≪ W s (ω) {\displaystyle W_{q}(\omega)\ll W_{s}(\omega)\,\!} имеем H (ω) ⇒ 1 {\displaystyle H(\omega)\Rightarrow 1\,\!} и фильтр воспроизводит входной (или заданный) сигнал без искажений (помех нет, исправлять нечего). Отметим также, что помеха, коррелированная с полезным сигналом, как это следует из (3.5) {\displaystyle \color {Maroon}(3.5)\,\!} , используется фильтром для повышения точности воспроизведения сигнала. С другой стороны, при W q (ω) ≫ W s (ω) {\displaystyle W_{q}(\omega)\gg W_{s}(\omega)\,\!} имеем H (ω) ⇒ 0 {\displaystyle H(\omega)\Rightarrow 0\,\!} и сигнал будет сильно искажен, но никакой другой фильтр лучшего результата обеспечить не сможет.

Пример: Расчет оптимального фильтра воспроизведения сигнала. Выполняется в среде Mathcad. Форма входного сигнала считается известной и близка к гауссовой. На входной сигнал наложен статистический шум с равномерным распределением мощности по всему частотному диапазону (белый шум), некоррелированный с сигналом, и функцию W z q {\displaystyle W_{zq}\,\!} принимаем равной нулю. Для наглядного просмотра связи параметров фильтра с параметрами сигнала задаем модели сигналов в двух вариантах: K:= 1000 k:= 0.. K A:= 50 {\displaystyle K:=1000k:=0..KA:=50\,\!} S 1 k:= A ⋅ e x p [ − 0.0005 ⋅ (k − 500) 2 ] s 2 k:= A ⋅ e x p [ − 0.00003 ⋅ (k − 500) 2 ] {\displaystyle s1_{k}:=A\cdot exp[-0.0005\cdot (k-500)^{2}]s2^{k}:=A\cdot exp[-0.00003\cdot (k-500)^{2}]\,\!} полезные сигналы Q:= 30 q k:= r n d (Q) = Q / 2 x 1 k:= s 1 k + q k x 2 k:= s 2 k + q k {\displaystyle ~Q:=30q_{k}:=rnd(Q)=Q/2x1_{k}:=s1_{k}+q_{k}x2_{k}:=s2_{k}+q_{k}} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} входные сигналы Рис. 3.1. Модельные сигналы. В качестве выходных сигналов задаем те же самые функции s 1 {\displaystyle s1\,\!} и s 2 {\displaystyle s2\,\!} . Быстрым преобразованием Фурье вычис-ляем спектры сигналов и формируем спектры плотности мощности: S 1:= C F F T (s 1) S 2:= C F F T (s 2) Q:= C F F T (q) {\displaystyle S1:=CFFT(s1)S2:=CFFT(s2)Q:=CFFT(q)\,\!} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} спектры мощности W s 1:= S 1 ⋅ S 1 ¯ S 2:= S 2 ⋅ S 2 ¯ W q:= Q ⋅ Q ¯ {\displaystyle Ws1:={\overline {S1\cdot S1}}S2:={\overline {S2\cdot S2}}W_{q}:={\overline {Q\cdot Q}}\,\!} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} спектры мощности D s 1:= v a r (s 1) D s 2:= v a r (s 2) D x 1:= v a r (x 1) D x 2:= v a r (x 2) D q:= v a r (q) {\displaystyle Ds_{1}:=var(s1)Ds_{2}:=var(s2)Dx_{1}:=var(x1)Dx_{2}:=var(x2)Dq:=var(q)\,\!} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} дисперсии D s 1 = 124.308 D s 2 = 310.264 D x 1 = 202.865 D x 2 = 386.78 D q = 79.038 {\displaystyle Ds1=124.308Ds2=310.264Dx1=202.865Dx2=386.78Dq=79.038\,\!} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} информация M e a n (W q) = 0.079 {\displaystyle mean(Wq)=0.079\,\!} W q 1:= (D x 1 − D s 1) / (K + 1) W q 1 = 0.078 {\displaystyle ~Wq1:=(Dx1-Ds1)/(K+1)Wq1=0.078} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} информация W q 2:= (D x 2 − D s 2) / (K + 1) W q 2 = 0.076 {\displaystyle ~Wq2:=(Dx2-Ds2)/(K+1)Wq2=0.076} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} информация W q k:= W q 1 {\displaystyle Wqk:=Wq1\,\!} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} замена на постоянное распределение Для воспроизведения сигналов вычисления функций W z s {\displaystyle W_{zs}\,\!} не требуется, т.к. W z s = W s {\displaystyle W_{zs}=W_{s}\,\!} . Вычисление W q {\displaystyle W_{q}\,\!} также имеет только информативный характер. Передаточные функции оптимальных фильтров (приведены на рис. 3.2): Рис. 3.2. Передаточные функции оптимальных фильтров в сопоставлении с нормированными модулями спектров сигналов. Как следует из рисунка 3.2, для плавных монотонных функций, основная энергия которых сосредоточена в низкочастотной части спектра, передаточные функции оптимальных фильтров, по существу, представляют собой низкочастотные сглаживающие фильтры с автоматической подстройкой граничной частоты пропускания под основные частоты входного сигнала. Операторы фильтров можно получить обратным преобразованием Фурье: H 1:= I C F F T (H 1) / (K + 1) h 2:= I C F F T (H 2) / (K + 1) {\displaystyle h1:=ICFFT(H1)/(K+1)h2:=ICFFT(H2)/(K+1)\,\!} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} обратное преобразование Фурье Рис. 3.3. Импульсные отклики фильтров. Оператор фильтра, в принципе, бесконечен. В данном случае, при использовании БПФ максимальное число отсчетов равно К/2 = 500. Усечение размеров оператора может выполняться по типовым методам с применением весовых функций (в расчете операторы нормируются к 1, весовые функции не применяются). N 1:= 160 n 1:= 0.. N 1 N 2 ; = 500 n 2:= 0.. N 2 {\displaystyle N1:=160n1:=0..N1N2;=500n2:=0..N2\,\!} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} размеры и нумерация операторов H m 1:= h 1 0 + 2 ∑ n 1 h 1 n 1 h m 1 = 0.988 h 1:= h 1 / h m 1 {\displaystyle hm1:=h1_{0}+2\sum _{n1}h1_{n1}hm1=0.988h1:=h1/hm1\,\!} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} нормировка H m 2:= h 2 0 + 2 ∑ n 1 h 2 n 2 h m 2 = 1.001 h 2:= h 2 / h m 2 {\displaystyle hm2:=h2_{0}+2\sum _{n1}h2_{n2}hm2=1.001h2:=h2/hm2\,\!} ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} нормировка Y k:= h 0 ⋅ x k + ∑ n = 1 N h n ⋅ i f (k − n < 0 , 0 , x k − n) + ∑ n = 1 N h n ⋅ i f (k + n < K , 0 , x k + n) {\displaystyle y_{k}:=h_{0}\cdot x_{k}+\sum _{n=1}^{N}h_{n}\cdot if(k-n<0,0,x_{k-n})+\sum _{n=1}^{N}h_{n}\cdot if(k+n ⇐ {\displaystyle \color {Maroon}\Leftarrow \,\!} свертка Рис. 3.4. Проверка действия оптимальных фильтров. Коэффициент усиления дисперсии помех K d:= (h 0) 2 + 2 ⋅ ∑ n = 1 N h n , K d 1 = 0.021 K d 2 = 0.0066 {\displaystyle Kd:=(h_{0})2+2\cdot \sum _{n=1}^{N}h_{n},Kd1=0.021Kd2=0.0066\,\!} Среднеквадратическое отклонение воспроизведения сигнала ⇒ {\displaystyle \color {Maroon}\Rightarrow \,\!} e:= 1 K + 1 ∑ k (s k − y k) 2 {\displaystyle e:={\sqrt {{\frac {1}{K+1}}\sum _{k}^{}(s_{k}-y_{k})^{2}}}\,\!} Проверка действия оператора фильтра приведена на рис. 3.4.

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

Рис. 3.5. Оптимальная фильтрация сигнала со сложным спектральным составом.

Рис. 3.6. Оптимальная фильтрация радиоимпульса.

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

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

Если в правой части уравнения (3.3) {\displaystyle \color {Maroon}(3.3)\,\!} желаемым сигналом задать входной сигнал со сдвигом на величину k Δ t {\displaystyle k\Delta t\,\!} , то при этом B (m) = R (m + k) {\displaystyle B(m)=R(m+k)\,\!} и уравнение принимает вид:

h (n) × R (m − n) = R (m + k) (3.8) {\displaystyle h(n) \times R(m-n) = R(m+k) \qquad \color{Maroon} (3.8) \,\!}

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

Оптимальные фильтры сжатия сигналов

Условие оптимальности. Фильтр сжатия сигнала x (t) {\displaystyle x(t)\,\!} , по существу, представляет собой фильтр формирования сигнала z (t) {\displaystyle z(t)\,\!} с эффективной шириной длительности, меньшей по сравнению с эффективной шириной длительности полезного сигнала s (t) {\displaystyle s(t)\,\!} во входном сигнале x (t) {\displaystyle x(t)\,\!} . Расчет оптимального фильтра сжатия может выполняться непосредственно по выражениям (3.3) {\displaystyle \color {Maroon}(3.3)\,\!} .

В предельном случае сжатия сигнала до импульса Кронекера положим, что z (k) = δ (k) {\displaystyle z(k)=\delta (k)\,\!} при статистической независимости сигнала и шума. Отсюда:

B s z (m) = δ (m) × s (k + m) = s (− m) {\displaystyle Bsz(m)=\delta (m)\times s(k+m)=s(-m)\,\!} h (n) × (R s (m − n) + R q (m − n)) = s (− m) (4.1) {\displaystyle h(n) \times (Rs(m-n)+Rq(m-n)) = s(-m) \qquad \color{Maroon} (4.1) \,\!} H (w) = S ∗ (w) | S (w) | 2 + W q (w) (4.2) {\displaystyle H(w) = \frac{S*(w)}{|S(w)|^2+Wq(w)} \qquad \color{Maroon} (4.2) \,\!}

При некоррелированной помехе с дисперсией система уравнений для определения значений коэффициентов h (n) {\displaystyle h(n)\,\!} :

h 0 (R (0) + σ 2) + h 1 R (1) + h 2 R (2) + h 3 R (3) + . . . + h M R (M) = s (0) (4.3) {\displaystyle h_0 (R(0) + \sigma^2) + h^1 R(1) + h_2 R(2) + h_3 R(3) + ... + h_M R(M) = s(0) \qquad \color{Maroon} (4.3) \,\!} h 0 R (1) + h 1 R (0) + h 2 R (1) + h 3 R (2) + . . . + h M R (M − 1) = 0 , {\displaystyle h_{0}R(1)+h_{1}R(0)+h_{2}R(1)+h_{3}R(2)+...+h_{M}R(M-1)=0,\,\!} h 0 R (2) + h 1 R (1) + h 2 R (0) + h 3 R (1) + . . . + h M R (M − 2) = 0 , {\displaystyle h_{0}R(2)+h_{1}R(1)+h_{2}R(0)+h_{3}R(1)+...+h_{M}R(M-2)=0,\,\!} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . {\displaystyle .................................................................\,\!} h 0 R (M) + h 1 R (M − 1) + h 2 R (M − 2) + . . . . . . . + h M R (0) = 0 {\displaystyle h_{0}R(M)+h_{1}R(M-1)+h_{2}R(M-2)+.......+h_{M}R(0)=0\,\!}

При расчете коэффициентов фильтра значение s (0) {\displaystyle s(0)\,\!} обычно принимается произвольным, чаще всего равным площади сигнала s (t) {\displaystyle s(t)\,\!} . Тем самым делается попытка полного сжатия площади сигнала до весовой функции Кронекера, что возможно только для сигналов со спектром в главном диапазоне до частоты Найквиста.

Рис. 4.1. Сжатие гладких сигналов с разным уровнем шумов.

Для гладких и монотонных функций со спектром в низкочастотной части главного диапазона сжатие до импульса Кронекера невозможно, и в зависимости от уровня шумов фильтр поднимает насколько возможно высокие частоты сигнала, учитывая значение уровня шумов. При этом нарушаются условия нормировки площади оператора фильтра к 1, о чем можно судить по значению передаточной функции H (ω) {\displaystyle H(\omega)\,\!} при ω {\displaystyle \omega \,\!} , которое становится меньше 1, и при обратном преобразовании H (ω) ⇒ h (m) {\displaystyle H(\omega)\Rightarrow h(m)\,\!} оператор h (m) {\displaystyle h(m)\,\!} требует нормировки к 1. Все эти факторы можно наглядно видеть на рис. 4.1.

На рисунке приведены три сигнала с одной и той же базовой функцией, на которую наложены шумы разного уровня. При малом уровне шумов (сигнал x 1 {\displaystyle x1\,\!} ) фильтр в максимальной степени использует высокие частоты сигнала ( | H 1 | ≫ 1 {\displaystyle |H1|\gg 1\,\!} на этих частотах), сохраняя устойчивость работы фильтра при достаточно удовлетворительном (хотя и больше 1) коэффициенте усиления дисперсии помех при максимально возможном сжатии сигнала. При повышении уровня шумов (сигналы x 2 {\displaystyle x2\,\!} и x 3 {\displaystyle x3\,\!} ) подъем высоких частот сигнала уменьшается, и сжатие сигнала соответственно также уменьшается, предпочтение отдается максимальному подавлению шумов.

Рис. 4.2. Сжатие сигнала с высокочастотным спектром.

На рис. 4.2. приведен пример сжатия сигнала, близкого к прямоугольному импульсу. Базовая функция сигнала s (k) {\displaystyle s(k)\,\!} имеет достаточно высокочастотный спектр мощности W s (ω) {\displaystyle W_{s}(\omega)\,\!} , и при задании формы выходного сигнала сжатия в виде гауссовой функции z (k) {\displaystyle z(k)\,\!} передаточная функция фильтра H (ω) {\displaystyle H(\omega)\,\!} обеспечивает достаточно уверенное сжатие сигнала (при уменьшении уровня шумов практически до заданной формы).

В пределе, при W q = 0 {\displaystyle W_{q}=0\,\!} фильтр сжатия превращается в обратный фильтр (фильтр деконволюции):

H (ω) = S ∗ (w) | S (w) | 2 (4.4) {\displaystyle H(\omega) = \frac{S*(w)}{|S(w)|^2} \qquad \color{Maroon} (4.4) \,\!}

На выходе такого фильтра имеем:

Y (ω) = H (ω) X (ω) → 1 , {\displaystyle Y(\omega)=H(\omega)X(\omega)\rightarrow 1,\,\!} при X (w) → S (w) (4.4) {\displaystyle X(w) \rightarrow S(w) \qquad \color{Maroon} (4.4) \,\!}

Реализация фильтра возможна только при условии S (ω) {\displaystyle S(\omega)\,\!} на всех частотах в главном частотном диапазоне. В противном случае, при S (ω i) → 0 {\displaystyle S(\omega _{i})\rightarrow 0\,\!} , H (ω i) → ∞ {\displaystyle H(\omega _{i})\rightarrow \infty \,\!} и фильтр становится неустойчивым. Для исключения возможности такого явления в фильтр (4.4) {\displaystyle \color {Maroon}(4.4)\,\!} вводится стабилизатор a:

H (ω) = S ∗ (w) | S (w) | 2 + a (4.5) {\displaystyle H(\omega) = \frac{S*(w)}{|S(w)|^2 + a} \qquad \color{Maroon} (4.5) \,\!} где | S (w) | 2 + a > 0 {\displaystyle |S(w)|^{2}+a>0\,\!} во всем частотном диапазоне.

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

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

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

Частотная характеристика

расчета фильтра требуется задать известную форму полезного сигнала s (k) ⟺ S (ω) {\displaystyle s(k)\iff S(\omega)\,\!} и функцию автокорреляции или спектр мощности помех R q (m) ⟺ W q (ω) {\displaystyle R_{q}(m)\iff W_{q}(\omega)\,\!} . Полный входной сигнал принимается по аддитивной модели: x (t) = s (t) + q (t) {\displaystyle x(t)=s(t)+q(t)\,\!} . На выходе проектируемого фильтра h (n) ⟺ H (ω) {\displaystyle h(n)\iff H(\omega)\,\!} для составляющих выходного сигнала имеем:

y (t) = 1 2 π ∫ − ∞ ∞ H (ω) S (ω) e j ω t d ω (5.1) {\displaystyle y(t) = \frac{1}{2 \pi} \int\limits_{- \infty}^{\infty} H(\omega) S(\omega) e^{j \omega t } d \omega \qquad \color{Maroon} (5.1) \,\!} σ 2 = 1 2 π ∫ − ∞ ∞ | H (ω) | 2 W q (ω) d ω (5.2) {\displaystyle \sigma^2 = \frac{1}{2 \pi} \int\limits_{- \infty}^{\infty} |H(\omega)|^2 W_q(\omega) d \omega \qquad \color{Maroon} (5.2) \,\!} где σ {\displaystyle \sigma \,\!} - средняя квадратическая амплитуда выходной помехи.

Оптимальным в задаче обнаружения одиночного сигнала конечной длительности является фильтр, обеспечивающий на выходе максимальное отношение пиковой мощности сигнала к мощности шума в момент окончания импульса. Значения (5.1 − 5.2) {\displaystyle \color {Maroon}(5.1-5.2)\,\!} используются для задания критерия максимума отношения сигнал/шум (2.3) {\displaystyle \color {Maroon}(2.3)\,\!} для произвольной точки :

ρ = [ y (t i) ] 2 δ 2 (5.3) {\displaystyle \rho = \frac{^2}{\delta^2} \qquad \color{Maroon} (5.3) \,\!}

Исследование функции (5.3) {\displaystyle \color {Maroon}(5.3)\,\!} на максимум показывает, что он достигается при частотной характеристике фильтра:

H (ω) = e − j ω t i | S ∗ (ω) | W q (ω) (5.4) {\displaystyle H(\omega) = \frac{e^{-j \omega t_i} |S*(\omega)|}{W_q(\omega)} \qquad \color{Maroon} (5.4) \,\!}

Для физически реализуемых фильтров в качестве точки t i {\displaystyle t_{i}\,\!} целесообразно использовать интервал длительности импульса τ {\displaystyle \tau \,\!} , при этом:

H (ω) = e − j ω τ | S ∗ (ω) | W q (ω) = e − j φ s − j ω τ | S (ω) | W q (ω) (5.4 ′) {\displaystyle H(\omega) = \frac{e^{-j \omega \tau} |S*(\omega)|}{W_q(\omega)} = \frac{e^{-j \varphi_s -j \omega \tau} |S(\omega)|}{W_q(\omega)} \qquad \color{Maroon} (5.4") \,\!}

Аргумент φ s {\displaystyle \varphi _{s}\,\!} в этом выражении компенсирует фазовые сдвиги составляющих спектра сигнала, а ω τ {\displaystyle \omega \tau \,\!} обеспечивает их задержку на время длительности сигнала. Таким образом, на концевой части сигнала фильтр выполняет синфазное суммирование всех полезных частотных составляющих входного сигнала с весами, пропорциональными отношению | S (ω) | W q (ω) {\displaystyle {\tfrac {|S(\omega)|}{W_{q}(\omega)}}\,\!} , что обеспечивает накопление амплитуды полезного сигнала на интервале всей длительности входного импульса и формирует максимум сигнала на момент его окончания. Вместе с тем фильтр ослабляет спектральные составляющие шума тем сильнее, чем меньше модуль | S (ω) | {\displaystyle |S(\omega)|\,\!} , и полная мощность шума на выходе фильтра оказывается меньшей, чем на входе.

Для получения линейных уравнений расчета коэффициентов фильтра без потери общности можно принять , при этом:

H (ω) = S ∗ (ω) W q (ω) = | S (ω) | e − j φ s (ω) W q (ω) (5.5) {\displaystyle H(\omega) = \frac{S*(\omega)}{W_q(\omega)} = \frac{ |S(\omega)| e^{-j \varphi_s (\omega)}}{W_q(\omega)} \qquad \color{Maroon} (5.5) \,\!}

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

H (ω) W q (ω) = S ∗ (ω) ⟺ h (n) × R q (n − m) = s (− m) (5.6) {\displaystyle H(\omega) W_q(\omega) = S*(\omega) \iff h(n) \times R_q(n-m) = s(-m) \qquad \color{Maroon} (5.6) \,\!}

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

h 0 R q (0) + h 1 R q (1) + h 2 R q (2) + h 3 R q (3) + . . . + h M R q (M) = S (− M) , {\displaystyle h_{0}R_{q}(0)+h_{1}R_{q}(1)+h_{2}R_{q}(2)+h_{3}R_{q}(3)+...+h_{M}R_{q}(M)=S(-M),\,\!} h 0 R q (1) + h 1 R q (0) + h 2 R q (1) + h 3 R q (2) + . . . + h M R q (M − 1) = S (− M + 1) , {\displaystyle h_{0}R_{q}(1)+h_{1}R_{q}(0)+h_{2}R_{q}(1)+h_{3}R_{q}(2)+...+h_{M}R_{q}(M-1)=S(-M+1),\,\!} h 0 R q (2) + h 1 R q (1) + h 2 R q (0) + h 3 R q (1) + . . . + h M R q (M − 2) = S (− M + 2) , {\displaystyle h_{0}R_{q}(2)+h_{1}R_{q}(1)+h_{2}R_{q}(0)+h_{3}R_{q}(1)+...+h_{M}R_{q}(M-2)=S(-M+2),\,\!} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . {\displaystyle .............................................................................\,\!} h 0 R q (M) + h 1 R q (M − 1) + h 2 R q (M − 2) + h 3 R q (M − 3) + . . . + h M R q (0) = S (0) , {\displaystyle h_{0}R_{q}(M)+h_{1}R_{q}(M-1)+h_{2}R_{q}(M-2)+h_{3}R_{q}(M-3)+...+h_{M}R_{q}(0)=S(0),\,\!}

При задании t i {\displaystyle t_{i}\,\!} по центру симметричных входных сигналов можно получить симметричные двусторонние фильтры, не изменяющие фазы сигнала, что удобно для цифровой обработки данных.

На рис. 5.1 приведен пример фильтрации фильтром обнаружения сигнала радиоимпульса (информационный сигнал) в сумме с шумами (входной сигнал) при отношении сигнал/шум по средним амплитудным значениям на входе фильтра порядка 1. Аналогичное отношение сигнал/шум на выходе фильтра повышается до 7 по интервалу полезного сигнала в целом, и превышает 8 в центральной части интервала сигнала.

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

Эффективность фильтра

Из выражения можно видеть, что фильтр имеет максимальный коэффициент передачи на частотах доминирования сигнала и минимальный коэффициент передачи на частотах доминирования помех. Синфазность суммирования всех частотных составляющих выходного сигнала обеспечивает максимальную амплитуду выходного сигнала в заданный момент времени t i {\displaystyle t_{i}\,\!} . Значение максимальной амплитуды можно оценить, приняв t i = 0 {\displaystyle t_{i}=0\,\!} , при этом выходной сигнал:

y (0) ⟺ S (ω) ⋅ H (ω) = 1 2 π ∫ − ∞ ∞ S (ω) S ∗ (ω) W q (ω) {\displaystyle y(0)\iff S(\omega)\cdot H(\omega)={\frac {1}{2\pi }}\int \limits _{-\infty }^{\infty }{\frac {S(\omega)S*(\omega)}{W_{q}(\omega)}}\,\!}

Коэффициент передачи фильтра прямо определяется спектром подлежащего обнаружению сигнала, его формой и длительностью. Для оценки эффективности фильтра зададим входной сигнал в виде прямоугольного импульса амплитудой u 0 {\displaystyle u_{0}\,\!} длительностью τ {\displaystyle \tau \,\!} на интервале 0 − τ {\displaystyle 0-\tau \,\!} . Спектральная плотность прямоугольного импульса при интегральном преобразовании Фурье:

Π (ω) = 1 − e − j ω τ j ω , Π ∗ (ω) = e j ω τ − 1 j ω {\displaystyle \Pi (\omega)={\frac {1-e^{-j\omega \tau }}{j\omega }},\Pi *(\omega)={\frac {e^{j\omega \tau }-1}{j\omega }}\,\!}

подстановке в (5.4 ′) {\displaystyle \color {Maroon}(5.4")\,\!} , принимая W q (ω) = c o n s t {\displaystyle W_{q}(\omega)=const\,\!} , коэффициент передачи фильтра:

H (ω) = α (1 − e j ω τ) e − j ω τ j ω = α 1 − e j ω τ j ω {\displaystyle H(\omega)=\alpha {\frac {(1-e^{j\omega \tau })e^{-j\omega \tau }}{j\omega }}=\alpha {\frac {1-e^{j\omega \tau }}{j\omega }}\,\!}

где - коэффициент пропорциональности с размерностью, обратной спектральной плотности, для получения безразмерных значений коэффициента H (ω) {\displaystyle H(\omega)\,\!} . При a = 1 {\displaystyle a=1\,\!} (нормировка оператора фильтра производится, как правило, по коэффициенту усиления постоянной составляющей входного сигнала) сигнал на выходе фильтра:

U o u t (t) = u 0 2 π ∫ − ∞ ∞ Π (ω) H (ω) d ω = ∫ − ∞ ∞ (1 − e − j ω τ) 2 ⋅ e j ω τ d ω {\displaystyle U_{out}(t)={\frac {u_{0}}{2\pi }}\int \limits _{-\infty }^{\infty }\Pi (\omega)H(\omega)d\omega =\int \limits _{-\infty }^{\infty }(1-e^{-j\omega \tau })^{2}\cdot e^{j\omega \tau }d\omega \,\!} U o u t (t) = U 0 { t | t > 0 − 2 (t − τ) | t > 0 + (t − 2 τ) | t > 0 } {\displaystyle U_{out}(t)=U_{0}{\big \{}t|_{t>0}-2(t-\tau)|_{t>0}+(t-2\tau)|_{t>0}{\big \}}\,\!}

Как можно видеть на рис 5.2, выходной сигнал для входного прямоугольного импульса представляет собой треугольный импульс длительностью 2t по основанию с максимальным значением амплитуды на концевой части входного импульса. Это определяется тем, что при W q (ω) = 1 {\displaystyle W_{q}(\omega)=1\,\!} оператор фильтра полностью повторяет форму входного сигнала (прямоугольного импульса), а выходной сигнал в отсутствие шумов представляет собой свертку двух одинаковых импульсов, максимальное значение которой достигается при полном входе сигнала в оператор фильтра ( t = τ {\displaystyle t=\tau \,\!} ) и равно полной энергии входного импульса:

U 0 (t) = ∫ 0 τ Π (t) h (t) d t = ∫ 0 τ Π (t) 2 d t = u 0 2 τ {\displaystyle U_{0}(t)=\int \limits _{0}^{\tau }\Pi (t)h(t)dt=\int \limits _{0}^{\tau }\Pi (t)^{2}dt=u_{0}^{2}\tau \,\!}

Значение U 0 {\displaystyle U_{0}\,\!} определяется нормировкой оператора фильтра α {\displaystyle \alpha \,\!} . Что касается усиления дисперсии (мощности) шумов, то, как известно, дисперсия шума на выходе фильтра равна входной дисперсии входных шумов σ 2 {\displaystyle \sigma ^{2}\,\!} , умноженной на интеграл квадрата импульсного отклика фильтра (для цифровых систем – сумма квадратов коэффициентов оператора фильтра):

σ o u t 2 = σ 2 ∫ 0 τ h 2 (t) d t = σ 2 2 π ∫ − ∞ ∞ | H (ω) | 2 d ω {\displaystyle \sigma _{out}^{2}=\sigma ^{2}\int \limits _{0}^{\tau }h^{2}(t)dt={\frac {\sigma ^{2}}{2\pi }}\int \limits _{-\infty }^{\infty }|H(\omega)|^{2}d\omega \,\!}

Для вычисления интеграла модуль передаточной функции фильтра для прямоугольного импульса может быть представлен в виде интегрального синуса:

∫ − ∞ ∞ | H (ω) | 2 d ω = 2 τ u 0 2 ∫ − ∞ ∞ s i n c (ω τ 2) d (ω τ 2) = 2 π u 0 2 τ {\displaystyle \int \limits _{-\infty }^{\infty }|H(\omega)|^{2}d\omega =2\tau u_{0}^{2}\int \limits _{-\infty }^{\infty }sinc\left({\frac {\omega \tau }{2}}\right)d\left({\frac {\omega \tau }{2}}\right)=2\pi u_{0}^{2}\tau \,\!}

Дисперсия шумов на выходе:

σ o u t 2 = σ 2 u 0 2 τ {\displaystyle \sigma _{out}^{2}=\sigma ^{2}u_{0}^{2}\tau \,\!}

С использованием этого выражения для отношения мощности сигнала к мощности шума для сигналов на входе и выходе фильтра имеем:

ρ i n = u 0 2 σ 2 , ρ o u t = u 0 4 τ 2 σ 2 u 0 2 τ = u 0 2 τ σ 2 {\displaystyle \rho _{in}={\frac {u_{0}^{2}}{\sigma ^{2}}},\rho _{out}={\frac {u_{0}^{4}\tau ^{2}}{\sigma ^{2}u_{0}^{2}\tau }}={\frac {u_{0}^{2}\tau }{\sigma ^{2}}}\,\!}

Соответственно, для отношения амплитудных значений сигнала к среднеквадратическим значениям шума:

ρ i n = u 0 σ , ρ o u t = u 0 τ τ {\displaystyle \rho _{in}={\frac {u_{0}}{\sigma }},\rho _{out}={\frac {u_{0}}{\tau }}{\sqrt {\tau }}\,\!}

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

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

Согласованный фильтр

При помехах типа белого шума W q (ω) = σ 2 {\displaystyle W_{q}(\omega)=\sigma ^{2}\,\!} и H (ω) = S ∗ (ω) σ 2 {\displaystyle H(\omega)={\frac {S*(\omega)}{\sigma ^{2}}}\,\!} . Постоянный множитель 1 σ 2 {\displaystyle {\frac {1}{\sigma ^{2}}}\,\!} может быть опущен. Частотная характеристика фильтра определяется только спектром сигнала, при этом:

h (n) = s (− n) (5.7) {\displaystyle h(n) = s(-n) \qquad \color{Maroon} (5.7) \,\!}

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

Обратный фильтр

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

W q = σ 2 | S (ω) | 2 {\displaystyle W_{q}=\sigma ^{2}|S(\omega)|^{2}\,\!}

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

H (ω) = S ∗ (ω) σ 2 | S (ω) | 2 = 1 σ 2 S (ω) (5.8) {\displaystyle H(\omega) = \frac{S*(\omega)}{\sigma^2 |S(\omega)|^2} = \frac{1}{\sigma^2 S(\omega)} \qquad \color{Maroon} (5.8) \,\!}

Выражение (5.8) {\displaystyle \color {Maroon}(5.8)\,\!} с точностью до постоянного множителя соответствует фильтру сжатия сигнала. Но если согласованный фильтр и фильтр сжатия рассматривать в качестве предельных случаев при полной неопределенности характеристики помех, то в качестве модели помех можно принять их суперпозицию:

W q = a 2 | S (w) | 2 + b 2 {\displaystyle Wq=a^{2}|S(w)|^{2}+b2\,\!}

Подставляя это выражение в (5.5) {\displaystyle \color {Maroon}(5.5)\,\!} , с точностью до множителя получаем:

H (ω) = S ∗ (ω) | S (ω) | 2 + γ 2 (5.8) {\displaystyle H(\omega) = \frac{S*(\omega)}{|S(\omega)|^2 + \gamma^2} \qquad \color{Maroon} (5.8) \,\!} где γ = b a {\displaystyle \gamma ={\tfrac {b}{a}}\,\!} - отношение дисперсий шума и сигнала.

Фильтр стремится к согласованному при больших γ {\displaystyle \gamma \,\!} , и к обратному (фильтру сжатия) при малых.

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

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

Критерий оптимальности

Энергия сигнала на выходе фильтра:

E s h = ∑ k s k 2 = ∑ k ⋅ (∑ n h n s k − n) 2 = ∑ k h k ∑ n h n R s (k − n) (6.1) {\displaystyle E_{sh} = \sum_{k} s_k^2 = \sum_{k} \cdot \left (\sum_{n} h_n s_{k-n} \right)^2 = \sum_{k} h_k \sum_{n} h_n R_s (k-n) \qquad \color{Maroon} (6.1) \,\!} где R s {\displaystyle R_{s}\,\!} - функция автокорреляции сигнала.

В векторной форме:

E s h = h T ¯ R s ¯ h ¯ (6.2) {\displaystyle E_{sh} = \overline{h_T} \overline{R_s} \overline{h} \qquad \color{Maroon} (6.2) \,\!}

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

E q h = ∑ k h k ∑ n h n R q (k − n) = h T ¯ R s ¯ h ¯ (6.3) {\displaystyle E_{qh} = \sum_{k} h_k \sum_{n} h_n R_q (k-n) = \overline{h_T} \overline{R_s} \overline{h} \qquad \color{Maroon} (6.3) \,\!} где R q {\displaystyle R_{q}\,\!} - функция автокорреляции помех.

При некоррелированной помехе E q h = σ 2 {\displaystyle E_{qh}=\sigma ^{2}\,\!} .

Подставим (6.2 − 6.3) {\displaystyle \color {Maroon}(6.2-6.3)\,\!} в выражение (2.4) {\displaystyle \color {Maroon}(2.4)\,\!} :

ρ = h T ¯ R s ¯ h ¯ h T ¯ R q ¯ h ¯ (6.4) {\displaystyle \rho = \frac{\overline{h_T} \overline{R_s} \overline{h}}{\overline{h_T} \overline{R_q} \overline{h}} \qquad \color{Maroon} (6.4) \,\!}

Расчет векторов операторов фильтров

Для определения значений вектора продифференцируем по h ¯ {\displaystyle {\overline {h}}\,\!} , и приравняем производную к нулю:

(h T ¯ R q ¯ h ¯) ⋅ (R s ¯ h ¯) − (h T ¯ R s ¯ h ¯) ⋅ (R q ¯ h ¯) (h T ¯ R q ¯ h ¯) 2 = 0 {\displaystyle {\frac {\left({\overline {h_{T}}}{\overline {R_{q}}}{\overline {h}}\right)\cdot \left({\overline {R_{s}}}{\overline {h}}\right)-\left({\overline {h_{T}}}{\overline {R_{s}}}{\overline {h}}\right)\cdot \left({\overline {R_{q}}}{\overline {h}}\right)}{\left({\overline {h_{T}}}{\overline {R_{q}}}{\overline {h}}\right)^{2}}}=0\,\!} (R s ¯ − ρ R q ¯) h ¯ = 0 (6.5) {\displaystyle \left (\overline{R_s} - \rho \overline{R_q} \right) \overline{h} = 0 \qquad \color{Maroon} (6.5) \,\!}

В системе уравнений неизвестны собственные значения ρ {\displaystyle \rho \,\!} матрицы и значения коэффициентов h n {\displaystyle h_{n}\,\!} , при этом система имеет N + 1 {\displaystyle N+1\,\!} ненулевых решений относительно значений ρ {\displaystyle \rho \,\!} и соответствующих этим значениям векторов h ¯ {\displaystyle {\overline {h}}\,\!} . Для определения коэффициентов фильтра приравнивается к нулю и решается относительно ρ {\displaystyle \rho \,\!} определитель матрицы (R s ¯ − ρ R q ¯) {\displaystyle \left({\overline {R_{s}}}-\rho {\overline {R_{q}}}\right)\,\!} , после чего максимальное значение ρ m a x {\displaystyle \rho _{max}\,\!} подставляется в (6.5) {\displaystyle \color {Maroon}(6.5)\,\!} и система уравнений решается относительно коэффициентов h i {\displaystyle h_{i}\,\!} вектора . При фильтрации сигнала вектор h 1 ¯ {\displaystyle {\overline {h_{1}}}\,\!} обеспечивает выделение первой по мощности главной компоненты сигнала, т.е. составляющей сигнала, которая имеет наибольшую энергию и отношение сигнал/шум. В сложных геофизических полях такая компонента, как правило, соответствует региональному фону.

В принципе, расчет может быть продолжен и для других значений ρ < ρ m a x {\displaystyle \rho <\rho _{max}\,\!} , и определены значения коэффициентов векторов , h 3 ¯ {\displaystyle {\overline {h_{3}}}\,\!} и т.д., с использованием которых могут выделяться вторая и прочие компоненты сигнала. Наиболее эффективно такой метод используется для разделения сигналов (полей) при некоррелированных помехах. В этом случае корреляционная матрица помех является единичной (единицы по диагонали, остальное - нули) и уравнение (6.5) {\displaystyle \color {Maroon}(6.5)\,\!} имеет вид:

(R s ¯ − ρ I ¯) h ¯ = 0 (6.6) {\displaystyle \left (\overline{R_s} - \rho \overline{I} \right) \overline{h} = 0 \qquad \color{Maroon} (6.6) \,\!}

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

h 0 (R s (0) − ρ) + h 1 R s (1) + h 2 R s (2) + h 3 R s (3) + . . . + h M R s (M) = 0 {\displaystyle h_{0}(R_{s}(0)-\rho)+h_{1}R_{s}(1)+h_{2}R_{s}(2)+h_{3}R_{s}(3)+...+h_{M}R_{s}(M)=0\,\!} h 0 R s (1) + h 1 (R s (0) − ρ) + h 2 R s (1) + h 3 R s (2) + . . . + h M R s (M − 1) = 0 {\displaystyle h_{0}R_{s}(1)+h_{1}(R_{s}(0)-\rho)+h_{2}R_{s}(1)+h_{3}R_{s}(2)+...+h_{M}R_{s}(M-1)=0\,\!} h 0 R s (2) + h 1 R s (1) + h 2 (R s (0) − ρ) + h 3 R s (1) + . . . + h M R s (M − 2) = 0 {\displaystyle h_{0}R_{s}(2)+h_{1}R_{s}(1)+h_{2}(R_{s}(0)-\rho)+h_{3}R_{s}(1)+...+h_{M}R_{s}(M-2)=0\,\!} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . {\displaystyle ..................................................................................\,\!} h 0 R s (M) + h 1 R s (M − 1) + h 2 R s (M − 2) + h 3 R s (M − 3) + . . . + h M (R s (0) − ρ) = 0 {\displaystyle h_{0}R_{s}(M)+h_{1}R_{s}(M-1)+h_{2}R_{s}(M-2)+h_{3}R_{s}(M-3)+...+h_{M}(R_{s}(0)-\rho)=0\,\!}

Выражение (6.5) {\displaystyle \color {Maroon}(6.5)\,\!} при малом уровне шумов позволяет вместо ФАК какого-либо определенного сигнала использовать ФАК непосредственно зарегистрированных данных (поля). Если при этом в зарегистрированных данных кроме помех присутствуют два (и более) сигналов, например, региональный фон и локальная составляющая (аномалия), то расчет векторов h i {\displaystyle h_{i}\,\!} приобретает конкретный практический смысл: после первой фильтрации оператором h 1 ¯ {\displaystyle {\overline {h_{1}}}\,\!} и выделения региональной составляющей, массив данных (исходный или с вычитанием из него региональной составляющей) может быть профильтрован повторно оператором h 2 ¯ {\displaystyle {\overline {h_{2}}}\,\!} , что позволит выделить и локальную аномалию (и т.д.). Разделение сигналов будет тем надежнее, чем сильнее они отличаются друг от друга по энергии и интервалу корреляции.