M. С. Денисов
Д. Б. Фиников
ООО "ГЕОТЕХСИСТЕМ", МОСКВА
ООО "ГЕОТЕХСИСТЕМ", МОСКВА
МЕТОДЫ ПОДАВЛЕНИЯ КРАТНЫХ ВОЛН
В СЕЙСМОРАЗВЕДКЕ. Часть 1
ВВЕДЕНИЕ. Волновое поле, получаемое в сейсморазведке, всегда представляет собой и тгерференцию полезных
сигналов и шумов. Они могут быть как некогерентными так и когерентными. К последним, в частности, относятся
кратные волны. Они, как правило, рассматриваются как помеха, подлежащая подавлению на самых ранних этапах
обработки данных. Действительно, зачастую кратные волны столь интенсивны, что они полностью маскируют одно-
кратные отражения, и это значительно затрудняет интерпретацию волнового поля. Разработке и анализу способов
подавления кратных волн посвящена обширная литература. Здесь мы дадим обзор хорошо зарекомендовавших себя
на практике подходов к решению задачи ослабления кратных волн, рассмотрим их теоретические предпосылки,
изучим область применимости каждого метода, а также проиллюстрируем их работоспособность на модельных и
реальных данных.
С другой стороны, существует и иной подход, в рамках которого утверждается, что в сейсморазведке отраженных
волн следует использовать все зарегистрированные сигналы, и информация о глубинном строении Земли, извле-
ченная из кратных волн, может успешно дополнять соответствующие построения, производимые только по полю
однократных отражений. Действительно, комбинирование информации, содержащейся в поле кратных и однократ-
ных волн, может обеспечить получение дополнительных сведений о свойствах исследуемой среды, но анализ таких
подходов выходит за рамки статьи.
ОДНОКРАТНЫЕ И КРАТНЫЕ ОТРАЖЕННЫЕ ВОЛНЫ
Сейсморазведка методом отраженных волн основана на регистрации колебаний, вызванных источником, располо
женным на свободной поверхности. Как правило, сейсмоприемники также расположены на поверхности земли
Волна, распространяясь от источника в нижнее полу-
пространство, достигает границ раздела слоев. На каж-
дой границе часть энергии волны проходит вниз, в сле-
дующий слой, а другая часть отражается вверх. Отра-
женные волны подходят снизу к свободной поверхнос-
ти, где они регистрируются сейсмоприемниками. Так
как и сама свободная поверхность является контраст-
ной границей раздела сред (земля-воздух), коэффици-
ент отражения от которой почти равен -I, то волны
вновь персотражаются в нижнее полупространство,
вновь распространяются в исследуемой среде и вновь
возвращаются к свободной поверхности, на этот раз в
виде кратных отражений. Этот принцип иллюстрирует-
ся на рис. 1, где условно показан ход луча однократной
и кратной волн, при этом источник колебаний распо-
ложен в точке й, сейсмоприемники - в точках <7, Ь и др.
Коэффициенты отражения от глубинных границ обо-
значены как /-|, г2, /3. Волна, претерпевшая один акт
отражения от какой-либо глубинной границы, имеет
Рис. 1. Ход луча отраженной волны:
луч, распространяющийся по траектории a-q (указаны ко-
ординаты выхода луча волны на свободную поверхность),
соответствует однократно-отраженной волне, а луч a-q-b -
кратной волне
амплитуду, пропорциональную rf, г2 или г3, и является однократно-отраженной (на рис. 1 ей соответствует ход
луча из а в q). Волна, претерпевшая более одного акта отражения от глубинных границ имеет амплитуду, пропор-
циональную старшим степеням коэффициентов отражения, и является кратной (на рис. 1 ей соответствует ход луча
из а в Ь). Глубинные границы, порождающие кратные волны, будем называть кратнообразуюшими горизонтами.
Ниже приводится классификация кратных волн по их динамическим (порядок кратности) и кинематическим
(виды кратных волн) параметрам.
ПОРЯДОК КРАТНОСТИ
Кратные волны первого порядка кратности претерпели два отражения от глубинных горизонтов, следовательно, их
амплитуда пропорциональна второй степени коэффициента отражения, например, rt2, г22, г^2, г^ и т. д. Волны
второго порядка кратности претерпели три отражения от глубинных горизонтов, их амплитуда пропорциональна
третей степени коэффициента отражения, например, г2, г^, Г|Г2г3, г\г2 11 т- д- Аналогично определяются волны более
высоких порядков.
Понятно, что наиболее существенными, конечно же, являются кратные первого порядка, так как в реальных
ситуациях всегда выполняется условие г, « 1, r2 << 1, г2 << 1. С другой стороны, интенсивность полезных отраже-
ний может столь сильно ослабевать с ростом времени прихода, что их могут маскировать волны старшего порядка
кратности.
ВИДЫ КРАТНЫХ ВОЛН
Все кратные волны можно разделить на три вида: полнократные, частично-кратные и внутренне-кратные. Для соот-
ветствующих видов волн в англоязычной литературе применяются термины - multiple, peg-leg, internal multiple. Ис-
пользуя трехслойную модель среды, представленную на рис. 1, дадим пример хода луча каждой из таких волн.
Полнократные волны соответствуют отражениям луча от свободной поверхности и только от одной кратнообра-
зующей границы. На рис. 2 показана полнократная волна первого порядка кратности. Как можно легко видеть, ее
амплитуда будет пропорциональна -г^, так как она имеет два отражения вверх от глубинной границы и одно отра-
жение вниз от свободной поверхности.
Частично-кратным волнам соответствует ход луча, имеющий как минимум два отражения от глубинных горизон-
тов, причем разных (например. и г3, - пример такой волны показан на рис. 1), и как минимум одно отражение в
нижнее полупространство от свободной поверхности.
Лучевая схема внутренне-кратных волн не имеет точек выхода луча на свободную поверхность, а переотражения
происходят между глубинными горизонтами (рис. 3).
Так как полнократные и частично-кратные волны связаны с переотраженисм в нижнее полупространство от
свободной поверхности, имеющей коэффициент отражения, приблизительно равный -1, то они являются наиболее
интенсивными. При этом внутренне-кратные волны обычно малоинтенсивны, так как они имеют как минимум три
отражения от кратнообразующих границ, коэффициенты отражения от которых, как правило, очень малы.
На рис. 4 показан разрез равных удалений, полученный выборкой трасс из данных, зарегистрированных при
морской сейсморазведке. Здесь четко выделяются основные виды волн, с которыми приходится иметь дело на
практике. Интенсивное отражение, обозначенное буквой А, - однократное отражение от морского дна, В - первая
полнократная волна, связанная с переотражением от морского дна, D - вторая полнократная волна. Укажем также
о	q	b
г2
г3 -----------------------------
Рис. 2. Ход луча a-q-b полнократной волны
первого порядка кратности:
волна имеет промежуточную точку, q, выхода луча
на свободную поверхность
Рис. 3. Ход луча внутренне-кратной волны:
лучевая схема не имеет промежуточной точки выхода луча
на свободную поверхность
6
Рис. 4. Разрез равных удалений (300 м), полученный
при морской сейсморазведке на Черном море
402
Рис. 5. Фрагмент исходной сейсмограммы ОГТ
на однократное отражение от глубокого горизон-
та - С, интерферирующее с первой полнократ-
ной волной В. Е - частично-кратная волна, пере-
отразившаяся в нижнее полупространство от сво-
бодной поверхности и дважды вверх: от дна моря
и от глубокого горизонта, которому соответству-
ет однократное отражение С.
Сейсмограмма ОГТ, полученная в области ин-
терференции однократной волны С и первой пол-
нократной В, представлена на рис. 5. На сейсмо-
грамме выделяется интенсивный цуг кратных волн,
обладающих малой скоростью Еогт, а также вы-
сокоскоростная однократная волна, выходящая из
области интерференции в области больших уда-
лений. Именно такая сейсмограмма была выбрана
неслучайно. При применении различных проце-
дур ослабления кратных волн удобно контроли-
ровать качество получаемого результата по крите-
рию сохранения динамики полезного отражения
после процедуры вычитания: амплитуда волны,
выделяемая из интерференции с помехой, долж-
на соответствовать динамике, прогнозируемой по
ее доступной асимптоте.
На рис. 6 дан другой пример волнового поля,
представляющего собой интерференцию однократ-
ных и кратных отражений. Материал был получен
приморской сейсморазведке. Изображенный фраг-
мент поля содержит однократно-отраженную от
морского дна волну (указана стрелкой А), первую
полнократную волну В, связанную с переотраже-
нием от морского дна, интенсивную однократ-
ную волну С, отраженную от неглубокого гори-
Рис. 6. Разрез равных удалений (200 м), полученный в
результате инженерных сейсмических наблюдений в
прибрежной мелководной части Каспийского моря
зонта, внутренне-кратную волну D, имеющую дополнительный пробег в слое, заключенном между морским дном и
упомянутым горизонтом, два цуга частично-кратных волн Е и F - соответственно первого и второго порядков
кратности, связанные с переотражениями от морского дна и горизонта С. Выделяется также интенсивное однократ-
ное отражение G и другие волны.
КЛАССИФИКАЦИЯ МЕТОДОВ ПОДАВЛЕНИЯ КРАТНЫХ ВОЛН
Все известные на сегодняшний день подходы к решению задачи подавления кратных можно условно разделить на
две большие группы (рис. 7): те, которые основаны на кинематическом различии полезных волн и помех (суммиро-
вание по ОГТ, веерная фильтрация, линейное и нелинейное преобразование Радона и т. д.) и те, которые явно или
неявно используют процедуру прямого продолжения волнового поля с целью прогнозирования кратных волн по
исходным данным (при этом прогнозирование может осуществляться как с учетом имеющейся оценки глубинно-
скоростнои модели среды, так и в условиях априорной неопределенности). Как бы промежуточное положение зани-
мает фильтр-маска: хотя он обнаруживает и вычитает кратные волны по кинематическому критерию, но для на-
стройки самого оператора преобразования поля требуется прогнозирование модели кратных волн.
Волновое поле, полученное в сложнопостроенной среде, целесообразно подвергать нескольким процедурам,
которые ослабляют фон кратных отражений, основываясь на различных принципах.
Наш опыт работы с достаточным количеством реальных данных показывает, что вначале всегда полезно приме-
нить двухшаговые алгоритмы (в рамках которых производится прогнозирование поля кратных волн по зарегистриро-
ванным данным с последующим его адаптивным вычитанием). При этом после адаптивного вычитания, как прави-
ло, фон кратных волн подавляется не до конца.
Остаточную энергию регулярных помех удобно ослабить уже кинематическими методами (пространственными
фильтрами, выделяющими и подчеркивающими однократные волны на фоне регулярных и нерегулярных помех). В
условиях нерегулярной геометрии 20-наблюдений или при большом расстоянии между соседними профилями в
ЗО-сейсморазведке необходимо особенное внимание уделять вопросам устойчивости используемых кинематических
фильтров, так как полезные однократные отражения, которые удалось выделить из их интерференции с кратными
волнами, могут быть вновь потеряны в их интерференции с артефактами многоканальных фильтров. Такие фильтры
оказываются эффективными в случаях, когда сигнал и помеха хорошо разделяются по своим кинематическим
параметрам. Сейсмограмма ОГТ, показанная на рис. 5, содержит цуг кратных отражений, интерферирующий с
полезной волной. Очевидно, что годографы сигнала и волн-помех имеют существенно различные производные в
значительном диапазоне удалений, хотя на ближних каналах они очень близки. Такие различия встречаются далеко
не всегда. Они свидетельствуют о большом контрасте пластовых скоростей в водном слое и подстилающих породах. В
подобных ситуациях критерием различения сигнала и помехи может являться производная годографа. Задачу разде-
Рис. 7. Условная классификация методов подавле-
ния кратных волн
8
ления волн удобно решать методами кинематической фильтрации. На основе априорной информации о кинемати-
ческих параметрах сигнала и регулярной помехи рассчитываются такие пространственные фильтры, которые выде-
ляют полезные отражения из их интерференции с шумом.
ДВУХШАГОВЫЕ МЕТОДЫ
История исследования проблемы ослабления кратных волн при помощи двухшаговых схем насчитывает уже не-
сколько десятилетий. В начале предлагались очень простые эвристические, чаще всего одноканальные, методы,
которые, хотя и находили свое применение, но показывали удовлетворительную работоспособность далеко не на
всех реальных данных. Позже они были улучшены, а также были разработаны новые, более совершенные и универ-
сальные подходы. Был выработан новый взгляд на их теоретические основы, а также были выявлены границы
области применимости методов, а именно, было показано, что самые простые подходы работают в рамках одномер-
ной ID-модели (распространение плоской волны в горизонтально-слоистой среде) сейсмического эксперимента,
более сложные - в полутаромерной 1,50-модели (точечный или линейный источник, возбуждающий упругие волны
в горизонтально-слоистой среде) и т. д. Поэтому в данном разделе, следуя за историческим ходом развития алгорит-
мической основы подходов к решению интересующей нас задачи, мы будем придерживаться логики анализа алго-
ритмов с единой точки зрения теории продолжения волновых полей, постепенно усложняя модель среды и последо-
вательно рассматривая случаи ID, 1,5D, 2D, 2,5D, 3D.
Здесь будет изучена ID-модель сейсмического эксперимента, т. е. скорее гипотетическая ситуация, но очень
полезная с точки зрения анализа работоспособности классических алгоритмов.
ID-МОДЕЛЬ СЕЙСМИЧЕСКОГО ЭКСПЕРИМЕНТА
Анализ ID-эксперимента мог бы быть совсем опущен, так как соответствующие алгоритмы уже не выдерживают
конкуренции с современными методами, ориентированными на обработку 2D- и ЗВ-материалов. Но все же мы
приведем краткое изложение особенностей одномерных подходов, так как они являются частными случаями много-
мерных, и на их примере удобно и наглядно анализировать некоторые динамические и кинематические аспекты
соответствующих преобразований. Будет рассмотрена одномерная сейсмическая задача, т. е. случай, когда регистри-
руемые колебания вызваны вертикальным падением плоской волны на систему плоско-параллельных слоев. Понят-
но, что в этом случае геометрическое расхождение волны отсутствует.
Предполагается, что свойства среды изменяются лишь с глубинной координатой г, но инвариантны относитель-
но х и у. Формой волны источника является неизвестный импульс s(t). Тогда зарегистрированное поле не будет
зависеть от координат л и у, и можно ограничиться рассмотрением лишь одной трассы.
ID-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН ПРИ НАЛИЧИИ МОДЕЛИ СРЕДЫ
Чаще всего прогнозирование поля кратных волн осуществляется следующим образом. На трассе p(t) выделяется
однократное отражение от интенсивного кратнообразующего горизонта, после чего в трассу' вводится сдвиг, равный
времени вступления г, этого сигнала, а результат рассматривается как модель кратных отражений l(t) = p(t - t}). Ее
преобразование Фурье будем обозначать через £(со). Так как в ID-случае первая полнокрагная волна имеет время
вступления 2/|, то, как легко видеть, в модели кратных
волн, полученной сдвигом на величину /(, ей будет соот-
ветствовать сигнал однократного отражения. Волнам стар-
шего порядка кратности в модели будут соответствовать
сигналы меньшего на единицу порядка. После этапа адап-
тации модель кратных волн вычитается из исходной трас-
сы, в чем и заключается традиционная реализация пред-
сказывающей деконволюции. Затем проводится интерпре-
тация результата и. возможно, выбирается ешс один крат-
нообразуюший горизонт. Тогда деконволюция повторяет-
ся с целью ослабления кратных волн, связанных с этим
горизонтом. Рассмотрим динамические особенности такой
процедуры последовательного вычитания кратных волн от
нескольких горизонтов.
Если коэффициент отражения кратнообразующего го-
ризонта равен /'|. то адаптация полнокрагных волн в трас-
се модели к соответствующим сигналам в исходной трас-
се осуществляется простым умножением первых на коэф-
фициент -/‘|. Но если, кроме полнокрагных в данных со-
Рис. 8. Лучевая схема частично-кратной волны со-
держит дополнительное переотражение со стороны
источника а и со стороны приемника b
9
держатся и частично-кратные волны, то такое соотношение амплитуд для адаптации не будет корректным. Рассмотрим
весь цуг волн, проходящих в среду ниже первой жесткой границы, отражающихся от нижерасположенных горизонтов
и возвращающихся к свободной поверхности. Каждому из таких отражений будут соответствовать две частично-крат-
ные волны, связанные с первой жесткой границей и имеющие отражения от нее со стороны источника и со стороны
приемника (рис. 8). Очевидно, что при сдвиге трассы на величину в полученной модели кратных волн этим двум
сигналам (в 1 D-случае времена их вступления совпадают) будет соответствовать лишь одна отраженная волна. Анало-
гичное несоответствие амплитуд будет наблюдаться и для волн старшего порядка кратности. Как уже отмечалось, в
исходной трассе присутствуют две частично-кратные волны первого порядка и, нетрудно видеть, три - второго поряд-
ка. После сдвига на /| первые две волны переместятся так, что время их вступления будет соответствовать последним
трем. Таким образом, в рамках ID-задачи не удается построить модель частично-кратных волн простым сдвигом на ф
Рассмотрим этот эффект на примере модели среды, состоящей из двух горизонтов и изображенной на рис. 9.
Пусть в момент времени г, зарегистрировано однократное отражение амплитуды Г| от первого горизонта и пусть
в момент времени Г2, t2 > г(, зарегистрировано отражение от второго горизонта. Его амплитуду обозначим г2 (где
учтено и влияние характеристики прохождения первого слоя). Легко показать, что в данном случае поле кратных
волн (вплоть до волн третьего порядка кратности включительно и без учета внутренне-кратных отражений) в частот-
ной области будет выглядеть следующим образом
£(со) = S(w)(-r,2e /<о2/| +	- rfe--’1"41'...) +	(1а)
+5(со)(-2г1г2е^'ш('1+'2) +3/-|2>2<?’-А,,(2'1+,2) -4^r2e-Ja(i,i+t2\..) +	(16)
+ 5(ш)(-~/-22е“^°2'2 +г23е~/и3'2 -г24е^ы4'2...) +
+5(со)(Зг22г1е_/'м(2'2+'|) - 4г23г)е“7Со(3,2+,|)...),
(1в)
(1г)
где .S(ctj) - спектр импульса s(f). В этом выражении, для наглядности, волны сгруппированы в строки, соответствую-
щие (1а) - полнократным волнам, связанным с первым горизонтом; (16) - частично-кратным волнам, связанным
с переотражением от первого горизонта; (1в) - полнократным волнам, связанным со вторым горизонтом; (1г) -
частично-кратным волнам, связанным с переотражением от второго горизонта.
На рис. 9 представлены лучевые схемы наиболее интенсивных отражений, входящих в выражение (1), - кратных
волн первого и второго порядков кратности. Показано, что существует одна полнократная волна первого порядка
а ь
гг-----
d
кратности и две частично-кратные волны (после-
дние формируют первый член строки (16)). В ситу-
ации 1D изображенные на рисунке частично-крат-
ные волны совпадают как по интенсивности (амп-
литуда каждой из них пропорциональна произве-
дению коэффициентов отражения Г[Т2), так и по
кинематике (время вступления каждой из них рав-
но Г| + т2). В случае 2D или 3D это свойство, вооб-
ще говоря, перестает выполняться. Использован-
ная при описании рисунка классификация (см. под-
рисуночную подпись): "‘переотражение со сторо-
ны приемника” и “переотражение со стороны ис-
точника” условна, так как одну и ту же волну мож-
но рассматривать с разных точек зрения. Напри-
мер, лучевая схема, показанная на рис. 9, Ь. соот-
ветствует однократной волне, отраженной от гра-
ницы с коэффициентом г2 и переотраженной со
стороны приемника h от границы с коэффициен-
том Г|. Она же может интерпретироваться как од-
Рис. 9. Лучевые схемы кратных волн в двухслойной
горизонтально-слоистой среде;
а - полнократная волна первого прядка кратности; b - ча-
стично-кратная волна первого порядка кратности, переот-
ражение со стороны приемника; с - частично-кратная волна
первого порядка кратности, переотражение со стороны ис-
точника; d - полнократная волна второго прядка кратнос-
ти; е - частично-кратная волна второго порядка кратности,
переотражение со стороны приемника; f - частично-крат-
ная волна второго порядка кратности, переотражения со
стороны источника и приемника; д - частично-кратная вол-
на второго порядка кратности, переотражение со стороны
источника
нократная волна, связанная с отражением от го-
ризонта с коэффициентом /у, и переотраженная
со стороны источника а от границы с коэффици-
ентом г2. Рис. 9 также иллюстрирует соотношение
амплитуд для волн второго порядка кратности: су-
ществует одна полнократная волна, отраженная от
границы с коэффициентом rt, и три связанных с
этой границей частично-кратных волны (последние
формируют второй член строки (16)).
Выберем первый горизонт как кратнообразую-
щую границу и получим модель кратных волн пу-
тем сдвига исходной трассы на величину г(:
10
5(oj)(/-1e'-ZOj2Zl -/12e“7°j3zi +z-3e^°j4,i -r14e“-/m5?l...) +
+5(ы)(-2/-1/-2е"/<0(2/| +Г2)+Зг|2г2е~/ю(3,|+/2)-4г|3г2е'-/“(4'|+,2)...) +
+5(w)(r2e~'т (rt+'^^r2e-A<1(t]+2l2)+r3e-j>O(t^3ti) _ r4e,-/<»(ri+4/2). ) +
+5(<o)(Зr1r22e'/“(2^l+2^2, -4r1r23e“/^0(2Z|+ЗZ2,...).
(2a)
(26)
(2b)
(2r)
На практике коэффициенты отражений от кратнообразующих горизонтов неизвестны, следовательно, нет воз-
можности учесть изменение амплитуды волны за счет отражения. Поэтому предсказанное поле кратных волн будет
получено с точностью до этого коэффициента, а в более сложных ситуациях - с точностью до характеристики
отражения от тонкослоистой пачки, рассматриваемой как одна отражающая граница. Таким образом, приходится
прибегать к процедуре адаптации модели кратных волн к исходному полю. Например, из (1а) и (2а) следует, что
модель полнократных волн переходит в искомую последовательность после простого умножения на коэффициент -rf.
В общем случае производится оценивание оптимального формирующего винеровского фильтра, “переводящего”
предсказанное поле кратных волн в исходную трассу. Здесь чаще всего используется критерий наименьших квадратов
отклонений, т. е. МНК (хотя, иногда, применяются и другие функции отклонений).
Попробуем адаптировать модель (2) к истинной последовательности (1). Очевидно, что удовлетворявшее нас
решение (умножение на коэффициент -г() для адаптации полнократных отражений теперь не преобразует модель в
поле реально зарегистрированных кратных волн, или, в более общей постановке задачи адаптации, не существует
такого фильтра, имеющего смысл отражательной характеристики тонкослоистого горизонта, который, будучи при-
менен к модели (2), обеспечил бы ее динамическое соответствие (1).
В этой ситуации можно предложить следующий подход к “улучшению” модели. Построим трассу, в которую будут
входить волны кратности два и выше. Такая трасса может быть получена дополнительным сдвигом (2) на величину /):
S(w)(/-1e“/'“3Z| - г12е'/ю4/| +	...) +	(За)
+S(co)(-2/-1r2e--/co(3zi+,2) +3r|2r2e''"f4/i+/’)...) +
+S(oj)(/-2c'/'"(2'i“/2> - r2e-j^2i\+2t2) _ } +
+5(to)(3r1r22e’/,,,(3/|+2Z2) -4rlr23e^'t’,(3f|+3'2)...).
(36)
(Зв)
(Зг)
Кроме того, выделим в модели (2) сигнал, соответствующий полнократному отражению первого порядка от
кратнообразующего горизонта, умножим его на 0,5 и получим
5(со)(0,5г|С >2z' -г12е--/“3'| +п3е->4,1 -/-4e'/<o5Z|...) +	(4а)
+5(w)(-2r|r2e“/fc,(2'|i/2) +3/'|2г2е-7’со(3'|+'2) -4rl3/-2e^‘”(4z'+'2)...) +	(46)
+5(w)(/-2e~/<D(Z|+'2) - г22е-/ш('1+2/2> + r3e->OJ<zi+3z2l _ r4e’>(zi +4z2>.„) +	(4в)
+S(w)(3r]r22e?ZM(2Z|4-2Z2) -4/-|r2ie^<u(2zl+3z2)...).	(4r)
Теперь уже можно подобрать такие операторы, после применения которых к выражениям (3) и (4) будет достиг-
нуто динамически правильное соотношение амплитуд модели и истинных волн. В нашем случае трассу в виде выра-
жения (4) следует умножить на коэффициент -2г|. трассу в виде выражения (3) на г2 и сложить результаты. После
вычитания этой суммы из исходной трассы получим следующий остаток
5((l>)(-r22e’7“2z2 +/-2 e“/<s3r2 - rfe уга4'2...) +
+.S(W)(rlr22e z<"'i ,2'2) - 2/|Г23е‘
что соответствует кратным волнам, связанным уже со вторым кратиообразующим горизонтом. Их подавление можно
производить аналогичным образом, т. е. при помощи описанной процедуры.
Если в распоряжении имеется информация о времени вступления однократного отражения от второго горизонта,
то ситуация значительно упрощается. А именно, если аналогичным образом получить и поле кратных волн, связан-
ных с отражением от второго горизонта
11
(5а)
(56)
(5в)
(5г)
5(оз)(г1е'/ш(,|+'2) -r|2e“/Q)(2Z1+Z2) + /-13е“7“(ЗГ|+/2) -/-4е~7“(4/1+'2)...) +
+5(оз)(-2г1Г2е^(о(г'+2,2) +ЗГ12г2е'/'“(2'|+2/2) -4г|3г2е^'(о(3'|+2?2)...) +
+..8(щ)(/'2е~7<',2'2 -/^e 7' -3'2 + ^е 7'04'2 -r24e-7w5f2...)-F
+5(со)(Зг1г22е"7“(,|+3'2) -4г]Г23е 7<о<,|+4'2)...).
то, как легко видеть, имея трассы (2) и (5), можно получить правильное соотношение амплитуд кратных волн
умножением (2) на коэффициент -Гр а (5) на -г2 (на практике подбор этих множителей осуществляется адаптацией).
Наконец, если после этого кратные волны вычесть из исходной записи, то получится трасса, состоящая только из
однократных отражений
5(го)(г1е“7“?| + r2c7'N'2).
Таким образом, вычитание полного набора кратных волн из исходного поля не сопряжено с принципиальными
трудностями.
При выводе общего случая следует ввести в рассмотрение вместо коэффициентов отражения соответствующие
операторы. Тем самым можно показать, что все приведенные рассуждения справедливы для произвольного числа
границ. Соответствующие формулы выписаны здесь лишь с целью продемонстрировать динамические особенности
прогнозирования кратных волн и пути совершенствования традиционных алгоритмов. Эти рассуждения справедливы
и для тех более универсальных подходов, ориентированных на обработку 2D- и ЗО-данных, которые будут рассмот-
рены в дальнейшем.
Ю-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН
БЕЗ ПРИВЛЕЧЕНИЯ АПРИОРНОЙ ИНФОРМАЦИИ О СТРУКТУРЕ СРЕДЫ
Зарегистрированная сейсмическая трасса содержит как отраженную волну, распространявшуюся по лучевой траек-
тории a-q, так и волну q-b (см. рис. 8, при этом в ID-случае координаты точек a, q, b совпадают). Тогда самый
непосредственный способ получить кратную волну a-q-b по зарегистрированным волнам a-q и q-b - это “сложить”
их. Такое “сложение” осуществляется автосверткой (иначе называемой ретрокорреляиия) трассы p(t) *p(f). при этом
не требуется применения процедур выделения волн a-q и q-b.
Иначе говоря, если волна подошла к поверхности и отразилась от нее в нижнее полупространство, то трассу,
зарегистрированную в точке отражения, можно рассматривать как функцию вторичного источника (принцип Тюй-
генса). Тогда, если запись в точке приема была реакцией на точечный импульсный источник, то, чтобы полудить
кратную волну, надо свернуть трассу с этой функцией источника. Отсюда приходим к алгоритму, основанному на
автосвертках трасс исходного волнового поля.
Для построения поля кратных волн /-го порядка кратности следует произвести /-кратную автосвертку рассматри-
ваемой нами трассы. Нетрудно показать, что в автосвертку первого порядка будут входить сигналы, представленные
выражением
52(ы)(г12е-7“2'1 - 2л-13е-Ло3,г’ +Зг14е>“7“4/1...) +
+52(td)(2ote'7“(Z1+f2) -6//2Z/e“7(0(2,|+r2) + 12c3te“7W(3'l+,2)...) +
+,8’2(ш)(/'22<?7<',2/2 -2/'23t>7<l3/2 + 3/'24<?7<i4/2...)t
+ 52(co)(-6r1/-22e"7“(Z|+2Z2) + 12/y23e~7“(Z|+3?2)...).
(6a)
(66)
(6b)
(6г)
Как видно из (6), полученные амплитуды и времена вступления волн первого порядка кратности точно (вернее,
с точностью до знака) соответствуют истинным (см. (1)), но такое соответствие не соблюдается для волн старшего
порядка кратности. Поэтому для получения более точного соотношения амплитуд производится дальнейший расчет
автосверток второго и, возможно, более высоких порядков.
Кроме того, отличие поля (6) от (1) заключается еще и в том, что после ретрокорреляции трассы вместо
импульса s(r) фигурирует его автосвертка	(спектральный эквивалент - 52(к>)). В силу того, что форма
импульса неизвестна, до этапа вычитания следует производить адаптацию. Из сопоставления (6) и (1) становится
понятным смысл оцениваемого оператора адаптации - это одноканальный стационарный фильтр .г'(f), такой, что
= 5(Г), где 6(/) - дельта-функция.
По аналогии с (6) можно получить и модели волн старшей кратности. Здесь мы их выписывать не будем, но еше
раз отметим, что в большинстве практически важных случаев можно ограничиться одним или, что бывает значи-
тельно реже, двумя наборами кратных волн. Дело в том, что коэффициенты отражения от реальных горизонтов
12
невелики, а поэтому волны третьего порядка кратности уже почти незаметны на фоне более интенсивных волн и
аддитивного шума. Как мы выяснили из сопоставления амплитуд волн первой кратности, содержащихся в модели
(6) и в исходной трассе (1). модель i = 1 обеспечивает их соответствие. При этом в модели / = 2, по построению, уже
не будут содержаться волны первого порядка кратности, а ее расчет позволит получить динамическое соответствие
и для волн второй кратности. Поэтому, если предположить, что волны третьей и более высокой кратности не вносят
заметного вклада в волновую картину, можно ограничиться использованием лишь двух моделей.
Наконец заметим, что, как и в предыдущем параграфе, если рассматривать вместо коэффициентов отражения
соответствующие операторы, то все выводы легко обобщаются на случай произвольного числа горизонтов.
ОСОБЕННОСТИ МЕТОДОВ АДАПТАЦИИ
Рассмотрев методы прогнозирования кратных волн, мы увидели, что этап адаптации является необходимым шагом
обработки, так как получаемая модель всегда отличается от зарегистрированных волн либо характеристикой отраже-
ния от кратнообразующего горизонта (предсказание кратных волн при наличии модели среды), либо тем, что в модель
дважды входит форма импульса источника (предсказание кратных волн без привлечения априорной информации о
структуре среды). В силу этого для построения динамически адекватного поля кратных волн всегда приходится прибе-
гать к процедуре адаптации, и здесь мы очень коротко остановимся на особенностях соответствующих алгоритмов.
Детальный разговор об особенностях адаптации мы отложим до рассмотрения случаев более сложных сред.
На данном этапе обработки осуществляется оценивание оптимального формирующего фильтра g(t), а физичес-
кий смысл этого фильтра будет различным при использовании различных методов расчета модели кратных волн.
Решается оптимизационная задача
g(t): [|/>(0 - l(t) * /(0|	min,
где * - символ свертки, а || || обозначает некоторую меру отклонения исходной трассы p(t) от результата адаптации
модели кратных волн £(/)*/(/). Как правило, используется квадратичная норма L2, т. е.
(/>(/)-g(r)W))2 min,
однако такой критерий не всегда приводит к надежным результатам.
Если прогнозирование кратных волн производилось в условиях априорной неопределенности относительно мо-
дели среды, то фильтр g(t) должен соответствовать оператору Иначе говоря, здесь мы имеем дело с операто-
ром обратной фильтрации, а следовательно, его оценивание может быть неустойчивым.
Если модель поля кратных волн была получена прогнозированием с использованием информации о строении
среды, то мы сталкиваемся с совсем другой задачей: требуется оценить либо коэффициент отражения от кратнооб-
разующего горизонта, либо оператор, соответствующий характеристике отражения от тонкослоистой пачки, т. е. нет
необходимости решать обратную задачу, а оцениваться будет короткий формирующий фильтр g(t), и сама процеду-
ра становится более устойчивой. Если, аналогично (2), модель кратных волн /(/) получена сдвигом исходной трассы
на величину l(t) = p(t - rjj, то Е2-адаптация производится при помощи минимизации функционала, соответству-
ющего дисперсии ошибки предсказания
£(0 : £(Р(0-g(0* р(г-Ц))2 -эти:.	(7)
t
Поставленная оптимизационная задача допускает аналитическое решение, основанное на теории оптимальных
фильтров Винера. Найденный таким образом оператор g(t) имеет смысл характеристики отражения от тонкослоис-
того кратнообразующего горизонта.
Функционал (7) имеет важную интерпретацию с точки зрения теории прогнозирования. Так как модель кратных
волн получена простым сдвигом исходной трассы Z(r) = pit - /|), а смысл выражения (7) заключается в том, что
требуется найти такой оператор, который бы, по критерию наименьших квадратов, обеспечил наилучшее соответ-
ствие самой трассы p(t) и /(/), то фактически здесь речь идет о предсказании (или, что то же самое, - прогнозиро-
вании) текущих значений трассы по предыдущим. Невязка e(i) = д(Г) - g(t)*p(t - Г0 называется ошибкой прогно-
зирования, а функционал, будучи суммой квадратов невязок, называется дисперсией ошибки прогнозирования.
Такой алгоритм получил название метод предсказывающей деконволюции, а способ оптимизации - минимизация дис-
персии ошибки предсказания.
КАКОЙ АЛГОРИТМ ПРОГНОЗИРОВАНИЯ КРАТНЫХ ВОЛН ПРИМЕНЯТЬ НА ПРАКТИКЕ?
Нами рассмотрены вопросы, связанные с прогнозированием цуга кратных волн, обусловленных переотражением в
нижнее полупространство от свободной поверхности. Хотя был изучен случай, соответствующий простейшей моде-
ли сейсмического эксперимента, динамические особенности процедуры расчета поля кратных волн окажутся спра-
ведливыми и для более общих постановок задачи.
13
Выбор между двумя алгоритмами прогнозирования кратных волн следует осуществлять отдельно в каждом конк-
ретном случае. Хотя, казалось бы, прогнозирование без модели среды является более привлекательным, но адаптация
к исходным трассам может быть проще и устойчивее, если в модели кратных отражений присутствует только один
набор волн, а не все отражения, связанные со свободной поверхностью. Геофизический смысл адаптивных фильтров
различен для моделей, получаемых различными методами. При прогнозировании без привлечения априорной инфор-
мации о структуре среды осуществляется расчет взаимных сверток трасс, что приводит к тому, что в результирующем
волновом поле импульс представлен в виде своей автосвертки. Таким образом, при адаптации модели к исходным
сейсмограммам фактически приходится оценивать обратный оператор на форму импульса, что негативно сказывается
на устойчивости решения. При прогнозировании с привлечением информации о строении среды нет возможности
учесть эффекты отражения от неизвестной тонкослоистой
пачки кратнообразующего горизонта. То есть при адапта-
ции будет оцениваться “прямой” оператор отражения. Сле-
довательно, в данном случае процесс расчета фильтра об-
ладает большей устойчивостью, а характеристика фильтра
становится более короткой.
Кроме того, в соответствии с динамическими особен-
ностями процедуры, не привлекающей информацию о мо-
дели среды, для корректного прогнозирования кратных
волн /-го порядка кратности необходимо получить / набо-
ров трасс (а в общем случае - сейсмограмм) кратных волн.
В реальных ситуациях, в особенности при обработке ма-
териалов 3D-сейсморазведки, вычислительные мощнос-
ти позволяют посчитать не более двух-трех наборов. Хотя
чаще всего этого бывает достаточно, но все же такое свой-
ство ограничивает возможности алгоритма. При этом по-
строение 7-го набора кратных волн основывается на ис-
пользовании (7 - 1)-го, т. е. соответствующая вычислитель-
ная схема представляет собой рекурсивную процедуру,
что ведет к быстрому накоплению погрешностей модели-
рования и шумов. В то же время при прогнозировании
поля кратных отражений! с использованием модели сре-
ды правильная динамика волн всех порядков достигается
расчетом всего двух моделей. Сразу скажем, что на прак-
тике мы не сталкивались с улучшением качества вычита-
ния при использовании более двух моделей. Поэтому дек-
ларируемое здесь преимущество приведено, скорее, для
полноты картины и является умозрительным.
Наглядное представление о динамических особеннос-
тях алгоритмов прогнозирования кратных волн могут дать
результаты обработки синтетических данных. Была полу-
чена сейсмическая трасса, соответствующая ID-наблю-
дениям в среде, состоящей из двух отражающих границ,
находящихся на глубинах 300 и 900 м. интервальные ско-
рости - 1500 и 2500 м/с. Были построены модели кратных
отражений при помощи обоих методов предсказания. По-
лученные результаты сопоставлены с исходной трассой
на рис. 10.
Очевидно, что прогнозирование с использованием
модели среды (были посчитаны волны, связанные с пере-
отражением от первого кратнообразуюшего горизонта, -
трасса С) обеспечило более точное по сравнению с про-
гнозированием без модели среды (трасса В) динамичес-
кое соответствие кратных волн реальным отражениям. Не-
избежные погрешности в амплитудах волн старшего по-
рядка кратности присутствуют как в трассе С, так и в
трассе В. Этот эффект обусловлен тем, что, как упомина-
лось выше, процедура прогнозирования без привлечения
модели среды является итерационной, и на каждом ее
шаге получаемый результат кинематически адекватен ре-
альным волнам, но динамически корректно прогнозиру-
ются волны до 7-го порядка кратности включительно, где
7 - номер итерации. Здесь была применена одна итерация,
Рис. 10. Сравнение методов расчета модели крат-
ных волн на синтетических данных:
результаты прогнозирования кратных волн при помощи
алгоритма, не привлекающего модель среды (В), с привле-
чением модели среды (С), исходная трасса (А); наиболее
интенсивные отражения: а - первая полнократная волна, от-
раженная от первого горизонта, время вступления - около
800 мс: b - однократная волна, отраженная от второго гори-
зонта, время вступления - около 880 мс; с - вторая полно-
кратная волна, отраженная от первого горизонта, время вступ-
ления - около 1200 мс; d - частично-кратная волна младшего
порядка кратности, отразилась по одному разу от первого и
от второго горизонтов, время вступления - около 1280 мс;
е - частично-кратная волна старшего порядка кратности,
имеет два отражения от первого горизонта и одно - от вто-
рого, время вступления - около 1680 мс; f - вторая полно-
кратная волна от второго горизонта, время вступления - око-
ло 1760 мс; однократная волна, отраженная от первого го-
ризонта, имеет время вступления около 400 мс, и на рисунке
не показана
14
т. е. следует ожидать, что корректно предсказаны волны первого порядка кратности. В данном случае это сигналы а
и f Как видно из рисунка, их амплитуды предсказаны с удовлетворительной точностью. С другой стороны, и одна
итерация алгоритма прогнозирования с учетом модели среды не обеспечивает динамически корректного результата.
Однако, как следует из рис. 10, несмотря на некоторую динамическую неточность, результат предсказания кратных
волн предпочтителен для их вычитания. Потенциальным недостатком такого способа является возможность расчета
цуга кратных волн, связанных с переотражением лишь от одного кратнообразующего горизонта. Для вычитания
остальных кратных волн расчеты приходится повторять. Так, например, трасса С не содержит первой полнократной
волны f отраженной от второго горизонта.
Если рассматривается задача прогнозирования с целью дальнейшего адаптивного вычитания всего цуга кратных
волн, связанных с переотражением в нижнее полупространство от свободной поверхности, по результатам морской
сейсморазведки, полученным в глубоководной части разреза, то во временном интервале регистрации кратные
волны старшего порядка кратности могут не вносить существенного вклада в волновое поле, и применение метода,
не привлекающего априорную информацию о модели среды, является предпочтительным. Кроме того, и на мелко-
водье двух наборов данных (что соответствует двум итерациям схемы прогнозирования) обычно бывает достаточно
для получения удовлетворительных результатов.
ЗАКЛЮЧЕНИЕ. В первой части статьи была рассмотрена одномерная сейсмическая задача, и было показано, что
нельзя получить модель кратных волн с правильным соотношением амплитуд сигналов простым сдвигом исходной
сейсмической трассы на величину, равную времени прихода однократно отраженной волны от выбранного кратно-
образующего горизонта (что соответствует традиционной реализации метода предсказывающей деконволюции в
задаче подавления ревербераций). Был описан альтернативный подход, основанный на привлечении дополнитель-
ной трассы, получаемой сдвигом на удвоенную величину. Кроме того, было показано, что поле волн первого поряд-
ка кратности с правильным соотношением амплитуд может быть получено автосверткой исходной трассы, а поле
волн старшего (/-го) порядка кратности - привлечением автосверток всех порядков кратности до /-го включительно.
Таким образом, теоретически обосновано использование обоих подходов к прогнозированию кратных волн. Мето-
дом автосверток можно получить все кратные волны, связанные с отражением в нижнее полупространство от сво-
бодной поверхности. Методом сдвига трассы можно получить кратные волны, связанные с отражением от свободной
поверхности и лишь от одного кратнообразующего горизонта.
При рассмотрении задачи 1D в явном виде информация о глубинно-скоростной модели среды не привлекалась,
а была осуществлена идентификация волн с целью выделения вступления однократных отражений и применения их
в качестве “сдвигающих” трассу операторов. При переходе к случаям 2D и 3D этот подход уже будет основываться на
расчетах с использованием оценки глубинно-скоростной модели среды, которая должна быть известна, по крайней
мере, в интервале глубин от свободной поверхности до кратнообразующего горизонта.
ПРИМЕЧАНИЯ
ВВЕДЕНИЕ, Одним из наиболее оригинальных способов
использования кратных волн с целью получения допол-
нительной информации о строении исследуемой среды
является “пересчет” кратных отражений в однократные.
Соответствующие исследования активно развиваются
группой исследователей из Дельфтского университета
(Нидерланды) под руководством проф. А. Дж. Беркхаута.
Предложенный ими алгоритм (см., например, [5|) осу-
ществляет прогнозирование кратных волн и их адапта-
цию к реально наблюденным кратным отражениям, со-
держащимся в исходных сейсмограммах. На следующем
этапе производится преобразование кратных волн в од-
нократные при помощи операторов миграционного типа,
полученных из анализа самих же исходных сейсмограмм.
Такая процедура позволяет получить информацию, от-
сутствующую в исходном волновом поле, но необходи-
мую для дальнейшей обработки. Рассмотрим типичную
схему возбуждения и регистрации волнового поля при
морской сейсморазведке (рис. 11). На косе, прикреплен-
ной к кораблю, расположены как источник колебаний
а, так и регистрирующие устройства Ь. При этом бли-
жайший к источнику приемник расположен, как пра-
вило, на значительном расстоянии. При формировании
сейсмограмм по результатам таких наблюдений прихо-
дится сталкиваться с эффектом отсутствия трасс ближ-
них удалений, что негативно сказывается на результатах
применения любых многоканальных процедур обработ-
ки, т. е. пространственных фильтров. В такой ситуации
обычно применяются различные процедуры интерполя-
ции, но все они имеют свои недостатки и ограничения.
При этом можно получить, например, сигнал однократ-
Рис. 11. Типичная схема возбуждения и регистра-
ции колебаний при морской сейсморазведке:
схематично показан корабль, буксирующий косу, на которой
расположены источник колебаний а и сейсмоприемники Ь;
для простоты изображена лишь одна глубинная граница -
дно моря
15
ного отражения от морского дна, который был бы заре-
гистрирован в точке с координатой q, преобразовав сиг-
нал первой полнократной волны, зарегистрированной
в точке с координатой Ь. Действительно, кратная волна
такого типа “содержит’’ в своей лучевой схеме искомую
однократную волну a-q, и применение специального опе-
ратора миграционного типа позволяет устранить из ее
лучевой схемы путы/-/?. Тем самым открывается возмож-
ность преобразования кратных волн в однократные, при
этом также реализуется и процедура экстраполяции вол-
нового поля на ближние удаления.
ВИДЫ КРАТНЫХ ВОЛН. При наземной сейсморазведке
проблема кратных волн зачастую не менее актуальна, чем
при морской сейсморазведке. Здесь показаны только мор-
ские данные, так как на таких разрезах и сейсмограммах
можно проводить более уверенную идентификацию от-
раженных волн, т. е. материалы морской сейсморазведки
предпочтительны для демонстрационных целей.
Отдельный важный вопрос, актуальный для данных,
полученных при наземной сейсморазведке, - это учет ста-
тических поправок и вообще влияния неоднородностей
верхней части разреза при прогнозировании кратных волн.
ДВУХШАГОВЫЕ МЕТОДЫ. Развитие подходов к пробле-
ме подавления кратных волн по принципу прогнозиро-
вания/вычитания шло по пути последовательного услож-
нения, отталкиваясь от алгоритмов, ориентированных
на простейшие модели сейсмического эксперимента. В
ранней работе Э. Робинсона [8] был предложен способ
подавления ревербераций по одной трассе - предсказы-
вающая деконволюция. Рассматривая этот алгоритм с
позиций теории продолжения волновых полей, можно
показать, что неявно используемый в методе этап пря-
мого продолжения волнового поля обеспечивает дина-
мически адекватные результаты только в случае 1D, т. е.
одномерной модели сейсмического эксперимента. По-
этому понятно, почему на практике алгоритм показы-
вал недостаточную эффективность и универсальность.
Если же иметь в виду, что при обработке приходится
иметь дело не только с трассами нулевого удаления, но
и с сейсмограммами, то следует учесть различные про-
странственные эффекты, что было предложено в работе
М. Тансра [9]. Соответствующий алгоритм применяется
в области разложения по плоским волнам (преобразо-
вание Радона).
Ю-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН ПРИ НАЛИЧИИ МО-
ДЕЛИ СРЕДЫ. На динамическое несоответствие модели
кратных отражений, полученной сдвигом исходной трас-
сы. реально зарегистрированным волнам указывалось,
в частности, в работе [7]. Метод построения динамичес-
ки адекватной модели кратных волн при помощи при-
влечения трассы с удвоенным сдвигом обсуждался в пуб-
ликации [1].
Ю-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН БЕЗ ПРИВЛЕЧЕНИЯ
АПРИОРНОЙ ИНФОРМАЦИИ О СТРУКТУРЕ СРЕДЫ. Наи-
более универсальным подходом к проблеме прогнози-
рования кратных волн является анализ задачи с пози-
ций теории продолжения волновых полей. Но при этом
он допускает более прозрачную интерпретацию в про-
стых случаях, например, если речь идет об одномерной!
модели сейсмического эксперимента, когда для прогно-
зирования кратных волн нет необходимости примене-
ния пространственных операторов, так как координаты
источника, промежуточной точки q выхода луча и при-
емника совпадают. Следовательно, процедура получения
кратных волн сводится к свертке трассы с самой же со-
бой или автосвертке - autoconvolution [4, 10], иначе на-
зываемой “ретрокорреляция” [2].
ОСОБЕННОСТИ МЕТОДОВ АДАПТАЦИИ. Авторы работы
[6] видят причину неустойчивости схем адаптации в ис-
пользовании квадратичного критерия L2, заменяя его
на норму L1 (метод наименьших модулей). С другой сто-
роны, как указывается в [3], такая замена влияет на
результат вычитания несущественно, но при этом зна-
чительно возрастают вычислительные затраты ввиду не-
обходимости решения системы нелинейных уравнений.
В последующих публикациях мы подробнее рассмотрим
эту проблему.
ЛИТЕРАТУРА
[.Денисов М. С., Дингман С. Л., Фиников Д. Б.. 2002. Экстрапо-
ляция волнового поля в задаче моделирования кратных волн
(с целью их подавления): Геофизика. 6.
2.	Мушин И. A., 1983, Конструирование алгоритмов и графов
обработки данных сейсморазведки: М., Недра.
3.	Abma R., Kabir N., Matson К., Michell S., Shaw S., VcLain B.,
2005. Comparison of adaptive subtraction methods for multiple
attenuation: The Leading Edge, 27, 3, 277 - 280.
4.	Berkhout A. J., 1999. Multiple removal based on the feedback
model: The Leading edge, 18, 127 - 131.
5.	Berkhout. A. J. and Verschuur D. J., 2006. Focal transformation,
an imaging concept for signal restoration and noise removal: Geo-
physics. 71. 6, A55 - A59.
B.Guitton A., Verschuur D.. 2004. Adaptive subtraction of multiples
using the LI norm: Geophysical Prospecting, 52, 1.
7.	Lokshtanov D., 1995, Multiple suppression by single channel
and multichannel deconvolution in the tau-p domain: 65th Ann.
Internal. Mtg. SEG.
8.	Robinson E. A., 1957. Predictive decomposition of seismic trace*
Geophysics, 22, 767 - 778.
9.	Taner M. T. 1980, Long-period sea-floor multiples and their!
suppression: Geophysical Prospecting, 28, 30 - 48.
10.	Tsai C. J.. 1985, Use of autoconvolution to suppress first-order,
long-period multiples: Geophysics, 50, 9, 1410 - 1425.
КОРОТКО ОБ АВТОРАХ
Михаил Сергеевич ДЕНИСОВ - ведущий математик ООО “Геотехсистем”, доктор физ.-мат. наук.
Дмитрий Борисович ФИНИКОВ - сотрудник ООО “Геотехсистем”, кандидат техн. наук.
16

M. С. Денисов
Д. Б. Фиников
ООО "ГЕОТЕХСИСТЕМ", МОСКВА
ООО "ГЕОТЕХСИСТЕМ", МОСКВА
МЕТОДЫ ПОДАВЛЕНИЯ КРАТНЫХ ВОЛН
В СЕЙСМОРАЗВЕДКЕ. Часть 2
АННОТАЦИЯ. В статье рассматриваются алгоритмы подавле-
ния кратных волн в рамках так называемой полуторамерной
модели сейсмического эксперимента. Анализируется волновое
поле вызванное линейным источником и распространяющее-
ся в среде, состоящей из плоско-параллельных границ. Иссле-
дуются алгоритмы, привлекающие априорную оценку глубин-
но-скоростной модели среды, и алгоритмы, работающие в ус-
ловиях отсутствия такой информации. Показано, что решение
проблемы подавления кратных волн может быть сведено к ре-
шению аналогичной одномерной задачи. Для этого требуется
разложить исходное волновое поле по базису плоских волн.
Приводятся результаты обработки модельных и реальных дан-
ных. В статье обобщаются и систематизируются известные и
применяемые на практике алгоритмы.
ABSTRACT. The paper considers the algorithms for multiple
attenuation in case of the so-called one and a half dimensional seismic
problem. A wavefield radiated by a point source in horizontally-layered
earth is analyzed. Studied are both the methods that utilize the a
priori macro model estimate and the data-only-driven ones. It is
shown that the solution of the multiple elimination problem might
be reduced to the one for the one-dimensional case. This requires
transformation of the input data into the plane-wave decomposition
domain. The results of synthetic and real data processing are given.
The paper generalizes the known and widely used in industrial
processing methods for multiple attenuation.
ВВЕДЕНИЕ. В первой части статьи1 была дана класси-
фикация известных на сегодняшний день методов по-
давления кратных волн. Также были предложены основ-
ные определения, необходимые для анализа природы
кратных отражений и разработки методов их ослабле-
ния. Было введено понятие двухшаговых алгоритмов, в
рамках которых осуществляются прогнозирование поля
кратных волн, его адаптация к реально зарегистриро-
ванным отражениям и последующее вычитание из ис-
ходных данных. При этом особенности двухшаговых ал-
горитмов рассматривались на примере простейшей мо-
дели среды, так называемой одномерной задачи - 1D,
т. е. был изучен случай вертикального падения плоской
волны на систему плоско-параллельных слоев, и отсут-
ствовало такое понятие, как удаление источник-приём-
ник. Однако выводы, которые были сделаны, обобща-
ются и на более сложные ситуации.
Во второй части статьи мы сделаем следующий шаг
по пути усложнения структуры сейсмического экспери-
мента и осуществим переход к так называемой полуто-
рамерной модели - 1,5D, которая подразумевает регис-
1 Денисов М. С., Фиников Д. Б., 2007, Методы подавления крат-
ных волн в сейсморазведке. Часть 1. Технологии сейсморазвед-
ки, № 1.
трацию волнового поля, вызванного точечным или ли-
нейным источником и распространяющегося в горизон-
тально-слоистой среде. Положение пункта регистрации
может отличаться от положения источника колебаний.
Здесь уже приходится учитывать такие параметры вол-
нового поля, как геометрическое расхождение, угол
выхода луча и др. Сейсмические трассы будут зависеть
от взаимного расположения источника колебаний и при-
ёмника, и, в отличие от ID-эксперимента, процесс
прямого продолжения волнового поля с целью прогно-
зирования кратных отражений не будет сводиться к про-
стому сдвигу зарегистрированной записи.
Несмотря на ограничения, накладываемые сделан-
ными предположениями, методы, которые будут рассмот-
рены ниже, часто применяются на практике. Основными
причинами их популярности являются наглядность, про-
стота реализации и вычислительная эффективность. Об-
ласть применимости таких алгоритмов можно несколь-
ко расширить, введя понятие “локальной 1,5D” зада-
чи, т. е. случая, когда среда устроена так, что наклоны и
кривизны отражающих горизонтов невелики на некото-
рой небольшой пространственной базе, сопоставимой
с длиной косы приёмников.
В дальнейшем мы будем широко использовать прин-
цип разложения поля упругих колебаний по базису плос-
ких волн вида <5(7 + ух), где 5() - дельта-функция, а
7
параметр у характеризует угол, образуемый фронтом
плоской волны и свободной поверхностью. Как след-
ствие этого, будем прибегать к рассмотрению частного
случая волнового поля в виде изолированной плоской
волны. Понятно, что с такой ситуацией мы никогда не
столкнемся на практике, но соответствующий анализ
алгоритмов и их “калибровка” по плоским волнам не
только добавляют ясности рассуждениям, но и делают
их более прозрачными и наглядными. С другой стороны,
так как любое волновое поле может быть представлено
в виде разложения по плоским волнам, то в силу ли-
нейности всех рассматриваемых операторов, изучив их
отклик на каждую плоскую волну, можно построить
реакцию на реально зарегистрированное поле. И тогда
такой анализ важен именно для практических нужд. Этот
способ аналогичен хорошо зарекомендовавшему себя в
теории линейных систем методу спектрального анализа
стационарных объектов, когда как сигнал входного воз-
действия, так и результат его преобразования системой
раскладываются по базису Фурье, состоящему из гар-
монических функций. Характеристика системы изучает-
ся отдельно для каждой гармоники, а общий результат
получается сложением результатов, полученных отдель-
но по каждой гармонике.
1,50-МОДЕЛЬ СЕЙСМИЧЕСКОГО ЭКСПЕРИМЕНТА
Как было оговорено выше, среда предполагается гори-
зонтально-слоистой. Расположим декартову систему ко-
ординат (х, у, z} так, чтобы ось z была направлена в
глубь среды и перпендикулярно к свободной поверхно-
сти. Тогда ввиду инвариантности строения среды отно-
сительно х и у можно выбрать произвольную ориента-
цию системы координат. Поэтому условимся считать,
что линия профиля сейсмических наблюдений всегда
имеет нулевую координату у, т. е. положение всех ис-
точников и приёмников описывается только коорди-
натой х, а если используется линейный источник, то он
ориентирован перпендикулярно к линии профиля.
Ситуации, с которыми приходится иметь дело при
анализе волновых полей, возбуждаемых точечным и ли-
нейным источниками, во многом похожи. В обоих слу-
чаях волновое поле имеет свойство симметрии (для то-
чечного источника - поле зависит только от удаления
источник-приёмник, для линейного источника - поле
не зависит от координаты у), которое учитывается в
соответствующих рассуждениях. При этом проще и на-
гляднее рассматривать ситуацию линейного источника,
тем более что для 1,5О-задачи разработана специальная
процедура динамически корректного преобразования
волнового поля точечного источника в поле линейного
источника. Поэтому здесь мы ограничимся анализом
именно волнового поля линейного источника.
1,50-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН
ПРИ НАЛИЧИИ МОДЕЛИ СРЕДЫ
В отличие от случая 1D, т. е. вертикального распро-
странения плоской волны, при котором для анализа всего
волнового поля достаточно было рассмотреть лишь одну
трассу, теперь приходится учитывать пространственные
эффекты, возникающие при распространении волн. В
рамках двухшаговых схем прогнозирование любой крат
ной волны производится по однократной волне (илипД
кратной меньшего порядка кратности). Такое преобра-
зование всегда осуществляется методами прямого про-
должения волнового поля, которые, вообще говоря, тре-
буют применения процедур пространственного накапли-
вания - операторов миграционного типа. Однако в рамках
ID-задачи анализ схем прогнозирования существенно
упрощался. Операторы миграционного типа вырождалиа
в операторы простого сдвига зарегистрированной трассь.
на величину, равную удвоенному времени распростране-
ния волны от свободной поверхности до выбранного
кратнообразующего горизонта, т. е. времени вступления
на зарегистрированной трассе сигнала однократного от-
ражения от кратнообразующего горизонта. Можно ска-
зать, что это и было равносильно привлечению априор-
ной информации о глубинно-скоростной модели среды.
Такое прогнозирование поля кратных отражений исполь-
зуется в методе предсказывающей деконволюции.
При изучении задачи 1,5D мы также будем прогно-
зировать кратные отражения по однократным волнам
(или по кратным меньшего порядка кратности). Но те-
перь пространственные операторы миграционного типа
уже не будут сводиться к одноканальным операторам
сдвига трассы, и нам придется учитывать параметр 6 -
угол выхода луча на свободную поверхность (напомним,
что луч соответствует нормали к фронту волны).
Действительно, если модель среды известна в интер-
вале глубин от свободной поверхности до выбранного
кратнообразующего горизонта, и, например, решается
задача прогнозирования кратных волн, связанных с до-
полнительным переотражением от верхнего горизонта,
то, как видно из рис. 1. Г, - время распространения вол-
ны в верхнем слое будет зависеть от угла 6, т. е. /| = /((0).
Рассмотрим этот эффект подробнее. Известно, что
на сейсмограмме общей глубинной точки (ОГТ) асимп-
тоты годографов регистрируемой отражённой волны и
порождаемой ею кратной волны смыкаются, т. е. раз-
ность времен вступления сигналов уменьшается с рос-
том удаления, или, что то же самое, с возрастанием 6
(рис. 2). Обратимся к рис. 3, а, где более детально пока-
зан ход луча в интересующем нас приповерхностном
слое, мощность которого /г, а соответствующая интер-
вальная скорость г. Пусть к свободной поверхности под-
ходит плоская волна. Путь вдоль луча, показанного на
рис. 3, а, т. е. подходящего к точке с координатой q,
Рис. 1. 1,5О-модель сейсмического эксперимента
8
Рис. 2. Фрагмент реальной сейсмограммы ОГТ, по-
лученной при морской сейсморазведке:
время вступления однократно отражённой от морского дна
волны на трассе ближнего удаления - около 1600 мс. За ней
следует цуг частично-кратных волн и полнократных волн.
Очевидно, что разность времён вступления отдельных сиг-
налов уменьшается с ростом удаления, т. е. годографы со-
ответствующих волн смыкаются
Как следует из рис. 3, а, при заданном угле 0 можно
получить следующее выражение для удаления: Дх = 2/г tgO.
Тогда окончательно приходим к
2й
Д/(0) = — cos 0.
V
Если для записи формулы использовать так называе-
мый лучевой параметр а, где а = sin0/v, и для наглядно-
сти изменить обозначения с Д/(0) на Д/(сс), то
Д/(«) = 2h ~ -а2
N V2
(1)
Нетрудно видеть, что выражение (1) соответствует
эллиптической траектории ДГ(а) как функции от а.
Если глубинный кратнообразующий горизонт не
верхний, то функция запаздывания кратной волны от-
носительно порождающей её волны выглядит сложнее,
но также зависит только от 0 и одинакова для всех волн,
подходящих под этим углом к свободной поверхности.
В самом деле, в наших рассуждениях непринципиальна
прямолинейность луча, соединяющего точку отраже-
ния q и точку/? выхода на свободную поверхность. Пусть
годограф однократно отражённой от кратнообразуюше-
го горизонта волны на сейсмограмме общего пункта
взрыва (ОПВ) с источником, расположенным, напри-
мер, в точке с нулевой координатой х = а = 0, описыва-
ется функцией g(x). Пусть х* такова, что
затем испытывающего переотражение в нижнее полу-
пространство от поверхности, вверх от кратнообразую-
щего горизонта и приходящего в точку регистрации с
координатой Ь, равен 1h/cosQ. Следовательно, дополни-
тельное время распространения кратной волны за счёт
переотражения
_ 1 2/т
*777	м ’
V cos 6
При этом в точке Ь также регистрируется и одно-
кратная волна. Обозначим время вступления волны,
порождающей кратное отражение и регистрируемой в точ-
ке q, через /(,. Запаздывание между моментами вступле-
ния кратной и однократной волн в точке с координатой
b определяется разностью времён /() + tm и /() + tp. где t -
время, требующееся для того, чтобы сигнал, принадле-
жащий рассматриваемой нами плоской волне, прошёл
расстояние, равное длине отрезка ОЬ в треугольнике qOb
(см. рис. 3, б). Очевидно, что Ob = Axsin6. Тогда
t „ = - Дх sin е,
v
Рис. 3. Вычисление запаздывания вступления крат-
ной волны относительно порождающего её однократ-
ного отражения:
а - расчёт времени распространения волны в приповерх-
ностном слое; б - расчёт относительного запаздывания
вступления кратной волны
а разность времен вычисляется как
или
Дг(б) = (г0 + г,„) - (Г0 + 9
дг(0) = 1
V
-------Дх sin 0
COS0
9
Если использовать понятие обратной функции g ()
к производной годографа g'(), то ** = g'-l(a). Из про-
зрачных геометрических соображений находим
A/(a) = g(g'_|(a)) -ag'-'(a).
(I)
Здесь предполагается, что g'_|(a) существует и что
для 1,5Э-задачи всегда справедливо. Очевидно, что крат-
ные волны старшего порядка кратности имеют те же
запаздывания относительно породивших их волн млад-
шего порядка, т. е. ситуация аналогична случаю, рас-
смотренному при анализе ID-модели сейсмического
эксперимента.
Этот принцип лег в основу алгоритма предсказыва-
ющей деконволюции, но реализуемой не по исходным
трассам, а по результату разделения волн в соответствии
с углом, образуемым фронтом волны и свободной по-
верхностью.
Пусть зарегистрировано волновое поле, вызванное
воздействием на среду источника, расположенного в
точке х = а = 0. Обозначим подборку трасс в виде сейс-
мограммы ОПВ как р(х, /). Предположим, чго нам уда-
лось представить это волновое поле в виде
p(x,/) = [ p(a, t + ax)da.
(2)
Здесь мы не конкретизируем способ получения тако-
го разложения по плоским волнам, описываемым функ-
циями р(а, /). При весьма необременительных условиях
(достаточно существования двумерного преобразования
Фурье от р(х. ()) можно показать, что разложение (2)
существует и единственно. В реальности мы можем осу-
ществить разложение по конечному набору а дискретно-
го по х и t пространственно временного поля, зарегист-
рированного на конечном пох и / объеме выборки. Здесь
надо было бы обсуждать множество деталей реализации
разложения (2), так как в этих условиях оно осуществ-
ляется неточно, и актуальной становится задача апп-
роксимации фрагмента исходного поля набором плос-
ких волн наилучшим в некотором смысле образом. Это
обширная тема, которой мы касаться нс будем, а отме-
тим лишь, что обычно для получения приблизительно-
го выполнения (2) используют направленное суммиро-
вание
p(a. t) ~ r(t) * У p(x. t - ax) ,
(3)
где х - дискретный индекс, т. е. номер трассы в сейсмо-
грамме, а интервал между трассами условно положен
равным единице Дх = 1. Символ * обозначает процедуру
свёртки по переменной /. Оператор фильтра /•(/) имеет
преобразование Фурье /?(со) = | w|/2л. При обработке сейс-
мических данных приходится вместо интегрального пре-
образования (2) применять суммирование с шагом Да.
при этом шаг выбирается подробным, а диапазон изме-
нения а достаточно широким.
Точным выражением, описывающим перевод сейс-
мограммы из области (х - г) в (а - г), является интеграль-
ный аналог формулы (3), рассчитанный на непрерыв-
ные по пространственной координате данные
p(a,t) = r(t)* [ p(x, t - ax)dx .
(3'1
Преобразования (2) и (3) носят название соответ-
ственно обратного и прямого преобразований Радона. Но.
сделав надлежащие оговорки, в настоящей работе это же
название закрепим и за парой выражений (2) и (3).
Пренебрегая пока неточностями, обусловленными
дискретностью данных, заметим, что если разложение
(2) реализовано, то задача прогнозирования и вычита-
ния кратных волн сводится к одномерной задаче, рас-
смотренной нами в первой части статьи, которая может
быть решена методом предсказывающей деконволюции.
Её аналог, применяемый к волновому полюр(а, /), вы-
глядит следующим образом. Пусть в качестве кратнооб-
разуюшего выбран верхний горизонт, расположенный
на глубине h, и требуется спрогнозировать кратные вол-
ны, связанные с переотражением от него. Сигнал соот-
ветствующего однократного отражения на сейсмограм-
ме ОГТ имеет годограф
/(х) = - J(2/!)2 +х2.
v
Как следует из (3), прямое преобразование Радона
представляет собой двухшаговую процедуру, на первом
шаге которой осуществляется ввод в сейсмограмму ли-
нейных кинематических поправоках, т. е. годограф пре-
образуется в
г(х) = г(х) - ах = —
v
- ax.
На втором шаге производится суммирование трасс
по пространственной координате х. Понятно, что при
этом сигнал однократного отражения будет синфазно
накоплен в окрестности точки х*, определяемой урав-
нением
ЭЦх)
дх
= 0.
Из этого условия легко вычислить значение искох
координаты
а
Г1 7
а"
v
x* = 2h
Подставив полученное значение в выражение для /(х)
приходим к
т(а) = 2й1— - а2,
V г“
где т(а) = 7(х*) - ах* - значение временной координа
ты, на которой будет накоплен сигнал. Таким образом
после применения преобразования Радона гиперболи
ческий годограф волны преобразуется в эллиптическую
траекторию т(а), что полностью согласуется с (1). Следо
вательно, прогнозирование кратных отражений для пред
сказывающей деконволюции, реализуемой в области раз
ложения Радона, производится сдвигом на величину т(а)
определяемую образом годографа однократно отражён
10
ной от кратнообразуюшего горизонта волны: получим
сейсмограмму 1\(а, г), /((а, t) = р(а, t - т(а)), которую в
дальнейшем будем адаптировать к полю Да, t).
Точно так же легко убедиться, что для любого друго-
го кратнообразующего горизонта
т(а) = ^'"'(а))	cig'’1 (а),
и это выражение согласуется с формулой (Г).
В соответствии с принципами, изложенными в пер-
вой части статьи, для получения динамически коррект-
ного результата прогнозирования следует также получить
дополнительное волновое поле кратных отражении стар-
шего порядка кратности /фа, г). Для этого к сейсмо-
грамме /|(а, 1) применим фазовые сдвиги т(а): /2(а, /) —
7|(а, г-т(а)). Кроме того, умножим сигнал, являющийся
образом первой полнократной волны на сейсмограмме
7|(а, t), на 0,5.
Затем осуществляется адаптация полученных двух
моделей кратных отражений к исходному волновому
полю. Этот этап может быть реализован также и после
пересчета моделей обратно в пространственно-времен-
ную область. Преимуществом такой адаптации и вычи-
тания будет высокое быстродействие, так как, как пра-
вило, в соответствии с требованиями однозначности
прямого и обратного преобразования Радона требуется
получить намного больше трасс-элементов разложения.
чем трасс, содержащихся в исходной сейс-
мограмме. С другой стороны, иногда при-
ходится применять вычитание в области
преобразования Радона, так как зависи-
мость отражающей характеристики кратно-
образующего горизонта от угла подхода луча
может играть существенную роль. Простей-
шим проявлением этой зависимости явля-
ется AVO-эффект. Кроме того, если глубин-
ный горизонт представляет собой тонкосло-
истую пачку, то оператор, описывающий
эффект отражения от неё, также зависит
от угла падения-отражения луча, так как
лучевая траектория распространения в си-
стеме тонких слоев будет различной для
различных углов.
Результат адаптивного вычитания крат-
ных отражений преобразуется обратно в
пространственно-временную область в со-
ответствии с формулой (2).
Проиллюстрируем описанные шаги на
примере обработки модельной сейсмограм-
мы ОГТ, показанной на рис. 4. Волновое
поле получено в двухслойной горизонталь-
но-слоистой среде и содержит два однократ-
ных отражения и полный цуг кратных волн
всех возможных типов. Времена вступления
однократных отражений на трассе ближне-
го удаления около 550 и 850 мс.
На первом этапе обработки применим
к сейсмограмме преобразование Радона в
соответствии с выражением (3). Его резуль-
тат, Да, т), показан на рис. 5, а. По гори-
зонтали отложены значения параметра а.
Первая (крайняя слева) трасса сейсмограм-
мы соответствует значению а = 0, г. е. плоской волне,
падающей вертикально. Последняя (крайняя справа)
трасса соответствует максимальному значению парамет-
ра а, которое реализуется на исходной сейсмограмме.
Это значение оценивается по асимптоте самого круто-
наклоненного годографа. На рис. 5, б изображён резуль-
тат прогнозирования кратных волн, связанных с пере-
отражением от верхнего горизонта и от свободной по-
верхности. Кривая фазовых сдвигов т(а) на сейсмограм-
ме в области разложения Радона совпадает с образом
однократно отражённой от кратнообразующего горизонта
волны. Как следует из рис. 5, о, все кратные отражения
заданного типа спрогнозированы кинематически кор-
ректно. Результат адаптивного вычитания в области пре-
образования Радона представлен на рис. 5, в. Выделяется
неослабленная кратная волна, время вступления кото-
рой на трассе а = 0 около 1700 мс. Это отражение не
было получено в результате прогнозирования, так как
является полнократной волной первого порядка крат-
ности, переотражённой от второго горизонта. Для ос-
лабления цуга кратных волн, связанных со вторым го-
ризоитом, следует повторить этапы прогнозирования и
вычитания. При этом фазовые сдвиги должны опреде-
ляться параметрами второго горизонта и вычисляться
по формуле (Г), иначе говоря - совпадать с образом
однократно отражённой от второго горизонта волны в
области преобразования Радона. Некоторые динамичес-
Рис. 4. Модельная сейсмограмма ОГТ, содержащая
два однократных отражения и весь набор полнократ-
ных, частично-кратных и внутренне-кратных волн:
а - исходная сейсмограмма; б - сейсмограмма, полученная
путем преобразования в пространственно-временную об-
ласть результата адаптивного вычитания
11
кие аспекты процедуры вычитания кратных волн от не-
скольких горизонтов подробно рассмотрены в первой
части статьи.
Наконец, на рис 4, б показан результат подавления
кратных волн, связанных только с первым горизонтом,
после его пересчета обратно в пространственно-времен-
ную область по формуле (2). Очевидно, что удалось не-
плохо ослабить энергию кратных волн в диапазоне ма-
лых и средних удалений, при этом на больших удалени-
ях заметны асимптоты годографов кратных отражений.
Причины, приводящие к такому эффекту, мы проана-
лизируем ниже. Здесь же отметим, что, как правило,
подавление кратных волн в области разложения Радона
сопровождается также применением процедуры мьютин-
га или обнуления некоторых фрагментов сейсмограммы
р(а, т) после её интерпретации с целью обнаружения
волн-сигналов и помех. Как правило, волны-сигналы и
регулярные помехи лучше локализуются и лучше разде-
ляются в пространстве преобразования Радона, чем на
исходных сейсмограч,мах. Такой приём обычно оказыва-
ется наиболее эффективным при ослаблении регуляр-
ной помехи на трассах больших удалений, а его резуль-
таты близки к результатам применения веерных режек-
торных фильтров. С другой стороны, здесь нам уже при-
ходится касаться одношаговых алгоритмов ослабления
регулярных помех, основанных на их кинематических
отличиях от сигналов, что пока преждевременно и вы-
ходит за рамки второй части статьи.
Эта же последовательность применения процедур
обработки демонстрируется на реальной сейсмограмме,
полученной в результате морской сейсморазведки на
Черном море (рис. 6, а). Преобразование Радона сейс-
мограммы показано на рис. 7, а. За кратнообразующую
границу было принято морское дно, отражение от ко-
торого имеет время вступления на трассе ближнего уда-
ления около 2900 мс. В соответствии с выражением (I)
была рассчитана траектория фазовых сдвигов т(а), ко-
торая изображена на рис. 7 зелёным цветом. Модель крат-
ных волн, полученная сдвигом трасс преобразования
Радона, представлена на рис. 7. о. Результат адаптивного
вычитания показан на рис. 7, в. Очевидно, что в области
интерференции первой полнократной волны и однократ-
ной волны (около 5800 мс) удалось ослабить энергию
кратного отражения и, тем самым, выделить полезный
Рис. 5. Прогнозирование и вычитание кратных волн
в области разложения Радона:
по горизонтальной оси отложен так называемый параметр мед-
ленности, измеряемый в отсчётах за канал, т. е. безразмерный
аналог величины, измеряемой в секундах за метр. Диапазон
изменения медленности - от нуля до 1/й, где v - минимальная
из кажущихся скоростей всех зарегистрированных на сейс-
мограмме ОГТ волн. В рамках задачи 1,5D отрицательные зна-
чения параметра не реализуются, а - результат преобразова-
ния Радона, полученного по исходной сейсмограмме; б - спро-
гнозированное поле кратных волн, связанных с переотраже-
нием от первого глубинного горизонта; в - результат адап-
тивного вычитания модели кратных волн
в
12
сигнал. Обработанная таким способом сейсмограмма
преобразуется обратно в пространственно-временную
область в соответствии с (2) (см. рис. 6, б). Очевидно,
что удалось достичь ослабления энергии помех, а одно-
кратное отражение, время вступления которого на трас-
се ближнего удаления около 5700 мс, может быть уве-
ренно идентифицировано и использовано при дальней-
шем анализе волнового поля. В то же время значитель-
ная часть энергии кратных отражений не была удалена.
Скорее всего, этот эффект не обусловлен несоответстви-
ем реальной среды сделанным нами предположениям о
1,5О-структуре сейсмического эксперимента. Действитель-
но, в данном случае регистрация осуществлялась на глу-
боководном участке при почти горизонтальной геомет-
рии морского дна. Это видно из рис. 8, где показан раз-
рез равных удалений, полученный по исходным сейс-
мограммам (в данной ситуации приходится надеяться
на отсутствие существенных вариаций геометрии дна в
направлении, перпендикулярном к линии профиля).
Теперь мы переходим к рассмотрению ограничений
метода, которые обусловлены несоответствием сделан-
ных нами выше предположений ситуациям, с которы-
ми приходится иметь дело на практике. Причём связаны
они будут не с отклонением структуры исследуемой сре-
ды от предполагаемой модели, но со спецификой ре-
альных условий наблюдений.
Первой особенностью большинства реальных систем
регистрации является ненулевое расстояние от источ-
ника колебаний до ближайшего приёмника. Иначе гово-
ря, на сейсмограммах почти всегда отсутствуют трассы
ближних удалений. Непосредственное применение к та-
кому волновому полю преобразования Радона приведёт
к неустойчивому результату, содержащему интенсивные
помехи, что существенно затруднит вычитание кратных
волн. Действительно, суммирование (3) ведётся не в
бесконечных симметричных пределах, а отсутствие в
сейсмограмме трасс, которые не могут быть учтены при
преобразовании, фактически эквивалентно их замене на
нулевые трассы. Поэтому до преобразования Радона сле-
дует применить процедуру экстраполяции, т. е. "пересчё-
та” трасс средних и больших удалений в трассы малых
удалений. С другой стороны, ограничения преобразова-
ний такого типа известны.
Кинематически и динамически корректное восста-
новление трасс ближних удалений важно и для приме-
нения алгоритмов прогнозирования кратных волн. Как
следует из рис. 1 и рис. 3. сигнал кратной волны на уда-
лении (Ь - а) прогнозируется по сигналу волны, зареги-
стрированной на удалении (q - а). Иначе говоря, крат-
ные волны на больших удалениях прогнозируются по
волнам, зарегистрированным на трассах малых удале-
ний. Как будет показано ниже, это свойство сохраняет-
а
750	2250	3750	5250 м
2500
3000
3500
4000
4500
5000
5500
6000
6500
7000
7500
8000
8500
мс
750	2250	3750	5250 м
3750
5250 м
2500
3000
3500
4000
4500
5000
5500
6000
6500
7000
7500
8000
8500
мс
В
750	2250
№1
'll
*?<. К:'41 Л



Рис. 6. Реальная сейсмограмма ОГТ, полученная при
морской сейсморазведке:
а - исходное волновое поле; б - результат адаптивного вы-
читания кратных волн, прогнозирование осуществлялось
сдвигом на величину (1); в - результат адаптивного вычита-
ния кратных волн, прогнозирование осуществлялось мето-
дом автосвёртки трасс в области преобразования Радона
13
2500
3000
3500
4000
4500
5000
5500
6000
6500
7000
7500
8000
8500
мс
a
О 0,5 1	1,5	2 2,5 3
Рис. 7. Волновое поле в области разложения по
плоским волнам:
зелёным цветом показана функция сдвигов, вычисленная по
формуле (1); за кратнообразующий горизонт принято мор-
ское дно; а - результат преобразования Радона, получен-
ный по исходной сейсмограмме; б - результат прогнозиро-
вания кратных волн, связанных с переотражением от морс-
кого дна; в - результат адаптивного вычитания кратных волн
в
0 0,5 1	1,5	2 2,5 3 а,о/к
ся и для алгоритмов расчета кратных волн без привлече-
ния информации о модели среды.
При подготовке данных к преобразованию Радона
также полезно воспользоваться свойством симметрии
годографов отражённых волн на сейсмограммах ОГТ,
т. е. получить трассы отрицательных удалений, симмет-
рично отобразив трассы из области положительных уда-
лений. Этот приём позволяет ослабить помехи, обуслов-
ленные невозможностью обеспечить в выражении (3)
бесконечные симметричные пределы суммирования. При
этом краевой эффект преобразования Радона за счет от-
сутствия трасс больших удалений (больше максималь-
ного удаления источник-приёмник в геометрии реаль-
ных наблюдений) остаётся. Он может быть ослаблен при
помощи пространственного сглаживания сейсмограммы:
взвешивания трасс таким образом, чтобы их энергия
плавно стремилась к нулю с ростом удаления.
Все эти преобразования в большей или меньшей сте-
пени приводят к несоответствию реальных данных сде-
ланным модельным предположениям.
Также важную роль играют искажения, обусловлен-
ные пространственной дискретностью результатов сейс-
мических наблюдений. Динамически корректный алго-
ритм прогнозирования поля кратных отражений осно-
ван на принципе выделения плоских волн и ввода фазо-
вого сдвига в каждый элемент разложения. Иначе гово-
ря, предполагается, что в результате направленного сум-
мирования (3) форма импульса, принадлежащего син-
фазности с криволинейным годографом, сохраняется,
а амплитуда сигнала на трассах в области преобразова-
ния Радона зависит от локальной кривизны годографа в
10000 15000 20000 25000 30000 35000 40000 45000 м
3000
3500
4000
4500
5000
5500
6000
6500
7000
7500
8000
8500
мс
Рис. 8. Разрез равных удалений (300 м), полученный
при морской сейсморазведке на Чёрном море
14
точке его касания прямолинейной траекторией сумми-
рования t(x) = ах. Причем считается, что указанная ам-
плитудная зависимость такова, что обеспечивается кор-
ректное динамическое соответствие спрогнозированных
кратных волн реально зарегистрированным отражениям.
То есть за счёт направленного суммирования происхо-
дит адекватный учёт геометрического расхождения при
прямом продолжении волнового поля.
Описанные теоретические заключения о динамике
волнового поля в области разложения Радона всегда
нарушаются ввиду того, что преобразование реализует-
ся не в виде (3'), но (3). То есть оно применяется к дан-
ным, дискретным по пространственной координате. Как
показывает наш опыт обработки реальных материалов,
при небольшом шаге дискретизации обеспечивается
вполне удовлетворительное соотношение амплитуд вол-
новых полей, достаточное, после этапа адаптации, для
уверенного вычитания кратных отражений. В противном
случае наблюдается динамическое несоответствие мо-
дели кратных волн реально зарегистрированным сигна-
лам, что может существенно затруднять последующий
этап адаптивного вычитания.
Иной вид динамических искажений, обусловленных
пространственной дискретностью волнового поля, свя-
зан с нарушением свёрточной модели в трассах, полу-
ченных в результате преобразования Радона. Предпола-
гается, что после разложения по плоским волнам трас-
сы представляют собой последовательность вступлений
импульсов, одинаковых по форме. Амплитуды их тако-
вы, что обеспечено корректное динамическое соответ-
ствие модели реально зарегистрированным волнам-по-
мехам. Пространственная дискретность сейсмограмм при-
водит к нарушению этого принципа. Действительно, при
направленном суммировании сигналов происходит на-
капливание импульсов с относительными запаздывани-
ями, определяемыми разностью траектории оператора
и годографа синфазности, которой эти сигналы при-
надлежат. Разумеется, результат накапливания будет за-
висеть от формы годографа каждой волны. Это приводит
к тому, что трассы в области разложении Радона не об-
ладают свёрточной структурой, но будут состоять из
импульсов различной формы. Прогнозирование модели
кратных волн осуществляется сдвигом на величину т(а).
Следовательно, на этапе адаптации не удастся достичь
динамического соответствия, подбирая один фильтр на
всю трассу.
Известной проблемой, возникающей при примене-
нии любых пространственных фильтров к дискретным
данным, является аляйсинг-эффект. Помехи такого рода
хорошо изучены, поэтому здесь мы не будем подробно
останавливаться на их анализе. Заметим только, что в за-
даче прогнозирования кратных волн имеется своя специ-
фика. а именно, при интерференции интенсивных крат-
ных и слабых однократных отражений аляйсинг- помеха,
вызванная направленным суммированием кратных волн,
может доминировать над однократными волнами, что
чрезвычайно затрудняет интерпретацию волнового поля
и контроль качества подавления помех.
Изучая 1,50-задачу прогнозирования кратных отра-
жений, мы убедились, что в области разложения по
юским волнам её решение выглядит особенно просто.
Если дополнительно учесть фактор с, отвечающий за
вклад в амплитуду сигнала коэффициентов прохожде-
ния и отражения глубинных границ раздела, встречен-
ных волной на пути от промежуточной точки выхода
луча q к приёмнику b (в т. ч. можно учесть инверсию
амплитуды за счёт отражения от свободной поверхнос-
ти), то набор волн 7|(а, f) рассчитывается при помощи
ввода сдвигов т(а)
/|(а, t) = ср (a, t - т(а)).
(4)
Искомую модель кратных отражений получим, пере-
ведя /|(а, /) в (/- х)-область.
Это же преобразование может быть реализовано и
непосредственно в пространственно-временной облас-
ти, где в рамках традиционно привлекаемых для ана-
лиза волнового поля асимптотик оно приобретает форму
преобразования миграционного типа, т. е. взвешенного
суммирования исходных трасс со сдвигами и после-
дующей одноканальной фильтрации. Пусть требуется по-
лучить модель кратных волн /| (/?, I), соответствующую
положению источника колебаний х = а = 0 и приёмника
х = Ь. Тогда
/1(Ь,О = с >'(/)* ^P(b~q,t~ g(b-q))w(b-q) , (5)
L ч	J.
где r(t) - оператор, имеющий спектральную характе-
ристику А(ю) = J — е 4 , весовые множители равны кор-
11 2л
ню из второй производной траектории суммирования
, s I() 2g(q) -г
w(q) = —. Также можно показать, что выражение
V dq2
(5) совпадает с дискретным аналогом интегрального
оператора, осуществляющего прямое продолжение вол-
нового поля по лучу отраженной от кратнообразующего
горизонта волны.
Здесь важно отметить, что оператором, осуществля-
ющим преобразование исходной сейсмограммы в поле
кратных отражений, является производная по нормали
к свободной поверхности от динамического годографа
отражённой от кратнообразующего горизонта волны,
вызванной линейным источником, пересекающим про-
филь наблюдений в точке с координатой^ = q. Действи-
тельно, спектральным эквивалентом соотношения (5)
является
£1 (Ь, го) ~ с
P(b - q,(£>)w(b- q)e J0)8(b 44
ч
а лучевое представление для динамического годографа:
С е
K{b-q)
где K{b - q) - фактор геометрического расхождения. С
точностью до оператора одноканальной фильтрации
асимптотическое выражение для производной по нор-
мали к свободной поверхности динамического годогра-
фа может быть записано в виде cw(b - q)ej'N^b~4\ В соот-
ветствии с (6) именно последнее выражение является
оператором, применяемым к наблюдённому полю.
15
Подводя итог описанию алгоритма подавления крат-
ных волн, реализуемому в области разложения по плос-
ким волнам, представим его в виде последовательности
процедур обработки исходного поля:
Шаг 1. Применение преобразования Радона к сейс-
мограмме ОГТ после её предварительной подготовки.
Тем самым осуществляется переход из пространствен-
но-временной области в область разложения по плос-
ким волнам.
Шаг 2. Выбор кратнообразующего горизонта и оце-
нивание кинематических параметров волны, однократ-
но отраженной от этого горизонта. Расчет функции сдви-
гов т(а).
Шаг 3. Прогнозирование кратных волн путем сдвига
трасс сейсмограммы в области разложения Радона.
Шаг 4. Получение дополнительной модели кратных
волн старшего порядка кратности при помощи сдвига
на удвоенную величину 2т(а), и умножение на 0,5 пер-
вой полнократной волны от выбранного кратнообразу-
ющего горизонта.
Шаг 5. Возможно применение мьютинга (обнуления
областей, занятых исключительно помехой), если пос-
ле преобразования Радона сигналы и помехи разделя-
ются.
Шаг 6. Адаптация и вычитание.
Шаг 7. Обратное преобразование Радона результата
адаптивного вычитания. Переход в пространственно-вре-
менную область.
Альтернативным способом прогнозирования цуга
кратных волн является использование преобразования
миграционного типа (5), реализуемого в пространствен-
но-временной области. В данном случае в качестве апри-
орной информации привлекается оценка годографа от-
раженной от кратнообразующего горизонта волны. Ос-
новным преимуществом такого расчёта является возмож-
ность корректного учёта краевых эффектов за счёт огра-
ничения апертуры суммирования по переменной q.
Понятно, что при неплоских границах раздела сред,
т. е. при анализе задач 2D и 3D уже нельзя осуществлять
прогнозирование путём введения постоянного фазового
сдвига, т. е. использовать принцип переноса плоских волн.
Действительно, наличие кривизны у отражающих гра-
ниц приводит к эффекту фокусировки и расфокусиров-
ки волн. Иначе говоря, если на систему слоёв падает
одна плоская волна, то на выходе из среды фронт вол-
ны перестаёт быть плоским, а регистрируется сложное
волновое поле, которое, в свою очередь, будучи разло-
жено на плоские волны, содержит целый их набор.
В заключение еше раз отметим важную особенность
алгоритма прогнозирования и вычитания кратных волн:
он является полным аналогом традиционного метода
предсказывающей деконволюции, но применяется к
трассам, преобразованным в область разложения по плос-
ким волнам. Сохраняются и все динамические особен-
ности процедуры, подробно рассмотренные в первой
части статьи. В частности, для получения модели крат-
ных волн с правильным соотношением амплитуд требу-
ется построение двух трасс. Первая получается сдвигом
исходной трассы на время, равное времени вступления
однократного отражения от выбранного глубинного крат-
нообразуюшего горизонта, вторая - введением удвоен-
ного сдвига. Аналогично случаю 1D, для вычитания цуга
кратных отражений, связанных с другим глубинным го-
ризонтом, процедуру приходится повторять.
1,50-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН
БЕЗ ПРИВЛЕЧЕНИЯ АПРИОРНОЙ ИНФОРМАЦИИ
О СТРУКТУРЕ СРЕДЫ
В горизонтально-слоистой среде плоская волна, подхо-
дящая к свободной поверхности под углом 6, переотра
зившись в нижнее полупространство, порождает краг
ные отражения, подходящие к поверхности под этим
же углом. Таким образом, в области разложения по плос-
ким волнам прогнозирование цуга кратных отражении
для случая 1,5D аналогично их моделированию в ID-
задаче. В самом деле, поле колебаний каждой плоской
волны, зарегистрированное на поверхности, может ин-
терпретироваться как источник вторичной волны, по-
сылаемой в среду. Тогда, чтобы спрогнозировать кра|
ную волну, надо произвести свёртку зарегистрирован
ного поля с функцией источника, т. е. посчитать авто-
свёртку трасс зарегистрированного поля в области пре-
образования Радона.
Для построения модели кратных отражений всех по-
рядков кратности нужно, вообще говоря, получить бес-
конечное число автосвёрток
+ ...,
где т(а, t) - полная модель кратных волн; дгДа, t) -
А-кратная автосвёртка. Однако ввиду быстрого затухания
амплитуды волны с ростом порядка её кратности, а так-
же с учетом того факта, что сейсмическая запись огра-
ничена по времени, и, следовательно, регистрируются
волны не всех старших порядков кратности, достаточно
бывает ограничиться расчётом автосвёрток второго
т । (a, t) = р(а, t) * р(а, t)
или третьего порядка
/й2(а, t) = р(а, г) * р(п, t) * р(а, t) .
тем более что уже результат первой автосвёртки содер-
жит полный набор волн, т. е. отражения всех порядков
кратности. Расчёт дополнительных автосвёрток необхо-
дим для достижения лучшего соответствия амплитуд
модели кратных волн реально зарегистрированным сиг-
налам. Соответствующие вопросы подробно обсуждалиа
в первой части статьи и будут проиллюстрированы ниже.
При применении процедуры прогнозирования поля
кратных волн без привлечения априорной информации
о структуре среды к реальным наблюдениям приходится
сталкиваться со всеми проблемами, которые обсужда
лись в предыдущем параграфе. Следует специальным
образом подготовить исходные сейсмограммы: проэкст
раполировать трассы ближних удалений, ослабить крае
вые эффекты и т. д.
Можно аналогично тому, как мы это сделали выше,
рассмотреть процесс расчёта модели кратных волн в
пространственно-временной области. Для этого полу-
чим двумерное преобразование Фурье исходного поля
p(x,t), которое обозначим Р(к, со). Если через Р(а, и)
16
обозначить преобразование Фурье р(а, I) по перемен-
ной /, то понятно, что
Р(а, со) = Р(асо, со).
Получим спектральный аналог автосвёртки т^а, t) ~
p(a,t)*p(a. Г) в виде Л/] (а, со) = Р(а, а>)Р(а, со) или
Л/|(а, ш) = Р(аы, со) Р(аа, со). Следовательно, для полу-
чения поля кратных волн в (х - /)-области достаточно
возвести в квадрат двумерный спектр исходного поля и
сделать обратное преобразование Фурье. Такая после-
довательность вычислений эквивалентна двумерной
автосвёртке сейсмограммы р(х, t)
т 1 (b, t) = 51X	т) P(b	т)-
q т
Однако в формуле (7) имеется одна неточность. Так
как в этом выражении в качестве оператора, осуществ-
ляющего прогнозирование кратных волн по исходному
полю, выступает само же исходное поле. При этом оно
может рассматриваться как суперпозиция динамических
годографов всех наблюденных волн. Как мы знаем из пре-
дыдущего раздела, динамически корректным оператором
прогнозирования цуга кратных отражений, связанных с
некоторым глубинным кратнообразующим горизонтом,
является не сам динамический годограф отражённой от
этого горизонта волны, но его производная по нормали
к свободной поверхности. Вычисление производной по
нормали оказывается асимптотически эквивалентным
, cos 6
умножению динамического годографа на---- и после-
v
дующей одноканальной фильтрации. Такое преобразо-
вание динамического годографа удобно осуществлять в
области его разложения по плоским волнам а - I. где
c	cos 6 б ~ 2л
множитель приобретает вид --= yl - (аг) и может
г
быть применён к каждой трассе-элементу разложения.
Олноканальную фильтрацию целесообразно производить
на этапе обратного преобразования в (х - г)-область (2),
совместив её оператор с фильтром г(Р). Обозначим ре-
зультирующее волновое поле как р(х, t) и перепишем
(7) в виде
= ^^p(q,T)p(b-	(8)
Q t
Двумерная свёртка сейсмограмм (8) может быть ре-
ализована в виде произведения их спектров в {к - оэ)-
области
Мх(к, то) = Р(к, а)Р(к, (о)
с последующим переводом результата обратно в про-
странственно-временную область. Такая реализация мо-
жет приводить к повышению быстродействия алгорит-
ма прогнозирования. С другой стороны, возникают слож-
ности с соблюдением пределов суммирования (8), а так-
же усиливаются краевые эффекты, приводящие к воз-
растанию уровня артефактов.
Еще одно отличие результата прогнозирования, по-
лучаемого без привлечения информации о модели сре-
ды, от (5) обусловлено эффектом автосвёртки сейсми-
ческого импульса. При моделировании в соответствии с
(5) оператор совпадал с производной от динамическо-
го годографа отражённой волны, вызванной линейным
источником с формой сигнала в виде дельта-функции. В
(8) используется реально наблюдённый годограф, что
приводит к эффекту автосвёртки импульса. Если форма
импульса источника стационарна (так бывает при морс-
кой сейсморазведке) или сейсмограммы предваритель-
но скорректированы за приповерхностные условия, то
можно избежать процедуры оценивания и коррекции
импульса и отложить компенсацию эффекта автосвёрт-
ки до этапа адаптации модели кратных отражений к ре-
ально зарегистрированным волнам.
Результат вычитания кратных волн с использовани-
ем прогнозирования методом автосвёрток трасс в обла-
сти преобразования Радона показан на рис. 6, в. Адап-
тивное вычитание модели кратных волн из исходной
сейсмограммы производилось в условиях, тождествен-
ных тем, при которых было получено волновое поле на
рис. 6, б. Очевидно, что энергия помех на рис. 6, в ниже,
чем на рис. 6, б. Следовательно, модель кратных отраже-
ний была получена с лучшими динамическими свой-
ствами. В данном случае эту особенность можно объяс-
нить тем, что в результате прогнозирования с привлече-
нием информации о модели среды мы получаем крат-
ные волны, связанные с переотражением только от мор-
ского дна. В то же время, как видно из рис. 8, цуг крат-
ных волн формируется тонкослоистой пачкой горизон-
тов, имеющей значительную протяженность по глубине.
В рамках алгоритмов, основанных на предварительном
построении модели среды, рассчитать кратные волны,
обусловленные переотражением от такой пачки, не пред-
ставляется возможным. В то же время для альтернатив-
ного подхода к прогнозированию обработка таких дан-
ных не вызывает затруднений. Таким образом, в первом
случае был получен не весь набор волн-помех, что на
этапе адаптации привело к получению смещённых оце-
нок фильтров и что обусловило меньшую степень ос-
лабления кратных волн.
Подводя итог описанию алгоритма подавления крат-
ных волн, которому была посвящена настоящая глава,
представим его в виде последовательности процедур
обработки исходного волнового поля:
Шаг 1. Применение преобразования Радона к сейс-
мограмме ОГТ после её предварительной подготовки.
Шаг 2. Прогнозирование кратных волн путём расчета
автосвёрток трасс сейсмограммы в области разложения
Радона (с целью повышения быстродействия расчёт ав-
тосвёрток целесообразно выполнять в спектральной об-
ласти). Для достижения лучшей динамической адекват-
ности модели кратных волн реально зарегистрирован-
ным помехам до k-vo порядка кратности включительно
следует получить дополнительные наборы кратных волн
путём расчёта автосвёрток до к-ro порядка. Обычно к < 3.
Шаг 3. Возможно применение мьютинга, или обну-
ления областей, занятых исключительно помехой, если
после преобразования Радона сигналы и помехи разде-
ляются.
Шаг 4. Адаптация и вычитание.
Шаг 5. Обратное преобразование Радона результат
адаптивного вычитания, т. е. переход в пространствен-
но-временную область.
17
Альтернативным способом прогнозирования цуга крат-
ных волн является использование преобразования (8),
реализуемого в пространственно-временной области. В
данном случае оператором прогнозирования является
специально подготовленная сейсмограмма р(х. t).
Реализация алгоритма расчета модели кратных волн
способом одноканальной автосвёртки трасс в области
разложения Радона аналогична рассмотренному в пер-
вой части статьи алгоритму ретрокорреляции и являет-
ся его обобщением на случай более сложной модели
среды. Их различие заключается в том, что для получе-
ния динамически адекватного результата одну из трасс,
полученную преобразованием Радона, следует умножить
на весовой коэффициент, пропорциональный косину-
су угла выхода луча, что требует некоторой, пусть и
минимальной, априорной информации, а именно -
скорости распространения волны в приповерхностном
слое. Алгоритм ректрокорреляции ориентирован на од-
номерную модель сейсмического эксперимента, при
которой косинус угла выхода луча любой волны равен
единице, поэтому удается построить динамически адек-
ватную модель кратных отражений самым простым спо-
собом.
КАКОЙ АЛГОРИТМ
ПРИМЕНЯТЬ НА ПРАКТИКЕ?
Несмотря на принципиальные ограничения описанных
во второй части статьи методов прогнозирования и вы-
читания кратных волн, они часто применяются на прак-
тике. Популярность подходов, основанных на разложе-
нии Радона, объясняется их простотой и быстродействи-
ем. Область наибольшей эффективности метода - “ло-
кальная 1,5О-задача”, т. е. случай, когда можно пренеб-
речь наклонами и кривизнами отражающих горизонтов
на некоторой небольшой пространственной базе, сопо-
ставимой с длиной косы приемников. Неточности, воз-
никающие при прогнозировании поля кратных отраже-
ний, компенсируются уже на этапе адаптации к исход-
ным данным. Высокая вычислительная эффективность
обусловлена возможностью применения процедур обра-
ботки к сейсмограммам ОГТ.
В силу того, что оба метода имеют практическую зна-
чимость, мы дадим достаточно подробное их сравнитель-
ное описание и вновь, как и в первой части статьи, за-
дадимся вопросом: какой из двух предложенных алго-
ритмов прогнозирования кратных волн предпочесть?
Наглядное представление о динамических особенно-
стях процедур могут дать результаты обработки модель-
ных данных. Был получен набор сейсмограмм, соответ-
ствующий сейсмическим наблюдениям в случае гори-
зонтально-слоистой среды, состоящей из двух отража-
ющих границ, находящихся на глубинах 300 и 900 м,
интервальные скорости - 1500 и 2500 м/с. Фрагмент сей-
смограммы ОГТ показан на рис. 9, а. Были построены
модели кратных отражений при помощи обоих методов
прогнозирования. Полученные результаты сопоставлены
с исходной записью на рис. 9, б. Изображены трассы,
соответствующие удалению источник-приёмник 100 м.
Наиболее интенсивные отражения обозначены буквами.
Очевидно, что прогнозирование с использованием мо-
дели среды (были посчитаны волны, связанные с пере-
отражением от первого кратнообразующего горизонта)
обеспечило более точное по сравнению с прогнозиро-
ванием без модели среды динамическое соответствие
кратных волн реальным отражениям.
Неизбежные погрешности в амплитудах волн стар-
шего порядка кратности присутствуют в обоих результа-
тах. Этот эффект обусловлен тем, что процедура про-
гнозирования без привлечения модели среды является
итерационной, и на каждом её шаге получаемое поле
кинематически адекватно реальным волнам, но дина-
мически корректно прогнозируются волны до к-го по-
рядка кратности включительно, где к - номер итерации.
Здесь была применена одна итерация, т. е. следует ожи-
дать, что корректно предсказаны волны первого поряд-
ка кратности. В данном случае это сигналы a vt f Как вид-
но из рисунка, их амплитуды предсказаны с удовлетво-
рительной точностью. С другой стороны, одна итерация
алгоритма прогнозирования с учётом модели среды также
не обеспечивает динамически корректного результата.
Однако, как видно из рис. 9, б, несмотря на некоторую
динамическую неточность, этот результат предсказания
кратных волн предпочтителен для их последующего вы-
читания. Потенциальным ограничением такого способа
прогнозирования является расчет цуга кратных волн,
связанных с переотраженисм лишь от одного кратнооб-
разующего горизонта. Для получения остальных кратных
волн расчёты приходится повторять. Так, например, трас-
са С не содержит первой полнократной волны f отра-
женной от второго горизонта.
Заметим, что оба метода прогнозирования правильно
учитывают эффект зависимости коэффициента отра-
жения от угла падения и отражения волны (AVO-эф-
фект). Если модель среды не привлекается, динамичес-
кая адекватность AVO-эффекта обусловливается тем
что оператор прямого продолжения волнового поля
(здесь - наблюденные годографы отражённых волн) за-
регистрирован уже “с учётом” этого фактора. Если опе-
ратор прямого продолжения вычисляется в рамках из-
вестной модели среды, то AVO-эффект может быть лег-
ко посчитан явным образом при получении динамичес-
кого годографа однократного отражения. На рис. 9, в срав-
ниваются фрагменты исходной сейсмограммы ОГТ и
поля кратных волн, спрогнозированного с использо-
ванием информации о модели среды. Хотя с ростом
удаления источник-приёмник в горизонтально-слоис-
той среде геометрическое расхождение волны также воз-
растает, и это приведи! к затуханию её амплитуды, но
за счёт AVO-эффскта затухание не только компенсиру-
ется, но для небольших удалений даже наблюдается рост
амплитуды. Разумеется, такое явление имеет место как
для однократных, так и для кратных отражений. Как
следует из рис. 9. в. модель кратных волн получена ди-
намически корректно, т. е. так, что AVO-эффекты в ис-
ходной сейсмограмме и модели кратных волн одинако-
вы. То, насколько существен учет такого рода эффек-
тов для последующего вычитания полученной модели
кратных волн, мы здесь не обсуждаем. По нашему мне-
нию, гладкие по времени искажения амплитуд могу:
быть скорректированы при адаптации модели кратных
волн к реально зарегистрированным сигналам. Также
отметим, что наличие дополнительных возможностей
18
для более точного моделирования может быть полез-
ным в разных сейсмических задачах.
Ещё раз обратим внимание на важное свойство пре-
образования (8). В отличие от (5), осуществляющего мо-
делирование волн, связанных с переотражением лишь
от одного глубинного кратнообразуюшего горизонта,
расчёт в соответствии с (8) приводит к получению поля
m^b, t) всех типов кратных волн, связанных с переот-
ражением в нижнее полупространство от свободной по-
верхности и вверх от всех глубинных границ.
ЗАКЛЮЧЕНИЕ. Прогнозирование кратных отражений
всегда производится методами прямого продолжения вол-
нового поля, независимо от того, используется ли при
этом оценка глубинно-скоростной модели среды. Вооб-
ще говоря, преобразование реализуется в виде взвешен-
ного суммирования со сдвигами и последующей однока-
нальной фильтрацией исходных трасс, т. е. применяются
операторы миграционного типа. Однако в рамках про-
стых моделей сейсмического эксперимента операторы вы-
рождаются и приобретают простую структуру.
О	1000	2000 м
450
600
750
900
1050
1200
1350
1500
1650
1800
1950
2100
2250
2400
2550
2700
2850
МС
Рис. 9. Сравнение работоспособности методов
прогнозирования кратных волн на синтетичес-
ких данных:
а - фрагмент исходной сейсмограммы ОГТ, получен-
ной для горизонтально-слоистой модели среды, содер-
жащей две отражающих границы; б - результаты про-
гнозирования кратных волн при помощи алгоритма, не
привлекающего модель среды (В), с привлечением мо-
дели среды (С), исходная трасса (А); наиболее интен-
сивные отражения: а - первая полнократная волна, от-
раженная от первого горизонта, время вступления около
800 мс; b - однократная волна, отраженная от второго
горизонта, время вступления около 880 мс; с - вторая
полнократная волна, отраженная от первого горизонта,
время вступления около 1200 мс; d - частично-кратная
волна младшего порядка кратности, отразилась по од-
ному разу от первого и от второго горизонтов, время
вступления около 1280 мс; е - частично-кратная волна
старшего порядка кратности, имеет два отражения от
первого горизонта и одно - от второго, время вступле-
ния около 1680 мс; f - вторая полнократная волна от
второго горизонта, время вступления около 1760 мс (од-
нократная волна, отраженная от первого горизонта,
имеет время вступления около 400 мс и на рисунке не
показана); в - фрагмент исходного волнового поля (сле-
ва) и результат прогнозирования кратных волн с ис-
пользованием модели среды (справа)
19
Если строение среды в интервале глубин от свобод-
ной поверхности до кратнообразуюшего горизонта из-
вестно, то в задаче 1D динамически адекватное прогно-
зирование цуга кратных отражений было достигнуто при
помощи внесения постоянного фазового сдвига в ис-
ходные трассы. В задаче 1,5D корректное прогнозирова-
ние может осуществляться при помощи внесения по-
стоянного фазового сдвига в результат разложения по
плоским волнам. Фазовый сдвиг определяется годогра
фом волны, отраженной от кратнообразующего гори-1
зонта. Таким образом, задача свелась к одномерной, но
реализуемой в области преобразования Радона.
ВВЕДЕНИЕ. Вопросы, касающиеся разложения поля
точечного источника по базису плоских волн, наиболее
подробно освещены в работах [2, 4J.
Мы будем прибегать к лучевым представлениям для
анализа волнового поля и получения некоторых анали-
тических результатов. Обоснование лучевого метода мож-
но найти в работе [3].
При выводе формул используется асимптотическая
оценка интегралов методом стационарной фазы. Описа-
ние метода содержится, например, в [8], а его обосно-
вание осуществляется в рамках той же высокочастотной
асимптотики, на которой базируется лучевой метод.
1,50-МОДЕЛЬ СЕЙСМИЧЕСКОГО ЭКСПЕРИМЕНТА
Специальный приём, позволяющий преобразовать ди-
намику волнового поля, возбужденного точечным ис-
точником и распространяющегося в горизонтально-сло-
истой среде, к динамике аналогичного поля, которое
было бы вызвано линейным источником, описан в ра-
боте [18].
1, SD-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН
ПРИ НАЛИЧИИ МОДЕЛИ СРЕДЫ В
В ранней работе Э. Робинсона [13] был предложен спо-
соб подавления ревербераций по одной трассе - пред-
сказывающая деконволюция. Отсутствие свойства пери-
одичности кратных отражений в трассах ненулевого уда-
ления и их периодичность после преобразования сейс-
мограммы в область разложения по плоским волнам были
замечены М. Танером [14]. При этом оказалось, что из-
за наличия как полнократных, так и частично-кратных
волн, корректного динамического соответствия прогно-
зируемой модели реально зарегистрированным сигна-
лам можно достичь лишь путем построения дополни-
тельной модели, содержащей волны старшего порядка
кратности [И].
Точное определение понятия динамического годо-
графа можно найти, например, в учебнике по сейсмо-
разведке [6].
Теория продолжения волновых полей изложена в ра-
боте [10], но во второй части статьи, как и в первой, в
своей общей формулировке она нам пока не понадобится.
Условия взаимной однозначности прямого и обрат-
ного преобразований Радона, применяемого к волново-
му полю, полученному на дискретной сети наблюде-
ний, тесно связаны с проблемой аляйсинг-эффекта. Этот
вопрос подробно изучен в работе [15]. Как правило, для
ПРИМЕЧАНИЯ
корректного отображения требуется получить значитель-
но большее число трасс сейсмограммы в области пре-
образования Радона, чем содержится в исходной сейс-
мограмме.
Дополнительное подавление энергии кратных волн
до обратного преобразования из области разложения
Радона в пространственно-временную область удобно
осуществлять обнулением зон на сейсмограмме, заня-
тых помехами. Можно показать (соответствующее обо-
снование дано в статье С. А. Нахамкина [9]), что такой
способ эквивалентен применению режекторной веерной
фильтрации с переменной по пространственной и врс
менной координатам характеристикой фильтра.
Комментарий к формуле (5)
Покажем, что с точностью до некоторой амплитуд-
ной поправки преобразование (4) эквивалентно двумер-
ной свертке исходного поля р(х, t) с оператором, вы-
числяемым по динамическому годографу волны, отра-
женной от кратнообразующего горизонта.
Получим модель кратных волн /](Ь. t), соответству-
ющую положению источника колебаний х = а = 0 и
приёмника х = Ь. Для этого подставим (4) в (2)
/| (b, t) = с J ~р{а, t - т(а) + ab)da.
Используя представление/? (a, t) в виде набора плос-
ких волн (3). перепишем последнее равенство
l[(b,t) =
= с /•(/) * J da У /?(</, t - т(а) + а(Ь - д))
L— I
(П1)
В спектральной области выражение (П1) преобразу-
ется в
Е|(7?,<д) = с А(со) J da У, Р(д, со)е7“( Т(а)+о(Ь 9»
L— ч
Пусть f(q, t) соответствует двумерному оператору,
преобразование Фурье по координате / которого опре-
деляется выражением
Д(</,£и) = J
(П2)
Тогда, используя введенное обозначение и меняя
порядок интегрирования и суммирования, (П1) можно
переписать как
20
I}(b,t) = c r(t)*
У p(Q-1) * f{b - q, t)
Q
Интеграл (П2) можно оценить методом стационар-
ной фазы. С учётом того, что сдвиг т(а) определяется
выражениями (1) или (Г), получим
1^с 7 4 l^2S(^) c-jwg(<j)
\Н v э^2
E(?,w)
Введя обозначение для спектра оператора однока-
нальной фильтрации
/?(о>) =
и. соответственно, для его временной характеристики -
?(/), получим выражение для прогнозирования кратных
волн в виде взвешенного пространственного суммиро-
вания трасс исходного поля со сдвигами и последую-
щей одноканальной фильтрации, т. е. преобразования
миграционного ти па
/|(й.г) = с /'(/)* У p(b-q,t-g(b-q))w(b- q) ,
L я	J,
(ПЗ)
где r(t) =	и, следовательно, Л(ю) = . -—-е 4
V 2п
а весовые множители определяются как w(q)	—-^2.
V а?2
Ниже мы покажем, что (ПЗ) соответствует дискрет-
ному аналогу интегрального оператора прямого продол-
жения волнового поля, где к качестве функции Грина
фигурирует волна, возбужденная линейным источником,
расположенным в заданном пункте взрыва, и отражен-
ная от кратнообразующей границы.
Вернёмся к рассмотрению процедуры прогнозирова-
ния цуга кратных волн. Изучим её специфику с позиций
теории продолжения волновых полей и покажем, что
соответствующие вычисления могут быть осуществлены
в области разложения по плоским волнам, а получае-
мый таким образом результат кинематически и динами-
чески адекватен искомому полю.
В современных графах обработки результатов сейс-
мических наблюдений часто применяются прямое и об-
ратное продолжения волнового поля. Соответствующие
процедуры хорошо знакомы геофизикам-обработчикам.
Методами прямого продолжения поля можно осуществ-
лять, например, коррекцию за рельеф. Обратное про-
должение поля применяется при миграции, “погруже-
нии” сейсмограмм и т. д. Здесь нас будут интересовать
вопросы построения модели кратных отражений мето-
дами прямого продолжения волнового поля.
Известно, что соответствующее преобразование до-
стигается взвешенным суммированием сейсмических
трасс. Пусть требуется осуществить пересчёт поля таким
образом, чтобы “переместить” сейсмоприёмники с ис-
ходного, возможно неплоского, уровня наблюдений
на новый уровень г2- Тогда траектория суммирования
трасс совпадает с годографом волны, возбужденной
источником, расположенным на уровне г2- 11 регистри-
руемой на уровне Z\- Такую траекторию часто называют
годографом дифрагированной волны, а сам оператор
интегрального преобразования (если скорость между Z\
и <2 постоянна) - интегралом Кирхгофа. Значения весо-
вых коэффициентов определяются через динамические
характеристики этой волны.
Если ?| и г2 разделены неоднородной средой, содер-
жащей отражающие горизонты, то построение точного
оператора наталкивается на существенные трудности.
Поэтому в большинстве практически важных случаев
считается оправданным использование высокочастотной
или, что то же самое, лучевой асимптотики при анали-
зе волнового поля, и это существенно упрощает расчёт
искомого оператора.
Преобразование вида прямого продолжения волно-
вого поля допускает обобщение, позволяющее обосно-
вать его использование при решении задачи прогнози-
рования кратных отражений по исходным сейсмограм-
мам. В соответствии с пространственным аналогом схе-
мы, показанной на рис. 1, процесс формирования сиг-
нала реально наблюденного кратного отражения про-
исходит следующим образом. Волна, вызванная источ-
ником, расположенным в точке а, отражается от глу-
бинных границ и подходит к свободной поверхности,
где регистрируется сейсмоприёмниками. Отметим, что
колебания наблюдаются во всех точках плоскости, но
при 2О-сейсморазведке регистрируются только на ли-
нии профиля (выше мы условились считать, что линия
профиля имеет нулевую координату у). Сигнал отража-
ется от свободной поверхности, вновь проходит через
среду до кратнообразующих горизонтов, где переотра-
жается и возвращается обратно к поверхности уже в
виде полнократной или частично-кратной волны. Ин-
тегральное преобразование прямого продолжения поля
призвано имитировать именно такую схему формиро-
вания волн-помех. В принципе, должна учитываться вся
выходящая на свободную поверхность энергия волны
первичного отражения путем суммирования по всем
значениям координатх и у. Тогда искомый сигнал крат-
ного отражения, наблюдаемый в точке регистрации Ь,
моделируется в результате сложения колебаний, заре-
гистрированных на каждом элементарном сегменте пло-
щади ДхДу (в пределе - бесконечно малой площади Ад-Ду).
Используя понятия, заимствованные из лучевого пред-
ставления волнового поля, можно сказать, что сумма
трасс, полученных для всех пунктов регистрации q,
после введения в них фазовых сдвигов, равных време-
ни распространения волны вдоль луча, связывающего
точки q и Ь, и весовых множителей, определяемых в
соответствии с затуханием энергии, обусловленным
геометрическим расхождением при распространении
вдоль этой траектории, представляет собой динамически
адекватную модель искомого кратного отражения. В со-
ответствии со схемой, показанной на рис. 1, такая тра-
ектория и такие веса определяются динамическим го-
дографом волны, вызванной точечным источником,
имеющим координату q, отраженной от некоторого
глубинного кратнообразуюшего горизонта и зарегист-
рированной в точке Ь.
21
Понятно, что при профильной схеме наблюдений мы
можем учесть вклад лишь от точек, лежащих на линии
профиля у = 0. Однако в рамках рассматриваемой нами
1,5Г)-залачи известное на линии профиля волновое поле
можно по аналитическим формулам пересчитать во все
точки плоскости наблюдений, или. что то же самое,
при интегрировании по двум пространственным пере-
менным х и у оценку интеграла по у можно получить
аналитическими методами.
Перейдём к математической формулировке проце-
дуры моделирования цуга кратных волн. Спектр Фурье
трассы кратных отражений будем искать как результат
действия интегрального оператора прямого продолже-
ния волнового поля в, вообще говоря, неоднородной
среде, т. е. преобразования вида пространственного сум-
мирования
= P(q,ofl~G(q,y, b,co)d^dy,	(П4)
Z7t J J	OZ
У 9
где G(q, у. b, cd) - динамический годограф волны, воз-
буждённой точечным источником, имеющим координа-
ты на плоскости наблюдений (q, у), и зарегистрирован-
ной сейсмоприёмником, расположенным на линии на-
блюдений у = 0 в точке х = Ь (функция Грина). P(q, со) -
спектр исходной трассы, наблюденной при положении
источника в точке с координатой х = а = 0 и приёмника
в точке х = q.
Лучевое представление динамического годографа от-
ражённой волны в неоднородной среде
С(?, у, Ь, CD) = - -С-А- е^я.у.ь\ (П5
k(q,y,b)	v >
где с - фактор, отвечающий за вклад в амплитуду сигнала
коэффициентов прохождения и отражения границ разде-
ла; ATq. у, Ь) - геометрическое расхождение; g[q, у, Ь) -
время распространения вдоль траектории луча. Отсюда
(в рамках традиционных и всегда справедливых в реаль-
ных условиях асимптотических допущений)
A G(q, у, b, cd) = ~Jco^ —С~—	(П6)
dz	v K(q,y,b)
где мы воспользовались тем свойством, что вертикаль-
ная производная годографа g(q, у, b) равна отношению
косинуса угла выхода луча на свободную поверхность к
скорости v в верхнем слое в окрестности точки выхода
cos 6	,
луча, т. е.--. Этот факт является следствием формулы
v
Бендорфа [6]. Понятно, что угол выхода луча зависит от
координаты точки регистрации, но для краткости обо-
значений в формулах эту зависимость мы будем лишь
подразумевать, но не указывать в явном виде.
Поле распространяющейся волны не зависит от ко-
ординаты у, так как мы полагаем, что оно вызвано ли-
нейным источником (или преобразовано так, что соот-
ветствует полю линейного источника), а строение среды
инвариантно относительно у. Тогда в интегральном выра-
жении (П4) от переменной у зависит только вертикаль-
ная производная динамического годографа, и интегри-
рование по у может быть произведено аналитически.
22
Подставим (П6) в (П4)
£](Ь.ш) =
-Juc
2лг
$РЫ
Q
г cos 6 e-jcag{q,y,b)^
\K(q,y,b)
(П?
Теперь воспользуемся формулой, связывающей про-
изводные годографа волны g{q, у, Ь) (в рассматривае
мой задаче такой годограф соответствует траектории
суммирования оператора продолжения поля) с геомет-
рическим расхождением [1, 17]
,.	1	cos 6 > cos 6 э
'.y.b) = - —--------L- у -2----,
v d2g(g>b^) d2g(g,yj)
у Эу2
(П8)
где через 6i и 62 соответственно обозначены угол вы-
хода луча из точечного источника с координатами (л fl
и угол подхода луча к приёмнику Ь. Очевидно, что для
среды, состоящей из плоско-параллельных границ,
,,	Э.е2(й,у,Ь) dg2(q,y,b)
cosO. = COS0-,, а --	< =	. То есть годе-
1	2 dqdb	dq2
графы волн на сейсмограммах ОПВ и ОПП (общегопун-
кта приема) совпадают. Учитывая эти свойства, упрос-
тим (П8)
. cos е 1	1
К(д, у, Ь) =-===== -===
V d2g(q, У,Ь) d2g(q,y,b)
\ dy2 \ dq2
и подставим результат в (П7). Сгруппировав сомножите-
ли, получим выражение для интеграла по у
д g(q, У, b)c ~jcag(q,y.b)^
ду2
и оценим его методом стационарной фазы
С учётом полученного результата выражение (П7)
сводится к
L|(b,OJ)
£(cd)F P(ftCD) i^L^e-P^Vdq
N dq2
где R(w) = . ——-e 4 - спектральная характеристика опе-
V 2л
ратора одноканальной фильтрации r(t), который уже
фигурировал в выражении (ПЗ), а переменная у в тра-
ектории суммирования g(q, у = 0, Ь) опущена, так как
теперь интегрирование ведётся только вдоль линии про-
филя. Наконец, привлекая оговоренное выше свойство
симметрии задачи 1,5D, согласно которому g(q,b) =
g(b - q) и переходя в пространственно-временную об-
ласть, получим
r(i)
J p(b - q., t - g(b - q))w(b-q)dq
ч
что является непрерывным аналогом дискретного сум-
мирования (ПЗ) и где, как и в (ПЗ), весовые множите-
ли w(q) определяются корнем из второй производной
траектории суммирования.
Заметим, что выражение для функции весов в виде
корня из второй производной траектории суммирова-
ния хорошо известно в теории сейсмической миграции
[5]. Смысл его заключается в том, что амплитуда резуль-
тата суммирования сигналов, принадлежащих годогра-
фу плоской волны, будет зависеть от кривизны g(q) в
точке их касания. Чем ближе g(q) к прямой, гем боль-
шая амплитуда суммарного сигнала накопится вследствие
почти синфазного суммирования. И наоборот, при не-
синфазном суммировании (большая кривизна g(q)) ам-
плитуда результата будет меньшей. Учет весовых коэф-
фициентов при суммировании призван скомпенсировать
такой эффект.
1,50-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН
БЕЗ ПРИВЛЕЧЕНИЯ АПРИОРНОЙ ИНФОРМАЦИИ
О СТРУКТУРЕ СРЕДЫ
Комментарий к формуле (8)
Иным способом получения интересующих нас соот-
ношений в пространственно-временной области явля-
ется подстановка в выражение для прямого продолже-
ния волнового поля (П4) сейсмограммы P(q, со) и ди-
намического годографа отражённой волны в виде их
разложения по плоским волнам в соответствии с (2)
(см. [7]). Это приведёт к несколько утомительным вык-
ладкам, результатом которых, после суммирования по
всем динамическим годографам волн, связанных с пе-
реотражением от всех глубинных кратнообразующих го-
ризонтов, станет выражение (8).
В отличие от динамически неточного (7), преобразо-
вание (8) содержит не только величины, явно присут-
ствующие в поле однократно отраженной волны. Отсут-
ствует информация об угле выхода луча 6. Поэтому, ради
простоты вычислений, на практике зачастую пренеб-
, ,	cos 6
регают весовым коэффициентом ---, а поле кратных
v
волн прогнозируют непосредственной пространственной
свёрткой (7) исходных сейсмограмм [12, 16]. Рассмот-
рим; какие динамические искажения вносятся при этом
a	q	b
Рис. 10. Модель среды, содержащая одну глубин-
ную отражающую границу
на примере простейшей ситуации. Пусть среда содержит
всего одну отражающую границу, расположенную на
глубине h, выше которой скорость постоянна и равна v
(рис. 10). Применение к полю однократно отраженной
волны процедуры прогнозирования путем аналитичес-
ких вычислений в соответствии с (8) приводит к ре-
зультату, имеющему спектральную характеристику
ЦЬ, со) = Гс'2-у2(()1) е -Jkr(a,b),	(П9)
V л(а, Ь)
которая, после применения фильтра ^'(со), кинемати-
чески и динамическими согласуется с характеристикой
а о	зооо	м
Рис. 11. Полнократная волна первого порядка:
а - результат прогнозирования, полученный взаимной свёр-
ткой исходных сейсмограмм; б - волна, полученная с пра-
вильным учетом геометрического расхождения
23
первой полнократной волны. Здесь г(а, Ь) - расстояние
от источника а до приёмника b вдоль траектории луча
кратной волны; с - коэффициент отражения от глубин-
ного горизонта; к = to/v. Аналитический расчет в соот-
ветствии с (7) приводит к иному результату
£(/>, и) = ~СУ(И) да е~^аЬ\ (п Ю)
4/г
который кинематически эквивалентен правильному, но
его динамические отличия уже не описываются однока-
нальным фильтром 5_|(ш). Действительно, амплитуда
волны ЦЬ, со) уже не затухает с ростом пройденного ей
расстояния, но увеличивается. Этот эффект иллюстри-
руется на рис. 11, где сопоставляются результаты чис-
ленных расчетов по формулам (7) и (8), полученные по
модельной сейсмограмме.
ЛИТЕРАТУРА
1.	Авербух А. Г., 1982, Изучение состава и своиств горных пород
при сейсморазведке: М., Недра.
2.	Аки К., Ричардс П., 1983, Количественная сейсмология. Т. 1, 2:
М., Мир.
3.	Алексеев А. С., Гельчинский Б. Я., 1959, О лучевом методе
вычисления полей волн в случае неоднородных сред с криво-
линейными границами раздела: Вопросы динамической тео-
рии распространения сейсмических волн. Ч. III: Л.. ЛГУ.
4.	БреховскихЛ. М., 1957, Волны в слоистых средах: М., Изд-во
АН СССР.
5.	Булатов М. Г., Телегин А. Н., Левый Н. В., 1980, Весовые
коэффициенты при дифракционном преобразовании сейсми-
ческих записей: Разведочная геофизика, 90.
6.	Гурвич И- И., Боганик Г. Н., 1980, Сейсмическая разведка:
М., Недра.
7.	Денисов М. С., 2006, Анализ метода прогнозирования крат-
ных волн без знания модели среды с позиций теории продол-
жения волновых полей: Геофизика, 1.
8.	Математическая энциклопедия, 1984: М., Советская энцик-
лопедия.
9.	Нахамкин С. А., 1969, О веерной фильтрации. Изв. АН СССР.
Физика земли, 11.
10.	Петрашень Г. И.. Нахамкин С. А., 1973, Продолжение вол-
новых полей в задачах сейсморазведки: Л., Наука.
1	VLokshtanov D.. 1995, Multiple suppression by single channel and
multichannel deconvolution in the tau-p domain: 65lh Ann. Internal.
Mtg. SEG.
12.	Matson К. H., Ahma R., 2005, Fast 3D surface-related multiple
elimination using azimuth moveout for multiples: 75th Ann.Internal.
Mtg. SEG.
13.	Robinson E. A., 1957, Predictive decomposition of seismic traces:
Geophysics, 22, 767 - 778.
14.	Taner M. T. 1980. Long-period sea-floor multiples and their
suppression: Geophysical Prospecting, 28, 30 - 48.
15.	Turner G., 1990, Aliasing in the tau-p transform and the removal
of spatially aliased coherent noise: Geophysics, 55, 1496 - 1503.
16.	Van Dedem E. J., Verschuur D. J., 2001, 3D Surface multiple
prediction using sparse inversion: 71th Ann. Internal. Mtg. SEG.
17.	Vidale J. E., Houston H., 1990, Rapid calculation of seismic
amplitudes: Geophysics, 55, 1504 - 1507.
18.	Wapenaar С. P. A., Verschuur D. J., Herrmann P., 1992, Amplit ide
preprocessing of single and multicomponent seismic data: Geophysics.
57. 1178 - 1188.
КОРОТКО ОБ АВТОРАХ
Михаил Сергеевич ДЕНИСОВ - доктор физ.-мат. наук, ведущий математик ООО “Геотехсистем”.
Дмитрий Борисович ФИНИКОВ - кандидат техн, наук, научный сотрудник ООО “Геотехсистем”.
24
М. С. Денисов,
Д. Б. Фиников
ООО "ГЕОТЕХСИСТЕМ", Москва
ООО "ГЕОТЕХСИСТЕМ", Москва
МЕТОДЫ ПОДАВЛЕНИЯ КРАТНЫХ ВОЛН
В СЕЙСМОРАЗВЕДКЕ. ЧАСТЬ 3
АННОТАЦИЯ. В третьей части статьи мы продолжаем рассмат-
ривать двухшаговые методы подавления кратных волн. На этот
раз изучаются 2D- и ЗО-схемы сейсмического эксперимента.
Для построения оператора прогнозирования цуга кратных от-
ражений мы используем анализ волнового поля с привлечени-
ем приближения Кирхгофа. Полученные результаты являются
обобщением описанных ранее подходов. Проводится сравни-
тельный анализ алгоритмов прогнозирования кратных отраже-
ний, выявляются границы их применимости.
ABSTRACT. The third part of the Tutorial paper is devoted to the
analysis of the two-step demultiple schemes. Now we consider 2D
and 3D situations. The Kirchhoff approximation is used for the design
of the multiple prediction operator. The results obtained are
generalization of the approaches described in the previous parts of
the paper.The methods for multiple prediction are compared and
their limitations are discovered.
ВВЕДЕНИЕ. В первой и второй частях статьи1 были под-
робно изучены двухшаговые алгоритмы ослабления крат-
ных отражений, в рамках которых осуществляются про-
гнозирование помех методами прямого продолжения
волнового поля и последующее их адаптивное вычита-
ние из исходных данных. Было показано, что предска-
зывающая деконволюция способна обеспечивать дина-
мически корректные результаты для одномерной моде-
ли сейсмического эксперимента (1D), т. е. для случая
вертикального падения плоской волны на систему плос-
ко-параллельных слоев. Вычитание кратных отражений
в области преобразования Радона способно “обслужи-
вать” и более общую ситуацию. Оно приводит к коррек-
тным результатам для полуторамерной сейсмической
модели (1,5D), т. е. волнового поля, вызванного точеч-
ным или линейным источником, и распространяюще-
гося в горизонтально-слоистой среде.
В третьей части статьи мы сделаем следующий шаг
по пути усложнения структуры сейсмического экспери-
мента и осуществим переход к двумерным (2D) и трех-
мерным (3D) схемам сейсмических наблюдений.
Методы, которые будут рассмотрены, являются обоб-
щениями алгоритмов, описанных в первых двух частях
1 Денисов М. С., Фиников Д. Б.. 2007. Методы подавления крат-
ных волн в сейсморазведке: Технологии сейсморазведки, № 1,
№ 2.
статьи. Их применение зачастую связано с большими
вычислительными сложностями (выдвигаются повышен-
ные требования к объёму памяти и быстродействию
вычислительной техники), что нередко приводит к по-
пыткам решить задачу упрощёнными средствами. Здесь
требуется понимание границ применимости тех или иных
подходов, чтобы обеспечить разумный компромисс меж-
ду затратами на вычисления и точностью получаемого
результата.
2D-, ЗЦ-МОДЕЛИ СЕЙСМИЧЕСКОГО ЭКСПЕРИМЕНТА
В рамках трехмерной сейсмической задачи рассматрива-
ется волновое поле, вызванное точечным источником,
распространяющееся в глубь среды, отражающееся от
границ раздела и регистрируемое сейсмоприёмниками.
Источник колебаний и приёмники расположены на сво-
бодной поверхности или в приповерхностном слое. Сре-
да состоит из границ, имеющих произвольную конфи-
гурацию.
В рамках двумерной задачи накладывается ограниче-
ние на геометрию отражающих горизонтов. Если задать
декартову систему координат (х, у, z) таким образом,
что ось z направлена в глубь среды, а плоскость < = 0
соответствует свободной поверхности, то можно подо-
брать такую ориентацию осей х и у, что строение среды
будет инвариантным относительно у. Условимся пола-
гать, что при профильной сейсморазведке линия наблю-
дений имеет нулевую координату у.
Как и в предыдущих частях статьи, здесь в первую
очередь нас будут интересовать кратные волны, связан-
ные с переотражением в нижнее полупространство от
свободной поверхности. Такие помехи оказываются наи-
более интенсивными. Они всегда доминируют над внут-
ренне-кратными отражениями, т. е. такими кратными
волнами, лучевая схема которых не имеет промежуточ-
ной точки выхода луча на поверхность.
2D-, ЗО-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН
ПРИ НАЛИЧИИ МОДЕЛИ СРЕДЫ
Рассмотрим процесс формирования кратной волны, свя-
занной с переотражением от некоторого глубинного крат-
нообразующего горизонта G, геометрия которого пред-
полагается известной. Кроме того, предположим, что мы
располагаем достаточно надежной оценкой глубинно-ско-
ростной модели среды в интервале глубин от свободной
поверхности до горизонта G. Пусть источник колебаний
расположен в точке с координатой (ха, уа, z = 0), для
которой мы будем использовать сокращенное обозначе-
ние а. Спрогнозируем модель кратной волны, зарегист-
рированной сейсмоприёмником, находящимся в точке Ь,
имеющей координаты (xb, уь, z = 0). Такая волна имеет
промежуточную точку выхода луча на свободную повер-
хность, которую обозначим q (рис. 1).
Здесь для нас важно, что при рассмотрении задачи
3D точка q может иметь произвольное положение на
плоскости наблюдений z = 0. В ситуации 2D она имеет
нулевую у-координату, что обусловлено свойствами сим-
метрии среды и соответствующим выбором ориентации
системы координат.
Как показано на рис. 1, волна, излученная источни-
ком с, распространяется в глубь среды, отражается от
глубинных границ и возвращается к свободной поверх-
ности, где регистрируется сейсмоприёмниками. При этом
Рис. 1. Схема хода луча кратной волны, вызванной
источником, расположенным в точке с координатой а,
прошедшей некоторый путь в среде, переотразившей-
ся от свободной поверхности в нижнее полупрост-
ранство в точке р, отражённой от глубинного гори-
зонта G и зарегистрированной приёмником, распо-
ложенным в точке с координатой b
она претерпевает переотражение в нижнее полупрост-
ранство от поверхности, вновь проходит через среду до
некоторого горизонта G, отражается от него вверх и воз-
вращается вновь к свободной поверхности, где регист-
рируется как кратное отражение.
Процедура прогнозирования кратных волн призвана
воспроизвести этот процесс распространения волны и
по зарегистрированному в точке q сигналу предсказать
его форму и момент вступления при наблюдении в точ-
ке Ь. Понятно, что такое прогнозирование реализуется
методами прямого продолжения волнового поля. В дан-
ной ситуации поле, зарегистрированное на сейсмограмме
общего пункта взрыва а (ОПВ(с)), подвергается действию
оператора миграционного типа, в результате чего полу-
чим сейсмограмму ОПВ(с), содержащую только крат-
ные отражения. Под оператором миграционного типа мы
понимаем пространственный линейный фильтр, осуще-
ствляющий взвешенное суммирование сейсмических
трасс со сдвигами и последующую одноканальную свер-
тку результата со стационарным фильтром. Построение
такого оператора требует привлечения оценки глубин-
но-скоростной модели среды и осуществляется метода-
ми решения прямой динамической задачи.
Следуя логике рассуждений, принятой во второй ча-
сти статьи, будем рассматривать зарегистрированное вол-
новое поле как суперпозицию плоских волн. Такой под-
ход не только удобен для дальнейших рассуждений, но
и придает им наглядность и простоту. Прямое продол-
жение поля осуществляется при помощи линейного пре-
образования. Таким образом, построив алгоритм прямо-
го продолжения для каждой плоской волны как элемен-
та разложения исходного поля, получим окончательный
результат прогнозирования кратных отражений в виде
суперпозиции результатов, полученных по каждой плос-
кой волне.
Кроме того, вопросы, касающиеся построения опе-
ратора прямого продолжения волнового поля, будем рас-
сматривать в рамках приближения Кирхгофа, т. е. будем
полагать, что радиусы кривизны отражающих горизон-
тов намного больше длины волны гармоники, домини-
рующей в сейсмическом импульсе. Следствием этого
является возможность изучения задачи прямого продол-
жения волнового поля в ее локально-плоской аппрок-
симации, а также привлечения лучевых представлений.
Иными словами, при продолжении плоской волны че-
рез среду, глубинно-скоростная модель которой пред-
полагается известной, можно производить прямое про-
должение зарегистрированного волнового поля вдоль
годографа однократно-отражённой волны. Если среда со-
держит криволинейные границы, то на каждой из них
описание акта отражения (преломления) осуществля-
ется при помощи выражений, полученных для плоских
волн. При этом сама граница аппроксимируется плоско-
стью, касательной к ней в точке отражения (преломле-
ния) луча. Как следствие этого, становится корректным
использование таких понятий, как луч,'годограф, ко-
эффициенты отражения и преломления и т. п. Известно,
что привлечение приближения Кирхгофа обусловливает
свойство локальности волновых продолжений, в соот-
ветствии с которым поле, отраженное (преломленное)
на некоторой границе, формируется в небольшой окре-
стности точки отражения (преломления).. К более под-
4
робному обсуждению этого важного свойства мы вер-
немся ниже. С другой стороны, при прямом или обрат-
ном продолжении поля в рамках такой асимптотики со-
храняются все эффекты фокусировки энергии волн.
В качестве обоснования предложенного подхода вы-
делим четыре основных, на наш взгляд, причины.
1.	В подавляющем большинстве реальных сейсмогео-
логических условий, в которых проводятся исследова-
ния, выполняется условие плавности границ, т. е. ради-
ус их кривизны много меньше длины волны гармони-
ки, доминирующей в импульсе.
2.	Даже если среда содержит негладкие элементы,
например, разломы или локальные включения, то, имея
в распоряжении только исходные сейсмические запи-
си, нет возможности восстановить их глубинно-скорос-
тную структуру с точностью, необходимой для коррек-
тного их учета при построении оператора прогнозиро-
вания волнового поля, учитывающего нелучевые эффек-
ты. Дело в том, что подобные объекты порождают отра-
жённые или дифрагированные волны, пространствен-
ный спектр которых насыщен высокими частотами. Тог-
да при сейсмических исследованиях выбор традицион-
ных на сегодняшний день интервалов между пунктами
взрыва и приёма приведёт к тому, что уже в исходных
данных высокие пространственные частоты сигналов
будут зарегистрированы в виде их аляйсинг-эквивален-
тов. Это сделает невозможным их корректный учёт при
дальнейшей обработке.
3.	При 2О-сейсморазведке колебания регистрируют-
ся только вдоль линии профиля наблюдений. Таким об-
разом, мы не располагаем информацией о поведении
волнового поля вне линии профиля. С другой стороны,
любые процедуры продолжения поля, зарегистрирован-
ного на свободной поверхности, требуют учёта всей вол-
новой картины. В простых ситуациях можно продолжить
зарегистрированные волны с линии профиля на всю
площадь. Однако и эти способы имеют свои ограниче-
ния. Казалось бы, площадные схемы наблюдений теоре-
тически должны снимать указанное ограничение. Но на
практике дискретность получаемых площадных сейсмо-
грамм такова, что обусловливаемые ею неточности про-
гнозирования делают бессмысленными все попытки
построить оператор, учитывающий “тонкие” нелучевые
эффекты.
4.	Реальное волновое поле содержит помехи различ-
ной природы, которые, как правило, интенсивны на-
столько, чтобы свести на нет все усилия, направленные
на построение “более точного” оператора.
Вначале рассмотрим 2О-задачу прогнозирования
кратных волн, а затем оговорим особенности обобще-
ния полученного решения на случай 3D. При выводе
оператора будем следовать логике, принятой нами во
второй части статьи, где был изучен процесс прогнози-
рования для среды, состоящей из плоско-параллельных
границ.
Тогда, в соответствии со схемой, показанной на
рис. 1, сформулируем задачу прогнозирования поля крат-
ных волн, связанных с переотражением от выбранного
кратнообразующего горизонта G (привлекая понятия,
заимствованные из лучевой теории), как задачу прямо-
го продолжения исходного волнового поля вдоль луча
однократно отраженной от G волны. Модель среды в
интервале глубин от свободной поверхности до грани-
цы G считается известной. Следовательно, при помощи
трассировки луча из точки с координатой q в точку b
(обе точки расположены на поверхности) можно полу-
чить годограф однократно-отражённой волны, который
мы обозначим через g(q, b).
Пусть зарегистрировано волновое поле, вызванное
воздействием на среду источника, расположенного в
точке а = 0. Обозначим подборку трасс в виде сейсмо-
граммы ОПВ(с) как р(х, t). Предположим, что нам уда-
лось представить это волновое поле в виде его разложе-
ния по плоским волнам
p(x,t) = J p(a,t + ax)da,
- оо
(1)
p(a.,t)~r(t)* £р(л,/-ат)
r(t) - оператор одноканальной компенсирующей филь-
трации. В последнем выражении х - дискретный индекс,
т. е. номер трассы в сейсмограмме, а интервал между
трассами условно положен равным единице Дх = 1. Сим-
вол * обозначает процедуру свёртки по переменной t.
Параметр а характеризует направление подхода плоской
волны к свободной поверхности, он же определяет на-
правленность суммирования (Г). Более подробное осве-
щение вопросов, связанных с разложением (1). (Г), дано
во второй части статьи.
Мы условились полагать, что прямое продолжение
волнового поля в точку приёма с координатой b может
осуществляться путем введения фазового сдвига т(а, Ь)
в каждый элемент разложения исходного поля по плос-
ким волнам. Полученный результат будет кинематичес-
ки адекватен реально наблюдённым помехам. Для их ди-
намического соответствия следует также учесть фактор
с, отвечающий за изменение амплитуды при продолже-
нии локально-плоской волны через среду за счёт актов
преломления и отражения. В более общем случае вместо
параметра с следует использовать фильтр. (Как правило,
этот фактор компенсируется на этапе адаптации спрог-
нозированной модели к реальным сейсмограммам.)
Тогда в области разложения по плоским волнам ал-
горитм предсказания кратных отражений младшего по-
рядка кратности /|(а, г), которые зарегистрированы в
точке приёма с координатой Ь, можно представить в
виде
/1(а,г) = ср(ос,г-т(а,Ь))	(2)
Во второй части статьи было дано правило пересчета
функции т() из g():
Дос, 6) = g(g'~\a.b),b) - ag'~1 (а,Ь),	(3)
где gl(a,b) - функция, обратная к g(q,b), g'~l(a. b) -
производная обратной функции.
5
Для получения динамически корректного результата
прогнозирования следует получить поля кратных волн,
связанных с переотражением со стороны источника и
со стороны приёмника, а также дополнительное поле
волн старшего порядка кратности. Две первых модели
рассчитываются применением оператора прямого про-
должения поля к исходным сейсмограммам в сортиров-
ках ОПВ и ОПП (общего пункта приёма) соответствен-
но. Третья модель вычисляется при помощи прямого про-
должения со стороны источника поля кратных волн
младшего порядка кратности, обусловленных переотра-
жением со стороны приёмника. Также необходимо при-
менить амплитудную коррекцию первой полнократной
волны, отражённой от кратнообразуюшего горизонта, в
соответствии с правилами, описанными во второй час-
ти статьи.
Процедура сохраняет многие особенности, присущие
её 1,5О-аналогу, которые были рассмотрены нами ра-
нее. В частности, реализация алгоритма в области разло-
жения по плоским волнам позволяет учесть зависимость
характеристики отражения от лучевого параметра а. К её
преимуществам следует также отнести и простоту реа-
лизации, так как прогнозирование осуществляется при
помощи ввода фазовых сдвигов (3) в каждую трассу-
элемент разложения исходной сейсмограммы по базису
плоских волн.
Альтернативным является преобразование, реализу-
емое в пространственно-временной области, т. е. приме-
няемое непосредственно к зарегистрированным сейсмо-
граммам (возможно, после их предварительной обработ-
ки). В рамках традиционно привлекаемых для анализа
волнового поля асимптотик оно приобретает форму пре-
образования миграционного типа. Пусть требуется полу-
чить модель кратных волн /|(Л, /), соответствующую по-
ложению источника колебаний а — 0 и приёмника Ь. Для
этого во второй части статьи была использована следу-
ющая последовательность рассуждений. Мы подставили
(2) в (1) и учли выражение (Г). Полученный интеграл
был оценен при помощи метода стационарной фазы. В
рамках локально-плоской аппроксимации задачи 2D
искомая характеристика фильтра в пространственно-
временной области получается воспроизведением этой
последовательности вычислений.
Мы не считаем целесообразным повторять описан-
ные во второй части статьи рассуждения и предложим
альтернативный и не менее наглядный способ расчёта
характеристики фильтра, который нам будет легко обоб-
щить на случай 3D. В соответствии с оговоренным при-
ближением будем искать оператор прогнозирования крат-
ных волн в виде пространственного суммирования трасс
со сдвигами и последующей одноканальной фильтраци-
ей. При этом, как следует из соображений о кинемати-
ческих свойствах преобразования, траектория суммиро-
вания (в случае 3D траектория суммирования переходит
в поверхность суммирования) совпадает с годографом
волны, отражённой от кратнообразуюшего горизонта.
Кроме того, как следует из приближения Кирхгофа,
действие оператора на локально-плоскую волну должно
приводить к её фазовому сдвигу, но не к изменению
амплитуды. Потеря элементарной волной энергии в ак-
тах отражения и преломления на границах среды учиты-
вается отдельно.
При расчёте кратной волны в пространственно-час-
тотной области действие оператора на плоскую волну,
имеющую спектральную характеристику	опи-
сывается выражением
=	^‘^'dq (4)
ч
Оператор преобразования осуществляет взвешенное про-
странственное суммирование исходного сигнала вдоль
траектории g(q, b) и последующую одноканальную
фильтрацию. Веса описываются функцией w(q, b), а спек-
тральная характеристика фильтра соответствует Л(со). Тре-
буется получить выражения для w(q. b) и Я(со). (Заме-
тим, что, в отличие от рассуждений во второй части
статьи, здесь мы пренебрегаем пространственной диск-
ретностью данных.)
Интегральное преобразование (4) соответствует дей-
ствию оператора на “эталонный сигнал” еЛд(*-?). При
этом нами сделаны такие предположения, при которых
известны свойства получаемого решения L^(b, со). Этот
принцип будет положен в основу алгоритма расчёта ве-
совых коэффициентов и характеристики фильтра.
Результат действия интегрального преобразования
может быть оценен при помощи метода стационарной
фазы, обоснование которого производится в той же асим-
птотике, в рамках которой мы осуществляем анализ за-
дачи. Тогда
Lx(b,со)® Л(со)и(q*,b)
ld2g(q,b)
V dq2
V ' ч=ч
 (5)
Здесь q* обозначает координату, удовлетворяющую ус-
ловию стационарности фазы Ф(<7, b) = y(b - q) - g(q, b),
т. е. получаемую как решение уравнения
.=-y-g(/,^ = 0.	(6)
'	10=0
Последнее условие, как нетрудно увидеть, соответству-
ет правилу пересчёта траектории суммирования опера-
тора в функцию фазовых сдвигов (3).
Требование сохранения амплитуды и формы сигнала
локально-плоской волны после преобразования (4) при-
водит к тому, что, как следует из (5), следует выбрать
У2л	\ дс1
Выражения согласуются с характеристикой фильтра и
весовыми коэффициентами, полученными во второй
части статьи при помощи более строгих рассуждений.
Причина этого, по-видимому, заключается в том, что
привлечение лучевых представлений, справедливых в
рамках высокочастотной асимптотики (и совпадающе-
го с ней по смыслу приближения Кирхгофа), позволя-
ет анализировать интегральный оператор продолжения
волнового поля при помощи метода стационарной фазы.
6
Иначе говоря, фактически осуществляется замена ре-
зультата интегрирования на значение подынтегральной
функции в одной точке (или нескольких изолирован-
ных точках), удовлетворяющей условию стационарно-
сти фазы. Анализируя такой подход с точки зрения те-
ории распространения упругих волн, можно дать сле-
дующую интерпретацию математических манипуляций,
которые были нами использованы. Если волновое поле
известно на некоторой замкнутой поверхности и тре-
буется продолжить его в точку, расположенную внутри
объёма, ограниченного этой поверхностью, то соот-
ветствующие расчёты производятся при помощи ин-
тегральных преобразований, аналогичных использован-
ным нами выше. То есть осуществляется взвешенное
интегрирование со сдвигами заданного поля по грани-
це области. Известно, что не все точки обеспечивают
одинаковый вклад в результат, но основной вклад обус-
ловлен стационарными точками и некоторой их окрест-
ностью. (Например, такое же накапливание осуществ-
ляется в процессе сейсмической миграции, а стацио-
нарные точки определяются как точки касания траек-
тории суммирования и годографа отражённой волны.)
Масштаб окрестности определяется размером зоны
Френеля, которая пропорциональна длине волны гар-
моники, доминирующей в спектре сейсмического им-
пульса. Мы условились полагать, что при анализе поля
справедливы высокочастотные приближения. Следова-
тельно, размер зоны Френеля стремится к нулю, и зна-
чение интеграла может быть получено как результат
“интегрирования” по бесконечно малой окрестности
стационарной точки. Понятно, что при этом эффекты,
обусловленные пространственной дискретностью, иг-
норируются.
Проиллюстрируем этот принцип на примере пря-
мого продолжения поля плоской волны от свободной
поверхности до кратнообразующего горизонта G и об-
ратно к свободной поверхности. На рис. 2 схематически
показан фронт плоской волны, распространяющейся
от поверхности в нижнее полупространство. Волна дос-
тигает границы G, отражается от неё и возвращается к
поверхности в виде кратного отражения. В точке b рас-
положен приёмник, для которого следует получить поле
кратной волны. Соответствующий расчёт подразумева-
ет интегрирование вдоль всей кривой G. Однако, если
кратнообразующая граница плавная, т. е. если в каждой
точке радиус её кривизны намного превышает длину
волны гармоники, доминирующей в спектре импуль-
са, то вкладом в результат интегрирования от всех уча-
стков, не обозначенных красным цветом, можно пре-
небречь. На этих локальных участках “работают” луче-
вые формулы, в том числе закон Снеллиуса, и можно
привлекать понятие коэффициента отражения плоской
волны, рассматривая процесс отражения от криволи-
нейной границы как от плоскости, касательной к ней
в точке отражения. В этом и заключается смысл исполь-
зуемого нами для вывода оператора прогнозирования
приближения Кирхгофа. Размер областей определяется
масштабом зоны Френеля, и в рамках высокочастот-
ной асимптотики стремится к нулю. Следовательно, и
сам интеграл в данной ситуации также определяется
как сумма значений подынтегральной функции в ста-
ционарных точках, и обоснование метода стационар-
Рис. 2. Принцип продолжения поля плоской волны
с целью прогнозирования кратных отражений, обус-
ловленных переотражением от свободной поверхно-
сти и кратнообразующего горизонта G:
для наглядности в качестве G здесь показан первый (верх-
ний) глубинный отражающий горизонт; координата b соот-
ветствует положению приёмника трассы, для которой про-
изводится расчёт; красным цветом условно показана об-
ласть границы G, вклад от которой в результат интегриро-
вания вдоль границы существен; её размер определяется
масштабом зоны Френеля и стремится к нулю, если при-
влекать высокочастотное приближение при анализе волно-
вого поля
ной фазы осуществляется в тех же предположениях, в
которых обосновано приближение Кирхгофа. В соответ-
ствии с геометрией задачи, показанной на рис. 2, та-
ких точек две. Иначе говоря, уравнение (6) имеет два
решения g*t и <7*2-
Рассмотренная нами ситуация также проясняет тот,
казалось бы, парадоксальный факт, что, несмотря на
пренебрежение кривизной границ и использование “ло-
кально-плоских” формул, мы правильно учитываем эф-
фекты фокусировки при обработке поля каждой отдель-
ной плоской волны. В данной ситуации фокусировка энер-
гии полечена за счёт учёта всех точек отражения от син-
клинального участка границы. Заметим также, что, хотя
на среду, состоящую из системы слоёв (в данном случае
одна, но криволинейная граница), падает только одна
плоская волна, отражённое поле, полученное при помо-
щи “локально-плоских” формул, содержит целый набор
плоских волн (на рис. 2 в точку регистрации приходят два
луча с различными лучевыми параметрами).
Получив выражения для одноканального компенси-
рующего фильтра и весовых множителей, мы можем
выписать оператор прогнозирования кратных волн по
исходной сейсмограмме ОПВ p(q, /) в пространственно-
временной области
Г
/](/;,/)-<? /(/)*
в
Ч=А
p(b-q,t~g(q,b))w(q,b)
(7)
7
Накапливание осуществляется в пределах апертуры [А. 5].
а весовые множители равны корню из второй производ-
ной траектории суммирования
d2g(g,b)
dq~
фильтр r(t) имеет спектральную характеристику
V 2л
Если производится обработка материалов 3D-ceiic-
моразведки. то в выражении (4) одномерное суммиро-
вание заменяется на суммирование по двум координа-
там, т. е. по поверхности наблюдений. ‘‘Эталонная” плос-
кая волна заменяется ее пространственным аналогом.
Тогда если под переменной интегрирования понимать
пространственную координату точки q = (qx. q^. то.
применяя к полученному двукратному интегралу дву-
мерное обобщение метода стационарной фазы, полу-
чим выражения для спектра фильтра r(f) в виде
/?(«)= ^ 4
2л
а для весовых множителей
w( д,Ь) =
где
-
dqxdqv 
Таким образом, веса равны корню из определителя мат-
рицы вторых производных (матрицы Гессе) поверхнос-
ти суммирования. Стационарная точка определяется из
условия равенства нулю градиента двумерной фазовой
функции
= 0,
ИЛИ
5Ф(с/,6)	() дФ(</,6)
5г/,.	*	’	5г/
л <7л=<7л
Можно аналогично тому, как мы это сделали во вто-
рой части статьи, показать, что выражение (7) соответ-
ствует дискретному аналогу интегрального оператора
прямого продолжения волнового поля (иногда также
называемому интегралом Кирхгофа). Однако при рас-
смотрении 2D- и 3D-задач это потребовало бы нс толь-
ко необоснованного увеличения объёма статьи, но и
насыщения её большим числом формул, что входит в
противоречие с заявленным нами вначале принципом
написания популярной работы.
Мы рассмотрели вопросы, связанные с предсказа-
нием цуга кратных волн, образованных со стороны при-
ёмника. Были получены алгоритмы прогнозирования
волн, связанных с переотражением от свободной по-
верхности и глубинного кратнообразующего горизонта,
применяемые к исходным сейсмограммам в сортировке
ОПВ. Реальное поле также содержит кратные волны,
переотраженные со стороны источника. Такой тип по-
мех прогнозируется при помощи аналогичных преобра-
зований, применяемых к сейсмограммам в сортировке
ОПП.
Подводя итог описанию алгоритма подавления крат-
ных волн с использованием модели среды, реализуемо-
го в области разложения по плоским волнам, предста-
вим его в виде последовательности процедур обработки
исходного поля:
Шаг 1. Применение преобразования Радона к сейс-
мограммам ОПВ и ОПП. Тем самым осуществляется пе-
реход из пространственно-временной области в область
разложения по плоским волнам.
Шаг 2. Выбор кратнообразующего горизонта и оце-
нивание кинематических параметров волны, однократ-
но отражённой от этого горизонта. Расчёт функции сдви-
гов т(а, Ь).
Шаг 3. Прогнозирование кратных волн путём сдвига
трасс сейсмограмм в области разложения Радона.
Шаг 4 Получение дополнительной модели кратных
волн старшего порядка кратности.
Шаг 5. Возможно применение мьютинга (обнуления
областей, занятых исключительно помехой), если пос-
ле преобразования Радона сигналы и помехи разделя-
ются.
Шаг 6. Адаптация и вычитание.
Шаг 7. Обратное преобразование Радона результата
адаптивного вычитания. Переход в пространственно-вре-
менную область.
Альтернативным способом прогнозирования цуга
кратных волн является использование преобразования
миграционного типа (7), реализуемого в пространствен-
но-временной области.
Также отметим, что если среда состоит из плоско-
параллельных границ, то выражение (7) переходит в его
упрощенный аналог, полученный нами для ситуации
1,5D. Если, кроме того, источник излучал бы плоскую
волну, что соответствовало бы задаче 1D, то (7) осуще-
ствляло бы одноканальное преобразование вида сдвига
трассы, т. е. прогнозирование кратных волн, реализуе-
мое в алгоритме предсказывающей деконволюции.
2D-, ЗО-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН
БЕЗ ПРИВЛЕЧЕНИЯ АПРИОРНОЙ ИНФОРМАЦИИ
О СТРУКТУРЕ СРЕДЫ
Одним из наиболее популярных и наглядных способов
получения трассы, содержащей как однократные, так и
многократные отражения в рамках ID-задачи является
процедура Баранова-Кюнетца. Алгоритм преобразования
допускает следующую интерпретацию. Пусть в соответ-
ствии с заданной априорной информацией о строении
среды мы смоделировали трассу, содержащую только
однократные отражения. Тогда модель кратных волн,
обусловленных переотражением в нижнее полупростран-
ство от свободной поверхности, первого порядка крат-
8
ности может быть получена по трассе однократных волн
уже без привлечения информации о модели среды, а
оператором построения последовательности кратных
отражений по однократным волнам является сама трас-
са однократных волн. Преобразование осуществляется
путём автосвёртки этой трассы.
Можно обобщить алгоритм Баранова-Кюнетца для
решения задач большей размерности - 2D и 3D. В соот-
ветствии с изложенными принципами построение мо-
дели кратных волн первого порядка кратности, связан-
ных с переотражением в нижнее полупространство от
свободной поверхности, может осуществляться по "сей-
смограммам однократных отражений”, при этом в ка-
честве оператора преобразования используются эти же
сейсмограммы. Вместо одноканальной свёртки приме-
няется пространственно-временная свёртка.
Сказанное верно и для расчёта поля кратных волн
старших порядков кратности: действуя на сейсмограм-
мы, содержащие однократные и кратные волны вплоть
до к-го порядка кратности оператором, представляю-
щим собой набор “сейсмограмм однократных волн”,
получим модель кратных волн вплоть до (к + 1)-го по-
рядка. Таким образом, динамические годографы одно-
кратных отражений являются операторами прогнозиро-
вания кратных волн.
Когда мы рассматривали способы моделирования
кратных отражений с привлечением априорной инфор-
мации о модели среды, то фактически мы решали пря-
мую задачу с целью получения динамического годо-
графа однократной волны. Результат расчёта использо-
вался при построении оператора прямого продолже-
ния поля. Теперь мы замечаем, что необходимую ин-
формацию можно заимствовать непосредственно из за-
регистрированных данных. Действительно, как следует
из рис. 1, лучевая схема любой кратной волны, имев-
шей хотя бы одно отражение в нижнее полупростран-
ство от свободной поверхности (в данном случае это
произошло в точке с координатой q), состоит из двух
сегментов - пути, пройденного волной из точки а в точку
q и из точки q в b. Этим сегментам соответствуют две
волны (однократные или кратные, на рис. 1 изображе-
ны однократные отражения), которые были зарегист-
рированы в исходных данных. Тогда самый непосред-
ственный способ получить кратную волну (с, q, b) из
волн (a, q) и (q. b) - это как бы “сложить” их, а такое
"сложение” осуществляется взаимной свёрткой трасс,
на которых зарегистрированы волны. То есть, если волна
подошла к поверхности и отразилась от неё в нижнее
полупространство, то трассу, зарегистрированную в
точке отражения, можно рассматривать как функцию
вторичного источника (принцип Гюйгенса). Тогда, если
запись в точке приёма была реакцией на точечный им-
пульсный источник, то, чтобы получить кратную вол-
ну, надо свернуть трассу с этой функцией источника.
Отсюда приходим к алгоритму, основанному на вза-
имных свёртках трасс исходного волнового поля. Но,
так как модель среды неизвестна, то неизвестна и точ-
ка q выхода луча на поверхность, благодаря отраже-
нию от которой образована и зарегистрированная крат-
ная волна. Тогда мы переберем все возможные точки
отражения в пределах некоторого интервала (аперту-
ры) и суммируем результаты взаимных свёрток. Так
осуществляется расчёт модели кратных волн без при-
влечения информации о глубинно-скоростном строе-
нии среды. Метод не требует привлечения процедур
обнаружения и идентификации волн на сейсмограм-
мах.
Разумеется, кратные волны рассматриваемого типа
являются реакцией среды не на совокупность вторич-
ных точечных источников с характеристиками, совпа-
дающими с зарегистрированными колебаниями, а на
всё, вообще говоря, непрерывное трехмерное поле, ко-
торое “принимается” и отражается всей земной поверх-
ностью. Для получения искомого корректного результа-
та - кратных волн всех порядков кратности, такую трёх-
мерную характеристику следовало бы сворачивать саму
с собой бесконечное число раз. Кроме того, зарегистри-
рованное поле должно было бы быть реакцией на иде-
альный импульсный точечный источник (типа дельта-
функции). Невыполнение этих условий приводит к не-
обходимости внесения поправок в алгоритм прогнози-
рования. Несомненным преимуществом метода является
то, что вид волнового уравнения, которое, фактичес-
ки, мы неявно решали в предыдущем параграфе, здесь
совершенно не важен. Оператор прогнозирования, ко-
торым по построению является зарегистрированное поле,
“всё знает” про эффекты рассеяния, поглощения, рас-
хождения, AVO и т. п. Однако для реализации этих по-
тенциальных возможностей требуется соблюдение мно-
гих, зачастую невыполнимых, условий.
Как показано во второй части статьи, упомянутая
нами выше амплитудная поправка, корректирующая тот
эффект, что свободная поверхность является не набо-
ром вторичных точечных источников, а "излучающей
площадкой”, равна косинусу угла выхода луча на по-
верхность. отнесенному к скорости распространения
волны в непосредственной близости от поверхности.
Тогда исходное волновое поле следует рассматривать
как два набора данных, первый из которых содержит
волны, подлежащие продолжению для прогнозирова-
ния кратных отражений. Второй набор представляет со-
бой оператор прямого продолжения. С этой целью ис-
ходные сейсмограммы удобно пересортировать по ОПП
и ОПВ. Затем к одному из наборов следует применить
предварительную обработку с целью учета амплитуд-
ной поправки. Для этого сейсмограммы надо перевести
из пространственно-временной области в область, где
волны разделены по углу выхода луча на поверхность
(например, область двумерного преобразования Фурье
или область разложения Радона). Введя поправку, про-
изводим переход обратно в пространственно-времен-
ную область.
Здесь, для удобства, мы несколько изменим приня-
тые выше обозначения. Пусть р(а, q, f) - исходная сей-
смограмма ОПВ(с). а р(а, q. t) - сейсмограмма, где
учтена поправка. Тогда модель кратных отражений млад-
шего порядка кратности получим как
/1(«,/?,/) = ^^(ц,б7,О*М9,/’Т),	(8)
ч
где p(q, b, I) - подборка трасс исходного поля в сорти-
ровке ОПП(й). Модель волн к-го порядка кратности вы-
9
числяется подстановкой результата прогнозирования
b, t) в (8).
Для построения динамически корректной модели крат-
ных отражений необходимо скомпенсировать эффект ав-
тосвёртки импульса, который в результате преобразова-
ния (8) дважды учитывается в результирующей трассе.
Если имеется оценка формы сигнала s(t), то, рассчитав
оператор обратной фильтрации можно применить
его к трассе /|(с, b, t). Если форма сигнала неизвестна, то
коррекция эффекта автосвёртки осуществляется на этапе
адаптации модели к реальным волнам-помехам путём
подбора формирующего фильтра. Отметим особенность
такой процедуры: получаемый фильтр может интерпре-
тироваться как оценка импульса, переменная по простран-
ственной координате (и, возможно, по временной коор-
динате). Тем самым этап адаптации модели кратных волн,
полученной без привлечения информации о модели сре-
ды, является одним из вариантов получения оценки фор-
мы импульса. При этом нет необходимости привлекать
традиционные для этой задачи статистические гипотезы,
например, описывать последовательность коэффициен-
тов отражения как реализацию случайного процесса типа
белого шума или делать предположения о фазовом спек-
тре сигнала. Однако форма импульса должна плавно ме-
няться по латерали, а требования к качеству исходных
данных повышаются.
Наши предыдущие работы, посвящённые отдельным
аспектам проблемы подавления кратных волн в сейсмо-
разведке, на которые мы ссылаемся в Примечаниях, со-
держат достаточное количество иллюстраций, демонст-
рирующих работоспособность различных способов. Поэто-
му здесь мы приведем небольшое число рисунков, пояс-
няющих наиболее характерные особенности алгоритмов
при их использовании для обработки материалов, полу-
ченных в 2D- и ЗО-ситуациях. Кроме того, наш опыт об-
работки реальных данных свидетельствует о том, что на
практике предпочтительным является применение про-
гнозирования, не привлекающего оценку модели среды.
Поэтому ниже мы проиллюстрируем возможности этого
подхода при обработке материалов 2D- и ЗО-сейсмичес-
ких наблюдений, а о его преимуществах мы скажем ниже,
когда, по традиции, будем проводить сравнительный ана-
лиз предложенных вычислительных схем.
На рис. 3 показан пример волнового поля, представ-
ляющего собой интерференцию однократных и кратных
отражений. Материал был получен при 2О-морской сей-
сморазведке и представляет особый для нас интерес, так
как исследуемая среда не может быть аппроксимирована
горизонтально-слоистой моделью, и при прогнозирова-
нии кратных отражений следует ожидать существенного
“2Е-эффекта”. Изображённый фрагмент содержит одно-
кратно-отраженную от морского дна волну (указана стрел-
кой А), первую полнократную волну В, связанную с пе-
реотражением от морского дна, интенсивную однократ-
ную волну С, отражённую от неглубокого горизонта,
внутренне-кратную волну D, имеющую дополнительный
пробег в слое, заключенном между морским дном и упо-
мянутым горизонтом, два цуга частично-кратных волн Е
и F - соответственно, первого и второго порядков крат-
ности, связанные с переотражениями от морского дна и
горизонта С. Выделяется также интенсивное однократное
отражение G и другие волны.
В результате применения процедуры прогнозирова-
ния было получено волновое поле, которое мы сопос-
тавляем с исходными данными на рис. 4.
На рис. 5 показан небольшой фрагмент разреза рав-
ных удалений, позволяющий детально проанализировать
волновую картину. В данном случае морское дно имеет
сложную геометрию. Как следствие этого, связанная с
ним первая полнократная волна обладает сложными
кинематическими характеристиками. На фрагменте вол-
нового поля, показанном на рис. 5, а, полнократная
Рис. 3. Разрез равных удалений (200 м), полученный
в результате инженерных сейсмических наблюдений
в прибрежной мелководной части Каспийского моря:
см. пояснения в тексте
Рис. 4. Фрагменты разреза равных удалений (200 м):
слева показано исходное волновое поле, справа - модель
кратных волн, полученная прогнозированием
10
a 19500	21000	22500 м
Рис. 5. Фрагменты разреза равных удалений (200 м):
а - исходное волновое поле; б - модель кратной волны,
полученная прогнозированием; в - результат адаптивного
вычитания модели из исходного поля
б
775
800
825
850
875
900
925
950
975
1000
1025
1050
1075
1100
1125
1150
1175
1200
мс
200	400	600	800 м
Рис. 6. Фрагменты сейсмограмм ОГТ:
а - исходное волновое поле; б - модель кратной волны,
полученная прогнозированием; в - результат адаптивного
вычитания модели из исходного поля
И
волна интерферирует с однократными отражениями. Как
следует из рис. 5, б, удалось хорошо предсказать кине-
матические и динамические особенности этой волны. На
рис. 5, в представлен результат вычитания модели крат-
ного отражения из исходного поля. Очевидно, что по-
меха существенно ослаблена, при этом полезные сигна-
лы сохранили свои амплитуды.
На рис. 6 представлен этот же набор данных в сорти-
ровке ОГТ. На фрагменте исходного волнового поля,
показанном на рис. 6, а, кратные волны интерферируют
с однократными отражениями. На рис. 6, б изображено
волновое поле, полученное прогнозированием, а на
рис. 6, в - результат вычитания.
Прогнозирование кратных волн по материалам, по-
лученным в результате ЗО-сейсмичсских исследований,
наталкивается на сложности, которые непреодолимы в
рамках используемых в настоящее время схем наблюде-
ний. С одной стороны, как отмечалось выше, для кор-
ректного прогнозирования кратных волн необходимо
учитывать всю энергию, переотразившуюся от поверх-
ности в нижнее полупространство. Следствием этого яв-
ляется неуниверсальность профильных методов сбора
сейсмической информации. С другой стороны, дискрет-
ность получаемых при современных способах наблюде-
ний площадных сейсмограмм такова, что обусловливае-
мые ею неточности прогнозирования делают бессмыс-
ленными все попытки построить более точную модель
кратных отражений. В некоторых случаях выигрыш в ка-
честве модели может быть получен за счет интерполя-
ции площадных сейсмограмм, хотя такой подход неиз-
бежно связан с повышенными и зачастую неоправдан-
ными вычислительными затратами.
Одним из примеров удачного применения такого под-
хода является ситуация, иллюстрируемая на рис. 7. Об-
щий вид волновой картины представлен на рис. 7, а, где
показан разрез равных удалений, полученный по одному
сечению площадных сейсмограмм. Материал получен в
глубоководной части района исследований, и время вступ-
ления однократной волны - около 2500 - 750 мс. Первая
полнократная волна имеет время вступления около 5000 -
5500 мс, за ней следует цуг частично-кратных отражений.
Специальные исследования показали, что морское дно
в о 5000	10000	15000	20000 и
5000
5200
5400
5600
5800
6000
6200
6400
6600
6800
МС
Рис. 7. Результаты обработки
реальных ЗО-данных, полученных
при морской сейсморазведке:
показаны фрагменты разреза равных
удалений (725 м); а - исходное вол-
новое поле; б - увеличенный фраг-
мент исходного волнового поля; в -
результат прогнозирования кратных
волн при помощи алгоритма 2D; г -
результат прогнозирования кратных
волн при помощи алгоритма 3D
12
обладает сложной геометрией и имеет существенный на-
клон и кривизну в направлении, ортогональном клиник
профиля. Результат прогнозирования кратных отражений,
полученный при помощи 2Б-алгоритма, т. е. на вход ко-
торого подавались лишь сечения в линии профиля пло-
щадных сейсмограмм, показан на рис. 7, в в увеличен-
ном масштабе. Этот же фрагмент исходного волнового поля
изображен на рис. 7, б. Очевидно, что волновые картины
различны по динамическим и кинематическим характе-
ристикам. и. как следствие этого, адаптивное вычитание
модели из исходных трасс может вызвать ослабление энер-
гии сигналов. Другая модель кратных волн была получена
непосредственно по исходным площадным сейсмограм-
мам при помощи ЗО-прогнозирования. Результат (здесь
не приводится) оказался хуже показанного на рис. 7, в.
Волновое поле содержало интенсивные артефакты, обус-
ловленные аляйсинг-эффектом. Они маскировали спрог-
нозированные кратные отражения. Кроме того, не наблю-
далось улучшения динамической и кинематической адек-
ватности смоделированных волн реальным помехам. Пос-
ле применения специальных способов интерполяции пло-
щадных сейсмограмм нами была получена модель крат-
ных отражений, которая представлена на рис. 7, г. Можно
отметить определенное улучшение кинематической и ди-
намической корректности спрогнозированных волн-по-
мех. При этом, правда, вычислительные затраты возрос-
ли в сотни раз по сравнению с 20-про1нозированием.
Подводя итог описанию алгоритма подавления крат-
ных волн без привлечения априорной информации о
модели среды, представим его в виде последовательно-
сти процедур обработки исходного поля:
Шаг 1. Пересортировка исходного набора трасс по
ОПВ и ОПП. Предварительная обработка одного набора
сейсмограмм с целью введения поправки за угол выхо-
да луча на поверхность.
Шаг 2. Суммирование взаимных свёрток трасс сейс-
мограмм ОПВ и ОПП.
ШагЗ. Получение дополнительных моделей кратных
волн старшего порядка кратности.
Шаг 4. Адаптация и вычитание.
Если сейсмические данные получены в горизонталь-
но-слоистой среде, то преобразование (8) переходит в
его 1,5О-аналог. реализуемый по сейсмограммам ОГТ.
Если волновое поле зарегистрировано в результате
1 D-сейсмического эксперимента, то алгоритм (8) пе-
реходит в его одномерный аналог - способ ретрокорре-
ляции, реализуемый по одной трассе.
КАКОЙ АЛГОРИТМ ПРИМЕНЯТЬ НА ПРАКТИКЕ?
Выбор между двумя методами прогнозирования крат-
ных волн производится исходя из специфики конкрет-
ного материала, а также особенностей сейсмогеологи-
ческих условий.
Несомненным преимуществом прогнозирования, не
привлекающего оценку модели среды, является то что,
как правило, кратные волны следуеКподавить до этапа
построения модели. Исходные сейсмограммы необходи-
мо очистить от нерегулярных и регулярных помех, а за-
тем, по результатам обработки, решать обратную задачу.
Кроме того, так как сейсмограммы содержат волны (ко-
торые в данном случае являются операторами прогно-
зирования), связанные с отражением от всех глубинных
границ, то за одну итерацию (8) мы получим модель
всех кратных волн всех порядков кратности. Эти сообра-
жения являются основными, по которым мы рекомен-
дуем этот подход как наиболее удобный и гибкий.
Кроме того, способ, привлекающий информацию о
модели среды, имеет следующий недостаток. Расчет опе-
ратора прогнозирования осуществляется методами ре-
шения прямой динамической задачи. При этом стано-
вится важным не только выбор вида волнового уравне-
ния, которое предстоит решить, но и та асимптотика, в
рамках которой это решение будет получено. Все это
может привести к неточным результатам, что, в свою
очередь, в дальнейшем отрицательно скажется как на
кинематических, так и на динамических характеристи-
ках модели поля кратных волн. Метод, не привлекаю-
щий информацию о модели среды, неявным образом
заимствует решение прямой динамической задачи не-
посредственно из самих же зарегистрированных сейс-
мограмм. Понятно, что здесь отсутствуют описанные
сложности. В то же время, если в глубинном интервале
от земной поверхности до кратнообразующей границы
среда “простая”, то расчёт прямой задачи не вызывает
принципиальных затруднений, и соответствующий ал-
горитм отличается хорошей работоспособностью.
В алгоритме, не использующем модель среды, опера-
тором продолжения волнового поля являются сами же
зарегистрированные сейсмограммы. Следовательно, они
должны быть “чистыми”. Например, для данных, полу-
ченных при морской сейсморазведке, оператором про-
должения поля через водный слой является однократно
отражённая от морского дна волна. Но если мощность
водного слоя незначительна, то эта волна будет интер-
ферировать с прямой волной, распространяющейся по
поверхности воды, и с головной волной, образующей-
ся на морском дне. Помехи такого типа должны быть
аккуратно подавлены, но на практике все соответству-
ющие процедуры так или иначе искажают динамику
однократной волны. Это свойство делает алгоритм про-
гнозирования кратных волн без модели среды малопри-
годным для сейсмограмм, полученных на мелководье
или при наземной сейсморазведке.
Проблема подавления кратных волн особенно важна
при обработке материалов морской сейсморазведки,
когда требуется ослабить цуг ревербераций, обусловлен-
ных переотражением волн от морского дна и свободной
поверхности. Однако кратные отражения могут быть столь
интенсивными, что будут маскировать интерферирую-
щие с ними сигналы и при наземной сейсморазведке.
Поэтому применение двухшаговых схем к наземным дан-
ным также актуально, а специфика обработки требует
отдельного рассмотрения. Дело в том, что сейсмограм-
мы всегда содержат весьма интенсивные поверхностные
волны, энергия которых локализована в области ближ-
них удалений. С помехами этого типа не приходится стал-
киваться при обработке глубоководных морских данных.
Цуг поверхностных волн не порождает кратных отраже-
ний, следовательно, он должен быть удален до этапа
прогнозирования.
Также необходимо устранить короткие и протяжен-
ные выбросы, т. е. такие участки трасс, энергия которых
13
существенно превышает средний уровень, иначе после
применения операторов пространственного суммирова-
ния они будут "размазаны”. Часто результаты регистра-
ции требуют интерполяции с целью получения данных
на регулярной и достаточно плотной сети наблюдений.
Важным подготовительным этапом является коррекция
за поверхностные условия, в частности - поверхностно-
согласованная деконволюция. Следует подавить регуляр-
ные помехи, под которыми здесь мы понимаем волны
различной природы, не порождающие цуга кратных от-
ражений. При такой предварительной обработке следует
избегать применения процедур, нарушающих баланс
амплитуд, типа выравнивания средней энергии в сколь-
зящем окне (АРУ) или коррекции расхождения. Как и
при применении любой процедуры, основанной на ис-
пользовании операторов миграционного типа, следует
рассмотреть возможность интерполяции сейсмограмм с
целью ослабления аляйсинг-помех.
Отличительной особенностью задачи моделирования
кратных отражений в 2D- и ЗО-средах является необхо-
димость формирования “полных” сейсмограмм ОПВ и
ОПП, т. е. содержащих трассы как положительных, так и
отрицательных удалений. Дело в том, что промежуточ-
ная точка выхода луча кратной волны q может находить-
ся в области отрицательных удалений, которые не реги-
стрируются при фланговой системе наблюдений. Такая
си туация иллюстрируется на рис. 8. Формирование трасс
отрицательных удалений сейсмограммы ОПВ(я) осуще-
ствляется путем подстановки на их место трасс положи-
тельных удалений сейсмограммы ОПП(а), что соответ-
ствует сейсмическому принципу взаимности.
Важным этапом подготовки данных также является
экстраполяция сейсмограмм с целью формирования
трасс ближних удалений.
При обработке материалов 2О-сейсморазведки на
динамическую адекватность получаемой модели реаль-
но зарегистрированным помехам может существенным
образом сказываться эффект рассеяния энергии волн вне
линии профиля, где она не регистрируется. В то же вре-
Рис. 8. Пример хода луча кратной волны, при кото-
ром промежуточная точка выхода на свободную по-
верхность для сейсмограммы ОПВ(а) находится в
области отрицательных удалений
мя для корректного прогнозирования её необходимо
учитывать при прямом продолжении поля до кратнооб-
разующей границы и обратно к поверхности. В лучшем
случае это приводит к несоответствию амплитуд смоде-
лированных помех реальным волнам, причём такое рас-
хождение носит нестационарный характер даже в пре-
делах одной трассы, что затрудняет последующий этап
вычитания и требует разработки специальных устойчи-
вых средств адаптации. Большие затруднения может выз-
вать кинематическое несоответствие модели реальным
волнам, которое будет наблюдаться при прогнозирова-
нии по результатам 2О-сейсморазведки в ЗО-среде. В ус-
ловиях интерференции помех и сигналов кинематичес-
кие погрешности могут приводить к эффекту подавле-
ния однократных отражений на этапе адаптивного вы-
читания.
Казалось бы, в подобной ситуации на помощь могут
прийти методы площадных наблюдений 3D. Но их пре-
имущества при ближайшем рассмотрении оказываются
не такими бесспорными. Использование соответствую-
щих интегральных преобразований требует разбивки
плоскости наблюдений на регулярные сегменты весьма
малого размера, причем необходимо наличие в каждом
таком сегменте как приемника колебаний, так и источ-
ника. На практике это требование, как правило, не удов-
летворяется: данные, полученные в результате назем-
ной сейсморазведки, крайне нерегулярны, имеют зна-
чительное количество пропусков как пунктов взрыва,
так и пунктов приёма, а морские данные фактически
являются набором независимо полученных линейных
профилей, и такая геометрия наблюдений не позволяет
производить адекватное прогнозирование кратных волн.
Отмеченная специфика материала приводит к тому, что
сводятся на нет почти все потенциальные преимуще-
ства пространственных наблюдений.
Следствием этого являются невозможность прогно-
зирования определенного набора кратных волн, а также
значительная зашумленность получаемых моделей. Дей-
ствительно, в силу того, что при продолжении поля
применяются операторы миграционного типа, форми-
рование сигнала в результирующей трассе происходит
за счет касания годографа зарегистрированной волны
линией суммирования оператора. Именно небольшая
окрестность точки касания и определяет как кинемати-
ку, так и динамику прогнозируемого сигнала. Если сфор-
мированные сейсмограммы ОПВ и ОПП существенно
нерегулярны, то велика вероятность того, что в иско-
мой точке касания будут отсутствовать трассы. Зашум-
лённость получаемых моделей кратных волн при обра-
ботке данных 3D обусловлена аляйсинг-эффектом, осо-
бенно сильно проявляющимся как раз в случае нерегу-
лярной и редкой сети наблюдений.
В такой ситуации единственным способом улучшения
качества прогнозируемого поля является интерполяция
исходных данных, причем интерполяция может присут-
ствовать неявно, т. е. производиться в процессе моделиро-
вания кратных волн. Однако при ЗО-сейсморазведке очень
существенно ограничение на объём материала, и всегда
приходится искать компромисс между плотностью ин-
терполируемых данных и точностью вычислений. Кроме
того, несмотря на обилие алгоритмов, в силу существен-
ной неопределенности и некорректности самой задачи
14
интерполяции, особенно в случае разведки сложнопост-
роенных сред, не существует вполне удовлетворительных
и универсальных способов ее решения. Таким образом,
мы приходим к выводу, что существуют принципиаль-
ные и пока непреодолимые погрешности при прогнози-
ровании кратных волн. Иначе говоря, если получить са-
мые точные формулы для расчета поля кратных волн,
то, будучи примененными к реальным данным с их ре-
альной геометрией наблюдений, они, как правило, всё
равно не обеспечат удовлетворительного результата. По-
этому, сделав всё возможное на этапе прогнозирования,
следует сосредоточить усилия на совершенствовании ал-
горитмов адаптивного вычитания модели кратных волн
из зарегистрированного поля. Эти вопросы мы осветим в
следующей части статьи.
Сочетая оба метода прогнозирования и используя
сильные стороны и преимущества каждого из них, можно
развить новый подход к решению задачи подавления
кратных волн. Пусть требуется получить кинематически
и динамически адекватную модель всего цуга кратных
волн, связанных с отражением в нижнее полупростран-
ство от земной поверхности, но специфика зарегистри-
рованных данных не позволяет это сделать с помощью
традиционных подходов. Как упоминалось выше, с та-
кой ситуацией приходится иметь дело, например, при
морской сейсморазведке на мелководье. В этом случае
однократно отражённая от дна моря волна сильно за-
шумлена, но именно она играет определяющую роль в
предсказании, фактически являясь оператором прогно-
зирования наиболее интенсивных волн, связанных с
переотражением от морского дна. Удовлетворительное
решение может быть получено при помощи комбиниро-
вания методов прогнозирования.
Пусть зарегистрированные данные pit) описываются
сумой однократных о(/) и кратных 1(f) волн p(f) = o(f) + 1(f).
Выделим из поля однократных отражений волну Oj(r),
которая не была уверенно зарегистрирована, а все ос-
тальные волны обозначим через о2(/): o(f) = О|(/) + о2(/).
Аналогично и для поля кратных волн: 1(f) = l\(f) + l2(f),
где /j(Z) - весь набор кратных волн (полнократные и ча-
стично-кратные), связанных с выбранным кратнообра-
зующим горизонтом, a l2(f) - остальные кратные волны.
Тогда, привлекая априорную оценку модели среды, ре-
шим прямую задачу и получим ‘‘чистый” (т. е. незашум-
лённый) динамический годограф однократного отраже-
ния О](/). Тем самым, фактически, имеем оператор про-
гнозирования волн, связанных с переотражением от это-
го горизонта. Используя методы прямого продолжения,
получим оценку / j(/) поля /](/). После адаптивного вы-
читания /](/) из p(f) приходим к p(t) = p(t) - I\(t) ~
0|(/) + o2(t) + l2(f). Теперь следует решить задачу про-
гнозирования и вычитания поля l2(t), но сделать это не-
посредственно по p(t) нельзя, так как последнее содер-
жит 0|(т). Можно “вырезать” эту волну простым мью-
тингом (обнулением некоторой области в окрестности
ее годографа) либо адаптивно вычесть из волнового поля
p(f) сигнал О](/). Тогда получим p(t) = p(t) - 5}(f) = o2(f)
+ l2(f). Очевидно, что p(f) - поле однократных и связан-
ных только с ними кратных волн. Применив кp(f) про-
цедуру прогнозирования кратных волн без привлечения
модели среды, получим кратные волны l2(f) и результат
подавления кратных волн p(f) = p(f) - l2(f). Таким обра-
зом, мы построили алгоритм прогнозирования и вычи-
тания всех кратных волн, связанных с отражением от
свободной поверхности, даже в случае отсутствия в ис-
ходном поле динамически адекватно зарегистрирован-
ных однократных отражений.
В первых двух частях статьи не была освещена пробле-
ма внутренне-кратных отражений. Помехи такого типа,
как правило, обладают существенно меньшей интенсив-
ностью по сравнению с кратными волнами, обусловлен-
ными переотражением от свободной поверхности. В боль-
шинстве практически важных случаев они не создают
проблем ни для обработки волнового поля, ни для его
интерпретации. С другой стороны, хотя и нечасто, но
приходится сталкиваться с ситуациями, когда внутрен-
не-кратная волна регистрируется в целевой зоне, маски-
руя полезные отражения. Подавление волн такого типа
можно осуществлять с использованием двухшаговых ал-
горитмов. При этом, однако, существенно увеличивают-
ся вычислительные затраты, что далеко не всегда оправ-
дано. Двухшаговые алгоритмы способны прогнозировать
и вычитать только те кратные волны, лучевая схема ко-
торых имеет промежуточную точку выхода луча на сво-
бодную поверхность, а внутренне-кратные волны такой
точки не имеют. Иначе говоря, у нас нет критерия для
различения внутренне-кратных и однократных волн, если
не привлекать априорную информацию о модели среды.
Тогда применяется следующая последовательность про-
цедур для подавления внутренне-кратных волн.
На первом этапе прогнозируются и вычитаются крат-
ные волны, связанные с переотражением от свободной
поверхности.
На втором этапе анализируется полученное волно-
вое поле, в нем идентифицируются внутренне-кратные
отражения и “угадываются” контрастные границы, по-
рождающие эти помехи. Производится построение глу-
бинно-скоростной модели среды в интервале от свобод-
ной поверхности др верхнего горизонта из выбранных
границ.
На третьем этапе производится обратное продол-
жение волнового поля с поверхности на глубинный уро-
вень выбранного горизонта. Тем самым лучевая схема
внутренне-кратных волн “приобретает” промежуточную
точку выхода луча на поверхность, с которой теперь
совпадает новый уровень приведения.
На четвертом этапе применяются двухшаговые схе-
мы для прогнозирования и вычитания волн-помех.
На пятом этапе полученный результат преобразует-
ся к исходному уровню наблюдения методами прямого
продолжения волнового поля.
Можно указать и на более экономичные алгоритмы,
направленные на подавление частично-кратных волн.
Наиболее благоприятной является ситуация, когда го-
дографы кратной волны и сигналов имеют существен-
ные различия, что дает возможность применения про-
цедур разделения волн по кинематическому признаку.
Для простых ситуаций предложены способы прогно-
зирования динамического годографа внутренне-кратной
волны по форме сигнала и годографам однократных от-
ражений от образующих её горизонтов.
ЗАКЛЮЧЕНИЕ. Мы рассмотрели методы прогнозиро-
вания кратных волн для 2D- и ЗО-моделей сейсмичес-
15
кого эксперимента. В большинстве ситуаций, с которы-
ми приходится иметь дело при промышленной обработ-
ке данных сейсморазведки, возможно использование
оператора, не привлекающего априорную информацию
о глубинно-скоростной модели среды, и применение
этого алгоритма предпочтительно.
Прогнозирование кратных отражений при помощи 3D-
алгоритмов по площадным сейсмограммам может быть
оправданным лишь в случаях чрезвычайно сложной гео-
метрии отражающих горизонтов, имеющей простран-
ственный характер. Такой обработке должна предшество-
вать процедура интерполяции сейсмограмм. В силу ука-
занных причин опыт промышленной обработки площад-
ных сейсмограмм при помощи ЗП-процедур прогнозиро-
вания скуден, соответствующие технологии пока не на-
работаны, а методические рекомендации отсутствуют.
Оба способа прогнозирования кратных волн при их
2О-реализации оказываются не настолько времяёмки-
ми, чтобы при обработке реальных данных даже в не-
сложных сейсмогеологических условиях сделать выбор в
пользу алгоритмов 1,5D, где прогнозирование осуще-
ствляется в области разложения Радона.
ПРИМЕЧАНИЯ
ВВЕДЕНИЕ. При выводе формул используется оценка
интегралов методом стационарной фазы (в т. ч. его мно-
гомерным вариантом). Описание метода содержится,
например, в [8], а его обоснование осуществляется в
рамках той же высокочастотной асимптотики, на кото-
рой базируется лучевое приближение, привлекаемое для
анализа волнового поля.
2D-, ЗО-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН ПРИ НАЛИЧИИ
МОДЕЛИ СРЕДЫ. Приближение Кирхгофа для анализа
акустических волновых полей исследовано и обоснова-
но в монографии [10], а также в современной книге,
написанной авторитетными западными исследователя-
ми [12].
Вопросы, связанные с использованием локального
преобразования Радона (т. е. направленного суммирова-
ния в пределах ограниченной пространственной базы)
вместо глобального преобразования (1) и (Г) освеще-
ны в работах [5, 6].
Определение зоны Френеля и зависимость ее разме-
ра от длины волны гармоники, доминирующей в спек-
тре импульса, можно найти, например, в учебнике по
сейсморазведке [11-
Процедура прямого продолжения поля по годографу
отраженной волны с целью прогнозирования кратных
волн тесно связана прямым продолжением поля с неко-
торого глубинного уровня Z\ на новый уровень z2 по го-
дографу дифрагированной волны. С другой стороны, она
не допускает непосредственной аналогии с традицион-
но используемыми алгоритмами и принципами, на ко-
торых они базируются. Оператор прямого продолжения
волнового поля осуществляет взвешенное суммирова-
ние исходных трасс, при этом весовая функция вычис-
ляется через функцию Грина. Функция Грина имеет
смысл волны, возбуждённой точеным источником, рас-
положенным на поверхности интегрирования, и регис-
трируемой на свободной поверхности. Такой оператор
также часто называют интегралом Кирхгофа. Классичес-
кая постановка задачи прямого продолжения поля от-
личается от рассматриваемой нами задачи прямого про-
должения поля с целью прогнозирования цуга кратных
отражений. Не вдаваясь здесь в подробности этих разли-
чий, отметим, что если классический оператор продол-
жения поля пересчитывается из динамического годог-
рафа дифрагированной волны, то нам требуется полу-
чить оператор, который осуществлял бы продолжение
поля по лучу отражённой волны. Действительно, исход-
ный уровень ?[, на котором расположены сейсмопри-
ёмники, совпадает с желаемым уровнем приведения z^.
причем Z] = ?2 = Такой подход требует отдельного
обоснования, которое пока отсутствует в литературе. Это
также является одной из причин, по которой здесь, в
отличие от второй части статьи, мы не проводим парал-
лели между полученным нами в рамках приближения
Кирхгофа оператора и его классической интегральной
формулировкой [9], основанной на вычислении и диф-
ференцировании функции Грина.
Традиционным способом ослабления помех, обус-
ловленных дискретностью данных (аляйсинг-эффект),
является интерполяция, применяемая явно или неявно.
Такой подход не всегда оказывается эффективным, тем
более при большом объёме материала, подлежащего
обработке. Альтернативный способ прогнозирования с
одновременным контролем уровня артефактов преобра-
зования описан в работах [6, 14].
Метод основан на предварительном расчете направ-
ленных сумм, по которым легко реализуется алгоритм
обнаружения аляйсинг-помех. Обнаружение осуществ-
ляется путём сопоставления амплитудных спектров трасс
до и после направленного суммирования.
2D-, ЗО-ПРЕДСКАЗАНИЕ КРАТНЫХ ВОЛН БЕЗ ПРИВЛЕ-
ЧЕНИЯ АПРИОРНОЙ ИНФОРМАЦИИ О СТРУКТУРЕ СРЕДЫ.
Логика рассуждений, основанная на обобщении пря-
мой одномерной задачи на случаи 2D и 3D с последую-
щей инверсией используемой вычислительной схемы,
предложена в работе [11].
Точное определение понятия динамического годог-
рафа можно найти, например, в учебнике по сейсмо-
разведке [1].
Преобразование вида (8), т. е. сумма взаимных свёр-
ток трасс сейсмограмм ОПВ и ОПП, может быть реали-
зовано после перевода этих сейсмограмм в область раз-
ложения Радона. Следует произвести взаимные свёртки
трасс-элементов разложения, соответствующих лучево-
му параметру а для сейсмограммы ОПВ и -а для ОПП. В
этом случае процедура прогнозирования имеет нагляд-
ный физический смысл. Действительно, в соответствии
с законом Снеллиуса, плоская волна, подходящая к
свободной поверхности под некоторым углом, переот-
16
ражается в нижнее полупространство под тем же углом.
Алгоритм прогнозирования в области разложения по
плоским волнам имитирует именно такой процесс фор-
мирования волн-помех. С другой стороны, как показано
в работе |2], если производить перебор по всем а, то
полученное волновое поле будет совпадать с результа-
том, вычисленным по формуле (8). Таким образом, про-
гнозирование в области разложения по плоским волнам
имеет смысл лишь в ситуациях, когда необходимо про-
изводить селекцию волн по их лучевому параметру.
Задача построения устойчивой модификации преоб-
разования (8), в которой осуществляется контроль уровня
артефактов вида аляйсинг-помех. решена в работе [7]. В
отличие от “глобального” суммирования (8) по всем
значениям q, метод основан на предварительном вы-
числении локальных сумм, которые легко поддаются
анализу с целью выявления областей, где артефакты
доминируют над сигналом, и последующей режекции
этих областей.
Результаты обработки данных, полученных при 3D-
морской сейсморазведке, заимствованы нами из рабо-
ты [4].
КАКОЙ АЛГОРИТМ ПРИМЕНЯТЬ НА ПРАКТИКЕ? Форми-
рование "полных” сейсмограмм осуществляется на ос-
новании сейсмического принципа взаимности, в соот-
ветствии с которым трасса, наблюдаемая при положе-
нии источника колебаний а и приемника Ь, совпадает с
трассой, которая была бы получена при положении ис-
точника b и приёмника а.
Вопросы реализации процедуры прогнозирования
кратных волн без привлечения информации о модели
среды для обработки реальных материалов, полученных
при морской ЗО-сейсморазведке, освещены в работе [4].
Если возникает необходимость удаления внутренне-
кратной волны, и при этом ее годограф существенным
образом отличается от годографов сигналов, то такую
процедуру предпочтительно осуществлять методами ки-
нематической селекции волн. К числу последних мы от-
носим и алгоритмы, осуществляющие выделение поме-
хи при одновременном подавлении сигнала средствами
веерной и режекторной фильтрации с последующим
адаптивным вычитанием результата из исходного поля.
Если кинематические различия не столь существенны,
чтобы достичь разделения сигнала и помехи, то наибо-
лее гибким подходом, хотя и требующим значительной
вычислительной нагрузки, является преобразование
внутренне-кратных волн в частично-кратные или пол-
нократные с последующим применением двухшаговых
процедур прогнозирования - вычитания [3|. Упрощён-
ный способ прогнозирования динамического годографа
внутренне-кратной волны по однократным волнам от
глубинных горизонтов, с переотражением от которых
она связана, предложен в работе 113].
* ЛИТЕРАТУРА
I.	Гурвич И. И Боганик Г. FL, 1980, Сейсмическая разведка: М.,
Недра.
2.	Денисов М. С.. 2006, Анализ метода прогнозирования крат-
ных волн без знания модели среды с позиций теории продол-
жения волновых полей: Геофизика, 1, 5 - 20.
3.	Денисов М. С., Кузнецов И. К., 2006, Пример использования
процедур прямого и обратного продолжений волнового поля
при решении задачи подавления кратных волн: Геофизика, 4.
5- 8.
4.	Денисов М. С., Курин Е. А., 2007, Способы прогнозирования
кратных волн поданным плошадных морских наблюдений: Тех-
нологии сейсморазведки, 2. 73 - 78.
5.	Денисов М. С., Фиников Д. Б., 2002, Использование локально-
го направленного суммирования для экстраполяции волнового
поля. Часть 1: Геофизика, 1, 59 - 67.
6.	Денисов М. С., Фиников Д. Б., 2002, Использование локально-
го направленного суммирования для экстраполяции волнового
поля. Часть 2: Геофизика, 2, 25 - 31.
7.	Денисов М. С., Фиников Д. Б., 2005, Способ подавления шумов
дискретизации при суммировании сейсмических трасс (на при-
мере моделирования кратных волн): Геофизика, 1, 12 - 16.
8.	Математическая энциклопедия, 1984: М., Советская энцик-
лопедия.
9.	Петрашень Г. // Нахамкин С. А., 1973. Продолжение волно-
вых полей в задачах сейсморазведки: Л., Наука.
10.	Рытое С. М., 1966, Введение в статистическую радиофизи-
ку: М., Наука.
11.	Berkhout A. J., Ferschuur D. J., 1997, Estimation of multiple
scattering by iterative inversion, Parti and 2: Geophysics, 62, 5,
1586 - 1595, 1596 - 1611.
12.	Bleisiein N., Cohen J. K., Stockwell J. IV., 2000, Mathematics of
multidimensional seismic imaging, migration, and inversion: Springer,
New York.
13.	Denisov M.. Finikov D., Langman S., Oberemchenko D., Spitz S.,
1997, Surface Related and Internal Multiple Attenuation.A Modeling
Approach: Internal. Geophysical Conference and Exhibition,
lstanbul’97.
14.	Lokshtanov D.,Denisov M..Finikov D., 2002, Multiple suppression
and datuming with an antialiased Radon transform for sea-floor
data: 64lh Ann. Internal. Mtg. EAGE.
КОРОТКО ОБ АВТОРАХ
Михаил Сергеевич Денисов - ведущий математик ООО "Геотехсистем”, доктор физ.-мат. наук.
Дмитрий Борисович Фиников - научный сотрудник ООО “Геотехсистем”, кандидат техн. наук.
17