Обратная фильтрация исходных сейсмических записей на основе кепстрального анализа Текст научной статьи по специальности «Математика»

УДК 550.834.05

ОБРАТНАЯ ФИЛЬТРАЦИЯ ИСХОДНЫХ СЕЙСМИЧЕСКИХ ЗАПИСЕЙ НА ОСНОВЕ КЕПСТРАЛЬНОГО АНАЛИЗА

© 2010 г. Ю.Д. Борисенко, Г.В. Калайдина

Кубанский государственный университет, Kuban State University,

ул. Ставропольская, 149, г. Краснодар, 355040, Stavropolskaya St., 149, Krasnodar, 355040,

[email protected] [email protected]

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

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

Suggested paper is devoted to the questions of inverse filtering (deconvolution) of .source seismic records on basis of the cepstral analysis. The canonical form of homomorphic system it is possible to represent as series connection of three systems: characteristic system, classic linear system and inverse characteristic system. The characteristic system transforms convolution vector space to additive vector space. This qualitatively new approach is based on the estimation of the seismic signal by means of the cepstral averaging and the design of shaping Wiener filter. The homomorphic deconvolution rises the time resolution of source seismic records.

Keywords: homomorphic systems, complex cepstrum, cepstral analysis, homomorphic deconvolution, shaping Wienerfilter.

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

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

- детерминистический подход М. Бакуса [1], основанный на применении операции линейной обратной фильтрации для подавления реверберации водного слоя;

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

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

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

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

Для принципа суперпозиции, применяемого к линейным системам L, требуется:

- свойство линейности между сигналами sj(t) и s2(t)

l[si (t)+со]=zh (t)]+l[s2 (t)]; (1)

- свойство однородности при умножении сигнала s(t) на скаляр c

L[c • s(t)]= c • L[s(t)], (2)

где t - время.

Рассмотрим законы комбинации входных сигналов □ и композиции А, которые служат для скалярного умножения входных сигналов, образующих векторное пространство G. Аналогично И и А - для выходных сигналов, образующих векторное пространство G'.

Тогда линейная система L преобразуется в систему Т, которая удовлетворяет обобщенным уравнениям (1), (2)

T[si (t) □ S2 (t)] = T[si (t)]0 T[s2 (t)] , (3)

T[c А s(t)] = c А T[s(t)]. (4)

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

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

Первая система А преобразует комбинацию сигналов □ в сумму и называется характеристической системой. Вторая система Ь - классическая линейная система. Третья система А"1 преобразует сумму сигналов в комбинацию 0 и называется обратной характеристической системой.

Если законы построения пространств тождественны (^=0, Д= А), то соответствующие им гомоморфные системы отличаются только линейной системой Ь. Системы А и А'1 определяются законами комбинации □ и композиции Д. Это основной результат, так как, если характеристические системы А и А"1 определены, то возвращаемся к вопросу о линейной фильтрации, в чем заключается гибкость данного класса.

□ + + + + □

s(t)eG А m 1 m А ï'(0 е G

Ч

m S(z) S(z) m

7 Ln 7-l

z

б

Z S\z) Exp S'( z) Z"1

Рис. 1. Блок-схемы систем: а - каноническая гомоморфная система; б - характеристическая система А; в - обратная характеристическая система А"1

Рассмотрим гомоморфную деконволюцию свертки двух сигналов. Поскольку мы имеем дело с зарегистрированными данными, то перейдем в дискретную область с шагом квантования, равным 1 Дt=1,

п=0,1,2,...), и представим входной сигнал s(n) виде свертки двух сигналов Sl(n) и S2(n)

s(n) = 2 (j)S2 (n - j) = Si (n) * s2 (n) .

(5)

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

Л^!(п),с1 *(п)'с2 ] = с • 5(п) + с • ^(п) . (6)

Здесь в1(п)*с1 = вг(п)*вг(п) ...*вг(п) - свертка 51 (п) с самим собой с раз.

Тогда характеристическая система А преобразует векторное пространство свертки в аддитивное векторное пространство. Рассмотрим прямое и обратное z-преобразования:

S(z) = 2 s(n)z—, S(z) = 2 S(n)z-n;

(7)

s(n) = — f 5(z)zn-ldz , S(n) = — f S(z)zn-ldz , (8) 2m с 2m с2

где С1 и С2 - замкнутые контуры в комплексной z-плоскости; i - мнимая единица. Так как входной сигнал s(n) представляет собой свертку двух сигналов sj(n) и s2(n), по теореме о свертке в комплексной z-плоскости имеем S(z) = S (z) • S2 (z), т.е. приходим к

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

S(z) = LnS(z) = ln | S(z) | +iArg[S(z)] (9)

Здесь возникает проблема многозначности комплексного логарифма LnS(z) [5], так как функция Arg\_S(z)]= arg[S(z)]+ 2mk, где -m<arg[S(z)]<m -главное значение аргумента комплексного числа S(z), определена с точностью до 2mk (k=0, +1, ±2, ...). Для того чтобы S(z) была аналитической на контуре С2, следует развернуть фазу z-преобразования Arg\S(z)] входного сигнала s(n) для получения непрерывной фазовой кривой в зависимости от круговой частоты а (в качестве контура интегрирования С2 берется

z = ae'a - окружность радиуса а в комплексной z-плоскости). Простые алгоритмы, основанные либо на исследовании дискретной фазы для определения скачков на 2 m либо на численном интегрировании фазовой производной, оказались в общем случае несостоятельными. Был предложен новый метод фазовой развертки путем адаптивного численного интегрирования [6] и реализован в данной работе.

Тогда характеристическую систему А можно представить, как показано на рис. 1б, где Z - прямое z-пре-образование; Ln - преобразование комплексного логарифма; Z"1 - обратное z-преобразование.

Выходные данные s(n) характеристической системы А будем называть комплексным кепстром входных данных s(n). Эта терминология была предложена в [7] для обнаружения отражений. Слово кепстр

(cepstrum) образовано парафразой слова «спектр» (spectrum). Термин «комплексный кепстр» подчеркивает, что использована фазовая информация в дополнение к амплитудной. Следует отметить, что комплексный кепстр представляет собой последовательность действительных чисел. Кепстральная область имеет размерность времени.

Затем производится линейное преобразование системой L: s(n) ^ s'(n).

Обратная характеристическая система A'1 преобразует аддитивное векторное пространство в векторное пространство свертки. Систему A"1 можно представить в виде, показанном на рис. 1в, где Exp - экспоненциальное преобразование.

При проведении обратной фильтрации на основе кепстрального анализа (гомоморфной деконволюции) необходимо сделать правильный выбор линейной системы L. Этот выбор тесно связан со свойствами комплексного кепстра:

1. Комплексный кепстр свертки двух или более сигналов равен сумме отдельных комплексных кеп-стров.

2. Комплексный кепстр минимально-фазовой последовательности равен нулю при n<0; а комплексный кепстр максимально-фазовой последовательности равен нулю при n>0.

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

4. Комплексный кепстр r(n) периодической последовательности импульсов r(n) также является периодической последовательностью импульсов с тем же периодом.

Заметим, что для получения комплексного кепстра сигнала s(n) надо сделать два дискретных преобразования Фурье. Здесь возникает проблема выбора шага квантования по времени At, чтобы избежать появления зеркальных частот. Для их подавления вводится экспоненциальное взвешивание входного сигнала, а также четырехкратное увеличение длительности входного сигнала за счет добавления нулей.

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

x(n) = w(n)* r(n) + na (n), (10)

где w(n) - сейсмический сигнал; r(n) - отражающий

ряд (совокупность коэффициентов отражения, которые определяются контрастностью акустических же-сткостей упругой среды); na(n) - аддитивный шум.

Если имеется набор из N сейсмических трасс (например, совокупность сейсмограмм ОСТ (общая средняя точка)), то с помощью кепстрального осреднения (10) можно получить оценку сейсмического сигнала w(n):

1 N

1 N

1 N

x(n) = — 2 Wj (n) +— 2 r j (n) +-2 naj (n), (11)

N j=i N j=i N j=i

где j (j = 1, 2, ..., N) - номер трассы.

+o>

j=—X)

n=—x

n=—x>

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

достаточно большом числе N второй и третий члены в правой части (11) стремятся к нулю. В таком случае X (п) = И (п) .

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

Для этого решается система уравнений Колмогорова-Винера, которая в матричной форме имеет вид [8]:

Bf = £ , (12)

где В - квадратная матрица, составленная из коэффициентов автокорреляции оценки сигнала И(п),

имеет тёплицеву форму; f - вектор-столбец, представляющий коэффициенты искомого фильтра; £ -

вектор-столбец, компонентами которого являются коэффициенты взаимной корреляции оценки сигнала И(п) и желаемого импульса ё(п).

Так как матрица В в системе линейных алгебраических уравнений (12) имеет тёплицеву форму, то система решается методом общей рекурсии Левинсо-на [9], что ведет к экономии машинной памяти и сокращению вычислительных затрат. Для повышения

р, мс/км 300 3 50 400

устойчивости численного решения системы (12) и для контроля за помехоустойчивостью добавляется «белый шум» порядка 1 % [10]. С помощью применения рассчитанного фильтра f к заданным трассам осуществляется гомоморфная деконволюция, приводящая к сжатию сейсмического сигнала.

На рис. 2а представлена исходная сейсмограмма ОСТ в 1аи-р области (р - лучевой параметр; =-рх -время двойной задержки плоской волны с лучевым параметром, где х - удаление) с интервалом дискретизации по времени Д!=2 мс.

Применив к 100 исходным сейсмограммам процедуру кепстрального осреднения, получим оценку сейсмического сигнала И(п) и ее спектра (рис. 3а). Оценка сигнала получена в интервале времен от 1300 -1600 мс, при кепстральном осреднении использовалось N = 400 трасс (из каждой сейсмограммы из трасс с номерами, принадлежащими отрезку [7; 16], случайным образом выбираются 4 трассы), коэффициент экспоненциального взвешивания был принят равным 0,965. Как показывают численные эксперименты, при кепстральном осреднении порядка нескольких сотен трасс значительно повышается статистическая устойчивость оценки сейсмического сигнала.

Максимум амплитудного спектра оценки сигнал равен Бтах=17,6 Гц, средняя частота Ртеап=50,9 Гц, медианная частота Рте<1 = 33,2 Гц (рис. 3а). Частота

Найквиста ^ равна 250 Гц ( г = —^ = —1— =

" 2Д 2 • 0,002

=250 Гц, где Д/=0,002 с - интервал дискретизации по времени).

р. мс/км

50 100 150 200 2 50 300 3 50 400

Рис. 2. Сейсмограммы ОСТ отражённых РР-волн в 1аи-р области: а - исходная сейсмограмма; б - сейсмограмма после гомоморфной деконволюции

б

Рис. 3. Оценки сейсмического сигнала и их амплитудные спектры: а - по исходным 1аи-р сейсмограммам; б - по 1аи-р сейсмограммам после гомоморфной деконволюции

В качестве желаемого импульса d(n) был выбран симметричный сигнал ограниченной длительности (пять дискретов в центральной части отличны от нуля). Полученный искомый фильтр был применен к 100 исходным сейсмограммам. На рис. 2б представлена сейсмограмма после гомоморфной деконволюции. На рис. 3б изображена оценка сигнала и её спектр после гомоморфной деконволюции. При этом характеристики амплитудного спектра сместились в область более высоких частот (Fmax=21,0 Гц, Fmean=55,4 Гц, Fmed=41,0 Гц).

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

Визуальное сравнение сейсмограмм ОСТ (рис. 2а и 2б) и оценок сейсмического сигнала и их амплитудных спектров (рис. 3 а и 3б) свидетельствует о повышении временной разрешенности исходных сейсмограмм после применения гомоморфной обратной фильтрации. Последнее обстоятельство обеспечивает получение более точных петрофизических характеристик упругой среды в результате решения обратной задачи сейсморазведки.

Литература

1. Backus M.M. Water reverberation - their nature and elimination // Geophysics. 1959. Vol. 24. Р. 233-261.

2. Robinson E.A. Predictive decomposition of seismic traces // Geophysics. 1957. Vol. 22. Р. 767-778.

3. Oppenheim A. V. Superposition in a class of non-linear systems // Tech. Rep. MIT. Res. Lab. оf Electr. 1965. № 432. Р. 62.

4. Oppenheim A., Schafer R. Digital processing. Engle-wood, 1975. 417 p.

5. Лаврентьев М.А., Шабат Б.В. Методы теории функций комплексного переменного : учеб. пособие. М., 1987. 688 с.

6. Tribolet J.M. A new phase unwrapping algorithm // IEEE Transactions on Acoustics, Speech and Signal Processing. 1977. Vol. 25, № 2. Р. 170-177.

7. Bogert B., Healy M., Tukey R. The frequency analysis of time series for echoes // Proc. Symp. On Time Series Analysis. N.Y., 1963. Ch. 15. Р. 209-243.

8. Применение цифровой обработки сигналов / под ред. Э. Оппенгейма : пер. с англ. М., 1980. 552 с.

9. Сильвиа М.Т., Робинсон Э.А. Обратная фильтрация геофизических временных рядов при разведке на нефть и газ : пер. с англ. М., 1983. 247 с.

10. Хаттон Л., Уэрдигтон М., Мейкин Дж. Обработка сейсмических данных. Теория и практика : пер. с англ. М., 1989. 214 с.

Поступила в редакцию

26 мая 2010 г