Файл: Цифровая обработка сейсмических данных..pdf

ВУЗ: Не указан

Категория: Не указан

Дисциплина: Не указана

Добавлен: 09.04.2024

Просмотров: 163

Скачиваний: 1

ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.

или

y(z)b(z) = y(z)a(z).

Перейдем к временным представлениям, заменив произведение z-трансформант сверткой соответствующих временных рядов:

у *Ъ = у *а.

Проделаем преобразования, соответствующие этой формуле, на конкретном примере. Выберем

 

 

 

a (z) = а0

+ axz

- j - a2z2,

 

 

 

 

 

 

 

b(z)^l

+ b1z+b2z2.

 

 

 

Действуя так, как показано на

рис. 80,

получим

 

Уо

 

 

 

 

= Уо^о,

 

 

 

 

УаЬ1 Jr

У1

 

 

= УоЯ] + УгЯо,

 

( 6

- 9 )

У»Ь2

+

уфх

+ у2

= у0аг

+ у^

+ У2а0,

 

 

 

 

У\h

+ УгЪ1 + Уз =

У\а<1 Л-У-г^-т Уз^о,

 

Мы видим,

что искомую последовательность у0,

уlt у2, • • •

на

выходе фильтра

можно найти из уравнений (6.9):

 

 

 

= Уоао,

 

 

 

 

 

 

 

 

VI =

Уоа1

+

У1ао—УоЬ1,

 

 

 

 

 

 

Уг = £/о«2

+

У1а1 + Угао — (УоЪг +

2/А)>

 

 

 

Уз =

У\аг

+

г«1 +

Уз"о ~ (?/А

+

J/26 i).

( 6 - 1 0 )

 

Ук = Ук-гЧ + yk-iai

+ Ука0

{Ук-Фг +

Ук-г^),

 

Или в общем виде

пm

Ук = Ъ"1Ук-1 - S

ь,-+1Ук-,--1.

(б/и)

г/=0

185


Процесс получения последовательности у0, уи у2, • • • с помощью (6.11) и есть рекурсивная фильтрация. Из (6.10) и (6.11) видно, что очередной отсчет выходного сигнала при рекурсивной фильтрации

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

с оператором а0 , а ь

. . .,

ап,

за вычетом свертки

предшествующих

отсчетов выходного

сиг­

нала

с оператором bи Ь2,

. . ., Ът. Использование предшествуютцих

значений входного сигнала и объясняет название «рекурсивная фильтрация».

 

Очевидно,

что

выходной сигнал

у0,

yt. . .

можно

получить

и

непосредственно

путем свертки входного сигнала yk с

оператором

4

= /0 , 1и

. . .,

1р.

Рекурсивная

фильтрация

может

быть

пред­

почтительной

лишь в

том случае, если

окажется, что т + п

<Ср.

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

ров с весьма «длинным» оператором lk.

Затраты машинного времени

в таких случаях могут быть в 5—10 раз меньше [115].

Проблема расчета весовых коэффициентов рекурсивных фильтров

по заданной частотной характеристике

рассмотрена в работах [92,

1151 и сводится к определению корней z полиномов числителя и зна­

менателя выражения (6.7)

С точностью до постоянного

множителя

эти корни, называемые для

числителя «нулями» и для

знаменателя

«полюсами», полностью характеризуют свойства фильтра. Действи­

тельно,

полином

всегда можно представить в

виде

(см. гл.

1)

 

 

a0

+ a x z + .

. . +

a „ 2 n = a 0 ( l

- r x z ) (1 -

r2z)

. . .

(1 -rnz),

(6.12)

где гi,

r2 , . . .

rp — корни

полинома.

(6.12) и

объединяя сла­

Произведя умножения в правой части

гаемые

с

одинаковой

степенью

z,

можем найти

непосредственно

значения

коэффициентов фильтра

й],

аг,

. . .

Множитель

а0

яв­

ляется масштабным и на форму

характеристики

не влияет. Если

мы

каким-либо способом установим положение

«полюсов» и ««нулей»

на

комплексной плоскости, то решение остальных проблем будет про­ стым.

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

мере фильтра

(6.8). Пусть

а0 = \, а{

= а2 = 0, bi =

—2,

Ъ2

-

0.

Тогда полюс

этого фильтра,

определяемый равенством

1 +

bz =

0,

будет г == —1/6 = 0,5. Импульсную

реакцию этого фильтра

опре­

делим из (6.8). Действуя по правилу деления многочленов, получим

У =

1, j / i =

2, уz =

4, у? == 8, г/4 = 16, уъ = 32 и т. д., т. е. опе­

ратор

фильтра

будет

расходящимся.

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

186


Фильтрация в области частот

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

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

Фурье

получают

ее комплексный

спектр 1

 

 

 

 

 

 

 

К-1

 

 

 

 

 

 

У (ге Дсо) =

Дг £

у (А: ДО е-'"

 

 

(6.13)

где

i — мнимая единица;

Дсо — интервал

дискретности

задания

спектра;

Л" = 0,1,

К — 1 отсчеты

трассы: п =

0,

1, . . .,

п—1 отсчеты спектра, причем К — N.

 

 

 

At

Выберем Дсо равным

2n/KAt.

Опуская

масштабные

множители

и Дсо,

имеем

 

 

 

 

 

 

 

 

 

 

К~1

 

 

 

 

 

 

У (п) = 2

У №) е-2 3 1 """* .

 

(6.14)

Затем находят комплексный спектр Y (п) выходной трассы путем перемножения Y (п) и заданной комплексной частотной характери­ стики L (п) фильтра:

Y (и) = Y (п) L (п).

После этого вычисляют выходную трассу yk путем обратного преобразования Фурье ее комплексного спектра:

 

IV-1

 

 

У* =

2 Y

(n)e^iknlN.

(6.15)

Следует подчеркнуть важность соблюдения соотношений

К = N

и Дсо = 2n/KAt при выполнении

прямого преобразования

Фурье.

Если выбрать Дсо ^>2n/KAt,

то возможны искажения типа эйлисинг-

эффекта, возникающего при декретировании непрерывных времен­ ных функций (см. гл. 1).

Фильтрации во временной и частотной области эквивалентны между собой. Выбор того или иного варианта расчетов определяется экономическими соображениями. Естественно, что перемножение спектров требует значительно меньше машинного времени, чем свертка (I.N операций). Однако вычисление прямого и обратного преобразования Фурье по обычным формулам (6.14) и (6.15) требует TV2 операций, что значительно превышает число операций при свертке.

В 1955 г. был найден новый алгоритм преобразования Фурье, получивший но имени авторов название «алгоритм Кули—Туки» или алгоритм быстрого преобразования Фурье (БПФ), который

1 Точнее, главный полупериод периодического комплексного спектра дискретной функции (см. гл. 1).

187


позволяет выполнить преобразование Фурье, используя 2NlgN операций. Это делает расчеты в частотной области более экономичными чем во времени, особенно для больших массивов чисел, характерных для задач обработки сейсмических данных.

Алгоритм быстрого преобразования Фурье. Предположим, что нам нужно рассчитать уравнение (6.14). Запишем его в виде матрицы,

заменив г/А наг/0

(Л) и обозначив e~2nikn/N

= W. Для простейшего при-

мера положим

N =

22 = 4

То гда

 

 

/ Г ( 0 ) \

/ 1

1

1

1

 

 

У(1)

 

Wl

W2

W3

(6.16)

 

^(2)

 

W2

И 7 4

W6

 

 

 

\ У ( 3 )

 

Ц/3

Ц/в

Ц/9

 

 

 

 

 

 

 

В силу периодичности функции Whn =

где

т = 1, 2, . . ., будем иметь W m J V

и

т. д.

 

Следовательно,

e~2nil'n/N

для всех пк=

mN

= W° =

1, W™;v+i _

ц/

1

1

1

 

И 7 1

И 7 2

W3

W2

И 7 0

И 7

(6.17)

2

W7 3

И 7 2

И 7

1

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

(6.18)

Введем обозначение

-УЛО)

 

Уг(1)

(6.19)

г/1 (2)

 

V i ( 3 )

Выражение (6.17) после подстановки в него (6.19) с учетом пере­ мены мест второй и третьей строки принимает вид

(6.20)

188


Обратим

внимание на

следующую важную особенность. Если

в матрице

У (п) аргумент

п заменить на его бинарный эквивалент

 

 

/ У

(0)

 

/

У

(1)

 

\ У

(3),

а затем изменить порядок бит в каждом аргументе на обратный, т. е.

00 00

01 -> 10

1001

1111,

то матрица [ Y (п) | примет вид

Следовательно, будет иметь именно тот порядок членов, который был принят для удобного разложения квадратной матрицы Ъ в выражении (6.17). Эта закономерность справедлива для любого N.

Продолжим основное решение. Из выражения (6.19), выполняя перемножение матриц, имеем

M 0 ) = ! / o ( 0 ) 4 - n o ( 2 ) ,

 

J/i(l) =

Z/od) +

^ V o (3),

(6.23)

 

 

 

.'/, О)

//„(•'-)

И'7/ н (3).

 

Аналогично из (6.20) получим:

У (0) =

Ух (0) +

W°yx (1),

 

У (2) =

^ ( 0 ) ^ - ^ ( 1 ) ,

(6.24)

^(1) = 2/х ( 2 ) 4 ^ ( 3 ) ,

 

y(3) = y i ( 2 ) +

WVi(3).

 

189