Оптимальная линейная динамическая фильтрация. Расчет векторов операторов фильтров

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

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.

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

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

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

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

где входной шум (или шум объекта) представляет собой белый шум с нулевым средним и ковариационной матрицей

Модель наблюдения или измерения задается линейным алгебраическим соотношением

. (7.3)

где шум измерения v представляет собой белый шум с нулевым средним и

. (7.4)

Ради простоты первоначальных выкладок предположим, что и некоррелированны, т. е.

Для всех , (7.5)

Начальное значение представляет случайную величину со средним значением и дисперсией , иначе говоря

; . (7.6)

Будем также полагать, что для всех .

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

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

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

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

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

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

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

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

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

Ортогональное проецирование. Теория ортогонального проецирования вкратце была рассмотрена в § 6.6. Здесь без доказательства будут представлены некоторые обобщения приведенных там результатов; они нам понадобятся в дальнейшем. Линейная оценка величины по критерию минимума дисперсии ошибки при заданном линейном пространстве наблюдений задается ортогональной проекцией на , т. е. .

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

. (7.8)

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

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

Легко показать, что вектор

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

.

Эти два выражения эквивалентны, так как содержится в пространстве наблюдений и, следовательно, не добавляется никакой дополнительной информации по сравнению с той, которая содержится в . Поскольку и ортогональны, можно воспользоваться ур-нием (7.9) и записать как . Так как , то это выражение можно представить в следующем виде:

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

Согласно определению и мы имеем

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

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

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

Если подставить ур-ние (7.12) в (7.11), то получим

Рассмотрим второе слагаемое в правой части этого уравнения. Используя ур-ние (7.8), можно записать в следующем виде:

Теперь исследуем отдельно каждый член, стоящий в правой части этого уравнения. Подставив (7.1) для , получаем для первого члена уравнения

Теперь, используя определения величин и [см. ур-ния (7.3) и (7.10)], можно записать в следующем виде:

где . Поэтому ур-ние (7.17) принимает вид

а после перемножения соответствующих членов преобразуется к виду

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

Полученное выражение можно еще более упростить, если учесть, что . При этом становится равным

Но первый член согласно лемме об ортогональном проецировании равен нулю. Поэтому ур-ние (7.18) можно записать в виде:

где Аналогичным образом можно показать, что

Если подставить уравнения (7.19), (7.20) и (7.10) в (7.16), то

Поэтому выражение для принимает вид

Этот результат можно представить в более удобной форме, если ввести обозначение

так что получаем окончательно

Величина называется коэффициентом усиления одношагового экстраполятора Калмана. Форма решения, представленного уравнениями (7.23) и (7.24), очень интересна и удобна с вычислительной точки зрения. Мы получили последовательный алгоритм вычисления по известной величине , вычисленной на предыдущем шаге, и новому наблюдению . Новая оценка здесь формируется как результат экстраполяции старой оценки и последующей коррекции при помощи взвешенного сигнала ошибки наблюдения Структурная схема экстраполятора Калмана показана на рис. 7.1б; для сравнения исходные модели сообщения и наблюдений показаны на рис. 7.1, а. Прежде чем воспользоваться полученным выше результатом, необходимо сначала найти выражение для , чтобы вычислить . Можно поступить иначе и найти . Для того чтобы определить , найдем сначала рекуррентное выражение для . Объединяя уравнения (7.1) и (7.24), получаем

Рис 7.1. Структурные схемы задачи одношагового предсказания: а) модели сообщения и наблюдений, б) устройство одношагового предсказания

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

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

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

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

Уравнение (7.26) совместно с (7.23) и (7.24) полностью определяют линейный последовательный одношаговый экстраполятор с минимальной дисперсией ошибки.

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

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

Все алгоритмы одношагового предсказания сведены в табл. 7.1.

Уравнение (7.26) можно переписать также в следующем виде:

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

Таблица 7.1. Дискретные алгоритмы одношагового предсказания

Модель сообщения

Модель наблюдений

Априорные данные

; ; ;

Алгоритм предсказания

Вычисление коэффициента усиления

Вычисление априорной дисперсии

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

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

Рис. 7.2. Структурная схема вычислений по алгоритмам предсказания

Пример 7.1 . Пусть модели сообщения и наблюдений заданы скалярными уравнениями:

; .

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

Для этого примера уравнение оценивания (7.24) принимает вид

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

Уравнение дисперсии имеет вид

Вычислим и в предположении, что у нас имеются наблюдения , . Вычисляем сначала коэффициент усиления , используя начальное условие :

; .

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

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

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

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

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

Если оценка , полученная как результат фильтрации, а именно , известна, то может быть получена как

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

Очевидно, что две оценки , основанные на наблюдении , должны быть эквивалентны. Следовательно, можно использовать ур-ние (7.28) для получения последовательного алгоритма оценивания из ур-ний (7.23), (7.24) и (7.26). Сначала подставим yp-н.ие (7.28) при в (7.24). В результате получим

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

Чтобы упростить полученное выражение, введем , определяемую как , или

если использовать для определения ур-ние (7.23). Поэтому записывается в виде

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

Это выражение можно еще более упростить, если ввести «обновляющую» величину , чтобы получить

Уравнения (7.29)-(7.31) или (7.32) совместно с ур-нием (7.26) полностью дают решение проблемы линейной фильтрации по критерию минимума среднеквадратической ошибки. Заданные начальные условия по , а именно и , используются для формирования начальных условий соответственно для и так же, как и в одношаговом экстраполяторе.

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

Если теперь подставить ур-ния (7.29) для и (7.19) и (7.20) для и в это выражение, то получим

Если воспользоваться ур-нием (7.29) для , то последнее выражение может быть переписано в виде

Согласно этому уравнению дисперсия ошибки фильтрации достаточно просто выражается через дисперсию ошибки одношагового предсказания. Использование величины позволяет также значительно упростить ур-ние (7.26). Перепишем его в виде

Воспользовавшись ф-лой (7.29) для , можно записать это выражение как

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

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

Уравнения (7.29), (7.30), (7.33) и (7.34) полностью определяют окончательный вариант дискретного фильтра Калмана. Эти уравнения сведены в табл. 7.2. Структурная схема вычислений согласно полученным алгоритмам приведена на рис. 7.3, а структурная схема дискретного фильтра Калмана - на рис. 7.4.

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

Таблица 7.2. Сводка дискретных алгоритмов фильтрации Калмана

Модель сообщения

Модель наблюдений

Априорные данные

Алгоритмы фильтрации

Вычисление коэффициента усиления

Вычисление априорной дисперсии

Уравнение для апостериорной дисперсии

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

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

Рис.7.3. Структурная схема вычислений по алгоритму фильтрации Калмана.

Рис. 7.4. Структурная схема дискретного фильтра Калмана.

Заметим, что структура фильтра Калмана, соответствующая ур-нию (7.30) и изображенная на рис. 7.4, очень напоминает структуру исходной модели сообщения, заданной ур-нием (7.1) и приведенной на рис. 7.1а. Алгоритм фильтрации строится на использовании «обновляющей» компоненты, которая содержит новую информацию, полученную в результате наблюдения.

Пример 7.2. Для иллюстрации применения алгоритма фильтрации Калмана рассмотрим двумерную модель сообщения, задаваемую уравнением

Наблюдение осуществляется согласно скалярной модели

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

Используя ур-ния (7.29) и (7.34), а также начальное условие , можно легко вычислить и , которые соответственно равны

Теперь с помощью ур-ния (7 23) можно вычислить апостериорную дисперсию

а также априорную дисперсию, которая изменяется для следующего шага согласно ур-нию (7.34) и становится равной

Рис. 7.5. Изменение коэффициентов усиления фильтра Калмана, рассмотренного в примере 7.2

Теперь можно вычислить и т.д. Компоненты вектора , при изменении от 1 до 10, показаны на рис 7.5. Отметим характерное увеличение коэффициента усиления для нечетных значений , в результате которого усиливаются относительно точные измерения. Можно заметить, что коэффициент усиления достигает своего установившегося периодически изменяющегося значения за несколько выборок. Вероятно, полезно вкратце и чисто качественно обсудить влияние соотношения величин и на , даже если трудно получить общие количественные результаты. Во-первых, здесь важны относительные значения, а не абсолютные. В частности, легко показать, что в том случае, когда , и умножаются на одну и ту же положительную скалярную постоянную, то не изменяется. Весьма приближенно можно лишь утверждать, что коэффициент усиления зависит от отношения сигнала к шуму . Элементы матрицы коэффициентов уменьшаются по мере уменьшения значений элементов матриц и [или только в ] или увеличения значений элементов матрицы . Этот результат представляется интуитивно вполне понятным, поскольку по мере уменьшения следует ожидать все меньших изменений в состоянии , а поэтому нет необходимости «отслеживать» наблюдения так точно. Аналогичным образом, если уменьшается, то повышается точность начальной оценки и потребность в информации, содержащейся в наблюдениях, снижается и, следовательно, коэффициент усиления уменьшается. С другой стороны, если будет возрастать, то коэффициент усиления снова уменьшается, препятствуя добавлению к оценке чрезмерного шума измерения. В пределе, когда стремится к нулю, как нетрудно показать, асимптотически приближается к нулю для больших значений . Когда стремится к нулю, дисперсии ошибок также стремятся к нулю и процедура оценивания становится не зависящей от наблюдения и входит в режим, известный под названием насыщения по входным данным. Этот режим может привести к серьезным проблемам расходимости. Методы коррекции расходимости подробно обсудим в разд. 8.5.

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

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

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

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

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

(7.35)

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

. (7.36)

при условии, что

(7.37)

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

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

Если рассматривать как объединение нового наблюдения и предыдущих наблюдений, то ур-ние (7.38) перепишется в виде

(7.39)

Рассмотрим числитель этого выражения. Применяя теорему умножения, можем записать

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

Применяя теорему умножения к знаменателю, полученное выражение запишем в виде

После сокращения на общую скалярную функцию вероятности получаем

(7.41)

Теперь можно определить условную плотность вероятности случайной величины при заданном путем вычисления каждого выражения для вероятности, стоящего в правой части ур-ния (7.41). Рассмотрим каждый член в отдельности, доказывая, что каждая плотность вероятности, входящая в (7.41), нормальная, и определяя первые два момента, характеризующие нормальное распределение. Исследуем сначала . Так как задается уравнением , a - нормальный случайный процесс, то плотность вероятности несомненно является нормальной, поскольку есть сумма нормального случайного процесса и постоянной величины . Среднее значение процесса равно

поскольку - случайный процесс с нулевым средним значением. Дисперсия случайного процесса равна по определению

а в данном случае

Отсюда плотность вероятности можно записать в следующем виде:

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

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

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

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

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

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

Следует отметить, что можно было бы непосредственно воспользоваться исходной формой записи плотности [выражением (7.52)], а не компактной формой (7.53). Такой подход представляется более привлекательным, так как в этом случае не требуется знание более компактной формы, которая не является достаточно простой и очевидной. Если воспользоваться выражением (7.52) для , то в результате преобразования ур-ния (7.57) имеем

Если теперь сгруппировать члены, включающие в себя , то получим

решение которого относительно приводит к следующему результату:

Хотя это решение для оптимальной оценки представлено не в такой удобной форме, как предыдущее, оно легко может быть приведено к (7.62), если воспользоваться леммой об обращении матриц либо непосредственно выражениями (7.55) и (7.56).

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

С ур-нием (7.70), получаем [которое представляет собой также оптимальную оценку действует на выходе системы, т. е. когда модель наблюдения имеет следующий вид:

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

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

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

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

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

Белый шум

Белый шум является стационарным случайным процессом , у которого автокорреляционная функция описывается дельта - функцией Дирака и, соответственно, спектральная плотность мощности не зависит от частоты и имеет постоянное значение 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}}}\,\!} , что позволит выделить и локальную аномалию (и т.д.). Разделение сигналов будет тем надежнее, чем сильнее они отличаются друг от друга по энергии и интервалу корреляции.

Как известно, сущность фильтрации состоит в не­прерывном оценивании изменяющихся во времени парамет­ров случайного процесса. Если сообщение является ска­лярным марковским процессом (для стационарного гауссовского процесса это означает, что ковариационная функ­ция имеет вид Aexp(-B|t-u|), то решение задачи может быть основано на следующих принципах, упрощающих дости­жение цели:

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

Линейную систему, генерирующую сообщение, следует описывать при помощи дифференциального уравнения, решением которого является искомое сообщение;

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

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

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

где -матричный коэффициент усиления опти­мального фильтра.

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

где - корректирующий сигнал вида:

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


Следовательно, алгоритм динамической фильтрации основан на классическом принципе регулирования по от­клонению с матричным коэффициентом усиления K(t), обеспечивающим минимальную среднюю квадратическую ошибку фильтрации. Корректирующий сигнал состоит из текущего сигнала наблюдения z(t) за сос­тоянием исходной системы, дополненного текущим сигна­лом состояния модели исходной системы. Сигнал является сигналом коррекции ошиб­ки фильтрации и характеризует дополнительную информацию между текущими измерениями z(t) и оценками состояния, полученными по результатам оценок , предшествующих текущим измерениям z(t). Матричная cxeма оптимального фильтра Калмана имеет вид, показанный на рис. 4.18. Эта схема реализует алгоритм динамической фильтрации, когда состояние исходной системы задается дифференциальными уравнениями, правая часть которых не зависит от наблюдения.

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

Рис. 4.17 . Матричная схема оптимального фильтра Калмана

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

Матричная схема дискретного фильтра Калмана по­казана на рис. 4.19 совместно с моделями исходной динамической системы и измерительной системы.

Рис. 4.18. Матричная схема дискретного фильтра Калмана

Основой для вывода уравнения фильтрации являются уравнения состояния динамической системы и уравнение наблюдения (измерения). Уравнение состояния линейной динамической систе­мы описывается системой разностных уравнений в векторно-матричной форме:

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

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




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

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


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

Подставив в это выражение (9.14), получаем уравнение дискретного фильтра Калмана в канонической форме:

О птимальный коэффициент передачи такого фильтра должен обеспечить минимум средней квадратической ошиб­ки фильтрации в соответствии с условием (4.152).

Контрольные вопросы к Главе 4

1. Какие критерии принятия решения применяются в ГАС НК?

2. В чём сходство и отличие критериев обнаружения «Идеального наблюдателя», «Неймана – Пирсона» и «Вальда»?

3. Какова физическая сущность вероятностей правильного обнаружения, правильного необнаружения, пропуска сигнала и ложной тревоги?

4. Как соотносится вероятность ложной тревоги «в точке» и многоканальной системы?

5. Как выбирается порог обнаружения при реализации критерия Неймана-Пирсона?

6. Как выбирается порог обнаружения при реализации критерия Котельникова-Зигерта?

7. Как выбирается порог обнаружения при реализации критерия обнаружения Вальда?

8. В чём адекватность и особенности корреляционного приёмника и согласованного фильтра?

9. В чём суть состоятельности оценки?

10. В чём суть эффективности оценки?

11. В чём суть несмещённости оценки?

12. Что представляет собой информационная матрица Фишера?

13. Как строится пеленгационная характеристика гидролокатора?

14. Как формируется словарь признаков и алфавит образов объектов гидролокации?

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

480 руб. | 150 грн. | 7,5 долл. ", MOUSEOFF, FGCOLOR, "#FFFFCC",BGCOLOR, "#393939");" onMouseOut="return nd();"> Диссертация - 480 руб., доставка 10 минут , круглосуточно, без выходных и праздников

Бирюков Руслан Сергеевич. Дискретное обобщенное H-оптимальное управление и фильтрация в линейных непрерывных объектах: диссертация... кандидата Физико-математических наук: 01.01.09 / Бирюков Руслан Сергеевич;[Место защиты: ФГАОУВО «Национальный исследовательский Нижегородский государственный университет им. Н.И. Лобачевского»], 2017

Введение

Глава 1. Обзор теории обобщенного -управления и фильтрации для линейных дискретных систем 8

1. Обобщенная -норма линейного объекта 8

2. Синтез обобщенного -управления 11

3. Синтез обобщенного -фильтра 13

Глава 2. Обобщенная -норма непрерывного объекта с дискретным целевым выходом 15

1. Уровень гашения возмущений в непрерывно-дискретном объекте 15

2. Наихудшие внешние возмущения и начальное состояние в непрерывно дискретном объекте 28

3. Уровень гашения возмущений в дискретно-дискретном объекте 32

4. Наихудшие внешние возмущения и начальное состояние в дискретно дискретном объекте 49

5. Уровень гашения возмущений в случае бесконечного горизонта 56

6. Характеризация уровня гашения возмущений в терминах LMI 61

7. Выводы 64

Глава 3. Дискретное обобщенное -оптимальное управление 66

1. Синтез оптимального управления по состоянию 66

2. Синтез оптимального управления по выходу 74

3. Управление электромагнитным подвесом 94

4. Выводы 101

Глава 4. Дискретная обобщенная -оптимальная фильтрация 102

1. Синтез оптимального фильтра 102

2. Фильтрация данных в задаче гашения колебаний зданий 108

3. Выводы 114

Заключение 115

Список литературы

Введение к работе

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

Для решения указанной задачи были предложены различные подходы. Одним из первых был подход, основанный на представлении исходной непрерывной системы с дискретным выходом как непрерывно-дискретной, поведение которой описывается совокупностью дифференциальных и разностных уравнений (Sun W., Nagpal K.M., Poolla K.R., Khargonekar P.P., Sagfors M.F., Toivonen H.T. и др.). В этом случае процедура синтеза дискретных 7^^-оптимальных регуляторов и фильтров основывалась на дифференциальных уравнениях Риккати, решения которых испытывают скачки в моменты времени, соответствующие наблюдениям. Практическая реализация предложенных алгоритмов синтеза наталкивается на ряд трудностей, связанных с решением нелинейной краевой задачи для дифференциальных уравнений Риккати.

Похожий подход использовался в работах Basar T. и Bernhard P., где задача дискретного ^^-оптимального управления непрерывным объектом рассматривалась с точки зрения теории игр. Условия существования %^-оптимальных регуляторов были сформулированы в случае измеряемого состояния объекта в терминах разностных уравнений Риккати, а процедура синтеза таких регуляторов также основана на решении нелинейной краевой задачи.

Другой подход основан на использовании метода лифтинга, в котором исходная непрерывная система преобразуется в эквивалентную дискретную (Bamieh B.A., Pearson J.B., Chen T., Francis B.A., Tadmor G., Sagfors M.F., Toivonen H.T., Lall S., Dullerud G. и др.). При этом, поскольку между моментами наблюдения внешнее возмущение, как и целевой выход исходного объекта, представляют собой кусочно-непрерывные функции, то возмущение и целевой выход эквивалентной дискретной системы уже принадлежат бесконечно-

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

Наконец, в работах Михеева Ю.В., Соболева В.А., Фридман Э.М., Shaked U., Suplin V. был предложен подход при котором задача синтеза дискретного -управления непрерывным объектом формально заменялась задачей синтеза -регулятора с запаздыванием. Условия существования -управления были сформулированы в форме достаточных условий в терминах линейных матричных неравенств.

Одним из существенных недостатков теории -управления является предположение о том, что в начальный момент времени объект находится в покое, то есть его начальное состояние нулевое. Если это требование не выполняется, то синтезированные регуляторы хорошо подавляют внешние возмущения, но не всегда адекватно справляются с задачей гашения начальных возмущений, порожденных ненулевыми начальными условиями. В этом случае в качестве единого критерия, учитывающего влияние как внешних, так и начальных возмущений, была предложена обобщенная -норма (Khargonekar P.P., Nagpal K.M. и Poolla K.R.). Эта норма совпадает с классической -нормой, если в начальный момент времени объект находится в покое, а когда начальное состояние объекта ненулевое, а внешнее возмущение отсутствует, то обобщенная -норма совпадает с 0 -нормой, определенной в работах Баландина Д.В. и Когана М.М. Для непрерывных объектов с непрерывным измеряемым выходом были синтезированы непрерывные законы управления и фильтрации в работах Khargonekar P.P., Nagpal K.M., Баландин Д.В., Коган М.М. и др. В случае непрерывного объекта с дискретным выходом известна работа Sun W., Nagpal K.M. и Khargonekar P.P., в которой решение задачи дискретного обобщенного -управления было получено для объекта на бесконечном горизонте. При этом сформулированные законы управления и фильтрации основаны на решении нелинейного дифференциального уравнения Риккати, что затрудняет их использование. Таким образом, дальнейшее развитие теории дискретного обобщенного -управления непрерывными системами является весьма актуальной задачей теории управления.

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

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

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

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

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

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

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

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

    Для линейных нестационарных объектов на конечном интервале времени получены необходимые и достаточные условия, а в случае неизмеряемого состояния только достаточные условия существования дискретных обобщенных -оптимальных законов управления. Эти законы управления синтезированы в классе линейных нестационарных обратных связей по состоянию и в классе линейных нестационарных динамических регуляторов по выходу (соответствует пункту 6 паспорта специальности 01.01.09).

    Для линейных стационарных объектов на бесконечном интервале времени получены необходимые и достаточные условия существования дискретных обобщенных -оптимальных законов управления. Эти законы управления синтезированы в классе линейных стационарных обратных связей по состоянию и в классе линейных стационарных динамических регуляторов по выходу (соответствует пункту 6 паспорта специальности 01.01.09).

    Для линейных нестационарных объектов на конечном (бесконечном) интервале времени получены необходимые и достаточные условия существования и осуществлен синтез нестационарных (стационарных) дискретных обобщенных "Н^ -оптимальных фильтров полного порядка в форме наблюдателя

    В качестве приложений синтезированы дискретные обобщенные Ti^ -оптимальные регуляторы в задаче управления телом в электромагнитном подвесом и дискретные обобщенные ^-оптимальные фильтры в задаче гашения колебаний высотных зданий и сооружений (соответствует пункту 6 паспорта специальности 01.01.09).

Соответствие шифру специальности. Работа соответствует формуле специальности 01.01.09 - Дискретная математика и математическая кибернетика и охватывает следующие области исследования, входящие в специальность 01.01.09: п. 6. Математическая теория оптимального управления.

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

Степень достоверности и апробация результатов исследования. Основные результаты диссертационной работы обсуждались на заседании Нижегородского научного семинара «Математическое моделирование динамики систем и процессов управления» в НИИ Прикладной математики и кибернетики, а также докладывались на следующих международных и всероссийских конференциях:

X Всероссийская научная конференция «Нелинейные колебания механических систем» им. Ю.И. Неймарка (Нижний Новгород, 2016);

XIII Международная конференция «Устойчивость и колебания нелинейных систем управления» (конференция Пятницкого) (Москва, 2016);

ХI Всероссийский съезд по фундаментальным проблемам теоретической и прикладной механики (Казань, 2015);

Международная конференция по математической теории управления и механике (Суздаль, 2015);

Шестая традиционная всероссийская молодежная летняя школа «Управление, информация и оптимизация» (Москва, 2014);

XII Всероссийское совещание по проблемам управления (Москва, 2014);

XIX Нижегородская сессия молодых ученых: Естественные, математические науки (Нижний Новгород, 2014).

В 2013-2014 гг. и 2014-2015 гг. исследования были поддержаны стипендией имени академика Г.А. Разуваева для аспирантов, а также стипендией Правительства Российской Федерации (2014-2015 гг).

Результаты первых трёх глав диссертации были получены при выполнении проекта № 14-01-31120 мол_а в 2014-2015 гг. (руководитель) и проектов № 12-01-31358 мол_а в 2012-2013 гг., № 14-01-00266 в 2014-2016 гг. (исполнитель), выполненных при финансовой поддержке Российского фонда фундаментальных исследований.

Результаты четвертой главы получены при финансовой поддержке Министерства образования и науки РФ в рамках Федеральной целевой программы «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2014-2020 годы» (соглашение 14.578.21.0110 от 27.10.2015, уникальный идентификатор RFMEFI57815X0110).

Публикации. Основные результаты по теме диссертации изложены в 10 печатных работах, включая 4 публикации в ведущих научных журналах, рекомендованных ВАК Минобрнауки РФ -], трудах двух международных конференций и четырех тезисах докладов региональных и Всероссийских конференций [-. В совместной работе ] автору принадлежат результаты численного моделирования.

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

Структура и объем работы. Диссертация состоит из введения, четырех глав, заключения и списка литературы. Работа изложена на 123 страницах, содержит 11 иллюстраций. Библиография включает 81 наименование.

Синтез обобщенного -управления

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

Приведем теперь основные факты, относящиеся к обобщенной Ноо-норме, при этом в изложении будем следовать работам . Для определенности рассмотрим линейный дискретный нестационарный объект вида Xk+i = Акхк + Bkvk, k = 0,...,N-l, zk = Ckxk + Dkvk, где х Є Ж1 - состояние, z Є Е"-2 - целевой выход иие Rnv - внешнее возмущение, N-l т ограниченное по 2-норме: vk vk oo. fc=0

Предположим, что в общем случае начальное состояние х0 ненулевое и неизвестно, а его влияние на динамику объекта интерпретируется как начальное возмущение.

Управляемый выход объекта для фиксированного начального состояния х0 и последовательности возмущений v0,... , vN_ і будем характеризовать значением функционала N-1 j(x0,v0,..., vN_ij = \\z\\i2 + xNSxN = У zk zk-\- xNSxN, (1.2) fc=0 где S = S 0 - весовая матрица, задающая приоритет между качеством переходного процесса и конечным состоянием объекта.

Сначала рассмотрим отдельно два крайних случая: на объект действует только начальное или только внешнее возмущение. Пусть объект в начальный момент времени находился в покое, что соответствует случаю, когда отсутствует начальное возмущение. Следуя , определим показатель влияния внешнего возмущений на целевой выход (1.1) - уровень гашения внешнего возмущения - как относительное значение функционала (1.2) в наихудшем случае: J(0,VO,...,VN_1) 2 = sup 2 0 2

Отметим, что если объект (1.1) является стационарным и рассматривается на бесконечном интервале времени, то, используя равенство Парсеваля, можно показать, что выражение (1.3) совпадает с 7 -нормой рассматриваемого объекта . Следующее утверждение характеризует уровень гашения внешнего возмущения в терминах решений линейных матричных неравенств .

Утверждение 1.1. Уровень гашения внешнего возмущения в системе (1.1) на конечном интервале времени удовлетворяет неравенству 7оо 7 тогда и только тогда, когда линейные матричные неравенства /AlXk+1Ak - Xk AjXk+lBk Ck\ BTkXk+lAk BjXk+1Bk--f2I Dj Ck Dk -i) 0, (1.4) разрешимы относительно матриц Xk = Xk 0, k = 0,..., N - 1, при XN = S.

Из утверждения следует, что уровень гашения внешнего возмущения 7оо находится как точная нижняя грань множества всех 7, для которых система линейных матричных неравенств (1.4) разрешима относительно матриц Хк = Хк 0 и 7 .

В случае, если внешнее возмущение отсутствует, то влияние начального возмуще ния на качество переходного процесса в системе (1.1) может быть охарактеризовано величиной 2 J(x0,0,...,0) 70 = sup 2 (1.5) х0ф0 \Х0\ которая называется уровнем гашения начального возмущения . В показано, что эта величина может быть найдена как решение оптимизационной задачи с ограничениями, заданными линейными матричными неравенствами.

Утверждение 1.2. Уровень гашения начального возмущения в системе (1.1) на конечном интервале времени удовлетворяет неравенству 70 7 тогда и только тогда, когда линейные матричные неравенства ATkXk+1Ak -Xk + ClCk О, Х0 -f2I, (1.6) разрешимы относительно матриц Хк = Хк О, к = 0,..., N - 1, при XN = S. Чтобы описать совместное влияние внешнего и начального возмущений на выход объекта (1.1), определим уровень гашения возмущений как своеобразную свертку двух рассматриваемых факторов : 7W = sup

Jx0,v0,. . . ,VN_1 =F , (1.7) где R = R 0 - весовая матрица, предназначенная для задания приоритета между внешним возмущением и компонентами начального состояния. Введенный таким образом показатель называется обобщенной 7 -нормой. Нетрудно видеть, что в крайних случаях выражение (1.7) превращается либо в (1.3), либо в (1.5), то есть, при х0 = 0 имеем 7w = 7оо, а при v = 0 получим 7«, = 70/ тах(-). Оказывается , что уровень гашения возмущений может быть выражен в терминах линейных матричных неравенств, для этого достаточно потребовать существование общего решения неравенств (1.4) и (1.6), характеризующих в отдельности уровень гашения внешнего возмущения и уровень гашения начального возмущения с учетом весового коэффициента.

Уровень гашения возмущений в системе (1.1) на конечном интервале времени удовлетворяет неравенству 7w 7 тогда и только тогда, когда линейные матричные неравенства (A.

Наихудшие внешние возмущения и начальное состояние в непрерывно дискретном объекте

Отметим, что согласно сформулированной теореме, уровень гашения возмущений 7С при помощи соотношения (2.45) выражается через значение матричной функции X(t). Однако, в силу уравнения (2.6a), величина X(t) неявно зависит от гус. Вследствие этого, для определения уровня гашения возмущений возникает нелинейная краевая задача для матричного дифференциального уравнения Риккати: найти решение уравнения (2.6a) с граничными условиями (2.6b) и (2.45), а также условием (2.6d).

Обратимся теперь к доказательству теоремы.

Доказательство теоремы 2.2. Нетрудно показать, что соотношение (2.4) эквивалентно выполнению равенства sup J(xo,v,w) = 0. (2.48) Иі!2 +ІНІ2 2 + 0Д 0=і Согласно формуле (2.39) функционал J(x0,v,w) может быть записан следующим образом: J(x0,v,w) = xUcJC0 + X(t0) - %R)X0 %\\v - v \\l2 + N-l + J2(wk w k)T(AjX(tk)Ak - %l) (wk - w k) + fc=l + wN - w N (АдгбАдг - 7C- (wN - w N где v и w k определены соотношениями (2.46). В силу справедливости неравенств (2.6b), (2.6c) и (2.6e) первое слагаемое является неположительно определенной, а оставшиеся - отрицательно определенными квадратичными формами, поэтому максимальное значение функционала J(x0, v, w) обращается в ноль при v = v и wk = w k, k = 1,... , N, и соответствующем выборе х0. Следовательно, возмущения v иwl являются наихудшими внешними возмущениями относительно критерия 7с. Подставим v и w k в соотношение (2.48), тогда: sup J(x0,v ,w)= sup xUx(t0) + CjC0--fcR)x0. \\v \\L+\\v \\2+x0R 0=l ll« llL+lh ll2+ ftr0 = l Теперь заметим, что и, и зависят от 0 и справедливы соотношения: v (t) = ъ1В (t)X(t) b(t,t0)x0, / -г- \ -1 -г w k = - (AjX(tk)Ak - 7c Л AjX(t (tk - 0, t0)x0, здесь Ф(Мо) - фундаментальная матрица решений замкнутой системы (2.115). Следовательно, ограничение есть квадратичная форма от х0: \\v \\l2 + \\w \\l+xUxo = x Qx0 = l, где tN Q = R + 1-2 Фт(т,і0)Х(т)В(т)Вт(т)Х(т)Ф(т,і0)(іт + «о N + fc=l J2 фТ(- 0, t0)X{tk)Ak(AlX(tk)Ak - 7сЛ \іХ(ік)Ф(ік - 0, t0). Таким образом, задача (2.48) свелась к следующей: sup x0 ПІКО =1 Xo(x(t0) + C0TC0 - 7ci?W

Для решения последней задачи воспользуемся правилом множителей Лагранжа: точка максимума х0 должна удовлетворять системе уравнений: (x(to) + СоТС0 - jcR\x0 +»Пх0 = 0 и ж аго = 1, (2.49) параметр /і есть множитель Лагранжа. Перепишем первое уравнение как (X(t0) + CQ С0 + /іП)х0 = lcRxo, откуда находим хо = «emax (R 1 \x(t0) + CjC0 + ц Ы V 7с = Атах (і?-1 [х(0) + С0ТСо + /х fil V значение а находится из второго уравнения (2.49). Подставим найденные значения в квадратичную форму и упростим: Xo(x(t0) + Со Со - 1CR)XQ = -iixoflx о = Iі Заметим, что по условию точная верхняя грань равна нулю, следовательно /і = 0. Подставляя найденное значение /і в выражение для х0 приходим к соотношениям (2.45) и (2.46c).

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

Следствие 2.5. В объекте (2.1), (2.2) и (2.3) уровень гашения начальных возмущений с = тах (J0 + (0)) (2.50) достигается при наихудшем начальном состоянии = max (J0 + (0)] , (2.51) где () - решение системы (2.41), найденное при с. Доказательство. Поскольку на объект не действует ни непрерывное, ни дискретное внешнее возмущение, то соотношение (2.51) получается из соотношений (2.46), если положить в последних = , () = 0 и к = 0, = 1,... , . Ш Следствие 2.6. В объекте (2.1), (2.2) и (2.3) уровень гашения непрерывных внешних возмущений ї = max (J0 + (0)) (2.52) достигается при наихудшем внешнем возмущении () = ") 1T()()(), (2.53) где () - решение системы (2.42), найденное при с.

Доказательство. Соотношение (2.53) получается из соотношений (2.46), если положить в последних к = 0, = 1,... , что равносильно тому, что на объект не действует дискретное внешнее возмущение, а в силу отсутствия начального возмущения необходимо отбросить условие (2.46c) и положить = в соотношении (2.45).

Следствие 2.7. В объекте (2.1), (2.2) и (2.3) уровень гашения дискретных внешних возмущений с = max (j 0 + (0)) (2.54) достигается при наихудшем внешнем возмущении / -г- \ -1 -г k = - (j(k)k - с") j(k)(k - 0), (2.55) где () - решение уравнения (2.43) с условиями (2.6b) и (2.6d), найденное при с" . Доказательство. Так как на объект не действует непрерывное внешнее возмущение, то соотношение (2.55) получается из соотношений (2.46), если положить в послед них В {і) = О, а в силу отсутствия начального возмущения необходимо отбросить усло вие (2.46c) и положить R = І в соотношении (2.45). Следствие 2.8. В объекте (2.1), (2.2) и (2.3) уровень гашения смешанных внешних возмущений lT = Amax (cJC0 + X(t0)) (2.56) достигается при наихудших внешних возмущениях (т- \ -1 -г AjX(tk)Ak - fc wl\ AjX(tk)x(tk-0), (2.57a) v (t) = (w) 1BT(t)X(t)x(t), (2.57b) где X{i) - решение системы (2.6a), (2.6b) и (2.6d), найденное при % w. Доказательство. Поскольку на объект не действует начальное возмущение, то, отбрасывая в соотношениях (2.46) условие (2.46c) и полагая R = І в формуле (2.45), получаем соотношения (2.57).

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

Синтез оптимального управления по выходу

Сгруппируем первое и второе слагаемые в (2.105) и упростим выражение для П2, для чего опять применим формулу Шермана-Моррисона-Вудбери, тогда: "Г 1 Г / 1 -г \-1 1 1 -г CTkGk+lWklx \l + Ek+lXk+l(І - Ек+1\к-ІгЕтк+1Хк+1) Ек+1]кЦ GTk+lCk = -г / -г \-1 -г = CTkGk+l(Wk+l - ETk+lXk+lEk+l) GTk+lCk и П2 = ATkXk+l l - Ek+l , формируется матрица S, например, по формуле

Теорема 3.4 также позволяет синтезировать и обобщенное ft -оптимальное управление по выходу на бесконечном интервале времени. Для этого достаточно найти решение задачи минимизации 7с() при ограничениях, задаваемых неравенствами (3.51), после чего оптимальный регулятор находится как решение (3.52).

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

Следствие 3.13. Для стационарного объекта (3.21), (3.22) при заданном 7 0 существует дискретное -управление по выходу на бесконечном интервале времени тогда и только тогда, когда линейные матричные неравенства Ah,XAh 0, С1 AhYAl Y C1YAl (Wc2 0 0 II МТ X I C1YCj WC 0 0 I M 0, (3.53a) (3.53b) x l Y 0, X yl, (3.53c) разрешимы относительно X = X 0, Y = У 0, при этом столбцы матриц Wr KJ 2 и M образуют базисы ядер матриц соответственно.

Следствие 3.14. Для стационарного объекта (3.21), (3.22) на бесконечном интервале времени существует дискретное И -управление по выходу, обеспечивающее гашение непрерывных внешних возмущений с заданным 7 0, тогда и только тогда, когда линейные матричные неравенства и первое неравенство (3.51c) разрешимы относительно X = X О, Y = У О, а столбцы матриц Wr и М образуют базисы пространств кет Со и кет [В.. D-,) соответственно.

Следствие 3.15. Для стационарного объекта (3.21), (3.22) на бесконечном интервале времени существует дискретное И -управление по выходу, обеспечивающее гашение дискретных внешних возмущений с заданным 7 О, тогда и только тогда, когда существуют матрицы X = X О, Y = У О, удовлетворяющие линейным матричным неравенствам и первому неравенству (3.51c), при этом столбцы матриц N и М образуют базисы пространств ker (С2 D2j и ker [Ви Dx) соответственно. Следствие 3.16. Для стационарного объекта (3.21), (3.22) на бесконечном интервале времени существует дискретное И -управление по выходу, обеспечивающее гашение смешанных внешних возмущений с заданным 7 0, тогда и только тогда, когда линейные матричные неравенства (3.51а), (3.51Ь) и первое неравенство (3.51с) разрешимы относительно матриц ХТ = Х 0 и Y = У О, при этом столбцы мат риц N и М образуют базисы пространств кег С2 О D2 0 и кег ВІ Dj О О соответственно.

Из замечания к теореме 2.8 следует, что существует такая конечная матрица R , что при любом весовом коэффициенте R R обобщенный Я -оптимальный регулятор по выходу на бесконечном интервале времени совпадает с // -оптимальным регулятором по выходу синтезированным по следствию 3.16 и обеспечивающим гашение смешан-ных внешних возмущений. Следовательно, для получения действительного компромисса при учете влияний как начального, так и внешних возмущений весовая матрица R должна удовлетворять условию Атах(Л_1Л) І. Численно граничное значение R весовой матрицы определяется следующим образом: brx y1 , где через Х обозначена матрица, удовлетворяющая неравенствам (3.51a), (3.51b) и первому неравенству (3.51c) при минимальном значении 7с.

Рассмотрим изображенную на рис. 3.3 механическую систему, состоящую из выве-шиваемого тела массы т и электромагнита . Левитация тела обеспечивается изменением магнитного поля, происходящим за счет изменения напряжения U, подаваемого на обмотку электромагнита. Динамика такого простейшего магнитного подвеса подчиняется двум уравнениям: ті) = F - та, (3.56) V + RI=U. Первое уравнение (3.56) выражает второй закон Ньютона и определяет изменение ко-ординаты s вывешиваемого тела под действием силы тяжести тд и силы F со стороны электромагнита, а второе - определяет изменение силы тока / в цепи электромагнита сопротивлением R при изменении подаваемого на него напряжения U и представляет собой закон Кирхгофа для электрической цепи электромагнита. Через Ф обозначено потокосцепление обмотки электромагнита, Ф = пФ, где Ф - магнитный поток, проходящий через один виток, а п - число витков в обмотке.

Потокосцепление Ф и сила тока / в цепи электромагнита связаны соотношением: = L(s)/, L{s) = , CL = /i0n2A/2, (3.57) где L(s) - индуктивность электромагнита, CL - конструктивный параметр и 6 - величина номинального зазора между электромагнитом и вывешиваемым телом. Если обозначить номинальную индуктивность как L0 = L(0), то С = L05, и тогда



Понравилась статья? Поделитесь с друзьями!