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

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

Эти способы позволяют при разумных затратах обработать большие объемы сейсмической информации; на практике ежедневно приходится решать не менее 5 млн. нормальных уравнений, многие из которых содержат до 100 и более переменных. При рассматриваемом подходе применяется особая форма автокорреляционной матрицы в скалярном случае называемая тёплицевой формой, а в матричном случае — блок-тёплицевой формой. Эту форму можно представить в виде матрицы

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

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

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

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

Прямая рекуррентная схема для скалярного процесса была впервые составлена Левинсоном [5]. Робинсон [15] обобщил ее на случай многоканальных записей и, наконец, Виггинс [16] распространил скалярную схему на случай многомерной информации. Рекуррентная косвенная схема предложена Симпсоном [17].

В случае одноканальных сигналов нормальные уравнения имеют вид

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

где — средняя квадратическая ошибка, а — дельта-функция Кронекера, по определению равная при при

Приведенный здесь вариант рекуррентного метода, относящийся к большим операторам, является доработанным вариантом метода Левинсона. Левинсон пользовался операторами предсказания, а не операторами ошибки предсказания. Кроме того, в первоначальном алгоритме Левинсона на каждом этапе необходимо вычислять пошарные скалярные произведения трех векторов. Одно из скалярных произведений используется для вычисления следующего значения дисперсии ошибки предсказания Поскольку ошибки предсказания могут становиться очень малыми, накопление ошибок округления может сделать вариант алгоритма, предложенный Левинсоном, неустойчивым. Как будет показано ниже, модификация алгоритма, предложенная Виггинсом и Робинсоном [18], позволяет обойти эту трудность при вычислении Эта модификация и ее значение рассматриваются также в работе Бёрга [13].

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

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

Величина и, равная

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

Чтобы найти множитель приравняем нулю сумму Тогда

Таким образом, равно отношению расхождения и к дисперсии ошибки предсказания у, взятому с обратным знаком. Новый оператор имеет вид

а новая дисперсия

Воспользуемся новым оператором ошибки предсказания для увеличения числа коэффициентов фильтра Опять в качестве первого приближения к добавим нуль в конец характеристики оператора Тогда

где

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

Выберем так, чтобы

Тогда характеристика нового фильтра будет иметь вид

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

— данный оператор ошибки предсказания,

— новый оператор ошибки предсказания,

— данный оператор ошибки вспоминания,

— новый оператор ошибки вспоминания,

отрезок автокорреляционной последовательности,

характеристика данного фильтра,

— характеристика нового фильтра.

Прежде всего вычисляется расхождение в виде скалярного произведения

а затем в виде дроби со знаком минус:

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

и находится дисперсия новой ошибки предсказания

(Данная операция характерна для прямого способа решения задачи и является существенным отличием от метода Левинсона, где величина вычислялась в виде скалярного произведения

Новый оператор ошибки вспоминания 6 получается путем расстановки коэффициентов нового оператора ошибки предсказания в обратном порядке. Затем вычисляется скалярное произведение

и постоянная

Наконец, получается характеристика нового фильтра

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

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

— данный оператор ошибки вспоминания,

— новый оператор ошибки вспоминания.

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

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

Сначала вычисляются скалярные произведения

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

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

Далее вычисляются весовые множители

Затем составляются новые операторы

и определяются новые значения дисперсий

После этого вычисляются скалярное произведение

и коэффициент

Наконец, получается новая характеристика фильтра

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