Содержание к диссертации
Введение
2. Устойчивые методы подавления кратных волн в сейсморазведке 11
2.1 Общая классификация алгоритмов 11
2.2 Двухшаговые методы 15
2.2.1 Прогнозирование волнового поля 15
2.2.2 Особенности реализации 34
2.2.3 Нестационарное адаптивное вычитание модели кратных волн из исходных сейсмограмм, примеры обработки реальных данных 38
2.2.3 О верификации результатов прямого и обратного продолжений волнового поля 53
2.2.4 Приложения к Главе 2 60
2.3 Кинематическая фильтрация как способ подавления регулярных помех 73
2.3.1 Реализация кинематических фильтров в t — X области 74
2.3.2 Коррекция амплитуд полезных отражений, искаженных кинематической фильтрацией 79
2.3.3 Фильтрация разрезов равных удалений, суммарных и мигрированных разрезов 82
2.3.4 Реализация кинематических фильтров в X — (0 области 83
2.3.5 Примеры обработки реальных данных 88
3. Построение устойчивого оператора продолжения волнового поля с использованием локального направленного суммирования 94
3.1 Алгоритмы продолжения поля 96
3.2 Расчет компенсирующих фильтров 101
3.3 Алгоритмы экстраполяции поля с сохранением амплитуд 107
3.4 Ослабление артефактов преобразования 109
3.5 Результаты опробования алгоритма на реальных данных 115
3.6 Особенности алгоритмов прямого продолжения поля в задаче прогнозирования кратных волн 120
3.7 Приложения к Главе 3 129
4. Способ подавления шумов дискретизации при суммировании сейсмических трасс 134
5. Параметрический оптимизационный способ оценивания амплитудного спектра сейсмического сигнала и декремента поглощения 146
5.1 Параметризация амплитудного спектра импульса и алгоритм амплитудной деконволюции 146
5.2 Амплитудная деконволюция с учетом частотно-зависимого поглощения 160
5.3 Приложение к Главе 5 177
6. Заключение.
- Нестационарное адаптивное вычитание модели кратных волн из исходных сейсмограмм, примеры обработки реальных данных
- Кинематическая фильтрация как способ подавления регулярных помех
- Алгоритмы экстраполяции поля с сохранением амплитуд
- Параметризация амплитудного спектра импульса и алгоритм амплитудной деконволюции
Введение к работе
Развитие алгоритмической базы обработки данных сейсморазведки идет параллельно с совершенствованием вычислительной техники. Внедрение в практику промышленной обработки материалов наблюдений мощных компьютеров позволило использовать новые устойчивые подходы, обеспечивающие повышение надежности и детальности построений.
Объектом исследования настоящей работы являются традиционные и новые методы обработки данных сейсморазведки, в частности, такие процедуры как подавление кратных волн, прямое и обратное продолжение волнового поля, спектральный анализ и деконволюция, оценивание и коррекция частотно-зависимого поглощения и т.д. При этом особое внимание уделяется повышению устойчивости алгоритмов к неточностям в оценке глубинно-скоростного строения среды, нерегулярности сети наблюдений, наличию в волновом поле факторов, не описываемых моделью, и другим помехам различной природы.
В работе рассматриваются традиционные и предлагаются новые подходы к проблеме ослабления регулярных шумов, в частности - кратных волн, исследуются границы применимости методов, производится их сопоставление на теоретическом уровне, а работоспособность соответствующих алгоритмов изучается на примерах обработки реальных данных. Анализ осуществляется с единых позиций, а именно с точки зрения теории продолжения волновых полей, что позволяет вскрыть недостатки некоторых методов и предложить пути улучшения их работоспособности и достичь универсальности.
Сразу же укажем, что наряду с отношением к кратным волнам, как к регулярному шуму, подлежащему ослаблению, следует развивать и альтернативный подход. Действительно, если зарегистрированные данные содержат как однократные, так и кратные отражения, то целесообразно использовать ту информацию, в первую очередь динамическую, которая содержится в последних, и только после этого осуществлять их подавление. Хотя такая задача здесь подробно не рассматривается, но в Заключении намечаются пути ее решения.
Важной спецификой прямого продолжения волнового поля в задаче прогнозирования кратных отражений является то, что оно верифицируемо. Иначе говоря, всегда имеется относительно надежный критерий проверки качества соответствующих алгоритмов. Действительно, в результате обработки зарегистрированных сейсмограмм прогнозируются кратные волны, которые, в свою очередь, тоже зарегистрированы в процессе записи. Поэтому точность, с которой математическое преобразование переводит исходные сейсмограммы в поле кратных волн, может быть оценена по степени кинематического и динамического соответствия результата реально зарегистрированным волнам. Но такая проверка может быть
осуществлена только после этапа адаптации. Этот вопрос подробно рассмотрен в работе и проиллюстрирован примерами обработки реальных данных.
Волновое поле, полученное в сложной среде, целесообразно подвергать нескольким процедурам, которые ослабляют фон кратных волн, основываясь на различных принципах. Так, на первом этапе удобно использовать двухшаговые алгоритмы, в рамках которых сперва производится прогнозирование поля кратных волн из зарегистрированных сейсмограмм, а затем - адаптивное вычитание. На втором этапе хорошую работоспособность показывают устойчивые кинематические фильтры, выделяющие и подчеркивающие однократные волны на фоне регулярных и нерегулярных помех. При этом неизбежно вносятся искажения формы полезных сигналов, и целью устойчивой фильтрации является учет и компенсация этих искажений для восстановления правильной динамики.
В работе также рассматриваются возможные пути улучшения и повышения устойчивости методов прямого и обратного продолжения волнового поля. Работоспособность предложенных алгоритмов иллюстрируется на примерах обработки реальных данных.
Важное место в любом графе обработки сейсмического материала занимает процедура пространственного суммирования трасс. Из алгоритмов, исследуемых в работе, такой этап содержится во всех методах продолжения волнового поля, в частности, при прогнозировании кратных волн, а также в схемах кинематической фильтрации. Наиболее распространенным же ее применением является построение временных разрезов из сейсмограмм ОГТ. Хорошо известно, что при этом, в силу дискретности данных по пространственной координате, результат суммирования обычно содержит так называемый аляйсинг - шум. Этот эффект особенно заметен при обработке материалов, полученных с редкой или нерегулярной сетью наблюдений и, в особенности, площадной сейсморазведки (3D). Повсеместное внедрение в практику промышленной обработки мощных компьютеров позволило с успехом применить нестандартный устойчивый способ улучшения оператора пространственного суммирования, включив в него процедуру предварительного накапливания локальных направленных сумм. Последние, в отличие от глобальной суммы, сразу преобразующейся в одну выходную трассу, легко поддаются анализу и дополнительной фильтрации.
Использование методов прямого продолжения волнового поля с целью прогнозирования кратных волн приводит к неизбежным искажениям формы импульса в получаемой модели. Такие искажения будут существенно ослаблены, если предварительно произведена коррекция импульса методами деконволюции. Дальнейшая компенсация производится уже на этапе адаптации модели кратных волн к исходным данным. Кроме того, коррекция формы импульса методами обратной фильтрации имеет и более общие приложения. Как правило, наиболее универсальными и эффективными являются статистические подходы к оцениванию оператора деконволюции. При этом сохраняется
особая актуальность такой задачи в условиях нестационарности формы импульса, обусловленной влиянием частотно-зависимого поглощения. Рассмотрены подходы, обеспечивающие устойчивость оценкам амплитудного спектра сигнала в условиях малой выборки, и предложена альтернатива традиционному методу многооконной деконволюции, обеспечивающая лучшие статистические характеристики получаемых параметров. Основой метода является компактная параметризация сейсмического сигнала при помощи разложения логарифма его амплитудного спектра по гладким базисным функциям. Это позволяет свести процедуру деконволюции (или, что в данном случае то же самое, - оценивание амплитудного спектра импульса) к поиску небольшого числа коэффициентов разложения, что придает ей устойчивость. Действительно, статистически надежных результатов оценивания амплитудного спектра сейсмического сигнала по малым выборкам можно добиться только при использовании компактной его параметризации, причем в ограниченной полосе частот. Используется оптимизационный подход, а в качестве критерия предложено применение хорошо зарекомендовавшей себя на практике дисперсии ошибки предсказания. При наличии возможной нестационарности импульса из-за влияния частотно-зависимого затухания, в компактную спектральную параметризацию добавляется всего один коэффициент. При этом можно так модифицировать вид оператора частотно - зависимого поглощения, что методом минимизации дисперсии ошибки предсказания удается получить как оценку спектра сигнала, так и декремент поглощения. Эти оценки могут быть использованы и с целью интерпретации, и для последующей деконволюции и коррекции поглощения. Важной чертой подхода является то, что все окна настройки включены в один функционал, и оценивание производится не отдельно по каждому окну, как в многооконной деконволюции (это значительно ухудшает статистические характеристики оценок), но функционал оптимизируется относительно одного вектора параметров для всех окон.
Изложение способов решения задачи подавления регулярных помех производится в следующей последовательности, В Главе 2 рассматриваются алгоритмы прогнозирования кратных волн при помощи методов прямого продолжения волнового поля и устойчивые способы их адаптивного вычитания из исходных данных. Там же описаны кинематические фильтры, применение которых целесообразно после адаптивного вычитания. Способы построения устойчивых операторов продолжения волнового поля рассмотрены в Главах 3 и 4. Область их применимости не ограничивается исключительно прогнозированием кратных волн, но использование соответствующих подходов оказывается полезным и при решении других задач, а разработанными алгоритмами подавления артефактов любого преобразования сейсмических данных, включающего пространственное суммирование трасс, полезно оснастить не только операторы продолжения волнового поля (как это сделано в Главах 3, 4), но и кинематические фильтры, описанные в Главе 2, причем соответствующие подходы
могут применяться при фильтрации как в пространственно-частотной, так и в пространственно-временной области.
С точки зрения проблемы подавления кратных волн алгоритм коррекции амплитудного спектра импульса (амплитудная деконволюция), предложенный в Главе 5 и основанный на статистически устойчивом оценивании амплитудного спектра сейсмического сигнала, является вспомогательным. Его место в графе обработки должно предшествовать этапу построения модели кратных волн, что приводит к получению более устойчивых результатов, в особенности при использовании алгоритма прогнозирования, не привлекающего информацию о глубинно-скоростной модели среды. Кроме того, надежные оценки амплитудного спектра импульса необходимы при применении методов продолжения волнового поля и кинематической фильтрации с подавлением артефактов этих преобразований вида аляйсинг- шума.
Амплитудную деконволюцию с одновременным оцениванием и коррекцией поглощения следует использовать на завершающем этапе, то есть после подавления регулярных помех, так как их наличие неизбежно приведет к смещению получаемых параметров. При этом само оценивание поглощения уже граничит с интерпретационными процедурами.
Актуальность работы
Повсеместное применение таких «тонких» процедур как AVO анализ [98], СВАН [99], интерпретация данных на основе вейвлет - разложения [36], [158], [218] и т.д., которыми оснащены все современные популярные пакеты обработки сейсмических данных, и которые являются неотъемлемыми этапами завершающих стадий обработки и интерпретации, требуют более тщательной предварительной обработки, в особенности применения динамически корректных устойчивых преобразований, не вносящих искажений в волновое поле. Например, методы подавления кратных волн должны не только эффективно ослаблять регулярный шум, но и при этом тщательно сохранять динамические и кинематические характеристики однократных отражений. Далеко не все даже самые современные алгоритмы обработки удовлетворяют этим требованиям, что не просто затрудняет последующее применение, например, AVO анализа, но и приводит к ложным интерпретационным заключениям (см., например, [115], [143]). Это же относится и к анализу данных на основе вейвлет- преобразования (см., например [184]).
Еще большие требования к качеству применяемых алгоритмов выдвигаются при обработке результатов площадной (3D) сейсморазведки. Неизбежно возникающие артефакты многоканальных преобразований (любое пространственное суммирование трасс порождает
аляйсинг- шум) могут существенно исказить динамику волновой картины, особенно в случае нерегулярной расстановки. В частности, использование продолжения волновых полей требует разбивки плоскости наблюдений на сегменты весьма малого размера, причем в каждом таком сегменте необходимо наличие как приемника колебаний, так и источника. На практике это требование, как правило, не удовлетворяется: данные, полученные в результате наземной сейсморазведки, имеют крайне нерегулярную геометрию со значительным количеством пропусков как пунктов взрыва, так и пунктов приема, а морские данные фактически являются набором независимо полученных линейных профилей. Такая специфика материала приводит к тому, что сводятся на нет почти все потенциальные преимущества площадных наблюдений.
Одной из рассматриваемых проблем является универсальный подход к ослаблению кратных волн. Сейсмическая запись всегда представляет собой интерференцию однократных и кратных отражений, причем последние обычно рассматриваются как помеха, подлежащая подавлению на самых ранних этапах обработки данных. Действительно, зачастую кратные волны столь интенсивны, что они полностью маскируют однократные отражения, и это значительно затрудняет интерпретацию волнового поля.
В последние годы все чаще применяется мониторинг уже разведанных месторождений (см., например, [146], [190]). Здесь изменение такого параметра как частотно-зависимое поглощение является критерием и средством контроля за процессом откачки нефти из продуктивного пласта. При этом традиционные статистические оценки, получаемые по малым окнам, оказываются ненадежными, то есть актуальной является разработка устойчивого алгоритма определения параметра затухания упругих отраженных волн по сейсмическим записям. Кроме того, всегда сохраняется и актуальность самой задачи нестационарной деконволюции с учетом искажения импульса за счет поглощения.
Цель работы
Целью работы является анализ сильных и слабых сторон традиционных алгоритмов подавления кратных волн, продолжения волнового поля, деконволюции, оценивания частотно-зависимого поглощения и др. с точки зрения новых требований, выдвигаемых к качеству применяемых процедур, а также разработка оригинальных методов и походов для помехоустойчивой обработки данных.
Задачи исследований
Получить метод прогнозирования поля кратных волн по зарегистрированным сейсмограммам в любых сейсмогеологических условиях. Выработать помехоустойчивый подход к прогнозированию с использованием оценки глубинно-скоростной модели среды. Исследовать возможность и получить алгоритм прогнозирования без привлечения информации о строении среды.
Разработать универсальный алгоритм помехоустойчивого в смысле минимизации артефактов преобразования (здесь - аляйсинг- шум) пространственного суммирования сейсмических трасс, полученных при редкой сети наблюдений, не прибегая к процедуре интерполяции.
Проанализировать специфику неизбежных ошибок и помех, возникающих при использовании операторов прогнозирования кратных волн как с привлечением информации о глубинно-скоростной модели среды, так и в условиях априорной неопределенности. Разработать алгоритм нестационарного адаптивного вычитания полученной модели кратных отражений из исходного поля, причем устойчивость соответствующего подхода должна проявляться в точном сохранении кинематических и динамических особенностей однократных отражений, присутствующих в зарегистрированном поле в интерференции с кратными волнами.
Исследовать возможности верификации алгоритмов прямого продолжения волнового поля в задаче прогнозирования кратных волн.
Получить алгоритм устойчивой кинематической фильтрации с подавлением артефактов преобразования в пространственно-временной t-x и в пространственно-частотной f-x областях.
Разработать алгоритм помехозащищенного прямого и обратного продолжения волнового поля.
Исследовать статистические свойства оценок амплитудного спектра сейсмического сигнала по малым выборкам. Предложить способ компактной параметризации импульса, позволяющий получать более надежные результаты. Определить вид функционала, обеспечивающего несмещенное оценивание параметров для сигнала с ограниченной полосой частот. Изучить возможность компактной параметризации и оценивания частотно-зависимого поглощения в рамках этого же алгоритма.
Разработать соответствующее программное обеспечение для предложенных алгоритмов и показать их эффективность на практических примерах.
Научная новизна и личный вклад
1. Разработан универсальный подход к задаче прогнозирования поля кратных волн по
исходным сейсмограммам. Рассмотрение ведется с точки зрения теории продолжения
волновых полей, что позволило не только точно очертить границы области применимости
каждого метода, но и провести параллель между различными алгоритмами на теоретическим
уровне, а также предложить улучшения с целью повышения их эффективности.
2. Предложен удобный способ расчета поля кратных волн с правильным
соотношением амплитуд без привлечения информации о глубинно-скоростном строении
среды. При наличии оценки модели среды возможно использование алгоритмов прямого
продолжения волнового поля, при этом гипотеза однородности среды не привлекается
(традиционные способы разработаны для однородных сред). Показано, что продолжение
волнового поля в задаче прогнозирования кратных волн корректнее осуществлять вдоль
семейства лучей ие дифрагированных, а отраженных волн.
3. Разработан метод параметризации вида нестационарности и оценки адаптивных
фильтров в больших окнах настройки для компенсации динамических погрешностей
моделирования поля кратных волн, которые носят нестационарный характер как по
пространственной, так и по временной координате.
4. Предложена процедура контроля уровня артефактов для алгоритмов
кинематической фильтрации, применяемых в пространственно- частотной области. Алгоритм
основан на методе масштабирования амплитудного спектра.
5. Разработаны способы коррекции динамических искажений полезных волн,
возникающих при нестационарной по временной и пространственной координатам
кинематической фильтрации, применяемой в / -х области.
Предложено преобразование волнового поля с целью его прямого и обратного продолжения, отличающееся повышенной устойчивостью за счет использования локального направленного суммирования с контролем помех дискретизации (аляйсинг- шума). Получены компенсирующие фильтры, корректирующие искажения динамики полезных отражений.
Предложен алгоритм выделения областей накапливания конструктивных сумм при пространственном суммировании сейсмических трасс. Это позволяет получать результаты, свободные от аляйсинг- помех.
8. Разработан способ компактной параметризации амплитудного спектра
сейсмического импульса в ограниченном диапазоне частот с целью устойчивого
спектрального анализа и деконволюции.
9. Предложена статистически устойчивая альтернатива методу многооконной
деконволюции, применяемому для коррекции сейсмических трасс с учетом нестационарности
сигнала, обусловленной частотно-зависимым поглощением. Важной особенностью такого подхода является возможность получения несмещенных оценок в ограниченной полосе частот.
10. Оценивание параметра частотно-зависимого поглощения производится оптимизационным способом, причем функционалом является дисперсия ошибки предсказания.
Предложенные алгоритмы и новые методы обработки сейсмических данных, обладающие научной новизной, получены лично автором или при его непосредственном участии.
Защищаемые положения
В диссертации защищаются следующие основные научные положения:
Предложены статистически устойчивые методы учета нестационарности сейсмических волновых полей для их обработки с целью построения сейсмо-геологической модели среды.
Разработаны устойчивые методы адаптивного вычитания кратных волн нестационарными фильтрами с сохранением динамики полезных отражений.
Предложен способ компенсации искажений динамики полезных отражений, возникающих при нестационарной кинематической фильтрации
Предложен метод ослабления шумов пространственной дискретизации (аляйсинг-помех) при прямом и обратном продолжении волнового поля.
5. Разработан алгоритм помехоустойчивого пространственного суммирования
сейсмических трасс при редкой сети наблюдений.
6. Предложен устойчивый метод оценивания амплитудного спектра сигнала и
коэффициента поглощения по малым выборкам.
Апробация
По теме диссертации опубликовано 24 научные работы. Результаты исследований были изложены в докладах на нескольких ежегодных международных конференциях SEG, EAGE, а также на международных конференциях Istanbul-97 и Геомодель-2005. Большое количество результатов обработки реальных данных представлено в отчетах ООО Геотехсистем (ген. директор А.А. Пудовкин, научный директор В.М. Глоговский). Все описанные алгоритмы обработки сейсмических материалов реализованы в системе VELINK (совместный продукт компаний Геотехсистем и CGG), ориентированной на обработку
сейсмических материалов 2D и в системе PRIME - обработка данных 3D сейсморазведки, используемыми как отечественными, так и зарубежными геофизиками. В самой диссертационной работе приводятся многочисленные иллюстрации, подтверждающие эффективность применения предложенных методов при обработке реальных данных, полученных, в первую очередь, в сложных сейсмогеологических условиях, а также при нерегулярной сети наблюдений. Производится сопоставление полученных результатов с традиционными методами обработки.
Кроме того, вопросы реализации алгоритмов и методологические аспекты приведены в работе [27] - Денисов М.С., Силаенков О.А., 2003, Пример использования процедур прямого и обращенного продолжения волнового поля в задаче построения глубинного разреза, Геофизика, N3. Там же даны разнообразные примеры обработки реальных данных, не приведенные в тексте диссертации.
Автор глубоко признателен проф. А.А. Никитину, фактически ставшему инициатором написания этой диссертации. Успешному проведению исследований способствовали консультации и поддержка В.М. Глоговского на протяжении более пятнадцати лет совместной работы. Значительная часть алгоритмических разработок была выполнена совместно с Д.Б. Финиковым, которого автор сердечно благодарит за многолетнее сотрудничество и который оказал решающее влияние на формирование научных взглядов соискателя. Автор выражает искреннюю признательность своим коллегам Ю.А. Харитонову, Д.М. Оберемченко, Е.А. Курину, А.Е. Фирсову, С.Л. Лангману, О.А. Силаенкову (ООО Геотехсистем) за предоставление соискателю возможности практической реализации научных разработок и за оказанную помощь при обработке реальных данных, анализе и интерпретации результатов. Автор благодарен Д. Локштанову (Norsk Hydro) за плодотворные обсуждения и обеспечение финансовой поддержки исследований.
Также хотелось бы отметить многолетнюю поддержку и доброжелательное отношение к работам автора со стороны редакционно-издательского центра ЕАГО с самого момента его основания и по настоящее время. На протяжении более десятилетия автора связывают теплые отношения с первым заместителем главного редактора журнала Геофизика Л.Д. Бовтом. Научный редактор изданий ЕАГО Л.Д. Овчининская неоднократно находила досадные опечатки не только в тексте рукописей, но даже и в формулах.
Нестационарное адаптивное вычитание модели кратных волн из исходных сейсмограмм, примеры обработки реальных данных
О качестве предсказанного поля можно судить по критерию его соответствия реально зарегистрированным волнам. Но предварительно следует применить процедуру адаптации. Действительно, необходимость такого этапа заложена уже в самих расчетных формулах. Преобразование (9), основанное на использовании функции Грина, осуществляет продолжение поля через среду известной структуры до кратнообразующего горизонта и обратно к свободной поверхности. Но при этом нам не известна характеристика отражения от этого горизонта, во всех практически важных случаях являющегося тонкослоистой пачкой. Вследствие этого модель кратных волн получается с точностью до такой характеристики отражения, и в данном случае целью адаптации является ее оценивание. Если прогнозирование производилось при помощи выражения (10), то фильтр должен соответствовать s l(t), иначе говоря, здесь мы имеем дело с оператором обратной фильтрации, а, следовательно, процедура его поиска может быть неустойчивой. Для повышения надежности вычислений полезно применить коррекцию формы сигнала при помощи метода, изложенного в Главе 5, до прогнозирования в соответствии с (10). Если модель получена при помощи (9), то при адаптации нет необходимости решать обратную задачу, а оцениваться будет короткий формирующий фильтр, и сама процедура становится более устойчивой.
Кроме подобных искажений, присущих самой структуре алгоритмов прогнозирования и устраняемых на этапе адаптации при помощи подбора стационарных фильтров, приходится иметь дело и с эффектами иного вида, приводящими к нестационарным искажениям модели.
Как уже упоминалось, специфика алгоритмов прогнозирования кратных волн состоит в том, что их кинематика обычно отображается вполне удовлетворительно, но амплитуды и форма колебаний существенно отличаются от тех, что зарегистрированы в исходной записи. Это связано со следующими факторами: 1) при моделировании кратных волн трудно, а зачастую невозможно правильно учесть геометрическое расхождение, 2) волны разного порядка кратности входят в смоделированное поле с различными множителями, 3) трудно учесть зависимость коэффициентов отражений кратнообразующих горизонтов от угла падения волн и пространственных координат, а также характеристики направленности групп приемников и источников. Из-за различия амплитуд волн различного порядка кратности приходится использовать для вычитания не одну, а несколько моделей кратных волн. Обычно достаточно двух моделей l(t) и l\i). Искажения формы волны, носящие сверточный характер, приводят к тому, что для вычитания приходится использовать многоканальные фильтры, адаптирующие поле кратных волн к исходной записи. Компенсация неточностей при учете геометрического расхождения будет рассмотрена ниже.
Стандартным методом адаптации (см. например, [177]) является подбор формирующего винеровского фильтра g(t) такого, что g(t) = argmin \(р(аМ - g(0 Ka,b,t))2dt. (12) № t
Учитывая все сказанное выше про динамическую и кинематическую адекватность модели кратных волн, становится понятным, что этап адаптации, помимо необходимой нагрузки, связанной с оцениванием характеристики отражения или s \t), является, пожалуй, ключевым звеном всей цепочки процедур подавления кратных волн, призванным, в первую очередь, скомпенсировать неизбежные погрешности прогнозирования. Закономерен вопрос: как создать алгоритм вычитания, который бы учитывал возможность несоответствия, в том числе нестационарного по t, полей p(a,b,t) и l(a,b,t), но, в то же время, не имел бы столько степеней свободы, чтобы скорректировать произвольные различия, так как эти различия, заключающиеся в том, что в p{a,b,t) присутствуют однократные отражения, a l{a,b,t) их не содержит, нас и интересуют, причем, как правило, однократные и кратные волны интерферируют. Таким образом, этап адаптации может не только, а иногда и не столько скомпенсировать эффект автосвертки импульса или неучета характеристики отражения, сколько иметь тенденцию "преобразовать" кратные волны в однократные, хотя с точки зрения критерия МНК будет получено наилучшее соответствие модели исходной трассе
На недостатки традиционной реализации адаптации неоднократно указывалось в различных публикациях. Действительно, одной из основных проблем является необходимость использования длинных формирующих фильтров, причем это физически обоснованно, так как иначе не удается скомпенсировать s l(t). Понятно, что такой режим обладает низкой помехоустойчивостью и зачастую имеет тенденцию ослабления полезных волн. Производились попытки повысить устойчивость, в первую очередь за счет снижения длины формирующего фильтра. Например, в [165] предложен алгоритм, в рамках которого полученная модель кратных волн подвергается нескольким одноканальным фильтрациям (преобразование Гильберта, дифференцирование и.т.д.), а затем такой набор трасс адаптивно вычитается из исходной записи при помощи подбора весовых коэффициентов. На недостаточную эффективность такого подхода указывалось в [209], где вместо весовых коэффициентов предложено использовать короткие фильтры. Такая оптимизационная задача приводит к необходимости решать систему линейных уравнений с блочно-Теплицевой матрицей, для чего может быть применен многоканальный аналог алгоритма Левинсона [44], [70]. Сразу же оговоримся, что подход, развиваемый ниже, чем-то напоминает такую процедуру, но позволяет принимать во внимание нестационарный характер динамических искажений модели кратных волн, что никак не учитывается в [165], [209].
При вычитании кратных волн по данным, полученным при наземной сейсморазведке, приходится сталкиваться с еще большими трудностями из-за большой зашумленности данных как регулярными, так и нерегулярными помехами [144].
Авторы работы [124] видит причину неустойчивости (12) в использовании квадратичного критерия L2, заменяя его на норму L1 (метод наименьших модулей). С другой стороны, как указывается в [77], такая замена влияет на результат вычитания несущественно, но при этом значительно возрастают вычислительные затраты ввиду необходимости решения системы нелинейных уравнений. Также с целью снижения риска подавления однократных предложено использовать раздельное вычитание из исходных сейсмограмм регулярной и нерегулярной компонент модели [139].
В серии исследований, инициированных работой С.Шпица [188] развивается альтернативный подход, основанный на совместном анализе исходного волнового поля и модели кратных волн в пространственно-частотной области, где они параметризуются как процесс авторегрессии, идентифицируются параметры, а на завершающей стадии применяется метод распознавания образов [111], [123]. Процедура адаптации и вычитания также может применяться в локальном режиме к выделенным синфазностям кратных отражений [ПО].
Авторы работ [129], [141] указывают на ненадежность традиционно используемых методов адаптации с вычитанием, а для подавления кратных волн применяют мьютинг в области нелинейного преобразования Радона. Для этого как зарегистрированные сейсмограммы, так и модель кратных волн переводятся в эту область при помощи накапливания сигналов вдоль линейной траектории, по полю кратных волн производится распознавание регулярных помех, которые затем подавляются (методом мьютинга) в исходном поле. Похожий подход, но использующий разложение по локальным линейным направлениям, рассмотрен и в настоящей работе (фильтр-маска - Раздел 2.3.1) уже не как альтернатива адаптации, а как дополняющая ее процедура.
Кинематическая фильтрация как способ подавления регулярных помех
Наряду с методами, основанными на адаптивном вычитании модели кратных волн из исходных данных, существует еще одна группа алгоритмов, осуществляющих выделение и подавление помех, основываясь на отличии их годографов от годографов полезных сигналов. Это так называемые кинематические фильтры, которые были упомянуты выше. Ввиду того, что подобные алгоритмы разделения регулярных волн и подавления помех основаны исключительно на их кинематических характеристиках, будем использовать общее наименование кинематическая фильтрация для всех соответствующих подходов. В отличие от методов, использующих прогнозирование с последующим вычитанием, эта группа содержит одношаговые алгоритмы, то есть соответствующий оператор применяется непосредственно к исходному волновому полю, а на выходе получается результат подавления регулярных помех. По сравнению с адаптивным вычитанием модели кратных волн такой подход имеет определенные преимущества, а именно: если двухшаговый алгоритм полностью игнорировал возможные кинематические отличия, то есть, как бы не использовал еще один возможный критерий разделения волн, то кинематический подход такую информацию привлекает. С другой стороны, он обладает и определенными недостатками: хотя очень часто волны различаются по кинематике на больших удалениях, но на ближних каналах почти всегда их годографы неотличимы (см. Рис. 8 и Рис. 12). Следовательно, и чисто кинематический критерий, как правило, недостаточен для уверенного выделения полезных волн из интерференции с регулярным шумом.
Этот же критерий разделения волн может быть использован и после миграционного преобразования волнового поля [221], но, понятно, что если сигнал и помеха близки по кинематике, то алгоритм теряет свою работоспособность.
Фильтр-маска занимает промежуточное положение между двухшаговыми алгоритмами и кинематическими фильтрами. Для работы фильтра-маски требуется прогнозирование поля кратных волн, по которому в дальнейшем будет настраиваться чисто кинематический алгоритм обнаружения сигналов.
Методы кинематической фильтрации хорошо исследованы и опробованы на модельных и реальных данных. Было предложено три основных способа применения фильтров: в пространственно-временной x области, в пространственно-частотной х-со области и в плоскости двумерного преобразования Фурье к -со. Совместный анализ этих алгоритмов дан в работе [58], а также в множестве последующих публикаций на эту тему, и здесь мы не будем касаться уже проработанного в литературе вопроса о сопоставлении способов реализации фильтрации и лишь заметим, что, как показывают соответствующие исследования, каждый из них имеет свои преимущества, и в современных системах промышленной обработки, как правило, присутствует весь их набор. Обзор традиционных методов кинематической фильтрации приведен, например, в монографии [42], алгоритмы веерной фильтрации описаны в [44], а оптимальные винеровские фильтры - в [12]. В последнее десятилетие тема кинематической фильтрации почти не развивается ни в отечественной, ни в зарубежной геофизике, но чаще используется обработка в области преобразования Радона, подразумевающего суммирование исходных сейсмограмм по линейному, параболическому или гиперболическому направлениям [223]. Реже применяется локальное описание поля набором отрезков годографов линейных синфазностей [187], но такой вариант фильтрации обладает низкой помехоустойчивостью из-за того, что он основан на разделении волн методом подгонки модели авторегрессии в пространственно- частотной области.
Динамические особенности этих, по выражению А.А.Никитина [62], "новых - старых" методов были подробно изучены еще С.А. Нахамкиным в его докторской диссертации [57], а также и в более ранних публикациях [58], [59], [61]. Там же были опробованы алгоритмы фильтрации в области параболических и гиперболических разложений сейсмограмм, аналогичные повсеместно демонстрируемым западными исследователями. Особенностью ситуации является то, что в настоящий момент доступны вычислительные мощности, которые отсутствовали в начале семидесятых годов, и современные геофизики находятся в лучшем положении, имея возможность интенсивного тестирования и совершенствования алгоритмов, вводя соответствующие процедуры в граф массовой обработки данных, тем самым, изучая сильные и слабые стороны таких фильтраций. И именно в этом и только в этом смысле можно говорить о научной новизне предложенных здесь кинематических алгоритмов. Их базовые формулы также выписаны в упомянутых работах С.А. Нахамкина, в ряде статей по методике регулируемого направленного приема (РНП) [73] и исследованиях С.А.Каца и его соавторов [38], [39], [40], [76], но внедрение в практику обработки новых компьютеров позволило, основываясь на этом базисе, продолжить исследования и предложить улучшения, сделав алгоритмы более универсальными, гибкими и устойчивыми. Например, новой является разработанная методика коррекции динамических искажений полезных волн, состоящая из трех различных преобразований. Ранее не изучались методы ослабления артефактов пространственной фильтрации, которые рассмотрены в настоящей работе и некоторые другие предложенные алгоритмы.
Пусть в результате предварительного анализа получены, по крайней мере приблизительно, годографы полезных волн. Фильтр должен пропустить (выделить) без искажения эти волны, а также все волны, наклоны годографов которых лежат в некоторой окрестности (веере направлений). Остальные волны, которые, возможно, присутствуют в сейсмограммах, рассматриваются как помехи и подлежат подавлению.
Такая задача может решаться при помощи преобразования волнового поля, близкого по своей структуре к РНП, то есть при помощи направленного пространственного суммирования исходных трасс в диапазоне направлений, принадлежащих вееру с выбранным раствором. Затем производится суммирование самих направленных сумм. При этом регулярные волны, наклон годографов которых принадлежит вееру, пропускаются, а волны, наклон годографов которых лежит за пределами веера, ослабляются. Эффективность ослабления зависит от параметров алгоритма, в первую очередь, от размера базы пространственного суммирования.
Выражение для преобразованиях\ осуществляемого кинематическим фильтром такого типа, можно записать в следующем виде: L da(x,t) = jp(x-y,t + ay)dy, (17а) -L b{x,t) p(x,t)= jda(x,t)da, (176) a(x,t) где p{x,t) - исходное волновое поле, da(x,t) - направленная сумма, ye(-L,L) пространственная база фильтра, a{x,t) и b{x,t) - границы перебора направлений суммирования (веер), p{x,t)- результат фильтрации. В данном режиме работы алгоритма всегда используется симметричный веер: a(x,t) = -b{x,t). Локальные направленные суммы могут быть получены при помощи эффективного и быстрого рекурсивного алгоритма, описанного в [28]. При фильтрации сейсмограмм ОГТ задается скоростной закон и диапазон возможных изменений скоростей, что и определяет вычисляемые границы веера в каждой точке.
Фильтр, компенсирующий искажение формы сигнала
При суммировании по методу (17) происходит искажение формы сигнала, зависящее от размера базы и раствора веера. Обычно база L выбирается как инвариантный параметр для всех трасс. Тогда, в случае стационарности раствора веера по времени (a(x,t) = а(х) и b(x,t) = b(x)), искажение формы импульса описывается его сверткой с некоторым оператором, который может быть рассчитан при помощи следующих рассуждений.
Алгоритмы экстраполяции поля с сохранением амплитуд
Одним из традиционных методов, используемых для ослабления аляйсинг- помех, является применение так называемого dip filter. На практике встречаются две основных модификации этого алгоритма: им оснащается оператор миграционного преобразования [121], или он встраивается в прямое разложение поля по плоским волнам. Основной принцип этого метода заключается в применении к исходным трассам полосового фильтра с верхней частотой среза, зависящей от производной линии суммирования (или параметра р в х-р преобразовании), при этом, в силу отсутствия априорной информации об углах наклона годографов полезных отражений, обычно используется гипотеза их горизонтальности (если кинематика известна, наиболее эффективным становится применение оператора с тонкой настройкой апертуры). Отметим, что, хотя такой подход иногда и помогает ослабить нежелательные шумы преобразования, он обедняет спектральный состав сигнала, принадлежащего крутым синфазностям [219] (этот эффект также проиллюстрирован на Рис. 4). Кроме того, его малоэффективность становится очевидной из следующего хорошо известного соображения: суммотрасса в х-р области для малого, пусть нулевого р, уже содержит аляйсинг- шумы, обычно накопленные с больших удалений в сортировке ОГТ и со всех удалений в сортировке ОПП или ОПВ, в то время как при этом верхняя частота среза полосового фильтра равна частоте Найквиста, то есть dip filter пропускает все частоты сигнала.
Таким образом, хотя реализация операторов продолжения волнового поля в х-р области улучшает условия накапливания сигнала, результат может оказаться неудовлетворительным в силу интенсивных артефактов преобразования, вызванных неограниченностью базы суммирования по пространственной координате.
Хорошие результаты показывает модификация алгоритма, использующая локальное суммирование сигналов в окрестности точек касания [186], но и такой подход имеет свои недостатки, в частности, требуется точная априорная информация о кинематике волн, что, вообще говоря, неприемлемо в задачах продолжения поля.
Преодоление проблемы аляйсинг- помех всегда связано либо с использованием априорной информации, которая обычно недоступна, либо с привлечением нелинейных алгоритмов анализа волнового поля. Нелинейные алгоритмы в миграционных преобразованиях имеют заслуженную репутацию весьма опасных и ненадежных, однако без них в задаче борьбы с шумами, вызванными дискретностью наблюдений, обойтись принципиально невозможно. Суть аляйсинг- эффекта, как известно, состоит в том, что в дискретном случае плоские монохромные волны определенных частот и наклонов становятся неразличимыми. Поэтому приходится либо привлекать какую-либо информацию об этих наклонах, либо распознавать их при помощи решающих правил на основе гипотез о локальном, например линейном, поведении синфазностей.
Оператор продолжения, основанный на локальном х-р преобразовании и рассмотренный выше, может быть оснащен двумя алгоритмами ослабления аляйсинг- шумов: во-первых, традиционным dip filter, обеспечивающим подавление высокочастотных компонент сигнала в зависимости от наклона локального суммирования и, во-вторых, схемой с адаптивным разделением сигнальной и шумовой компонент волнового поля. В основу последнего подхода положен очень простой и наглядный принцип, уже упоминавшийся в Главе 2 при выводе выражения для устойчивого кинематического фильтра вх-и области: максимумы характеристики суммирования линейной синфазности вдоль прямолинейного направления образуют монотонно затухающую последовательность для данных, непрерывных по пространственной координате, но периодическую в дискретном случае. По положению побочных максимумов последней можно судить о тех частотных компонентах сигнала, которые осложнены пространственным аляйсинг- эффектом. Таким образом, в результате направленного суммирования синфазности, годограф которой в пределах базы близок к линейному, оценка соотношения энергии входных трасс к суммотрассе, полученная в некотором низкочастотном диапазоне, не может превосходить аналогичного соотношения для высоких частот. Иными словами, зная, что аляйсинг- эффект проявляется на суммотрассе в виде локальных разрастаний ее энергетического спектра в области высоких частот и сделав предположение, что некоторый низкочастотный диапазон спектра сигнала не порождает аляйсинг- шума, можно контролировать и редактировать соотношение амплитуд спектральных компонент трасс до и после направленного суммирования.
Теперь перейдем к формальному описанию метода. Для оценивания амплитудного спектра сигнала после суммирования потребуется вспомогательный алгоритм разделения волнового поля на шумовую и сигнальную составляющие. Дело в том, что данные до суммирования обычно содержат весьма интенсивную аддитивную помеху. Даже синфазное накапливание полезного сигнала на ограниченной базе, ослабляющее нерегулярный шум, может приводить к результату с низким соотношением сигнал/помеха, тем самым оценка спектра, производимая по суммотрассе, может сильно отличаться от искомого амплитудного спектра сигнала.
Параметризация амплитудного спектра импульса и алгоритм амплитудной деконволюции
Обсуждается обычная постановка задачи, когда сейсмическая трасса p{t) описывается сверточной моделью y{t) = s(t) Z,(t) на фоне аддитивного шума «(/), то есть / (/) = Я0 + «(0 = (0 4(0 + «(0, (47) где t - индекс дискретного времени, s{t) - сигнал, (/) - импульсная трасса, n(t) - помеха, некоррелированная с y(t). Сигнал s{t) может включать в себя многие компоненты [72].
Импульсная трасса ,(() - реализация стационарного случайного процесса типа белого шума, то есть M{ (/ + x)} = a26(x). Задачи спектрального анализа и обратной фильтрации сейсмических записей так или иначе связаны с оцениванием или формы сигнала s(t) или его амплитудного спектра s(co), где s(&) - преобразование Фурье от s(t). Известно [71], что спектр мощности процесса (47) описывается выражением: Р(со) = а2ф )2+л(со)2, (48) где л(со) - спектр мощности помехи.
Многие трудности спектрального анализа сейсмических записей вызваны негауссовостью процесса ,(/) [24], [47], [74]. Это проявляется в так называемом эффекте "влияния импульсной трассы" на оценки спектра мощности [208]. Он обычно вызывает ложное расщепление или сильную изрезашюсть спектральных оценок, не поддающуюся сглаживанию [19]. Преодоление этих сложностей возможно либо при учете статистических свойств ,((), либо при помощи привлечения некоторых дополнительных сведений о возможных свойствах искомого спектра мощности. Первое из этих направлений неизбежно связано с необходимостью оценивания и фазового спектра сейсмического сигнала [164] , что представляет собой весьма сложную задачу, которую по ряду причин желательно рассматривать отдельно [52].
Ниже будет исследован подход к привлечению некоторых самых общих представлений о характере спектра мощности сейсмической трассы для спектрального анализа и обратной фильтрации.
Будем пока рассматривать задачу в отсутствие помехи n(t). Предположим, что на интервале соє(і,сй2) N s(co) = сс0 exp[a,\/,-(о)], (49) (=1 л(со-соі) H/,(CU) = COS[/ — -]. (49 ) (о2-а і Понятно, что при сколь угодно больших N формулой (49) можно описать практически любой амплитудный спектр [46], а само выражение соответствует разложению в ряд Фурье для симметричной функции, каковой и является логарифм амплитудного спектра вещественного сигнала. Ограничение числа N - количества базисных функций в разложении логарифма амплитудного спектра на интервале (сО], )» предполагает гладкость, то есть отсутствие резких выбросов и провалов этой функции, что сужает область возможных решений. Здесь и в дальнейшем спектр фильтра на интервале (-(й2,-щ) легко получить, воспользовавшись свойствами симметрии преобразования Фурье.
" Это обусловлено тем, что функция распределения случайного процесса (47) зависит от фазового спектра сигнала s(t).
Использование такой параметризации позволяет существенно упростить процесс оценивания и сделать его более устойчивым, так как количество искомых коэффициентов становится малым (обычно 5-7). На такой же подход к разработке алгоритмов устойчивого оценивания сейсмического импульса указывалось, например в исследованиях [163] и [168], а также в работе [172], посвященной методам оценивания неминимально-фазовых сигналов, но там рассматривалась иная задача - оценивание формы импульса. Использовалась его параметризация при помощи модели авторегрессии - скользящего среднего (АРСС). Но так как параметризация касалась самой формы сигнала, то есть как его амплитудного, так и фазового спектра, то она становилась менее компактной и не позволяла достичь существенного повышения устойчивости оценивания. Важной особенностью предлагаемого здесь подхода является независимое рассмотрение амплитудного и фазового спектров: их параметризация и последующее оценивание. Это позволяет получить более устойчивые алгоритмы спектрального анализа и деконволюции. Кроме того, еще более существенным отличием нового подхода к компактной параметризации сейсмического сигнала является то, что он, в отличие от, например, АРСС и других известных моделей [53], аппроксимирует спектр в ограниченном диапазоне частот, что принципиально в задачах обратной фильтрации. Тем самым оценивание оператора производится только на тех частотах, где энергия сигнала значительна, и сама процедура деконволюции становится устойчивой, чего невозможно добиться в рамках традиционных моделей.
