7.7. Двумерные формирующие фильтры

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

Сейсмограмма состоит из набора записей, описывающих колебания, происходящие в Земле, как функции времени. Это могут быть записи выходных сигналов отдельных датчиков или комбинации сигналов от группы датчиков, определенным образом размещенных в пространстве. Если точное положение этих датчиков в явном виде учитывать не нужно, то единственной независимой переменной является время. Если же пространственные координаты должны фигурировать в явном виде, то число независимых переменных может дойти до четырех: это время и, возможно, три пространственных координаты. Будем считать, что размерность процесса равна числу независимых переменных, а его порядок — числу зависимых переменных, описывающих процесс в каждой точке пространства. Так, например, совокупность сигналов, снимаемых с линейной решетки сейсмометров, измерявших колебания в Земле по трем осям, образует двумерный процесс третьего порядка. Более подробно эти вопросы рассмотрены в работе Виггинса [16].

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

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

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

или

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

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

Поэтому

и

Однако можно также записать

что

или

где

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

или после транспонирования обеих частей имеем

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

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

имеет размеры .

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

1) (-канального сигнала длительностью

2) одномерного -канального фильтра с характеристикой длиной

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

т. е. является полиномиальной с размерами причем из всех ее элементов только независимы, а

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

что приводит к эквивалентной многоканальной системе

где

Если массивы А и В имеют размеры соответственно, то в рассматриваемом случае эквивалентная многоканальная система состоит из

1) -канального входного сигнала длительностью ;

2) -канального одномерного фильтра с характеристикой длиной т.

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

причем только полиномов являются линейно-независимыми, а

элементов равны нулю.

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

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

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

Если же воспользоваться уравнением (7.8), то задача формулируется следующим образом:

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

где — векторы, образованные столбцами требуемого и фактического С выходных массивов соответственно.

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

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

Необходимо получить фильтр А с размером характеристики (), чтобы выходной массив С аппроксимировал заданный выходной массив размером ():

с минимальной квадратической ошибкой.

Если воспользоваться формулировкой (7.9), то

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

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

а взаимно-корреляционная матрица

Таким образом, искомые матричные коэффициенты автокорреляции размерами () равны

а искомые коэффициенты взаимной корреляции размерами () равны

Соответствующие нормальные уравнения, которые имеют вид

можно решить с помощью многоканального блок-тёплицева алгоритма, что дает

Фактический выходной сигнал

имеет нормированную среднюю квадратическую ошибку, равную

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

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

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

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

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

В последнее время в литературе уделялось достаточно много внимания конструированию двумерных рекурсивных фильтров (см., например, [19 и 20]), поэтому здесь нецелесообразно детально повторять эти работы. Однако один из важнейших аспектов двумерных рекурсивных фильтров, а именно вопрос об их устойчивости, имеет отношение к рассмотренному выше двумерному фильтру, дающему наименьшую среднюю квадратическую ошибку. Так же, как и в одномерном случае, неустойчивые двумерные рекурсивные фильтры создают выходные сигналы, не ограниченные по величине. Многие из существующих способов расчета приводят к созданию таких фильтров, которые нуждаются в дополнительной стабилизации. Это означает, что знаменатель дроби (7.11) [полином обратной связи] не может обращаться в нуль, когда одновременно равны или меньше единицы, а из этого условия следует, что функция является минимально-фазовой.

Разработан ряд методов стабилизации фильтров. В одном из них используется двумерный комплексный кепстр [21, 22], другой основан на вычислении двумерного преобразования Гильберта, связанного с логарифмом модуля спектра В обоих методах стараются превратить знаменатель в минимальнофазовую функцию, но в обоих случаях при реализации фильтра могут встретиться трудности.

Еще один метод основывается на постулате [24], который выполняется во многих случаях (хотя недавно найден [25] пример, в котором он несправедлив): если имеется произвольный конечный массив (матрица) X, составленный из действительных чисел, то двумерный массив (матрица), который с минимальной средней

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

Чтобы обеспечить устойчивость фильтра, необходимо выполнить следующее:

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

где — матрица с единичным импульсом, имеющая вид

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

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

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

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

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

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

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

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

откуда следует, что

Таким образом, имеет нулевую фазовую характеристику, а ее амплитудно-частотная характеристика равна квадрату амплитудной характеристики Ряд других способов синтеза фильтров с нулевыми фазовыми характеристиками описан в работе [20].