Бесплатный автореферат и диссертация по наукам о земле на тему
Устойчивые методы подавления регулярных и нерегулярных помех в сейсморазведке
ВАК РФ 25.00.10, Геофизика, геофизические методы поисков полезных ископаемых

Автореферат диссертации по теме "Устойчивые методы подавления регулярных и нерегулярных помех в сейсморазведке"

На правах рукописи

УДК 55(084.3):553.98:550.834

Денисов Михаил Сергеевич

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

Специальность 25.00.10 - геофизика, геофизические методы поисков полезных ископаемых

Автореферат диссертации на соискание ученой степени доктора физико-

математических наук

Москва, 2006

Работа выполнена в ООО Гготехсистем.

Официальные оппоненты: д.ф-м.н. О.К. Кондратьев д.ф-м.н. Ю.П. Ампилов д.т.н. Е.А. Козлов

Ведущая организация: Российский Государственный университет нефти и газа.

Защита состоится 18 мая 2006 г. в 15 часов на заседании Диссертационного совета Д212.121.07 при Московском Государственном геологоразведочном университете (МГГРУ) по адресу:

117997, г. Москва, ГСП-7, ул. Миклухо-Маклая, д.23, МГГРУ, ауд. 6-38

С диссертацией можно ознакомиться в библиотеке МГГРУ. Отзывы заверенные печатью учреждения, в двух экземплярах просим направлять пс адресу: 117997, г. Москва, ГСП-7, ул. Миклухо-Маклая, д.23, МГГРУ Ученому секретарю Диссертационного Совета.

Автореферат разослан 1 апреля 2006 г.

Ученый секретарь Диссертационного Совета кандидат технических наук, профессор

Г.Н. Богани

Общая характеристика работы

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

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

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

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

этапа адаптации. Этот вопрос подробно рассмотрен в работе и проиллюстрирован примерами обработки реальных данных.

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

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

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

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

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

Изложение способов решения задачи подавления регулярных помех производится в следующей последовательности. В Главе 2 рассматриваются алгоритмы прогнозирования кратных волн при помощи методов прямого продолжения волнового поля и устойчивые способы их адаптивного вычитания из исходных данных. Там же описаны кинематические фильтры, применение которых целесообразно после адаптивного вычитания. Способы построения устойчивых операторов продолжения волнового поля рассмотрены в Главах 3 и 4. Область их применения не ограничивается исключительно прогнозированием кратных волн, но использование соответствующих подходов оказывается полезным и при решении других задач, а разработанными алгоритмами подавления артефактов любого преобразования сейсмических данных, включающего пространственное суммирование трасс?, полезно оснастить не только операторы продолжения волнового поля (как это сделано в Главах 3 и 4), но и кинематические фильтры, описанные в Главе 2, причем соответствующие подходы могут применяться при фильтрации как в пространственно-частотной, так и в пространственно-временной области.

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

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

Актуальность работы

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

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

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

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

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

Цель работы

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

Задачи исследований

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

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

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

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

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

Получить алгоритм устойчивой нестационарной кинематической фильтрации с компенсацией возможных искажений формы полезных отражений в пространственно-временной 1-х и в пространственно-частотной / -х областях. Разработать и применить схему подавления артефактов преобразования.

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

поля.

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

Разработать соответствующее программное обеспечение и показать их эффективность методов на практических примерах.

Научная новизна и личный вклад

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

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

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

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

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

5. Разработаны способы коррекции динамических искажений полезных волн, возникающих при нестационарной по временной и пространственной координатам кинематической фильтрации, применяемой в 1-х области.

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

7. Предложен алгоритм выделения областей накапливания конструктивных сумм при пространственном суммировании сейсмических трасс. Это позволяет получать результаты, свободные от аляйсинг- помех.

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

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

10. Оценивание параметра частотно-зависимого поглощения производится оптимизационным способом, причем функционалом является дисперсия ошибки предсказания.

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

Защищаемые положения В диссертации защищаются следующие основные научные положения:

1. Предложены статистически устойчивые методы учета нестационарности сейсмических волновых полей для их обработки с целью построения сейсмо-геологической модели среды.

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

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

4. Предложен метод ослабления шумов пространственной дискретизации (аляйсинг- помех) при прямом и обратном продолжении волнового поля.

5. Разработан алгоритм помехоустойчивого пространственного суммирования сейсмических трасс при редкой сети наблюдений.

6. Предложен устойчивый метод оценивания амплитудного спектра сигнала и коэффициента поглощения по малым выборкам.

Апробация

По теме диссертации опубликовано 24 научные работы. Результаты исследований были изложены в докладах на нескольких ежегодных международных конференциях SEG, EAGE, а также на международных конференциях Istanbul-97 и Геомодель-2005. Большое количество результатов обработки реальных данных представлено в отчетах ООО Геотехсистем. (ген. директор А.А. Пудовкин, научный директор В.М. Глоговский). Все описанные алгоритмы обработки сейсмических материалов реализованы в системе VELINK (совместный продукт компаний Геотехсистем и CGG), ориентированной на обработку сейсмических материалов 2D, и в системе PRIME — обработка 3D, используемыми как отечественными, так и зарубежными геофизиками. В самой диссертационной работе приводятся многочисленные иллюстрации, подтверждающие эффективность применения предложенных процедур при обработке реальных данных, полученных, в первую

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

Кроме того, вопросы реализации алгоритмов и методологические аспекты применения процедур обработки приведены в работе Денисов М.С., Силаенков O.A., 2003, Пример использования процедур прямого и обращенного продолжения волнового поля в задаче построения глубинного разреза. Геофизика, N3. Там же даны разнообразные примеры обработки реальных данных, не приведенные в тексте диссертации.

Автор глубоко признателен проф. A.A. Никитину, фактически ставшему инициатором написания этой диссертации. Успешному проведению исследований способствовали консультации и поддержка В.М. Глоговского на протяжении более пятнадцати лет совместной работы. Значительная часть алгоритмических разработок была выполнена совместно с Д.Б. Финиковым, которого автор сердечно благодарит за многолетнее сотрудничество и который оказал решающее влияние на формирование научных взглядов соискателя. Автор выражает искреннюю признательность своим коллегам Ю.А. Харитонову, Д.М. Оберемченко, Е.А. Курину, А.Е. Фирсову, С.Л. Лангману, O.A. Силаенкову (ООО Геотехсистем) за предоставление соискателю возможности практической реализации научных разработок и за оказанную помощь при обработке реальных данных, анализе и интерпретации результатов. Автор благодарен Д. Локштанову (Norsk Hydro) за плодотворные обсуждения и обеспечение финансовой поддержки исследований.

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

Объем и структура диссертации

Диссертация состоит из пяти глав, общий объем - 198 страниц, включая 63 рисунка. Список литературы - 223 наименования.

Содержание работы

Глава 1. Введение

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

Глава 2. Устойчивые методы подавления кратных волн в сейсморазведке

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

С другой стороны, существует и иной подход, в рамках которого утверждается; что в сейсморазведке, основанной на методе отраженных волн, следует использовать все зарегистрированные сигналы, и информация о глубинном строении земли, извлеченная из кратных волн, способна успешно дополнять соответствующие построения, производимые только по полю однократных отражений (Z Jiang, J. Yu, G. Schuster, В. Hornby, 2005). Действительно, комбинирование информации, содержащейся в поле кратных и однократных волн, может обеспечить получение более подробных сведений о свойствах исследуемой среды. Подавление кратных отражений целесообразно производить после такого анализа. В диссертационной работе основное внимание сосредоточено на анализе свойств наиболее перспективных подходов к подавлению кратных волн и разработке новых устойчивых методов.

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

которые явно, а чаще всего неявно, заложены в каждом из них), а также их сильные и слабые стороны.

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

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

Прогнозирование волнового поля Прогнозирование кратных волн всегда соответствует преобразованию вида прямого продолжения волнового поля, даже если речь идет об одномерной задаче сейсморазведки и о предсказывающей деконволюции, где получение модели кратных волн сводится к простому сдвигу зарегистрированной трассы. Согласно теории, (Г.И. Петрашень, С. А. Нахамкин, 1973, N. Bleistein, J.K. Cohen, J.W. Stockwell, 2000) решение задачи определяется интегральным преобразованием

L(a, b, со) - J jP(a. q.coсо)dxdy, (1)

*y

где предполагается, что волновое поле, наблюдаемое на свободной поверхности z = 0 в точках с координатами q-{x,y), вызвано источником, расположенным в точке а = (ха.уа). Тогда подборка трасс p(a.q.t) будет соответствовать сейсмограмме

общего пункта взрыва (ОПВ(а)), a P(a,q.t>) - ее преобразованию Фурье по t. Ь -некоторая фиксированная координата пункта регистрации, для которого требуется получить модель кратных волн £(а,6,ю), G(q,b, а>) - функция Грина.

Оператор —G(c?,6,6i), который в (1) является весовой функцией, в нашем случае

oz

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

Так как в рамках лучевой асимптотики

dz K{q, b) dz dz v

L(a, b, to) = F(co) J \P{a, q, со) cos QP(q, b, a)dxdy, (2)

ХУ

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

8 - угол выхода луча, F{ со) = у— - одноканальный фильтр.

v

Преимущество преобразования (2) заключается в том, что оно не требует оценки модели среды. Приведенные выше рассуждения можно повторить для всех зарегистрированных отраженных волн (см., например, A.J. Berkhout, D.J. Verschuur, 1997), что приводит к выводу о принципиальной возможности построения поля всех кратных волн, связанных с переотражением от свободной поверхности, без использования глубинно- скоростной модели. Само выражение (2) соответствует взвешенному суммированию взаимных сверток трасс сейсмограмм ОПВ(д) и 01111(¿). Однако для учета множителя cos 9 потребуется специальное преобразование. Его удобно производить в области разложения по базису плоских волн (преобразование

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

является более строгим аналогом весовой функции — С(д,Ь,оз). То есть из (2)

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

Ца,Ь,<й) = Р(ю)ЦР(а,д,а)Р{д,Ь,а)<1хс1у. (3)

ху

Одномерными аналогами выражений (1) и (3) соответственно являются сдвиг трассы, осуществляемый в предсказывающей деконволюции, и метод ретрокорреляции (см., например, Мушин И.А., 1983). На практике выбор между подходами (1) и (3) базируется на специфике зарегистрированных данных и особенностях структуры исследуемой среды.

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

Рис.1 иллюстрирует переход от асимптотически точного алгоритма прогнозирования (1), основанного на расчете функции Грина, к выражению (3). Сигнал, излучаемый источником а. регистрируется как минимум дважды: в первый раз приемником Щ как однократное отражение, во Рис.1 второй раз - приемником Ь как кратная волна. Выражение

(1) преобразует однократное отражение в кратное при помощи расчета динамического годографа волны, распространяющейся из д в Ь, а (3) как бы заимствует этот оператор непосредственно из зарегистрированных сейсмограмм. Назовем точку с координатой д точкой выхода промежуточного луча.

Адаптивное вычитание модели кратных волн из исходных сейсмограмм

Оба подхода к прогнозированию поля кратных волн имеют и общие для них недостатки, приводящие к неизбежным динамическим искажениям получаемой модели. Во-первых, каждое из преобразований получает волновое поле с точностью до сверточной компоненты: (1) - не учитывает характеристику отражения от выбранного кратнообразующего горизонта, которая неизвестна; (3) - из-за взаимных сверток трасс форма импульса претерпевает автосвертку, поэтому следует производить коррекцию за форму сигнала, которая неизвестна (для повышения надежности вычислений полезно применить коррекцию сигнала при помощи метода, изложенного в Главе 5, до прогнозирования в соответствии с (3)). Таким образом, необходим этап адаптации модели к зарегистрированной трассе перед вычитанием, при этом упомянутые искажения описываются стационарным по времени формирующим фильтром.

Иной вид динамических искажений вызван тем, что при прогнозировании в соответствии как с (1), так и с (3), некоторые типы волн моделируются дважды. В результате разные сигналы кратных волн будут отличаться от реально зарегистрированных отражений разными множителями. При использовании (1) следует построить и использовать на этапе адаптации дополнительную модель волн старшего порядка кратности /*(а,Ь,/). Если применяется алгоритм прогнозирования (3), то необходим рекурсивный расчет модели кратных волн, при котором волновое поле, полученное на предыдущей итерации, подается на вход последующей итерации.

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

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

колебаний, так и источника. На практике это требование, как правило, не удовлетворяется: данные, полученные в результате наземной сейсморазведки, имеют крайне нерегулярную геометрию со значительным количеством пропусков как пунктов взрыва, так и пунктов приема, а данные морской сейсморазведки фактически являются набором независимо полученных линейных профилей. Такая специфика сводит на нет почти все потенциальные преимущества площадных наблюдений (R.G. Borselen, М.А. Schonewille, R.F. Hegge, 2005, R.G. Van Borselen, D.J. Verschuur, 2003). Следствием этого является невозможность прогнозирования определенного набора кратных волн, а также значительная зашумленность получаемых моделей, обусловленная аляйсинг- эффектом, особенно сильно проявляющимся как раз в случае нерегулярной и редкой сети наблюдений.

В этой ситуации стандартным подходом к ослаблению артефактов преобразования вида пространственного суммирования является предварительная интерполяция исходных сейсмограмм, но при 3D сейсморазведке очень существенно ограничение на объем данных, и всегда приходится искать компромисс между плотностью интерполируемых данных п точностью вычислений. Кроме того, несмотря на обилие алгоритмов, в силу существенной неопределенности самой задачи интерполяции, особенно в случае разведки сложнопостроснных сред, нет вполне удовлетворительных и универсальных способов ее решения. В настоящей работе предложено использовать альтернативный подход к проблеме шумов дискретизации, не прибегающий к интерполяции, то есть позволяющий вести обработку без увеличения объема исходных данных. В силу универсальности его приложений в сейсморазведке, метод рассмотрен в отдельных главах - Глава 3 н Глава 4. Глава 3 посвящена разработке устойчивой модификации алгоритма продолжения волнового поля с привлечением априорной информации о модели среды, то есть прогнозирования кратных волн в соответствии с выражением (1), так как расчет оператора соответствующего преобразования основан на привлечении явного выражения для его линии суммирования. Особенности реализации вычислений в соответствии с (2) обсуждаются в Главе 4.

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

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

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

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

где Т - длина трассы, а и'ДО - ступенчатая функция (именно она и определяет скачкообразное поведение нестационарности, которое допускается при адаптации), §, й - матрицы, содержащие набор фильтров £,(?) и §,(()■ Использование такой параметризации, не накладывающей ограничений на глобальное поведение формирующего фильтра и основанной на независимом его оценивании в пределах каждого окна, как раз и приводит к опасности ослабления всех, в том числе однократных, волн.

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

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

полиномы

«,-(/)-/' (5)

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

Для более уверенной адаптации, особенно на фоне шумов, полезно привлекать и соседние трассы, например, 1(а,Ь + ДхЛ), 1(а,Ь- Лх,() и т.д. Тем самым размерность фильтра увеличивается, но структура решаемой системы уравнений не изменяется. Тогда задача адаптации модели кратных волн к исходным записям формулируется следующим образом: требуется найти ¿¿(/,0 и ¿'к (/,/)» минимизирующие функционал т( мы, ' \2

2,8;' 04 А=-А/ /=0 )

где /,(а,Ь,г) = /(а,6,г)м>Дг)> 1Ца,Ь,1) = /'(а, 6,/)■№,(/), а выражение записано с учетом дополнительно посчитанной модели кратных волн.

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

Т(а. Ъ + кАх, 0=1 (;, 0 * I, (а, Ь + Ш,1) + Щ1,1)* //(а, Ь + кАх,0). ;=о

Затем оценим коэффициенты с* такие, что

ск = arg min J(c),

т( м _ Y2

J(c)= fi pia,b,t)- YckUa.b + k&xj) \ dt.

Процедуру повторяем несколько раз, подменяя в последнем выражении p{a,b,t) на результат вычитания на предыдущей итерации, пока не достигнем минимума J{c), пользуясь обычным критерием сходимости. Так как исходный функционал выпуклый, нетрудно доказать сходимость описанной процедуры. Поскольку размерность предложенной итерационной схемы существенно ниже, легче справиться с проблемами обусловленности.

Очень кратко проиллюстрируем свойства алгоритма прогнозирования (3) и многооконной адаптации на реальных данных, полученных при морской сейсморазведке. На Рис.2 показан фрагмент исходной сейсмограммы ОГТ (слева) а также результат адаптивного вычитания модели кратных волн (справа), построенной без привлечения информации о структуре среды.

Рис. 2. Вычитание модели кратных волн. Слева - фрагмент исходной сейсмограммы ОГТ, справа - результат вычитания при помощи устойчивого алгоритма адаптации.

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

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

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

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

На Рис.3 показан пример использования процедуры прогнозирования кратных волн без привлечения априорной информации о модели среды с последующим адаптивным вычитанием по данным, полученным при морской сейсморазведке. Слева изображен - " 20 - "

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

2.3 Одношаговые методы, кинематическая фильтрация

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

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

Устойчивый t-x фильтр

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

г(лг,0 М

p(,x,t)= J da(x,t)da, где da(x,t) = \p{x-yj + a.y)dy, (6)

f(x,l) -M

где p{x,t) - исходное волновое поле, da(x,t) - локальная направленная сумма, у е (~М,М) пространственная база фильтра, г(х,1) и r(x,t) - границы перебора направлений суммирования (веер), р(х,1) - результат фильтрации.

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

-иЛ и а

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

p(x,t) = P(x,l)*f(t). (7)

Пусть наряду с годографами полезных волн получена оценка кинематики регулярного шума, и при этом следует разработать алгоритм, подавляющий такой шум. Задача решается при помощи модификации преобразования (6), подразумевающей соответствующий выбор границ веера r(x,t) и r(x,t), при этом, как правило, они уже не стационарны по времени. Фильтрация вносит такие искажения в форму полезного сигнала, которые не могут быть скомпенсированы при помощи описанного выше одноканального алгоритма, и требуется иной, более сложный подход. Кроме того, в волновом поле после фильтрации могут появляться области, где вообще отсутствуют амплитуды: такое происходит, когда годографы сигнала и помехи имеют близкие производные, а веер на выделение сигнала не шире веера режекцни.

Коррекция за кривизну годографа Пусть полезный сигнал имеет годограф Г(х). Тогда действие оператора кинематической фильтрации (6) на такой сигнал может быть представлено в виде

м

Da(x,v) = S(ш)

Р(х,со) = ^Da(x,(o)da, где S(cl>) - спектр импульса.

Для простоты изложения здесь и ниже формулы преобразования выписаны в спектральной области. Интегрирование может быть произведено с использованием высокочастотной асимптотики, то есть по методу стационарной фазы {Математическая энциклопедия, 1984). Тогда, опуская все промежуточные выкладки, получим, что амплитуда результата фильтрации будет зависеть от удаления х, причем

эта зависимость пропорциональна величине ) ' Если сигнал имеет

гиперболический годограф Г(х) = J/q + ~ , то -/Г"(х) = -— -, что и будет

являться коэффициентом, корректирующим амплитуду. результата . фильтрации за кривизну годографа. К получаемым коэффициентам необходимо применить нормировку. Например, можно принять за единицу коэффициент, рассчитанный для нулевого удаления х = 0 на времени t — 1000ms .

Коррекция методом пространственного сглаживания

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

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

5( со) = 1,сое (со,,ш2), S(a>) = O.cog (И|,со2),

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

~ я ям

Р{х = 0,со) = }оа(х = 0,со)<&х = | ¡е^е/уаа, О 0-М

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

ю0 а

а амплитуда импульса после кинематической фильтрации, равна со2 ~

1(0|,=0 = / Р(х — 0,оэ)<Л». Так как нас интересует искажение амплитуды сигнала «1

относительно сигнала на выходе алгоритма суммирования с полным раствором веера (обозначим ее через 5(О|,_0 ), то необходимо применить нормировку

(8)

...............'

Величина может быть рассчитана теоретически по той же самой формуле при

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

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

Ат*\ .\=х

А= 1Р(.Г,?0 + Г(Д:)),

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

л- А

-тта\ 2>00

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

. p(x,t0 + Г(дг)> = А(1 - w(x)) + р{х,/0 + Г(.х))и-(х) для всех значений пространственной координаты х. Из этого выражения видно, что если полезный сигнал не был искажен кинематическим фильтром, то есть w(i) = 1, то преобразование не изменит его амплитуду. А если приходится иметь дело с областью, где произошло полное подавление энергии волнового поля, то есть w(x) = 0, то преобразование восстановит в этих точках амплитуду, равную средневзвешенному значению А .

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

P(.x,t0 + Г(х)) = + r(x))vv(x), где н(л) = 1/ w(x), либо с применением регуляризации, которая может достигаться

использованием гладких функций вида vv(x) =-W --. Здесь у - параметр

w"+4*)+Y

регуляризации.

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

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

Устойчивый / - х фильтр

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

М _

Р(х,со)= | £>аО,ш)с/а, где Оа(х,а)= \Р{х-у,<а)е^у(1у, и Р{х,т) = Г(а)Р(х,<о).

Пх) -М

F(<й) - фильтр, компенсирующий искажения формы импульса, вызванные суммированием.

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

\Оа(х,оз) = 4йа(х,со)],' где /$[•] обозначает применение алгоритма ослабления аляйсинг- эффекта, а пересчитанная таким образом локальная сумма далее подставляется в тот же интеграл. Этот алгоритм основан на принципе масштабирования амплитудных спектров локальных направленных сумм. В основе подхода лежит гипотеза, согласно которой при направленном суммировании трасс исходной сейсмограммы низкочастотная часть сигнала накапливается лучше, чем высокочастотная. Иначе говоря, отношение энергии спектров результата направленного суммирования к исходным данным, взятое для некоторого высокочастотного интервала, не может превосходить соответствующее отношение в назначенном эталонном низкочастотном интервале. Противная ситуация, встречаемая при локальных разрастаниях энергии спектра направленной суммы, всегда сигнализирует о наличии аляйсинг- помех на этих частотах. Алгоритм редактирует (масштабирует) амплитудный спектр локальных направленных сумм так, что условие, накладываемое на поведение амплитудного спектра, будет соблюдаться (более подробное изложение метода дано в Главе 3).

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

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

Глава 3. Построение устойчивого оператора продолжения волнового поля с использованием локального направленного суммирования

Этот и следующий разделы посвящены разработке новых методов построения устойчивых операторов прямого и обратного продолжения волнового поля. Задача пересчета данных, зарегистрированных в точках дневной поверхности, на другой уровень является одной из наиболее актуальных и, в то же время, хорошо изученных в сейсморазведке (Е.А. Козлов, 1986, Э.А. Робинсон, 1988, A.J. Berkhout, 1981). В последние годы проблема точности глубинных построений приобретает особую значимость в связи с появившейся возможностью использования новых мощных вычислительных средств для реализации времяемких алгоритмов обработки данных. В Главе 3 разработана устойчивая модификация метода продолжения волнового поля с привлечением априорной информации о модели среды, то есть прогнозирования кратных волн в соответствии с выражением (1).

Одним из способов улучшения характеристик операторов продолжения волнового поля является реализация соответствующих преобразований с использованием локального суммирования, что позволяет улучшать условия накапливания сигналов, принадлежащих синфазностям небольшой кривизны, а также повышает устойчивость алгоритма к нерегулярному шуму {<P.,\f. Гольциан, 1964, U. AI bertin, С. Stork, D. ringst, W. Chang. R. Fletcher, P. Kitchenside, 2004, N.R. Hill, R.T.

Langem. Т. Nemelh, M. Zhao, K.P. Bube, 2002). В работе изложен способ расчета оператора с использованием локальных направленных суммотрасс. Подход описан в самом общем виде, и может быть применен в любых задачах продолжения волновых полей. Его основные преимущества заключаются в лучшем (близком к синфазному) суммированию полезного сигнала, а также в возможности применения эффективного алгоритма ослабления артефактов преобразования. При разработке специальных приемов подавления аляйсинг- шумов был применен нетрадиционный в этой задаче способ разделения спектров сигналов и помех.

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

= р х=х*{р)

дх

то есть касательной в этой точке. В сумме участвуют трассы, находящиеся в пределах базы [х*(р)~М,х*(р) +М]. Таким образом, локальная сумма может быть выражена как

м

С1(х\р),0= 5>(х*(р)+ (9)

х=-М

здесь обозначение соответствует исходному волновому полю. Сигнал,

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

преобразование, осуществляемое на втором шаге:

"(р)

w(b,t) = d{x\p),t)*f^w(t)dp (10)

-Vi

где - и д2 соответствуют минимальной и максимальной производной Лйнии суммирования оператора экстраполяции.

Очевидно, что при М —► 0 алгоритм (9), (10) переходит в традиционную схему типа (1), а при Мсо в оператор, основанный на преобразовании Радона. Таким образом, способ занимает промежуточное положение между двумя хорошо изученными и применяемыми на практике алгоритмами.

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

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

<Г(дг*0Д/)= |>(* *(р)-х,г+-А-(д))(-1)\ (И)

х=-М

результат которого принимается за зависящую от направления оценку шума. Разность квадратов амплитудных спектров (9) и (11), jGp(©)j , соответствует несмещенной

оценке энергетического спектра сигнала. Тогда расчет фильтра осуществляется в соответствии с описанным ниже алгоритмом.

Выбрав некоторый низкочастотный интервал (ш, ,со2 ) в области, не

подверженной аляйсинг- эффекту, по характеристике |С(ш)| получим оценку энергии

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

,2

] Лю

= -.

I |{/(<в)|2С&) • ш,

|СР(»)|2

Если на какой-либо частоте, большей еъ, соотношение Кр(и>) =!--- превосходит

■ Р(<й)?

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

Fp(<a) = 1, юе[0,о2], где Лр(а) - амплитудный спектр суммотрассы (9).

Рис.4(а) Рис.4(б) Рис.4(в)

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

Разумеется, более «чистая» модель кратных волн предпочтительна для последующего адаптивного вычитания.

Глава 4. Способ подавления шумов дискретизации при суммировании сейсмических

трасс

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

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

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

Будем рассматривать вычисления, производимые в соответствии с выражением (3). Тогда получим

¿(я A to) = F(03)2 Р(а, а + /Дх, ш)Р(а + г'Дх, ю). (12)

i

Диапазон изменения индекса / определяет базу суммирования (апертуру). Дискретность данных по пространственной координате неизбежно приведет к тому, что в поле L(a,b, со) наряду с сигналом накопится и аляйсинг- шум. Основная идея алгоритма помехоустойчивого преобразования заключается в том, что в глобальном варианте суммирования (12),невозможно выделить те области, которые дадут вклад в формирование сигнала от тех, которые сформируют лишь аляйсинг- шум. Но если посчитать L(a,b,a) не непосредственно, а через предварительный расчет локальных сумм в пределах небольшой скользящей базы, а затём произвести их анализ, обработку и суммирование, то такое выделение станет возможным. Если говорить о пространственном суммировании (12) в кинематических терминах, то следует отделить область накопления конструктивной суммы (или, что то же самое, область касания годографа волны на сейсмограмме ОПВ Р(а,х, со) годографом волны на сейсмограмме ОПП Р(х,Ь,а)) от области, где условия касания не выполнены.

Введем понятие взвешенной локальной суммы, полученной на базе (2М +1) канатов, с весами и(/)

м _

D(a,b,iAx, со)= 2^P(a,a + (i + I)Ax,n)P(a + (i + l)Ax,b,<i>)u(0 (13)

и отметим, что при и{1) = 1 глобальная сумма (12) может быть получена как непосредственным суммированием в пределах апертуры, так и сложением локальных сумм:

Ца, Ь, со) = F(o)Y.D(a,b,iAx,a).

I

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

Рис. 5(а)

Рис. 5(6)

Рис. 5(в)

На Рис. 5(а) показана модельная сейсмограмма, содержащая три волны. Результат расчета локальных сумм с!(а, Ь.г'Лх,/) представлен на Рис. 5(6), горизонтальная ось соответствует координатам ¡Ах центров баз суммирования. Буквами обозначены Л -пример области, дающей вклад в конструктивное синфазное суммирование, В - пример области, дающей вклад в артефакты преобразования вида аляйсинг- шума.

Алгоритм обнаружения конструктивных областей применяется во временной области, то есть к последовательностям с1{а,Ь,¡Ах,(). Анализируются локальные суммы и производится их взвешивание: с1(а,Ь,1Ах,1) = м>{1Ах,1)с1(а,Ь,1Ах,1), причем уфАг,;) = 1- если в локальной сумме накоплен полезный сигнал, н^/Дх.О = 0 - если сумма содержит только аляйсинг- шум. Преобразованные таким образом результаты суммируются, а на выходе получается искомая трасса 1(а.Ь,1) с ослабленными артефактами.

Расчет весов производится в скользящей базе путем сопоставления энергии Е * (1,1 Ах) и Е~(иДх) результатов локального суммирования, полученных по формуле (13) соответственно при к(/) = 1 и аналогичного знакопеременного суммирования при «(/) = (-1)': А(1,х) = Е~(/,х)/Е + 0,х). Понятно, что случай Е~(1,х)/Е+(1,х)« 1 соответствует синфазному суммированию, а Е~(1,х)/Е + (1,х) = 1 - аляйсинг~эффекту).

Рис. 6. Выбор критерия для анализа областей накапливания сигнала и аляйсинг- шума. Обозначения Ь+ и Ь~ использованы- соответственно для знакопостоянного и знакопеременного суммирования.

Рис. 6 иллюстрирует процесс получения локальных сумм для различных случаев (области типа А и В). В левом верхнем углу изображен фрагмент из пяти взаимных сверток, которые будут суммироваться синфазно, здесь 2М + 1 = 5. Хотя все сигналы одинаковы, амплитуда отсчетов последовательности Ь~ ненулевая, так как количество каналов нечетное (в принципе, алгоритм допускает использование и несимметричных четных баз для получения локальных сумм). Если энергия измеряется во всем временном интервале, изображенном на рисунке, то для случая синфазного суммирования получим А(1,д)<< 1. Трассы, показанные во втором ряду, соответствуют ситуации, когда локальная сумма будет содержать лишь аляйсинг- шум. Легко видеть, что в этом случае энергии результатов суммирования с единичными весами и знакопеременными весами будут близки: А(1, д) « 1.

Тогда можно выбрать, например, следующую весовую функцию

*>{*) =---г>

1 + (А/р)"

где п - параметр, отвечающий за крутизну спада от интервала м'(Л) = 1 к интервалу 1у(/4) = 0 , р - параметр (пороговое значение).

Результат применения алгоритма взвешивания локальных сумм показан на Рис.5(в).

Рис.7 Результаты прогнозирования поля кратных волн по модельным данным.

На Рис.7 и Рис.8 сопоставляются результаты моделирования кратных волн но алгоритму (3) слева — без ослабления аляйсинг- шумов, справа - с ослаблением. Разумеется, более «чистая» модель кратных волн предпочтительна для последующего адаптивного вычитания.

Рис.8 Результаты прогнозирования поля кратных волн по реальным данным

Глава 5. Параметрический оптимизационный способ оценивания амплитудного спектра сейсмического сигнала и декремента поглощения Оценивание амплитудного спектра импульса - одна из задач, решаемых в различных областях применения цифровой обработки сигналов. Это и одна из самых старых проблем сейсморазведки, решению которой посвящено огромное количество работ. Многие алгоритмы, созданные специально для геофизики, получили широкое распространение в совершенно иных приложениях (например, метод максимальной энтропии Дж. Бурга - J. P. Btirg, 1975), но оказались недостаточно эффективны в сейсморазведке. Причины этого связаны со спецификой сейсмических записей. В диссертационной работе предложен алгоритм спектрального анализа и естественно вытекающий из него способ коррекции амплитудного спектра сигнала, который опирается на несколько известных и хорошо зарекомендовавших себя на практике идей и позволяет учесть важные особенности сейсмических.записей. При этом в структуру алгоритма легко включается фактор, отвечающий за нестационарность формы сигнала, обусловленную частотно-зависимым поглощением.

Параметризация амплитудного спектра импульса и алгоритм амплитудной

деконволюции

Обсуждается обычная постановка задачи, когда сейсмическая трасса p{t) описывается сверточпой моделью

Р(') = Х'>+"(') = *(')* 4(0+«(0,

где / - индекс дискретного времени, s(t) - сигнал, £(/) - импульсная трасса, n(t) -помеха, некоррелированная с y(t). Предположим, что на интервале о е (W|,c£>2)

LS(o)| - а0 exp[¿a,M/,(co)], у,(со) = cos[/n^°

,=i co2-cd,

Понятно, что при сколь угодно больших N этой формулой можно описать практически

любой амплитудный спектр. Ограничение числа N - количества базисных функций в

разложении'логарифма амплитудного спектра на интервале (а>|,ш2), предполагает

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

область решений. Другое важное ограничение на поведение возможных амплитудных

спектров сейсмического импульса описывается формулой

N ■}

Х(а,02<р. 1=1

i i N

Пусть /(/, a) - фильтр, амплитудный спектр которого |/г(со, a)| = exp[£a,v|;l(<o)]. Будем

искать вектор параметров a : {a |, a2,..., aN}, минимизирующий функционал

Ъ

S(/('.5)V(0)2 (И)

Л' ~ 2

при ограничении X(a,i) < р. Здесь (Г,,Г2) - окно настройки. Выбор фазового спектра 1=1

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

Амплитудная деконволюция с учетом частотно-зависимого поглощения Частотно-зависимое поглощение - один из основных факторов, приводящих к изменению формы регистрируемого сейсмического сигнала отраженной волны во времени. Большинство алгоритмов обратной фильтрации опираются на гипотезу стационарности формы импульса в окне настройки оператора. Поэтому в граф обработки требуется включать процедуру коррекции частотно-зависимого поглощения. Но применению любого способа коррекции должно предшествовать оценивание коэффициента поглощения, что представляет собой самостоятельную и весьма нетривиальную задачу (М.Б. Рапопорт, 1969, А.Г. Авербух, 1982, Э. Данквардт, У. Патцер, 1982, O.K. Кондратьев, 1986, Ю.П. Алтичов, 1992).

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

нестационарной сейсмической трассы использовать выборки значительно большего объема.

Аналогично (14), в многооконном режиме ; решение задачи сводится к минимизации функционала

г! к г*

где Г* и - времена начала и конца к-го окна настройки. Если стоит и задача

оценивания поглощения, то базис у, (со) дополняется функцией £>(&) = ю - Ю| * с

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

г*

•/р(а,0) = ££(/Н'^е)**(О)2. (15)

* г*

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

N V

(со, а, 9)| =

ехр

9*е(а)+Хч>,(<0) а,-<=1

1,о> еЦ,«^)

,ае(а1,в2)

" ~ 2

при прежнем ограничении 2(а/') -Р и с учетом условия > . Удобно

/=1

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

И2 -у ( ы Л

к ш, V 1 = 1 )

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

На Рис. 9 показан фрагмент разреза. Выбран небольшой интервал по латерали, чтобы можно было рассмотреть особенности волнового поля, но протяженный по времени, чтобы можно было визуально оценить эффект выравнивания частотного состава записи вследствие амплитудной деконволюции с коррекцией поглощения. Этот пример похож на результат работы обычной многооконной деконволюции, но получен оцениванием шести параметров амплитудного спектра в диапазоне (7,60)Гц при частоте Найквиста 250Гц и двух параметров поглощения на интервалах (1440,2050)гпз и (2050,2830)гп5 относительно интервала (580,1440)тз, при этом использовалось шесть коэффициентов разложения при параметризации амплитудного спектра. Понятно, что восемь параметров по таким временным интервалам оцениваются сравнительно надежно. В данном случае были получены значения 02 = 0.85 при декременте поглощения равном 0.015 и 93 = 5.43 при декременте 0.096.

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

минимально-фазовым операторами коррекции поглощения обычно несущественен, так как последний имеет практически линейную характеристику в рабочем диапазоне частот (Д.Б. Фиников, 1989, L.J. Gelius, 1987). Другими словами, если импульс имеет невыровненный амплитудный спектр, то эффект частотно-зависимого поглощения может несущественно отразиться на его длительности, а коррекция поглощения не приведет к заметному изменению его формы. Кроме того, если импульс имеет сложный фазовый спектр на тех частотах, которые подавлены поглощением или ослаблены в силу других факторов, то после коррекции амплитудного спектра и поглощения разрешенность трассы может даже ухудшиться. Таким образом, задача повышения разрешенности решается лишь совместным оцениванием и последующей коррекцией как амплитудного спектра и параметров поглощения, так и фазовой характеристики сигнала. Результат коррекции фазового спектра показан на Рис.9(в).

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

работах

1. Денисов М.С., 1995, О рекурсивном НЧ фильтре для обработки гравимагнитных данных, Физика Земли, N11.

2. Денисов М.С. 1998, Новые возможности спектрального оценивания по методу максимальной энтропии, Физика Земли, N3.

3. Денисов М.С., 2005, Динамические особенности современных алгоритмов прогнозирования и вычитания кратных волн, Технологии сейсморазведки, N3.

4. Денисов М.С., 2006, Анализ метода прогнозирования кратных волн без знания модели среды с позиций теории продолжения волновых полей, Геофизика, N1.

5. Денисов М.С., 1995, Как интерпретировать алгоритм предсказывающей деконволюции, Ггофизика, N3.

6. Денисов М.С., Бусыгин Г.В., 1995, Повышение временной разрешенное™ сейсмических записей методом робастной деконволюции, Геофизика, N5.

7. Денисов М.С., Фиников Д.Б., 1999, Критерий минимума дисперсии ошибки предсказания в задаче оценивания относительного декремента поглощения по сейсмическим записям, Геология иГеофгаика, N2.

8. Денисов М.С., Фиников Д.Б., .1997, Способ оценивания амплитудного спектра сейсмического импульса и алгоритм "амплитудной деконволюции", Геофизика, N2.

9. Денисов М.С., Силаенков O.A., 2003, Пример использования процедур прямого и обращенного продолжения волнового поля в задаче построения глубинного разреза, Геофизика, N3.

10. Денисов М.С., Силаенков O.A., Фиников Д.Б., 2003, Кинематическая фильтрация как способ подавления помех в системе VELINK: Геофизика, Спецвыпуск: Технологии сейсморазведки, N2.

11. Денисов М.С., Фиников Д.Б., 2005, Адаптивная нестационарная коррекция амплитуд при вычитании кратных волн, Технологии сейсморазведки, N1.

12. Денисов М.С., Лангман СЛ., Фиников Д.Б., 2002, Экстраполяция волнового поля в задаче моделирования кратных волн (с целью их подавления), Геофизика, N6.

13. Денисов М.С., Оберемченко Д.М., Фиников Д.Б., 1999, Амплитудная дсконволкшия сейсмических записей с учетом частотно-зависимого поглощения. Геофизика, N4.

14. Денисов М.С., Фиников Д.Б., 2005, Способ подавления шумов дискретизации при суммировании сейсмических трасс (на примере моделирования кратных волн), Геофизика, N1.

15. Денисов М.С., Фиников Д.Б., 2002, Использование локального направленного суммирования для экстраполяции волнового поля. Часть 1., Геофизика, N1.

16. Денисов М.С., Фиников Д.Б., 2002, Использование локального направленного суммирования для экстраполяции волнового поля. Часть 2., Геофизика, N2.

17. Денисов М.С., Фиников Д.Б., 2005, Современные методы подавления кратных волн в сейсморазведке: теория и опыт применения. Международная научно-практическая конференция Геомодель-2005. Сборник тезисов.

18. Denisov М., and Finikov D., 1997, Use of Wavelet Amplitude Spectrum Paiametrization for Deconvolution - Examples from Western Siberia, International Geophysical Conference and Exhibition, Istanbul'97.

19. Denisov, M.7 and Finikov, D., 2001, Wavefield extrapolation based on the iocai slant stacking: 63th Ann. Intemat. Mtg. EAGE.

20. Denisov, M., and Finikov, D., 2002, An alias protection scheme for Radon transform: 64th Ann. Internat. Mtg. EAGE.

21. Lokshtanov, D., Denisov, M., and Finikov, D., 2002, Multiple suppression and datuming with an antialiased Radon transform for sea-floor data: 64lh Ann. Internat. Mtg. EAGE.

22. Denisov, M., and Finikov, D., 2004, Efficient multiple attenuation with a dealiased velocity filtering in f-x domain: 66th Ann. Internat. Mtg. EAGE.

23. Denisov, M., and Finikov, D. 2001, A compact non-stationary wavelet parametrization for deconvolution and Q estimation: 71th Ann. Internat. Mtg. SEG.

24. Denisov, M., Finikov, D., Langman, S., Oberemchenko, D., and Spitz, S., Surface Related and Internal Multiple Attenuation. A Modeling Approach.: International Geophysical Conference and Exhibition, Istanbul'97.

■ V"

i с"

Подписано в печать 2о. оЗ 2006 г. Объем 2-5 пл. Тираж Л 00 экз. Заказ

Редакциокно-издательскнй отдел РГГРУ Москва, ул. Миклухо-Маклая, 23

Содержание диссертации, доктора физико-математических наук, Денисов, Михаил Сергеевич

1. ВВЕДЕНИЕ.

2. УСТОЙЧИВЫЕ МЕТОДЫ ПОДАВЛЕНИЯ КРАТНЫХ ВОЛН В СЕЙСМОРАЗВЕДКЕ.

2.1 Общая классификация алгоритмов.

2.2 Двухшаговые методы.

2.2.1 Прогнозирование волнового поля.

2.2.2 Особенности реализации.

2.2.3 Нестационарное адаптивное вычитание модели кратных волн из исходных сейсмограмм, примеры обработки реальных данных.

2.2.3 О верификации результатов прямого и обратного продолжений волнового поля.

Введение Диссертация по наукам о земле, на тему "Устойчивые методы подавления регулярных и нерегулярных помех в сейсморазведке"

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

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

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

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

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

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

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

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

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

Изложение способов решения задачи подавления регулярных помех производится в следующей последовательности, В Главе 2 рассматриваются алгоритмы прогнозирования кратных волн при помощи методов прямого продолжения волнового поля и устойчивые способы их адаптивного вычитания из исходных данных. Там же описаны кинематические фильтры, применение которых целесообразно после адаптивного вычитания. Способы построения устойчивых операторов продолжения волнового поля рассмотрены в Главах 3 и 4. Область их применимости не ограничивается исключительно прогнозированием кратных волн, но использование соответствующих подходов оказывается полезным и при решении других задач, а разработанными алгоритмами подавления артефактов любого преобразования сейсмических данных, включающего пространственное суммирование трасс, полезно оснастить не только операторы продолжения волнового поля (как это сделано в Главах 3, 4), но и кинематические фильтры, описанные в Главе 2, причем соответствующие подходы могут применяться при фильтрации как в пространственно-частотной, так и в пространственно-временной области.

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

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

Актуальность работы

Повсеместное применение таких «тонких» процедур как AVO анализ [98], СВАН [99], интерпретация данных на основе вейвлет - разложения [36], [158], [218] и т.д., которыми оснащены все современные популярные пакеты обработки сейсмических данных, и которые являются неотъемлемыми этапами завершающих стадий обработки и интерпретации, требуют более тщательной предварительной обработки, в особенности применения динамически корректных устойчивых преобразований, не вносящих искажений в волновое поле. Например, методы подавления кратных волн должны не только эффективно ослаблять регулярный шум, но и при этом тщательно сохранять динамические и кинематические характеристики однократных отражений. Далеко не все даже самые современные алгоритмы обработки удовлетворяют этим требованиям, что не просто затрудняет последующее применение, например, AVO анализа, но и приводит к ложным интерпретационным заключениям (см., например, [115], [143]). Это же относится и к анализу данных на основе вейвлет- преобразования (см., например [184]).

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

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

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

Цель работы

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

Задачи исследований

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

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

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

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

Получить алгоритм устойчивой кинематической фильтрации с подавлением артефактов преобразования в пространственно-временной t-x и в пространственно-частотной / - х областях.

Разработать алгоритм помехозащищенного прямого и обратного продолжения волнового поля.

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

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

Научная новизна и личный вклад

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

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

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

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

5. Разработаны способы коррекции динамических искажений полезных волн, возникающих при нестационарной по временной и пространственной координатам кинематической фильтрации, применяемой в t-x области.

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

7. Предложен алгоритм выделения областей накапливания конструктивных сумм при пространственном суммировании сейсмических трасс. Это позволяет получать результаты, свободные от аляйсинг- помех.

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

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

10. Оценивание параметра частотно-зависимого поглощения производится оптимизационным способом, причем функционалом является дисперсия ошибки предсказания.

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

Защищаемые положения

В диссертации защищаются следующие основные научные положения:

1. Предложены статистически устойчивые методы учета нестационарности сейсмических волновых полей для их обработки с целью построения сейсмо-геологической модели среды.

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

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

4. Предложен метод ослабления шумов пространственной дискретизации (аляйсинг-помех) при прямом и обратном продолжении волнового поля.

5. Разработан алгоритм помехоустойчивого пространственного суммирования сейсмических трасс при редкой сети наблюдений.

6. Предложен устойчивый метод оценивания амплитудного спектра сигнала и коэффициента поглощения по малым выборкам.

Апробация

По теме диссертации опубликовано 24 научные работы. Результаты исследований были изложены в докладах на нескольких ежегодных международных конференциях SEG, EAGE, а также па международных конференциях Istanbul-97 и Геомодель-2005. Большое количество результатов обработки реальных данных представлено в отчетах ООО Геотехсистем (ген. директор А.А. Пудовкин, научный директор В.М. Глоговский). Все описанные алгоритмы обработки сейсмических материалов реализованы в системе VELINK (совместный продукт компаний Геотехсистем и CGG), ориентированной на обработку сейсмических материалов 2D и в системе PRIME - обработка данных 3D сейсморазведки, используемыми как отечественными, так и зарубежными геофизиками. В самой диссертационной работе приводятся многочисленные иллюстрации, подтверждающие эффективность применения предложенных методов при обработке реальных данных, полученных, в первую очередь, в сложных сейсмогеологических условиях, а также при нерегулярной сети наблюдений. Производится сопоставление полученных результатов с традиционными методами обработки.

Кроме того, вопросы реализации алгоритмов и методологические аспекты приведены в работе [27] - Денисов М.С., Силаенков О.А., 2003, Пример использования процедур прямого и обращенного продолжения волнового поля в задаче построения глубинного разреза, Геофизика, N3. Там же даны разнообразные примеры обработки реальных данных, не приведенные в тексте диссертации.

Автор глубоко признателен проф. А.А. Никитину, фактически ставшему инициатором написания этой диссертации. Успешному проведению исследований способствовали консультации и поддержка В.М. Глоговского на протяжении более пятнадцати лет совместной работы. Значительная часть алгоритмических разработок была выполнена совместно с Д.Б. Финиковым, которого автор сердечно благодарит за многолетнее сотрудничество и который оказал решающее влияние на формирование научных взглядов соискателя. Автор выражает искреннюю признательность своим коллегам Ю.А. Харитонову, Д.М. Оберемченко, Е.А. Курину, А.Е. Фирсову, C.J1. Лангману, О.А. Силаенкову (ООО Геотехсистем) за предоставление соискателю возможности практической реализации научных разработок и за оказанную помощь при обработке реальных данных, анализе и интерпретации результатов. Автор благодарен Д. Локштанову (Norsk Hydro) за плодотворные обсуждения и обеспечение финансовой поддержки исследований.

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

Заключение Диссертация по теме "Геофизика, геофизические методы поисков полезных ископаемых", Денисов, Михаил Сергеевич

Результаты работы алгоритма демонстрируются на двух фрагментах реального временного разреза из Западной Сибири. Данный пример был выбран из тех соображений, что на Западно - Сибирских материалах верхняя часть разреза часто сильно отличается от нижней по частотному составу, и приходится выбирать сравнительно узкие окна настройки. Эти отличия можно видеть по оценкам спектра мощности, приведенным на Рис. 50. Отметим, что когда спектр оценивается в диапазоне частот, за пределами диапазона оценка дополняется константой единичной амплитуды (так же, как и на Рис. 47), что соответствует логике алгоритма, корректирующего спектр в заданном интервале частот и оставляющего его без изменения на остальных частотах. Напомним также, что оценивается спектр смеси сигнала и помехи.

Рис. 51(a) Рис. 51(6)

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

Результаты обработки разреза представлены на Рис. 51 и Рис. 52, где показаны: (а) -фрагменты исходного разреза, (6) - результаты предсказывающей деконволюции, (в) -результаты амплитудной деконволюции оператором, рассчитанным в точке глобального (то есть без учета ограничения (54)) минимума функционала (53), (г) - результаты амплитудной деконволюции оператором, рассчитанным с использованием процедур выбора р и количества базисных функций. Этим экспериментам соответствуют оценки спектров, показанные на Рис. 50.

Рис. 51(b) Рис. 51 (г)

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

Рис. 52(a) Рис. 52(6)

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

Можно видеть, что волновая картина на Рис. 51 (в), (г) и Рис. 52(в), (г) существенно проще, чем на Рис. 51(6) и Рис. 52(6) соответственно. Отличия фрагментов разрезов, изображенных на Рис. 51 (в) и Рис. 52(в) от Рис. 51 (г) и Рис. 52(г) менее существенны. Можно, однако, видеть, что записи на Рис. 51 (г) и Рис. 52(г) выглядят более регулярными, чем на Рис. 51(в) и Рис. 52(b). Некоторые отличия отмечены на разрезах Рис. 51 стрелками А и В. Иначе выглядит синфазность, проявившаяся на Рис. 52(г) на времени «730мс. Понятно, что такого рода изменения волновой картины не могут быть оценены с точки зрения повышения информативности данных. Для этого требуется более содержательная интерпретация волнового поля. Здесь же можно лишь утверждать, что наблюдаемые эффекты вполне объясняются теми теоретическими соображениями, на основе которых построен алгоритм, и которые изложены выше.

Рис. 52(b) Рис. 52(г)

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

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

5.2 Амплитудная деконволюция с учетом частотно-зависимого поглощения

Частотно-зависимое поглощение - один из основных факторов, приводящих к изменению формы регистрируемого сигнала отраженной волны во времени. Хорошо известно, что большинство алгоритмов обратной фильтрации опираются на гипотезу стационарности формы импульса в окне настройки оператора. Поэтому в граф обработки требуется включать процедуру коррекции частотно-зависимого поглощения. Такие алгоритмы известны довольно давно [88], [118], и ими оснащены практически все системы обработки сейсмических записей. Опыт свидетельствует, что применяются такие процедуры намного реже, чем это представляется необходимым или, во всяком случае, оправданным [128]. Дело в том, что применению любого способа коррекции поглощения должно предшествовать оценивание коэффициента поглощения (возможно, как функции от времени и пространственной координаты), что представляет собой самостоятельную и весьма нетривиальную задачу [1], [17].

Обзор подходов к решению задачи оценивания поглощения по отраженным волнам содержится в [45] и в более поздней книге [3]. Здесь мы не будем повторять рассмотрение традиционных методов, кратко упомянув лишь те из них, на основе которых в дальнейшем были развиты отдельные направления. Чуть подробнее будут освещены более современные алгоритмы. Как показано в [3], все известные методы основаны на совместном анализе амплитудных спектров (как вариант - лог- амплитудных спектров или кепстров) двух фрагментов сейсмической записи, при делении которых остаток соответствует оператору поглощения, по которому и вычисляется искомый параметр (например, [68], [220]). Здесь начинает существенно сказываться влияние специфики реализации импульсной трассы в каждом из окон. Ослабить этот эффект призван алгоритм, полностью игнорирующий все отсчеты спектра, кроме частоты, на которой реализуется его максимум [173], а о величине поглощения судят по смещению максимума во второй базе относительно первой. Понятно, что такой подход не только не имеет достаточных теоретических обоснований, но и, как показала практика, обладает еще меньшей устойчивостью.

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

В работе [125] фактически подводится своеобразный итог многолетним исследованиям и указывается на существенную неустойчивость традиционных статистических оценок, получаемых в многооконном режиме путем сопоставления лог-амплитудпых спектров фрагментов. Причиной этого, по мнению авторов, является влияние реализации импульсной трассы, которая пусть даже и является белым шумом, но короткие ее реализации могут иметь, вообще говоря, произвольные спектральные характеристики. Поэтому авторы предлагают привлекать оценки импульсной трассы из скважинных данных. Заметим, что сама идея совместного использования данных ГИС и сейсморазведки для повышения устойчивости алгоритмов определения поглощения была высказана значительно раньше [174]. Этот же подход развивают авторы исследования [191], также указывая на неустойчивость оценок, получаемых только по результатам регистрации отраженных волн.

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

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

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

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

Тогда нестационарная модель трассы p(t) может быть представлена в виде n{t) = p(t)=s(t) * /(/)+n{t) = - м)

Q \ T )

55) T где s(t) - стационарная форма сигнала, /(/) - нестационарная отражательная характеристика среды, 7(t, т)- изменяющаяся во времени форма сигнала, A(t,q) - линейный оператор, описывающий эффект частотно-зависимого поглощения, £(/) - импульсная трасса, n(t) -аддитивный шум. Традиционно используется оператор вида

A((0,q) = ехр(- c(q)qP((в) + jaq), (56) здесь A((o,q) - преобразование Фурье оператора A(t,q) по первой координате, a(q) -функция декремента поглощения, Р(со) = |со| + у'Я(со), причем #(со) и |со| связаны преобразованием Гильберта. Такой вид функции Р(со) соответствует минимально-фазовой модели поглощения [45], [161]. Величину "iq=a(q)q будем называть коэффициентом поглощения.

Пусть сейсмическую трассу можно разбить на несколько интервалов, в пределах которых форма сигнала постоянна. Тогда на каждом из этих интервалов фрагмент трассы pk(t) может быть описан сверточной моделью рк(0 = ч(0*Ш+ч( 0, (57)

162 где к - номер интервала, sk(t) - форма импульса с "включенным" в нее поглощением, которому соответствует коэффициент ук.

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

Модель (57) может быть проиллюстрирована реальными данными из Западной Сибири, представленными на Рис. 53. Справа от разреза (Рис. 53(6)) показаны формы сигналов, полученные в коротких, неперекрывающихся по времени окнах, причем четыре первых импульса оценивались в окрестности соответствующих интенсивных пачек. Отвлекаясь от способов получения таких оценок, отметим, что эффект частотно-зависимого поглощения существенно сказывается на форме импульса, а деконволюция требует выбора очень коротких окон настройки.

Надо сказать, что (57) - обычная модель, оправдывающая применение упомянутой выше многооконной деконволюции [72] или способов оценивания декремента поглощения пугем сопоставления амплитудных спектров интервалов сейсмических трасс [68].

В Разделе 5.1 (см. также [31], [105], [109]) предложен способ параметризации, позволяющий заменить процедуру прямого вычисления логарифма амплитудного спектра сигнала оцениванием коэффициентов его разложения по тригонометрическому базису в заданном диапазоне частот (©1,(02), что приводит к получению более устойчивых результатов при малых окнах настройки. Решение задачи сводится к минимизации функционала (53), который в многооконном режиме может быть записан в виде

Логарифм амплитудного спектра фильтра f(t, а) представлен разложением по тригонометрическому базису у, (со) (49), (51), (52), причем f(t,а) - минимально-фазовый фильтр, и /(0,а) = 1, а минимизация производится с учетом ограничения (54), где р -параметр, определяющий, наряду с N, "сложность" оцениваемого амплитудного спектра [26]. Как видно из (58), здесь производится поиск одного фильтра для всех окон настройки. Небольшое количество параметров а,- позволяет получать устойчивые оценки при сравнительно малых объемах выборок. Понятно, тем не менее, что, хотя амплитудная деконволюция более устойчива по отношению к влиянию реализации импульсной трассы, чем обычные методы, существуют ограничения, которые могут быть преодолены лишь увеличением окна настройки оператора. Надо сказать, что эти недостатки свойственны всем методам, опирающимся на статистическую модель данных и применяемых к малым фрагментам. Рассмотрим пример, приведенный на Рис. 54. Слева (Рис. 54(a)) показаны фрагменты трасс в трех небольших окнах, отличающиеся частотно-зависимым поглощением. Результаты деконволюции по каждому из окон отдельно (многооконный режим деконволюции) показаны на Рис. 54(6). Видно, что в последнем фрагменте результат совершенно неудовлетворителен, и сложность формы записи после деконволюции целиком обусловлена влиянием импульсной трассы. Идеальный результат в диапазоне частот показан на Рис. 54(г).

58) где и - времена начала и конца к -го окна настройки.

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

Если решается и задача оценивания поглощения, то базис ц/Дсо) дополняется функцией Q(со) с искомым коэффициентом 0^ перед ней, а функционал (58) преобразуется в jk

Jp(a,0) = XtOi(^a,O)*p(O)2, (59) к т{ но количество параметров, подлежащих оцениванию, увеличивается лишь на число временных окон настройки, причем 0! = 0. Функция Q(со) и коэффициенты 0^ связаны с Р(а) и ук: так как постановка задачи не подразумевает влияния фазового спектра оператора fk(t,а) (функционал (59) квадратичен), а также требует центрированности базиса в диапазоне частот [25], то вместо A((o,q) используется оператор, имеющий модуль амплитудного спектра

Гехр(0А0(со)),со е Ц,со2) д^ = ш сР|+со2 1,соеЦ,со2) ' 2

Такая замена позволяет применить критерий минимума дисперсии. Более подробно этот вопрос рассмотрен в Приложении к Главе 5.

Таким образом, при решении задачи (59) оцениваются параметры оператора few,а,9)| = N Л ехр 6*0(ю) +j]V/(<0)oi< ,сое(со,,со2)

V /=1 J со^сог) при прежнем ограничении (54) и с учетом условия

61)

6. Заключение

В работе рассмотрены устойчивые методы ослабления регулярных и нерегулярных помех в сейсморазведке, при этом особенно подробно изучалась задача подавления кратных волн. Были проанализированы традиционные методы, и были предложены новые подходы. В Главе 2 изложены теоретические основы прогнозирования поля кратных волн, при этом детальное рассмотрение соответствующих преобразований вида продолжения волнового не осуществлялось, а все выражения, относящиеся к пространственному накапливанию трасс, выписаны пе как суммы, а как интегралы. Подробности реализации таких преобразований при обработке реальных данных с учетом дискретности по пространственной и временной координатам изложены в Главах 3 и 4. Глава 3 посвящена разработке устойчивой модификации алгоритма продолжения волнового поля с привлечением априорной информации о модели среды, то есть прогнозирования поля кратных волн в соответствии с выражением (9). В Главе 4 изучены способы повышения устойчивости преобразования вида (10), не требующего информации о модели среды для получения модели кратных волн.

Многие из предложенных алгоритмов требуют оценки амплитудного спектра сейсмического импульса. Статистически устойчивый способ получения таких оценок описан в Главе 5. Разработанные алгоритмы могут применяться не только в упомянутых задачах, но область их применения значительно шире. Например, методы, описанные в Главе 4, могут быть использованы для помехоустойчивого суммирования трасс, методы Главы 3 - для любых процедур, подразумевающих продолжение волнового поля, а алгоритм амплитудной деконволюции (Глава 5) - для оценивания и коррекции частотно-зависимого поглощения.

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

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

Продолжение исследований особенно актуально в задаче анализа кратных волн, но уже не в русле совершенствования алгоритмов их прогнозирования или адаптивного вычитания, но разработки алгоритмов, которые будут применяться либо до этапа подавления, либо после, но будут использовать информацию, содержащуюся в оценке поля кратных отражений. Выше уже упоминалось, что существует подход, в рамках которого утверждается, что в сейсморазведке, основанной на методе отраженных волн, следует использовать все зарегистрированные сигналы, и информация о глубинном строении земли, извлеченная из кратных волн, может успешно дополнять соответствующие построения, производимые только по полю однократных отражений (можно, например, мигрировать кратные волны [136], [183], [93]). Действительно, комбинирование информации, содержащейся в поле кратных и однократных волн, может обеспечить получение дополнительных сведений о динамических свойствах исследуемой среды. Рассмотрение этого подхода остается за рамками работы, но здесь, намечая план дальнейших исследований, кратко укажем на его потенциальные возможности. В параграфе, посвященном адаптации поля кратных волн к исходным трассам, говорилось, что физический смысл получаемых по критерию МНК оптимальных фильтров это либо оператор деконволюции s~l(t) (при прогнозировании без модели среды), либо характеристика отражения от кратнообразующего горизонта (предсказание кратных волн в рамках известной модели).

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

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

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

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

Сочетая оба метода прогнозирования и используя сильные стороны и преимущества каждого из них, можно развить новый подход к решению задачи подавления кратных волн, который пока не применялся на практике и также является новым, но, скорее, с методологической точки зрения. Пусть требуется получить кинематически и динамически адекватную модель всего цуга кратных волн, связанных с отражением от дневной поверхности, но зарегистрированные данные не позволяют это сделать. Так бывает, например, при морской сейсморазведке на мелководье (соответствующий пример дан в Разделе 2.2.3, см. Рис. 16 - Рис. 20). В этом случае отсутствует (или сильно зашумлена) однократно отраженная волна, которая и играет определяющую роль при предсказании, фактически являясь функцией Грина для прогнозирования наиболее интенсивных волн, связанных с переотражением от морского дна. Выход может быть найдем путем комбинирования тех особенностей и преимуществ, которые имеет каждый из методов. Пусть зарегистрированные данные Р(со) описываются сумой однократных 0(a) и кратных К(со) волн Р(со) = 0((й) + К((й). Выделим из поля однократных волну 0\ (со), которая не была уверенно зарегистрирована, а все остальные волны обозначим через 02(&): 0(a) = 0\ (о) + 02(со). Аналогично и для поля кратных волн: А'(со) = (и) + К2 (со), где А] (со) - весь набор кратных волн (полнократные и частично-кратные), связанных с выбранным кратнообразующим горизонтом, а К2(со) - остальные кратные волны. Тогда, решив прямую задачу, получим "чистый" (то есть незашумлениый) динамический годограф однократного отражения Oj (со), тем самым, фактически имеем функцию Грина для решения задачи прогнозирования волн, связанных с переотражением от этого горизонта. Используя метод прямого продолжения (9), получим оценку /^(о) поля К] (со). После адаптивного вычитания АГ](со) из Р{со) получим Р(со) = Р(со) - ^i(co)« Oi(co) + 02(<х>) + К2(а). Теперь стоит задача прогнозирования и вычитания поля К2((£>), но сделать это непосредственно по Р(со) нельзя, так как последнее содержит О) (со). Можно «вырезать» эту волну либо простым мьютингом (обнулением некоторой области вокруг ее годографа), либо адаптивно вычесть из нее сигнал Ot (со). Тогда получим ?(со) = ?(со) - О, (со) = 02 (со) + К2 (ю). Очевидно, что Р{со)

- поле однократных и связанных только с ними кратных волн. Тогда, подставив Р{со) в (10) получим прогноз кратных 02(&). Таким образом, нам удалось построить алгоритм прогнозирования всех кратных волн, связанных с отражением от свободной поверхности, даже в случае отсутствия в исходном поле адекватно зарегистрированных однократных отражений.

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

Так, одним из приложений может быть поверхностно-согласованная деконволюция [15] [97] или иная коррекция за поверхностные условия [193]. Спецификой задачи является то, что уравнения, которые приходится решать, нелинейны, но линеаризуются в лог-спектральной области. Переход из временной области в лог-спектральиую настолько неустойчив [147], что сразу же теряются все потенциальные преимущества такого метода деконволюции. Действительно, соответствующие линеаризованные уравнения могут быть записаны только в отсутствие аддитивного шума, а его наличие всегда вносит непредсказуемое смещение оценок. Кроме того, при явном вычислении логарифма трудно контролировать негативное влияние тех частот, где амплитудный спектр близок к нулю, а это очень часто наблюдается и внутри полезного диапазона. С этими и другими сложностями позволяет справиться предложенный метод, использующий компактное разложение и оптимизацию по критерию минимума дисперсии ошибки предсказания. Параметризация осуществляется в области лог-спектров, так что линеаризованные уравнения просто-напросто запишутся относительно коэффициентов разложения. Но, что принципиально важно, для оценивания этих параметров нет необходимости перевода данных в лог-спектры, а вычисления можно проводить либо во временной, либо в спектральной областях. Это позволяет избежать неустойчивости на стадии логарифмирования. Оптимизационный подход к оцениванию позволяет контролировать влияние аддитивного шума на результаты.

Интересным представляется также модификация соответствующих алгоритмов с целью их применения к данным 4D сейсморазведки. На интересующей площади с некоторой периодичностью проводятся сейсмические наблюдения с одной и той же геометрией расстановки. Одним из критериев контроля качества разработки месторождения, находящего на этой площади, является изменение параметра частотно-зависимого поглощения. В настоящее время такие работы и оценивание поглощения проводят в традиционном режиме, сопоставляя спектры (чаще всего опять-таки в лог-спектралыюй области, где параметр поглощения определяется методами линейной регрессии) в двух временных интервалах, при этом спектральный анализ производится независимо по каждому интервалу. С целью повышения устойчивости и, как следствие, достоверности оценок, целесообразно объединить в один функционал (как это было сделано в (63)) все трассы, относящиеся к некоторой пространственной области и, согласно модели, оценивать один набор параметров разложения для всей группы трасс и один коэффициент относительного поглощения. Такое осреднение целесообразно осуществлять при помощи устойчивых методов суммирования [18].

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

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

Библиография Диссертация по наукам о земле, доктора физико-математических наук, Денисов, Михаил Сергеевич, Москва

1. Авербух А.Г., 1982, Изучение состава и свойств горных пород при сейсморазведке, М., Недра.

2. Ампилов Ю.П., 1992, Поглощение и рассеяние сейсмических волн в неоднородных средах. М, Недра.

3. Андерсон Т., 1976, Статистический анализ временных рядов. М., Мир.

4. Беркаут А., 1986, Миграция сейсмограмм и подавление кратных волн как решение обратной задачи сейсморазведки, ТИИЭР, Том 74, N3.

5. Бреховских Л.М., 1957, Волны в слоистых средах, М., Изд. Академии наук СССР.

6. Булатов М.Г., Телегин А.Н., Левый Н.В., 1980, Весовые коэффициенты при дифракционном преобразовании сейсмических записей, Разведочная геофизика, N90, М., Недра.

7. Варакин Л.Е., 1970, Теория сложных сигналов, М., Советское радио.

8. Глоговский В.М., Лангман С.Л., Фиников Д.Б., 1998, Погружение волнового поля -альтернатива миграции до суммирования, Нефтегаз, 165-171.

9. Глоговский В.М., Фиников Д.Б, 1987, Кинематические фильтры миграционного преобразования реальных сейсмических наблюдений, Сборник докладов третьего научного симпозиума стран-членов СЭВ по нефтяной геофизике, М.

10. Голд Б., Рейдер Ч., 1973, Цифровая обработка сигналов, М., Советское радио.

11. Гольдин С.В., 1974, Линейные преобразования сейсмических сигналов, М., Недра.

12. Гольдии С.В, 1987, Динамический анализ изображений в сейсмике, Геология и Геофизика, N2.

13. Гольдин С.В, 1985, Интегральные продолжения волновых полей, Геология и Геофизика, N3

14. Гольдин С.В., Митрофанов Г.М., 1975, Спектрально-статистический метод учета поверхностных неоднородностей в системах многократного прослеживания отраженных волн, Геология и Геофизика, N 6.

15. Гольцман Ф.М., 1964, Основы теории интерференционного приема регулярных волн, М., Наука.

16. Данквардт Э., Патцер У., 1982, Комбинированное определение поглощения с помощью кепстралыюго анализа. Сборник докладов второго научного семинара стран-членов СЭВ по нефтяной геофизике. Т.1, М.

17. Денисов М.С., 1995, О рекурсивном НЧ фильтре для обработки гравимагнитных данных, Физика Земли, N11.

18. Денисов М.С. 1998, Новые возможности спектрального оценивания по методу максимальной энтропии, Физика Земли, N3.

19. Денисов М.С., 2005, Динамические особенности современных алгоритмов прогнозирования и вычитания кратных волн, Технологии сейсморазведки, N3.

20. Денисов М.С., 2006, Анализ метода прогнозирования кратных волн без знания модели среды с позиций теории продолжения волновых полей, Геофизика, N1.

21. Денисов М.С., 1995, Как интерпретировать алгоритм предсказывающей деконволюции, Геофизика, N3.

22. Денисов М.С., 1992, Эффективные и робастные оценки корреляционных функций: Геология, геофизика и разработка нефтяных месторождений, М., ВНИОЭНГ, N 12.

23. Денисов М.С., Бусыгин Г.В., 1995, Повышение временной разрешенности сейсмических записей методом робастной деконволюции, Геофизика, N5.

24. Денисов М.С., Фиников Д.Б., 1999, Критерий минимума дисперсии ошибки предсказания в задаче оценивания относительного декремента поглощения по сейсмическим записям, Геология и Геофизика, N2.

25. Денисов М.С., Фиников Д.Б., 1997, Способ оценивания амплитудного спектра сейсмического импульса и алгоритм "амплитудной деконволюции", Геофизика, N2.

26. Денисов М.С., Силаенков О.А., 2003, Пример использования процедур прямого и обращенного продолжения волнового поля в задаче построения глубинного разреза, Геофизика, N3.

27. Денисов М.С., Силаенков О.А., Фиников Д.Б., 2003, Кинематическая фильтрация как способ подавления помех в системе VELINK, Геофизика, Спецвыпуск: Технологии сейсморазведки, N2.

28. Денисов М.С., Фиников Д.Б., 2005, Адаптивная нестационарная коррекция амплитуд при вычитании кратных волн, Технологии сейсморазведки, N1.

29. Денисов М.С., Лангман С.Л., Фиников Д.Б., 2002, Экстраполяция волнового поля в задаче моделирования кратных волн (с целью их подавления), Геофизика, N6.

30. Денисов М.С., Оберемчепко Д.М., Фиников Д.Б., 1999, Амплитудная деконволюция сейсмических записей с учетом частотно-зависимого поглощения, Геофизика, N4.

31. Денисов М.С., Фиников Д.Б., 2005, Способ подавления шумов дискретизации при суммировании сейсмических трасс (на примере моделирования кратных волн), Геофизика, N1.

32. Денисов М.С., Фиников Д.Б., 2002, Использование локального направленного суммирования для экстраполяции волнового поля. Часть 1., Геофизика, N1.

33. Денисов М.С., Фиников Д.Б., 2002, Использование локального направленного суммирования для экстраполяции волнового поля. Часть 2., Геофизика, N2.

34. Денисов М.С., Фиников Д.Б., 2005, Современные методы подавления кратных волн в сейсморазведке: теория и опыт применения. Международная научно-практическая конференция Геомодель-2005. Сборник тезисов.

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

36. Кауфман А.А., Левшин А.Л., Ларнер К.Л., 2003, Введение в теорию геофизических методов, акустические и упругие волновые поля, Часть 4., М., Недра.

37. Кац С.А., Киселевич В.Л., Шубик Б.М., 1973, Комплекс методов обнаружения сейсмических волн: В кн. Сейсмические волны в тонкослоистых средах. М., Наука.

38. Кац С.А., Птецов С.Н., 1978, Спектральный анализ поля регулярных сейсмических сигналов и помех, Изв. АН СССР, Сер. Физика земли, N 1.

39. Кац С.А., Шубик Б.М., 1977, Адаптивные веерные фильтры, Физика Земли, N8.

40. Клаербоут Д.Ф., 1989, Сейсмическое изображение земных недр, М. Недра.

41. Козлов Е.А., 1982, Распознавание и подавление многократных волн в сейсморазведке, М. Недра.

42. Козлов Е.А., 1986, Миграционные преобразования в сейсморазведке, М., Недра.

43. Кондратьев И.К., 1976, Линейные обрабатывающие системы в сейсморазведке. М. Недра.

44. Кондратьев O.K., 1986, Сейсмические волны в поглощающих средах, М. Недра.

45. Корн Г., Корн Т., 1977, Справочник по математике, М., Наука.

46. Малкин А.Л., 1989, Негауссовская статистическая модель сейсмической записи, Геология и геофизика, N 1.

47. Малкин А.Л., Сорин АЛ., Фиников Д.Б., 1985, Селективная предсказывающая деконволюция сейсмических записей., НТИС Сер. Нефтегазовая геология, геофизика и бурение. Вып. 10, М., ВНИИОЭНГ.

48. Малкин А.Л., Фиников Д.Б., 1988, Оптимизационный способ коррекции фазового спектра сейсмического сигнала, Геология и геофизика, N 6.

49. Малкин А.Л., Фиников Д.Б., 1989, Параметризация фазового спектра сейсмического сигнала, Геофизический журнал, т.11, N 3.

50. Малкин А.Л., Фиников Д.Б., 1988, Фазовая деконволюция. Вопросы реализации, Геология и геофизика, N 4.

51. Малкин А.Л., Фиников Д.Б., 1988, Фазовая деконволюция. Теоретический аспект, Геология и геофизика, N 3.

52. Марпл С.JI., 1990, Цифровой спектральный анализ и его приложения. М., Мир.

53. Математическая энциклопедия, 1984, М., Советская энциклопедия.

54. Мушип И.А., 1983, Конструирование алгоритмов и графов обработки данных сейсморазведки, М. Недра.

55. Наттерер Ф., 1990, Математические аспекты компьютерной томографии, М., Мир.

56. Нахамкин С.А., 1977, Интерференционные преобразования сейсмических полей, В кн. Вопросы динамической теории распространения сейсмических волн, Вып. XVII, Л., Наука.

57. Нахамкин С.А., 1969,0 веерной фильтрации, Физика земли, N11.

58. Нахамкин С.А., 1973, Параболическая фильтрация сейсмограмм, В кп. Вопросы динамической теории распространения сейсмических волн, Вып. XIII, Л., Наука.

59. Нахамкин С.А, Владимиров Ю.М., Решетников В.В., 1977, О методике обращенного волнового продолжения: В кн. Вопросы динамической теории распространения сейсмических волн. Вып. XVII, Л., Наука.

60. Нахамкин С.А., Рудаков А.Г., 1972, Обобщенные операторы некоторых интерференционных преобразований и их двумерные спектральные аналоги, Физика Земли, N10.

61. Никитин А.А., 2005. О новых- старых алгоритмах обработки геофизических данных, Международная научно-практическая конференция Геомодель-2005. Сборник тезисов.

62. Оппенгейм А.В., Шафер Р.В., 1977, Цифровая обработка сигналов. М., Связь.

63. Петрашень Г.И., Нахамкин С.А., 1973, Продолжение волновых полей в задачах сейсморазведки: Л., Наука.

64. Полшков М.К., Козлов Е.А., Мешбей В.И. и др., 1984, Системы регистрации и обработки данных сейсморазведки, М. Недра.

65. Полиа Г., Сеге Г., 1978, Задачи и теоремы из анализа. Часть1. М, Наука.

66. Рабинер Л., Гоулд, Б., 1978, Теория и применение цифровой обработки сигналов: М., Мир.

67. Рапопорт М.Б., 1969, О некоторых сейсморазведочных приложениях корреляционной теории. В кп. Прикладная геофизика, Вып. 56., М., Недра.

68. Робинсон Э.А., 1988, Метод миграции в сейсморазведке, М., Недра.

69. Робинсон Е., Трейтел С., 1980, Цифровая обработка сигналов в геофизике, В кн. Применение цифровой обработки сигналов, М., Мир.

70. Рытов С.М., 1976, Введение в статистическую радиофизику. Часть 1: случайные процессы, М., Наука.

71. Сейсморазведка: Справочник геофизика. Т2. Под ред. Номоконова В.П. М, Недра, 1990.

72. Теория и практика сейсмического метода РНП. 1962. Под редакцией J1.A. Рябинкина Гостоптехиздат.

73. ТИИЭР. Спектральное оценивание, тематический выпуск. 1982, т.70, N9.

74. Фиников Д.Б., 1989, Коррекция нестационарности сейсмических трасс, вызванной частотно-зависимым поглощением: Вопросы обработки и комплексной интерпретации в сейсморазведке, М.: ВНИИОЭНГ.

75. Шубик Б.М., 1980, Адаптивная фильтрация сейсмограмм общего пункта взрыва (ОПВ), Разведочная геофизика, N 90, М., Недра.

76. Abma R., Kabir N., Matson К., Michell S., Shaw S., VcLain В., 2005, Comparison of adaptive subtraction methods for multiple attenuation, The Leading Edge, Vol.27, N 3, pp. 277-280.

77. Al'Ai R., Verschuur D.J., Reagan R., 2003, Land data prestack surface multiple elimination strategy, case study in north Africa, 73rd Ann. Intemat. Mtg. SEG.

78. Albertin U„ Stork C., Yingst D., Chang W., Fletcher R., Kitchenside P., 2004, Amplitude behavior of Kirchhoff wavefield extrapolation, and beam migration in areas of poor illumination, 66th Ann. Internat. Mtg. EAGE.

79. Artman В., Shragge J., Biondi В., 2003, Operator aliasing in wavefield continuation migration, 73rd Ann. Internat. Mtg. SEG.

80. Baina R., Nguyen S., Noble M., Lambre G., 2003, Optimal anti-aliasing for ray-based Kirchhoff depth migration, 73rd Ann. Internat. Mtg. SEG.

81. Berkhout, A.J., 1981, Wave field extrapolation techniques in seismic migration, a tutorial: Geophysics, Vol. 46, pp. 1638-1656.

82. Berkhout A.J., 1999, Multiple removal based on the feedback model, The leading edge, Vol. 18, pp. 127-131.

83. Berkhout A. J., Verschuur D.J., 1997, Estimation of multiple scattering by iterative inversion, Parti and 2, Geophysics, Vol. 62, N5, pp. 1586-1595,1596-1611.

84. Berryhill J.R., 1979, Wave-equation datuming, Geophysics, Vol. 44, pp. 1329-1344.

85. Berryhill J.R., 1984, Wave-equation datuming before stack (short note), Geophysics, Vol. 49, pp. 2064-2066.

86. Berryhill J.R. Kim Y.C., 1986, Deep-Water Peg Legs and Multiples: Emulation and Suppression: Geophysics, Vol. 51, pp. 2177-2184.

87. Bickel S.H., Natarajan R.R., 1985, Plane-Wave Q-Deconvolution. Geophysics, Vol.50, N9.

88. Bleistein N., Cohen J.K., Stockwell J.W., 2000, Mathematics of multidimensional seismic imaging, migration, and inversion, Springer, New York.

89. Bleistein N., Jaramilo H., 1998, A platform for data mapping in scalar models of data acquisition, 68th Ann. Intemat. Mtg. SEG.

90. Borselen R.G., Schonewille M.A., Hegge R.F., 2005, 3D surface-related multiple elimination: acquisition and processing solutions, The Leading Edge, Vol. 24, N3, pp. 260-268.

91. Brown M., 2004, Least-squares joint imaging of multiples and primaries, 74th Ann. Internat. Mtg. SEG.

92. Burg, J. P., 1975, Maximum entropy spectral analysis: Ph.D. dissertation, Stanford Univ.

93. Burg, J. P., 1972, The relationship between maximum entropy spectra and maximum likelihood spectra, Geophysics, Vol. 37, pp. 375-376.

94. Brysk, H., McCowan, Douglas W., 1986, A slant-stack procedure for point-source data, Geophysics, Vol. 51, pp. 1370-1386.

95. Cambois G., Stoffa P., 1992, Surface-consistent deconvolution in the log-Fourier domain, Geophysics, Vol. 57, N 6, pp. 823-840.

96. Castagna J.P., Backus M.M. ed., 1995, Offset-dependent reflectivity theory and practice of AVO analysis, SEG, Tulsa.

97. Castanga J.P., Sun S., Siegfried R., 2003, Instanteneous spectral analysis: Detection of low-frequency shadows associated with hydrocarbons, The Leading Edge, Vol. 22, N2, pp. 120-127.

98. Chauris H., Salomons В., 2004, Data-driven Kirchhoff migration or how to take advantage of the slant stack domain, 66th Ann. Internat. Mtg. EAGE.

99. Carcione J., Herman G., Kroode A.P.E., 2002, Seismic modeling, Geophysics, Vol. 74, N4.

100. Clark R.A., Carter A.J., Nevill P.C., Benson P.M., 2001, Attenuation measurmnets from surface seismic data Azimuthal variation and time-lapse case studies, 63rd Ann. Internat. Mtg. EAGE.

101. Claerbout J.F., 1972, Downward continuation of moveout corrected seismograms, Geophysics, Vol. 37, pp. 741-768.

102. Dasqupta R., Clark R., 1998, Estimation of Q from surface seismic reflection data, Geophysics, N6, pp. 2120-2128.

103. Denisov M., Finikov D., 1997, Use of Wavelet Amplitude Spectrum Parametrization for Deconvolution Examples from Western Siberia, International Geophysical Conference and Exhibition, Istanbul'97.

104. Denisov M., Finikov D., 2001, Wavefield extrapolation based on the local slant stacking, 63rd Ann. Internat. Mtg. EAGE.

105. Denisov M., Finikov D., 2002, An alias protection scheme for Radon transform, 64th Ann. Internat. Mtg. EAGE.

106. Denisov, М., Finikov, D., 2004, Efficient multiple attenuation with a dealiased velocity filtering in f-x domain, 66th Ann. Internat. Mtg. EAGE.

107. Denisov M., Finikov D., 2001, A compact non-stationary wavelet parameterization for deconvolution and Q estimation, 71st Ann. Internat. Mtg. SEG.

108. Denisov M., Finikov D., Langman S., Oberemchenko D., and Spitz S., Surface Related and Internal Multiple Attenuation. A Modeling Approach., Internat. Geophysical Conference and Exhibition, Istanbul'97.

109. Doicin D., Spits S., 1991, Multichannel attenuation of high amplitude peg-legs: Examples from the North Sea, 53rd Ann. Internat. Mtg. EAGE.

110. Dragoset В., 1992, Surface multiple attenuation theory, practical issues, and examples, 54th Ann. Internat. Mtg. EAGE.

111. Dragoset W.H., Jericevic Z., 1998, Some remarks on surface multiple attenuation, Geophysics, Vol. 63, pp. 772-789.

112. Dragoset В., Weglein A., 1999, An introduction: The new world of multiple attenuation, The Leading Edge, Vol. 18, N1.

113. DuBose J., 2003, A modification of the parabolic Radon transform for the preservation of AVO effects, 73rd Ann. Internat. Mtg. SEG.

114. Ferreira C., Cruz J., 2004, Modified Kirchhoff prestack depth migration using the Gaussian beam operator as Green's function, 66th Ann. Internat. Mtg. EAGE.

115. Foster D.J., Mosher C.C., 1992, Suppression of multiple reflection using the Radon transform, Geophysics, Vol. 57, pp. 386-395.

116. Gelius L.J., 1987, Inverse Q-filtering. A spectral balancing technique, Geophysical Prospecting, Vol.35, N2.

117. Glogovsky V.M., Gogonenkov G.N., 1988, Study of Methods for Determining Velocity and Depth Parameters in Layered Realistic Media: Geophysical transactions, Vol. 33, N3-4, pp.157-173.

118. Gray, S., 1997, True-amplitude seismic migration: A comparison of three approaches, Geophysics, Vol. 62, pp. 929-936.

119. Gray S.H., 1992, Frequency-selective design of the Kirchhoff migration operator, Geophysical Prospecting, Vol. 40, pp. 565-572.

120. Guitton A., 2005, Sparse Radon transforms with bound-constrained optimization, 67 Ann. Internat. Mtg. EAGE.

121. Guitton A., Cambois G., 1999, Multiple elimination using a pattern-recognition technique, The Leding Edge, Vol.1, pp. 92-98.

122. Guitton A., Verschuur D., 2004, Adaptive subtraction of multiples using the LI norm, Geophysical Prospecting, Vol. 52, N 1.

123. Hackert С., Parra J., 2004, Improving Q estimates from seismic reflection data using well-log-based spectral correction, 74th Ann. Internat. Mtg. SEG.

124. Hampson D., 1988, Inverse velocity stacking for multiple elimination, J. Can. Soc. Expl. Goephys., Vol. 22, pp. 49-55.

125. Hargreaves N., Cooper N., 2001, High-resolution Radon demultiple, 74th Arm. Internat. Mtg. SEG.

126. Hargreaves N., Roberts G., Wombell R., 2005, Q-guided wavelet-domain amplitude correction, 67th Ann. Internat. Mtg. EAGE.

127. Hargreaves N., verWest В., Wombell R., Trad D., 2003, Multiple attenuation using the apex-shifted Radon transform, 73rd Ann. Internat. Mtg. SEG.

128. He R., 2004, Waveform prediction of water-layer related multiples, 74th Ann. Internat. Mtg. SEG.

129. Herrmann P., Mojesky Т., Magesan M., Hugonnet P., 2000, De-aliased, high-resolution Radon transforms, 70th Ann. Internat. Mtg. SEG.

130. Hicks G., 2001, Removing NMO stretch using the Radon and Fourier-Radon transforms: 63th Ann. Internat. Mtg. EAGE.

131. Hill N.R., Langan R.T., Nemeth Т., Zhao M., Bube K.P., 2002, Beam methods for predictive suppression of seismic multiples in deep water, 72nd Ann. Internat. Mtg. SEG

132. Hindriks C.O.H., Duijndam A.J.W., 1998, Radon domain reconstruction of 3D irregularly sampled VSP data, 68th Ann. Internat. Mtg. SEG.

133. Hugonnet P., Herrmann P., Ribeiro C., 2001, High resolution Radon: a review, 63th Ann. Internat. Mtg. EAGE.

134. Jiang Z., Yu J., Schuster G., Hornby В., 2005, Migration of multiples, The Leading Edge, Vol. 27, N3, pp. 315-318.

135. Jones I., 2003, A review of 3D PreSDM model building techniques, First break, Vol.21, N.3, pp. 45-60.

136. Jones I., Fruen K., 2003, Factors affecting frequency content in preSDM imaging, The Leading Edge, Vol. 22, N2, pp. 128-134.

137. Kabir N., Abma R., 2003, Weighted subtraction for diffracted multiple attenuation, 73rd Ann. Internat. Mtg. SEG.

138. Kabir N., Abma R., Ganyuan X., 2004, 3D wavefield extrapolation based demultiple in Ormen Lange, 74th Ann. Internat. Mtg. SEG.

139. Kabir N., Marfurt K.J., 1999, Toward true amplitude multiple removal, The Leading Edge, Vol.1, pp. 66-73.

140. Kabir M.M., Verschuur D.J., 1995, Restoration of missing offsets by parabolic Radon transform, Geophysical prospecting, Vol. 43, pp. 347-368.

141. Kapoor S.J., O'Brien M., Stork C., Woodward M., 2003, Integrating complementary tools fro improved depth imaging, 73rd Ann. Internat. Mtg. SEG.

142. Kelamis P.G., Verschuur D.J., 2000, Surface-related multiple elimination on land seismic data Strategies via case studies, Geophysics, Vol. 63, pp. 719-734.

143. Kostner U., Buske S., 1999, Computing geometrical spreading from traveltimes in 3D, 69th Ann. Internat. Mtg. SEG.

144. Lecerf D., Reiser C., 2004, High-resolution processing for time-lapse seismic, 74th Ann. Internat. Mtg. SEG.

145. Levin S., 1989, Surface-consistent deconvolution, Geophysics, Vol. 54, pp. 1123-1133

146. Levin S., 2002, Prestack poststack 3D multiple prediction, 72nd Ann. Internat. Mtg. SEG.

147. Levin F.K., Shah P. M., 1977, Peg-Leg Multiples and Dipping Reflectors, Geophysics, Vol. 42, pp. 957-981.

148. Levy S., Oldenburg D.W., 1987, Automatic phase correction of common-midpoint stacked data, Geophysics, Vol. 52, pp. 51-59.

149. Lin D., Young J., Huang Y., 2004, 3D SRME application in the Gulf of Mexico, 66th Ann. Internat. Mtg. EAGE.

150. Lin D., Young J., Huang Y., 2005, 3D SRME practice for better imaging, 67th Ann. Internat. Mtg. EAGE.

151. Lokshtanov D., 1995, Multiple suppression by single channel and multichannel deconvolution in the tau-p domain, 65th Ann. Internat. Mtg. SEG.

152. Lokshtanov D., 2001, Suppression of water-layer multiples and peg-legs by wave-equation approach, 63th Ann. Internat. Mtg. EAGE.

153. Lokshtanov D., 1999, Multiple suppression by data-consistent deconvolution, The Leading Edge, pp. 115-119.

154. Lokshtanov D., Denisov M., Finikov D., 2002, Multiple suppression and datuming with an antialiased Radon transform for sea-floor data, 64th Ann. Internat. Mtg. EAGE.

155. Malkin A.L., Finikov D.B., 1986, A statistical approach to correction of the seismic wavelet phase (phase deconvolution): Proceedings of the 31 international geophysical symposium, Gdansk, Vol. 1, pp. 130-139.

156. Mallat S., 1989, A theory for multiresolution signal decomposition: The wavelet representation, IEEE Trans. Pattern. Anal. Machine Intell., Vol. 11. pp. 674-693.

157. Martini F., Bean C., 2002, Application of pre-stack wave equation datuming to remove interface scattering in sub-salt imaging, The First Break, Vol. 20, pp.395-403

158. Matson K.H., Abma R., 2005, Fast 3D surface-related multiple elimination using azimuth moveout for multiples, 75th Ann. Internat. Mtg. SEG.

159. McCarley L.A., 1985, An Autoregressive Filter Models for Constant Q-Attenuation. Geophysics, Vol.50, N5, pp. 749-758.

160. McGlynn J.D., Ioup G.E., 1985, Phase-coherency filtering of reflection seismic data, Geophysics Vol. 50, N 9, pp. 1505- 1509.

161. Mendel J.M., 1990, Maximum likelihood deconvolution, Springer-Verlag, New York.

162. Mendel J.M., 1991, Tutorial on higher-order statistics (spectra) in signal processing and system theory: theoretical results and some applications, Proceedings of the IEEE, Vol. 79, pp. 278305.

163. Monk D.J., 1993, Wave-equation multiple suppression using constrained cross-equalization, Geophysical Prospecting, Vol. 41, pp. 725-736.

164. Moore I., Bisley R., 2005, 3D surface-related multiple prediction (SRMP): a case history, The Lading Edge, Vol.27, N.3, pp. 270-284.

165. Moore I., Dragoset В., 2004, Practical, 2D surface-related multiple prediction (SRMP), 72nd Ann. Internat. Mtg. SEG.

166. Neidell N. S., 1991, Could the processed seismic wavelet be simpler than we think?, Geophysics, Vol. 56, pp. 681-690.

167. Neidell N.S., Taner M.T., 1971, Semblance and other coherency measures for multichannel data, Geophysivs, Vol. 36, pp. 482-497.

168. Novotny M., 1990, Trace interpolation by slant-stack migration: Geophysical prospecting, Vol.38, pp. 833-851

169. Pica A., Poulain G., David В., Magesan M., Baldock S., Weisser Т., Hugonnet P., Herrmann Ph., 2005,3D surface related multiple modeling, The Leading Edge, Vol.24, N3, pp. 292-296.

170. Porsani M.J., Ursin В., 1998, Mixed-phase deconvolution, Geophysics, Vol. 63, pp. 637-647.

171. Quan Y., Harris J.M., 1995, Seismic attenuation tomography using the frequency shift method, Geophysics, Vol. 62, pp. 895-905.

172. Raikes S.A., White R.E., 1984, Measurements of earth attenuation from downhole and surface seismic recordings, Geophysical Prospecting, Vol. 32, pp. 892-919.

173. Riley D.C., Claerbout J.F., 1972,2D multiple reflections, Geophysics, Vol. 41, pp. 592-620.

174. Robinson E. A., 1957, Predictive decomposition of seismic traces: Geophysics, Vol. 22, pp. 767-778.

175. Robinson E.A., Treitel S., 1967, Principles of digital Wiener filtering, Geophysical Prospecting, Vol. 15, N3, pp. 311-333.

176. Sacchi M., Porsani M., 1999, Fast high resolution parabolic Radon transform, 69th Ann. Internat. Mtg. SEG.

177. Sacchi M.D., Ulrych T.J., 1998, Recovery of near offsets using a F-X gap filtering algorithm, 68th Ann. Internat. Mtg. SEG.

178. Saggaf M.M., Robinson E.A., 2000, A unified framework for the deconvolution of traces of nonwhite reflectivity, Geophysics, Vol.65, pp. 1660-1676.

179. Santos L., Schleicher J., Tygel M., 2000, Modeling, migration, and demigration, The Leading Edge. Vol. 22, pp. 712-715

180. Schultz, P.S., Claerbout J.F., 1978, Velocity estimation and downward continuation by wavefront synthesis, Geophysics, Vol.43, pp. 691-714.

181. Schuster G., 2003, Imaging the most bounce out of multiples, 65th Ann. Internat. Mtg. EAGE.

182. Shang В., Caldwell D., 2003, Seismic frequency enhancement through spectral borrowing, 73rd Ann. Internat. Mtg. SEG.

183. Shoenwille M., Hegge В., VanBorselen R., 2005, A comparison of sparse techniques for 3D SRME, 67th Ann. Internat. Mtg. EAGE.

184. Schultz, P.S. and Claerbout J.F., 1978, Velocity estimation and downward continuation by wavefront synthesis: Geophysics, Vol.43, pp.691-714.

185. Spitz, S., 1991, Seismic trace interpolation in the f-x domain: Geophysics, Vol. 56, pp. 785794.

186. Spitz, S., 1999, Pattern recognition, spatial predictability, and subtraction of multiple events, The Leading Edge, Vol.21.

187. Shtivelman V., Canning A., 1988, Datum correction by wave equation extrapolation, Geophysics, Vol. 53, pp. 1311-1322.

188. Spetzler J., Kvam O., 2003, Time-lapse monitoring in the prestack domain, 73rd Ann. Internat. Mtg. SEG.

189. Stainsby S.D., Worthington M.H., 1985, Q estimation form vertical seismic profile data and anomalous variations in the central North Sea, Geophysics, Vol.35, pp. 387-403.

190. Taner M.T., 1980, Long-period sea-floor multiples and their suppression, Geophysical Prospecting, Vol. 28, pp. 30-48.

191. Taner M.T., Koehler F., 1981, Surface consistent corrections, Geophysics, Vol. 46, N 1, pp. 7-22.

192. Tegtmeier S., Gislof A., Verschuur D.J., 2004, 3D sparse-data Kirchhoff redatuming, Geophysical Prospecting, Vol. 24, pp. 509-521.

193. Tonn R., 1991, The determination of seismic quality factor Q from VSP data A comparison of different computational methods, Geophysical Prospecting, Vol. 39, pp. 1-28.

194. Turner G., 1990, Aliasing in the tau-p transform and the removal of spatially aliased coherent noise, Geophysics, Vol. 55, pp. 1496-1503.

195. Trad D., 2002, Interpolation with migration operators, 72nd Ann. Internat. Mtg. SEG.

196. Trad D., Ulrych Т., Sacchi M., 2003, Latest views on the sparse Radon transform, Geophysics, Vol.68, pp. 386-399.

197. Trad D., Ulrych Т., Sacchi M., 2002, Accurate interpolation with high-resolution time-variant Radon transforms, Geophysics, Vol.67, N 2, pp. 644-656.

198. Treitel S., Shanks J.L., Frasier C.W., 1967, Some aspects of fan filtering, Geophysics, Vol. 32, pp. 789-800.

199. Tribolet J.M., 1978, Application of short-time homomorthic signal analysis to seismic wavelet estimation: Geoexploration, Vol. 16, pp. 75-96.

200. Tsai C.J., 1985, Use of autoconvolution to suppress first-order, long-period multiples, Geophysics, Vol. 50, N9, pp. 1410-1425.

201. Tygel M., Schleicher J., Hubral P., Santos L.T., 1998, 2.5-D true-amplitude Kirchhoff migration to zero offset in laterally inhomogeneous media, Geophysics, Vol. 63, pp. 557-573.

202. Van Borselen R.G., Verschuur D.J., 2003, Optimization of marine data acquisition for the application of 3D SRME, 73rd Ann. Internat. Mtg. SEG.

203. Van Dedem E.J., Verschuur D.J., 2001, 3D Surface multiple prediction using sparse inversion, 71th Ann. Internat. Mtg. SEG.

204. Verschuur D.J., Berkhout A.J., Wapenaar C.P.A., 1992, Adaptive surface-related multiple elimination, Geophysics, Vol.57, pp. 1166-1177.

205. Vidale J.E., Houston H., 1990, Rapid calculation of seismic amplitudes, Geophysics, Vol. 55, pp. 1504-1507.

206. Walden A.T., White R.E., 1992, Some fads and fallacies in seismic data analysis, Geophysical prospecting., V.40, N2

207. Wang Y., 2003, Multiple subtraction using an expanded multichannel filter, Geophysics, Vol. 68, Nl,pp. 346-354.

208. Wang Y., 2002, Quantifying the effectiveness of stabilized inverse filtering, Geophysics, Vol. 68, Nl,pp. 337-345.

209. West В., 2002, Suppressing peg-leg multiples with parabolic Radon demultiple, 64th Ann. Internat. Mtg. EAGE.

210. Wiggins W., 1999, Multiple attenuation by explicit wave extrapolation to an interpreted horizon, The Leading Edge, Vol. 18, pp. 46-54.

211. Wiggins J.W., 1984, Kirchhoff integral extrapolation and migration of nonplanar data, Geophysics, Vol.49, pp. 1239-1248.

212. Wiggins J.W., 1988, Attenuation of complex water-bottom multiples by wave-equation-based prediction and subtraction, Geophysics, Vol. 53, pp. 1527-1539.

213. Wiggins R., 1978, Minimum entropy deconvolution, Geoexploration, Vol. 16, pp. 21-35.

214. Wapenaar C.P.A., Peels G.L., Budejicky V., Berkhout A.J., 1989, Inverse extrapolation of primary seismic waves, Geophysics, Vol. 54, pp. 853-863.

215. Wapenaar C.P.A., Verschuur D.J., Herrmann P., 1992, Amplitude preprocessing of single and multicomponent seismic data, Geophysics, Vol. 57, pp. 1178-1188.

216. Yu Z., McMechan G.A., Ferguson J.F., Anno P.D., 2002, Adaptive wavelet filtering of seismic data in the wavelet transform domain, J. of seismic exploration, Vol. 11, pp.223-246.

217. Zhang Y., Gray S., Sun J., Notfors C., 2001, Theory of migration anti-aliasing, 71st Ann. Internat. Mtg. SEG.

218. Zhang C., Ulrych Т., 2002, Estimation of quality factors from CMP data, Geophysics, Vol. 67, pp. 1542-1547.

219. Zhang R., Ulrych Т., 2003, Multiple suppression based on the migration operator and a hyperbolic median filter, 73rd Ann. Internat. Mtg. SEG.

220. Zhao M., Hill N.R., Nemeth Т., Shank R.D., 2005, Application of a model-based multiple attenuation method to a Gulf of Mexico deepwater dataset, The Leading Edge, Vol.27, N3, pp. 285290.

221. Zhou В., Greenhalgh, S.A., 1994, Linear and parabolic x-p transforms revisited: Geophysics, Vol. 59, pp. 1133-1149.