УДК 550.834

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

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

Paper contains consideration of applicability of seismic reflection tomography in regard to reconstruction of absorption-scattering factor distribution above the reflecting boundaries. Paper presents results of tomographic reconstruction of mathematical and physical modeling data as well as real field seismic data. Efficacy of different algorithms is evaluated.

Keywords: tomography, seismic, reflection waves, kimberlite pipe

СЕЙСМИЧЕСКАЯ ТОМОГРАФИЯ С

ИСПОЛЬЗОВАНИЕМ АМПЛИТУД ОТРАЖЕННЫХ ВОЛН

А.Г. Колонин

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

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

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

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

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

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

Также зададимся условием независимости величины коэффициента отражения в фиксированной точке границы от различных углов падения на нее волны (соответствующих различным парам и — удалениям для конкретной общей глубинной точки. Последнее близко к действительности, так как при докритических углах падения и контрастности перепада акустической жесткости на границе до 50% изменения амплитуд отражений не превышают 15%, а при обычных реально существующих перепадах акустической жесткости — и того меньше [7].

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

Если и изменяются около некоторого среднего значения, то, при достаточно большом числе наблюдений с различными и фиксированном (сейсмограмма общего пункта возбуждения), можно записать на основании (2):

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

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

После этого можно перейти от (1) к выражению

получив известную [1,4] интегральную формулу.

Восстановление по логарифмам ослабления может вестись по широкому набору алгоритмов лучевой томографии: обратного проецирования, свертки обратной проекции, итерационного восстановления и т.д. [1,3].

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

1. Метод группирования сигналов [2], или осреднения амплитуд [3], основанный на вычислении средней энергии или амплитуды волн, прошедших по лучам через некоторую область среды, окружающую заданную точку , с компенсацией их геометрического расхождения:

2. Метод обратного проецирования [1], позволяющий оценить непосредственно значения коэффициента поглощения-рассеяния в среде.

3. Метод суммирования разностей, предложенный ранее автором [8], дающий оценку по "остаточному" логарифму ослабления на отрезке за вычетом "нормального" — на остальной части луча :

Как в этих трех способах, так и в соответствующих формулах других алгоритмов [1,4] величина играет роль своеобразного весового коэффициента, характеризующего степень зависимости логарифма ослабления конкретного луча от величины в точке и может быть, соответственно, заменена на .

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

С целью сопоставления результатов различных способов обработки были рассчитаны амплитуды отраженных волн для модели, указанной на рис.1,а. Расстояние между точками задания исходных значений коэффициента поглощения-рассеяния и размеры окружающих их областей среды брались равными 50 м, глубина горизонтально залегающей границы (низ разреза) от поверхности наблюдений (верх разреза) — 1000 м. Шаг источников — 200 м, шаг приемников — 50 м при наименьшем удалении 200 м и наибольшем — 1350 м. То есть, при 24-х канальных наблюдениях кратность перекрытия равнялась 6. При расчетах по (1) дифракция не учитывалась, условия наблюдений считались постоянными ().

По результатам обработки (рис.1,б) видно, что разрез средних амплитуд, полученный по методике простого группирования сигналов с осреднением амплитудных характеристик, содержит в верхней части периодические флуктуации, а "изображения" объектов (области ослабления амплитуд) — размазаны из-за отсутствия учета в (7) существующей квази-экспоненциальнй (а не линейной) зависимости U от . Различаются по разрешенности разрезы , полученные с помощью обратного проецирования и суммирования разностей (рис.1,в-г). Очевидная "размазанность" наиболее хорошо восстановленной (способ суммирования разностей) модели по вертикали объясняется преобладанием субвертикальных направлений имевшихся в распоряжении лучевых траекторий. Также велась обработка по алгоритму одновременного итерационного восстановления [1] с различным числом итераций, но даже при большом их числе результат был сопоставим с разрезом обратного проецирования, при существенных затратах машинного времени из-за большого объема данных.

Результаты обработки полевых данных

В целях апробации указанных алгоритмов на реальном материале, были обработаны данные метода общей глубинной точки, полученные с помощью вибросейсмического комплекса по профилю, проходящему через известную кимберлитовую трубку в Якутии (рис.2,а-б). Кратность перекрытия, изменяющаяся по средним точкам, достигает 16, шаг приемников — 20 м, источников — 80 м, диапазон удалений 80 — 1260 м при центральной системе наблюдений. Предобработка проводилась на ЭВМ ЕС-1055 в системе СЦС-3 и включала в себя препроцессинг, регулировку амплитуд и подготовку данных для ввода. От предобработки требовалось максимальное сохранение динамической информации, необходимой для количественного восстановления поглощающе-рассеивающих свойств.

Компенсация затухания амплитуд со временем была реализована статистическим способом, т.к. большая крутизна кривых затухания амплитуд со временем, вызванная наличием высокочастотных компонентов при зондирующем сигнале 35-140 Гц, а также наличие существенной дифференциации, как по степени кривизны, особенно для малых времен, так и по уровню этих кривых в зависимости от расстояния ПП - ПВ, не позволяют без значительных погрешностей обойтись одной компенсационной кривой для трасс всех удалений, что делается при детерминированном способе. Но прежде, чем решить задачу компенсации затухания амплитуд с временем, была проведена компенсация изменения усиления записи станции по профилю применением подпрограмма GAIN(Х). Затем посредством программы ADGIN было получено семейство кривых затухания амплитуд со временем в зависимости от удаления, после анализа которых проделана компенсация этого затухания путем применения подпрограмм GAIN(). Таким образом, длиннопериодные вариации коэффициента поглощения — рассеяния по профилю были сглажены компенсацией, и в дальнейшем возможно было лишь определение некоторого "остаточного" его значения, способного принимать как отрицательные, так и положительные значения, на неизвестном "сглаженном" фоне.

Коррекция амплитудных статистических поправок за неидентичность условий возбуждения и приема не проводилось, т.к. кимберлитовая трубка — локальная неоднородность, размеры которой сопоставимы с периодом изменения среднепериодных статистических амплитудных поправок (100 + 300 м), а короткопериодные поправки были сглажены в поле за счет размеров баз группирования вибраторов 50 м и приемников 33 м.

Для обработки использовались средние модули сейсмических трасс, вычисленные программой TRAM во временных окнах, содержащих отражение от "опорной", прорываемой трубкой, границы (нулевое время 300 — 400 мс). Эти величины рассматривались непосредственно как амплитуды сейсмических волн, от нее отраженных.

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

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

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

С целью определения граничных условий применимости лучевого подхода проводилось физическое моделирование. Была промоделирована сложная сейсмогеологическая ситуация, когда искомые объекты — зоны дробления и гидротермальной минерализации — располагаются над тектонически нарушенной границей раздела рыхлых терригенных отложений и кристаллического фундамента, причем сверху и снизу экранируются пластовыми телами плотных базальтов или конгломератов. Вмещающей средой служила вода со скоростью продольных волн =1500 м/с, роль подстилающей границы выполнял пенопласт с вырезанными в нем "тектоническими ступенями", пластовых экранирующих тел — оргстекло (=2200 м/с). Объектами являлись гетерогенные зоны из различным образом ориентированных в воде пенопластовых мелкомасштабных элементов, собранных в "грозди". При пересчете модели (рис.3,а) на "реальные условия" в соответствии с [5], можно получить, при сохранении скоростных параметров: преобладающая длина волны =100 м, полоса частот импульсного сигнала 40-120 Гц, глубина до "реперной" границы 1000-1100 м, амплитуда ее нарушений 100 м, диаметры гетерогенных зон — порядка 200 м. Геометрия наблюдений соответствовала описанной для математической модели. Если рассчитать диаметр первой зоны Френеля:

где и — расстояния от объекта до источника и приемника (по лучу), то для минимального удаления получим /=0,625, максимального — 0,5. Как показывает большой объем полученных экспериментальных данных [5,9], дифракция на самом объекте в этих случаях будет достаточно сильна — амплитуда интерференционной волны в зоне "геометрической тени" может составлять до половины от амплитуды "нормальной" волны. Сильное ослабление отраженных волн (рис.3,б) связано с "разрушительной" интерференцией на "ступенях" подстилающей границы и экранирующих слоев.

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

Результаты показывают, что по средним амплитудам выделяются, как области пониженных , все тектонические "ступени". Сами гетерогенные зоны из-за дифракции выделяются неясно, "сливаются" с расположенными под ними "нарушениями". Сглаженные картины получаются по результатам обработки с помощью обратного проецирования и итерационного восстановления. Имеющие место на всех разрезах противоположные изменения и в целом по профилю связаны с неучтенным ослаблением амплитуд за счет отражения на экранирующих пластах. Наиболее разрешенное восстановление коэффициента поглощения-рассеяния получено на разрезах суммирования разностей, причем заметно, что предварительный учет вариаций "коэффициента отражения" (рис.3,в) позволяет подавить влияние дифракционного ослабления волн на "ступенях" границы. Несмотря на это, явно неполный учет разрушения волнового поля на "ступенях" границ и экранирующих пластов не позволяет однозначно локализовать расположенные справа и слева объекты — они проявляются на фоне сильных побочных аномалий, связанных с "тектоникой". Даже центральная гетерогенная зона искажена странным образом и проявляется совместно с "разломом", под ней расположенным. Хотя ясно видна область повышенного поглощения-рассеяния (рис.3,в) в средней части разреза, в центре ее расположен локальный минимум, связанный, по-видимому, с отсутствием учета дифракции. То есть — по результатам обработки данных метода отраженных волн с 6-ти кратным перекрытием способом суммирования разностей, с учетом переменного коэффициента отражения, возможна лишь локализация структур типа "разлом — гетерогенная зона над разломом", в случае отсутствия экранирующих рассеивающих волны пластов-помех.

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

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

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

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

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

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

1. Бельфер И.К., Непомнящих И.А. Сейсмическая томография. — М., 1988. — 70 с: Ил. — (Разведочная геофизика: Обзор/ВНИИ экон. минер. сырья и геол. — развед. работ (ВИЭМС). — Библиогр. с. 66 — 70 (78 назв.)

2. Гольцман Ф.М., Калинин Д.Ф. Группирование сигналов в задачах обнаружения нарушений среды //СО АН СССР. Геология и геофизика. — 1986. — N° 5 — c. 85 — 94.

3. Колонин А.Г. Возможности использования результатов сейсмического просвечивания для обнаружения локальных неоднородностей //СО АН СССР. Геология и геофизика. — 1989.—N 3.

4. Кьяртанссон Э. Анализ вариации амплитуд и времен пробега в зависимости от удаления и положения срединной точки // Численные методы в сейсмических исследованиях. — Новосибирск: Наука, 1983, — с. 221 — 233.

5. Гик Л.Д. Сейсмическое моделирование сложнопостроенных структур. — Новосибирск: Наука, 1983 — 118 с.

6. Файзуллин И.С., Шапиро С.А. Рассеяние и зависимость затухания сейсмических волн от длины базы наблюдений.1. Элементы теории. // Изв. АН СССР. Физика Земли. 1988. N° 2. — с. 20 — 30.

7. Tooley R.D., Spenser T.W. and Sagosi H.F. Reflection and transmission of plane compressional waves: Geophysics, 30, 552-70, 1965.

8. Колонин А.Г. Автоматическая обработка данных просвечивания // Передовой науч.-произв. опыт / ВИЭМС. - М., 1989. - Вып.7. - С. 3-6.

9. Колонин А.Г. Дифракция на трехмерных изометричных неоднородностях. // Изв. АН СССР. Физика Земли.1991. N 7. С. 112-119.

Рис.1. Исходный разрез коэффициента поглощения-рассеяния и результаты обработки данных математического моделирования. а — исходный разрез; б — разрез средних амплитуд; в — разрез обратного проецирования; г — разрез суммирования разностей.

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

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