Содержание к диссертации
Введение
Глава 1. Численное моделирование сейсмических полей методом конечно-разностных схем 19
1. Математическая постановка задачи лля уравнении теории упругости 20
2. Конечно-разностные схемы 22
3. Порядок аппроксимации, сходимость и устойчивость разностной схемы 26
4. Сеточная дисперсия 29
5. Задание граничных условий, диссипативные граничные условия 30
Глава 2. Математическая постановка задачи миграции сейсмических полей 41
1. Формулировка задачи миграции. 42
2. Оптимизационная постановка задачи продолжения сейсмических полей
3. Особенности и ограничения процедуры экстраполяции волнового поля 46
Глава 3. Векторная инверсия сейсмических полей 53
1. Постановка задачи инверсии на основе лучевой селекции 53
2. Определение кинематических параметров среды и нормали к границе
3. Определение отношения плотностей и вычисление векторных коэффициентов отражений
4. Регуляризация процедуры инверсии 64
5. Метод селекции волновых полей ВСП 70
6. Построение изображения среды в векторных коэффициентах отражения J2
Глава 4. Алгоритмы решения, результаты численного моделирования и обработки реальных данных 76
1. Алгоритм вычисления продолженного поля 76
2. Алгоритм векторной инверсии, исследование устойчивости 77
3. Модельный эксперимент, описание модели, результаты, оценка эффективности метода
4. Применение алгоритмов миграции и инверсии для обработки реальных данных «4
Заключение 102
Итоги 102
Перспективы 103
Благодарности 103
Литература
- Конечно-разностные схемы
- Оптимизационная постановка задачи продолжения сейсмических полей
- Определение кинематических параметров среды и нормали к границе
- Алгоритм векторной инверсии, исследование устойчивости
Введение к работе
Метод вертикального сейсмического профилировании
Вертикальное сейсмическое профилирование [13] — метод сейсмической разведки, при котором регистрация колебаний, возбужденных источником на поверхности, осуществляется во внутренних точках среды на вертикальном профиле (а точнее, вдоль скважины). Существуют также и модификации метода, когда источник возбуждения помещается в скважину, а приемники -на земную поверхность.
ВСП выросло из сейсмического каротажа, призванного обеспечить глубинное преобразование наблюдений сейсморазведки на поверхности, получаемых исходно в масштабе времен. Использование последующих вступлений создало возможности для полноценного изучения волновых полей в скважине и дальнейших попыток исследования околоскважинного пространства с использованием выносных пунктов возбуждения, профилей возбуждения (метод обращенного годографа - МОГ) и двумерных систем возбуждения (3D ВСП).
К принципиальным особенностям метода ВСП относятся следующие [14]: изучается не только волновое поле, наблюдаемое па дневной поверхности, но и сам процесс его формирования. Одновременно выделяются, прослеживаются и изучаются волны разных типов и классов (прямые, отраженные, кратные, преломленные, продольные, поперечные и обменные), возбуждаемые в источнике и образовавшиеся на неоднородностях среды;
в отличие от большинства традиционных геофизических (электрических, акустических и др.) измерений в скважинах, изучающих разрез только в ближайшей окрестности ствола скважины, ВСП позволяет исследовать околоскважинное и межскважинное пространство на значительных расстояниях от скважины и в большом диапазоне геологических условий, не только в интервале глубин, вскрытых скважиной, но и глубже забоя;
- ВСП, являясь одновременно сейсмическим и скважинным методом, находится на стыке наземной сейсморазведки и ГИС - геофизических исследований скважин. ВСП позволяет, с одной стороны, достичь разрешающей способности, близкой к ГИС, при изучении разреза вдоль ствола скважины и, с другой — распространить эти данные на окрестности скважины.
Среди решаемых при помощи ВСП задач основными на сегодняшний день являются:
- изучение упруго-плотностных характеристик среды на сейсмических частотах;
- динамическая привязка отражений, регистрируемых на поверхности, к литологическому разрезу;
- детальное изучение околоскважинного пространства с использованием более широкого спектра частот, чем при сейсморазведке на поверхности;
- прогнозирование геологического разреза ниже забоя скважины. Очевидными преимуществами ВСП являются возможности деконволюции по форме прямой волны, а также достоверной оценки интервальной скоростной модели среды по временам первых вступлений с нескольких пунктов возбуждения [31, 34]. Кроме того, система наблюдений ВСП обеспечивает полноценную трехкомпонентную регистрацию колебаний среды, что делает возможным качественное разделение и изучение различных типов волн и условий их формирования, а использование полного вектора смещений при построении сейсмических изображений делает их динамически наиболее адекватными реальной среде. Наконец, только наличие информации о полном векторе колебаний является основой для построения алгоритмов инверсии поля упругих смещений с целью восстановления таких характеристик среды, как импедансы, соответствующие продольным и поперечным волнам, и/или векторные (вдоль нормали к границе раздела жесткостей) коэффициенты отражения.
Математическая модель упругой среды
При обработке и интерпретации данных сейсмической разведки одним из основных объектов изучения является модель физических параметров среды. Распределение физических характеристик горных пород, восстанавливаемое по данным сейсмических наблюдений, используется при построении геологической модели месторождения.
Математически модель упругой среды описывается тензором 4-го ранга, состоящим из 36 независимых материальных параметров. Этот тензор входит в закон Гука, определяющий связь между напряжением а и деформациями є материала:
& = Сє, или &J = CijUeu. (1)
Закон Гука, а также второй закон Ньютона ри„ = diva-, (2)
где р = p(x,y,z) - плотность и diver = V ® т, задают уравнение относительно вектора смещений u = u(x, ,z) для процесса распространения упругих колебаний. Коэффициентами этого уравнения являются компоненты тензора С материальных параметров среды.
В изотропном случае тензор С зависит только от двух параметров А и //, которые называются параметрами Ламэ и связаны со скоростями распространения продольных и поперечных волн а и Ъ соотношениями: Л = (а2-2Ь2)р, (3) М = Ъ Р-При этом уравнения (1) и (2) принимают вид:
a = Aldivu + 2/.іє, (4)
/лі„ = grad(2di vu) + 2div(//), (5)
где I — единичный тензор.
Коэффициенты уравнения (5) в предположении о слабом изменении параметров среды могут быть выражены через скорости продольных и поперечных волн, а само уравнение в этом случае расщепляется на независимые волновые уравнения относительно дилатации и ротора смещения, описывающие распространение соответственно продольных и поперечных волн [24].
На практике вместо параметров Ламэ для описания модели упругой среды часто используют другие величины: скорости продольных и поперечных волн и плотности, коэффициент Пуассона и модуль Юнга, импедансы (произведение скорости и плотности) или коэффициенты отражения для продольных и поперечных волн. Это обусловлено, прежде всего, требованиями конкретной задачи и соображениями удобства. Так, в процессе обработки данных сейсморазведки (в том числе и при миграционном преобразовании), а также при математическом моделировании в качестве опорной используется главным образом скоростная модель среды, в то время как для задач геологической интерпретации физически более значимыми являются импедансы или коэффициенты отражения.
Наибольший разведочный интерес представляют такие параметры реальных сред, как пористость, проницаемость, трещиноватость, флюидосодержание и т.п. Эти параметры довольно многочисленны и не могут быть ни напрямую учтены при описании математической модели среды, ни, соответственно, восстановлены по результатам сейсмических наблюдений. Однако об их изменении можно судить по характеру распределения скоростей, плотности, коэффициентов анизотропии и других величин, характеризующих свойства упругих сред.
Используемые на практике модели различаются по степени сложности распределения параметров среды. Выбор того или иного описания модели определяется главным образом условиями поставленной задачи, а также имеющимися в распоряжении вычислительными мощностями. Модели реальных сред можно разделить на одномерные, двумерные и трехмерные.
При рассмотрении одномерной модели предполагается, что параметры разреза слабо изменяются по латерали и этими вариациями можно пренебречь. В случае существенной латеральной неоднородности необходимо переходить к двумерным или трехмерным моделям среды.
Реальные среды в нефтегазовой геологии и геофизике часто описываются слоистыми моделями. При этом предполагается, что внутри слоев параметры разреза меняются гладко, а на границах раздела терпят разрыв. Границы слоев задаются в виде гладких или кусочно-гладких кривых (для трехмерных моделей — поверхностей).
Моделирование сейсмических полей, лучевое приближение, полное динамическое решение
Математическое моделирование получило широкое применение при решении различных задач геофизики и геологии. Моделирование сейсмических данных помогает решить многие вопросы методики проведения полевых работ, обработки волновых полей и интерпретации ее результатов. Под математическим моделированием в сейсморазведке подразумевают процедуры построения сейсмогеологической модели среды, математическое описание ее параметров, расчет синтетических волновых полей для заданных конфигураций систем наблюдения. Интерпретация реальных волновых полей на основе математического моделирования предполагает последовательное уточнение априорной модели путем сравнения теоретических сейсмограмм с реальными данными и корректировки параметров модели.
Методы расчета волновых полей можно разбить на две группы. К первой группе относятся асимптотические методы, основанные на лучевом приближении, их обзор см. в работе [2]. Ко второй группе относятся методы расчета полного волнового поля, основанные либо на точном решении уравнения распространения упругих или акустических колебаний с теми или иными ограничениями на рассматриваемую модель среды [38, 50, 68], либо на численном решении этого уравнения с помощью конечно-разностных схем [42, 46, 50].
Лучевой метод [1, 27, 28] позволяет рассчитывать волновое поле путем его представления в виде бесконечного ряда, каждый член которого описывает вклад в суммарное поле так называемой элементарной (обобщенной) волны и называется лучевым приближением. Каждая такая волна соответствует какому-либо допускаемому геометрической сейсмикой лучу, выходящему из источника и приходящему к приемнику, имеет определенную природу и характеризуется амплитудой, формой и временем прихода в точки наблюдения. Сейсмограмма образуется путем суммирования вкладов всех элементарных волн, причем в каждый момент времени число ненулевых вкладов конечно.
Описание волнового поля в рамках лучевого метода справедливо, когда поле времен, определяющее фронт волны, является аналитической функцией. Там, где это условие нарушается, происходит либо разрыв волновых фронтов, либо разрыв на фронте амплитуды волны, определяемой лучевым методом. Вблизи таких особенностей значительная часть энергии будет переноситься в направлениях, отличных от направления распространения волны и в которых законы лучевого распространения теряют силу. Явления распространения упругой энергии, не соответствующие лучевой схеме, относят к явлениям дифракции. Таким образом, с помощью лучевого метода невозможно получить полное решение прямой задачи при наличии границ, являющихся источниками дифрагированных волн: границ с изломами, выпуклых границ, при переходе через которые скорость терпит разрыв.
Точное решение прямой динамической задачи можно получить для вертикально неоднородных сред методом полного разделения переменных. При этом решение уравнения упругости, описывающего колебания среды, находится в виде двойного интеграла, где один интеграл берется по частоте, а второй — по параметру луча. Несмотря на то, что формальное решение и записывается в явном виде, для практических расчетов волновых полей необходимо вычислять интегралы с бесконечными пределами, что является самостоятельной и нетривиальной задачей.
В основе конечно-разностного подхода лежит аппроксимация производных конечными разностями при решении уравнений распространения колебаний. Конечно-разностные методы условно можно разделить на три группы [61]: 1) методы экстраполяции; 2) комбинированные методы, типа «предиктор-корректор»; 3) подстановочные методы, типа метода Рунге-Кутта. В силу исторических причин, а также благодаря своей простоте и универсальности наиболее широкое распространение в сейсмическом моделировании получили методы первой группы.
Являясь наиболее естественным подходом к математическому моделированию волновых полей, использование конечно-разностных схем потенциально свободно от ограничений на сложность используемых моделей. При этом в результате получаются все типы волн, которые могут образоваться в изучаемой среде при заданной модели распространения.
Заметим, что при использовании конечно-разностных методов следует помнить о таких проблемах, как точность и устойчивость используемой разностной схемы, выполнение законов сохранения, граничные условия, эффекты рассеяния на пространственной сетке и др.
Метод миграции волновых полей
Миграция сейсмических волновых полей - один из эффективных методов, применяемых в настоящее время при обработке данных сейсморазведки. Некоторые из теоретических и методических основ метода миграции были заложены еще в 50-60-х годах прошлого века до превращения сейсморазведки в промышленную технологию [57].
Идея миграционного преобразования записей сейсмических волновых полей впервые была четко сформулирована Ю.В. Тимошиным в 1960 г. Им было введено представление о сейсмическом поле как о суперпозиции дифрагированных волн и о среде как о совокупности точек дифракции, а также даны основы теории дифракционных преобразований [35]. Г.И. Петрашенем и С.А. Нахамкиным введено представление об обращенном продолжении сейсмического поля как об инструменте восстановления поля в среде в момент его возникновения в точках дифракции, обозначена роль фундаментальных законов теории распространения волн как теоретической основы продолжения полей [29].
За рубежом всплеск в развитии теории миграционных преобразований произошел в 70-е годы и был вызван появлением серии основополагающих работ Клаербоута и его коллег из Стэнфордского университета [51, 52, 53]. Спустя некоторое время были представлены методы решения задачи миграции на основе интегральной формулы Кирхгоффа, аналогичные дифракционным преобразованиям [62, 66, 69], в пространстве частота-волновое число [12, 55, 71], а также работы по миграции в обратном времени [44,48,49,63,74].
Основной задачей миграции является перенесение и фокусирование сейсмической энергии, рассеянной на неоднородностях среды и зарегистрированной на земной поверхности или в скважине, в реальные местоположения рассейвателей, к которым она относится. В результате такой процедуры формируется структурное изображение среды, которое в общем случае динамически должно отражать распределение упруго-плотностных параметров (например, коэффициентов отражения). Метод миграции логически состоит из нескольких этапов. Вначале производится продолжение зарегистрированных приемниками сейсмических полей в геологическую среду, описываемую опорной моделью. Затем к продолженному полю применяется условие отображения, сформулированное Клаербоутом [52]: отражатели находятся в тех точках среды, где пересекаются поля времен падающей и отраженной волн. И, наконец, проводится инверсия полученного изображения среды.
Процедуры миграции волнового поля могут выполняться на основе лучевого подхода, в частотной области либо продолжением сейсмического поля с помощью метода конечных разностей.
Применение лучевого подхода и расчет продолженного поля в частотной области, как и при решении задачи моделирования, накладывают определенные ограничения на сложность опорной модели. Кроме того, при использовании лучевого метода возникают трудности, связанные с необходимостью учета множественных траекторий. С другой стороны, в случае применения конечно-разностного подхода приходится решать типичные для разностных схем проблемы выбора граничных условий, устойчивости, рассеяния на сетке и др.
Для последующей инверсии принципиален вопрос о сохранении динамических характеристик волнового поля в процессе миграции. В этом отношении более предпочтительными являются конечно-разностные методы, так как вся динамика содержится в уравнениях для упругих колебаний, на основе которых и строятся разностные схемы, в то время как корректный учет всех факторов распространения (преломление, рассеяние, обмен, геометрическое расхождение и т.п.) в рамках лучевого приближения требует усложнения вычислений.
По способу экстраполяции волнового поля во внутренние точки среды методы миграции делятся на интегральные, алгоритмы продолжения с поверхности вдоль вертикальной оси, основанные на решении параксиальных аппроксимаций волновых уравнений, и методы продолжения путем решения нестационарных краевых задач распространения колебаний с обратным временем. Наиболее естественным и точным [56] методом экстраполяции сейсмического поля является последний из перечисленных, однако на практике (особенно при миграции данных наземной сейсморазведки) часто используются более простые и экономичные алгоритмы продолжения в нижнее полупространство вдоль оси глубин или интегральные методы. В то же время при обработке данных ВСП в силу специфики геометрии и ограниченности апертуры наблюдений экстраполяция в пространстве вдоль заданной оси теряет смысл, а применение процедур миграции, основанных на продолжении полей с использованием формул Кирхгоффа и Грина, является теоретически некорректным, а практически не всегда адекватным.
Результатом алгоритмов экстраполяции является волновое поле, первоначально заданное в точках приема, а затем продолженное в геологическую среду. Мигрированное волновое поле качественно отражает изучаемый геологический разрез в том смысле, что наибольшее изменение поля соответствует искомым рассеивающим неоднородностям среды, таким как границы слоев, нарушения границ и т. д.
Заключительным этапом процедуры миграции является построение изображения среды, количественно описывающее ее характеристики, то есть, по сути, решение задачи инверсии. Практически все известные подходы в этой области основаны на эвристическом положении, высказанном в первых работах по миграции [52, 19, 20]. Суть его состоит в том, что коэффициент отражения равен отношению амплитуд мигрированного и первично падающего поля. Ясно, что такой подход не применим при инверсии векторных полей упругих смещений для определения векторного коэффициента отражения, то есть коэффициента, соответствующего направлению вдоль нормали к отражающей границе. Кроме того, для получения такого коэффициента является необходимым использование информации обо всех типах волн (т.е. о полном векторе смещений), образующихся в каждой точке среды. Это требование редко выполняется на практике, так как в силу устоявшейся традиции большинство методов миграции базируются на скалярном волновом уравнении. При этом существующие алгоритмы упругой инверсии [65, 73], имеют ряд недостатков, среди которых недостаточная устойчивость, обусловленная объединением процессов миграции и инверсии в одну оптимизационную процедуру, и необходимость серьезных вычислительных затрат для достижения приемлемой точности.
Комплекс вышеизложенных фактов с учетом специфики и преимуществ метода вертикального сейсмического профилирования определяет актуальность разработки алгоритмов инверсии с целью восстановления векторных коэффициентов отражения и для формирования идеологии миграции и построения динамически представительных сейсмических изображений в истинных коэффициентах отражения по данным ВСП.
Цели, задачи, новизна, защищаемые положения
Цель работы
Разработка алгоритмов и программ моделирования и миграции волновых полей ВСП на основе конечно-разностного подхода, инверсии волновых полей для решения задачи восстановления векторного коэффициента отражения для продольных и поперечных волн и построения изображения среды в виде вертикальных профилей с количественной оценкой коэффициентов отражения границ раздела.
Основные задачи исследования
1. Разработка программы расчета сейсмических полей методом конечных разностей
2. Исследование свойств различных типов диссипативных граничных условий в задаче конечно-разностного моделирования сейсмических волновых полей
3. Разработка алгоритма и программы продолжения волнового поля ВСП в обратном времени
4. Разработка алгоритма и программы векторной инверсии волнового поля ВСП
5. Оценка эффективности разработанных алгоритмов и программ на модельных данных
6. Опробование разработанных алгоритмов миграции для обработки реальных данных ВСП
Научная новизна
1. Получено решение задачи продолжения сейсмического поля в оптимизационной постановке
2. Получено решение задачи восстановления векторного коэффициента отражения по данным ВСП (векторная лучевая инверсия)
3. На основе комбинации процедур миграции, селекции и инверсии волнового поля ВСП разработан метод построения сейсмического изображения среды в терминах нормальных коэффициентов отражений
Практическая ценность
1. Реализована программа конечно-разностного моделирования волновых полей
2. Составлен алгоритм и написана программа продолжения волновых полей ВСП в среду на основе конечно-разностного подхода
3. Разработан алгоритм и программа векторной лучевой инверсии данных ВСП
4. Разработана методика построения изображения среды в терминах векторного коэффициента отражения
Защищаемые положения
1. В оптимизационной постановке задачи миграции данных ВСП экстраполированное сейсмическое поле вычисляется методом градиентного спуска, направление которого определяется в результате решения нестационарной краевой задачи распространения упругих колебаний с обратным временем
2. Метод векторной лучевой инверсии на основе информации о векторах смещений всех типов волн позволяет эффективно определять коэффициенты отражения вдоль нормали к границе для продольных и поперечных волн, а также параметры залегания пластов
3. Сейсмическое изображение среды в истинных (по нормали к границе) коэффициентах отражения продольных и поперечных волн может быть построено как результат селекции и последующей векторной инверсии продолженного волнового поля на множестве (вертикальных) профилей, пересекающих исследуемую область
Личный вклад автора
1. Разработка алгоритма экстраполяции волнового поля ВСП на основе решения задачи распространения упругих колебаний методом конечных разностей
2. Разработка алгоритма динамической векторной инверсии данных ВСП
3. Исследование устойчивости и регуляризация процедуры инверсии, выбор оптимальных значений параметров регуляризации
4. Исследование эффективности алгоритмов экстраполяции и инверсии волнового поля ВСП на модельных данных
5. Опробование предложенного метода миграции на реальных материалах ВСП
Конечно-разностные схемы
Рассмотрим функцию целочисленного аргумента y(i), / = 0, ±1, ±2,..., и образуем в точке / разности: Ay, =y(i + l) y(i) - правую, Vy, = y{i) — y(i — \) - левую. Эти выражения можно считать формальными аналогами первой производной функции y(i). Обозначим yi = y(i) и введем разность второго порядка А2Уі = Д(Ду,) = A(yi+] -yt) = (yi+2 - .у;.+1) - (ум -у,) = yi+2 - 2yi+l + у,.
Заметим, что Ду._, = Vy,, то есть применение оператора левой разности эквивалентно применению оператора правой разности в точке, смещенной на единицу влево. Следовательно AVy, = А2 .., = ум - 2yt + у._,.
Аналогично второй разности определим разности более высоких порядков: А ,.=А(А .). Каждое последовательное применение оператора А захватывает еще одну точку вправо, так что после т применении оператор A yf содержит значения функции у., yj+], ..., у1+т в точках/, / + 1, ..., i + m.
Можно записать уравнение, содержащее конечные разности различных порядков: a0Amyi + alA "-lyi +... + ат_хАУі + атУі = f,, где а0 ах, ..., ат - некоторые коэффициенты. Подставляя выражения для разностей Ayt, ..., Amyj, получим ад+« + а\У т-\ + - + п-\Ум + атУ, = f, Если коэффициенты а0 Ф 0, ат Ф 0, то это уравнение называется разностным уравнением m-го порядка относительно искомой функции целочисленного аргумента уг Это разностное уравнение - формальный аналог дифференциального уравнения т -го порядка:
Подобно тому, как коэффициенты дифференциального уравнения могут являться функциями от х, коэффициенты разностного уравнения тоже могут быть переменными, то есть ак = cck{i), к = 0, 2,..., т.
Перейдем от дифференциальной задачи к разностной постановке. Пусть требуется приближенно вычислить решение и дифференциальной краевой задачи Lu = f, поставленной в некоторой области D с границей 3D.
Для этого необходимо выполнить следующие действия: 1) выбрать дискретное множество точек (сетку) DhczD + dD (h — шаг сетки); 2) ввести линейное нормированное пространство функций U/;, определенных на сетке Dh, и установить соответствие между решением и и функцией йн є ил, которая будет считаться искомой таблицей значений решения и; 3) определить оператор разностной задачи Lh, заменив в L производные на соответствующие конечно-разностные выражения; 4) ввести также сеточные аналоги для всех дополнительных условий: источника /, а также начальных и краевых условий.
Совокупность разностных уравнений, аппроксимирующих основное дифференциальное уравнение и дополнительные (краевые и начальные) условия, называют разностной схемой.
Уравнения (1.1) и (1.2) приводят, соответственно, к разностным задачам второго и первого порядка по времени. Конечно-разностное выражение для производной по времени может связывать различное число соседних точек временной сетки, которые принято называть слоями. По этой терминологии, как правило, различают двухслойные и трехслойные разностные схемы.
На практике расчет волнового поля осуществляется путем последовательного обновления поля смещений на пространственной сетке при переходе с одного временного слоя на другой.
Широкое распространение вычислительных машин и популярность конечно-разностных методов обозначили необходимость разработки универсальных разностных схем и алгоритмов, обеспечивающих решение целых классов задач, определяемых типом дифференциального уравнения, начальных и краевых условий, а также функциональным пространством, в котором заданы коэффициенты уравнения. Это требование приводит к понятию однородных разностных схем, вид которых не зависит ни от выбора конкретной задачи из класса, ни от выбора дискретной сетки. Во всех узлах сетки разностные уравнения однородной схемы имеют один и тот же вид, а их коэффициенты определяются как функционалы коэффициентов соответствующего дифференциального уравнения.
Особый интерес представляет построение однородных схем «сквозного» счета [32], пригодных для решения уравнений с разрывными коэффициентами без явного выделения точек или линий разрыва параметров среды (ситуация, типичная при моделировании сейсмических наблюдений — реальные среды чаще всего описываются слоистой моделью, в которой материальные параметры терпят разрыв первого рода на границах раздела слоев). Это значит, что схема в окрестности разрывов не меняется и вычисления во всех узлах производятся по тем же формулам, что и в случае непрерывных коэффициентов1. При конструировании таких разностных схем важно обеспечить выполнение в точках разрыва условий сопряжения, вытекающих из физической природы моделируемых процессов, а именно из законов сохранения.
Оптимизационная постановка задачи продолжения сейсмических полей
Рассматриваемый метод, как и любые процедуры экстраполяции волнового поля, использует априорную информацию о скоростной модели изучаемой среды. Очевидно, что если опорная модель в точности соответствует реальному распределению физических параметров, то и процедура продолжения поля (и построения изображения) будет давать точный результат как в смысле положения отражающих горизонтов, так и в смысле их динамических характеристик. Но так же очевидно и то, что построить абсолютно адекватную модель среды невозможно. Поэтому существенным является вопрос о том, насколько тот или иной метод миграции (экстраполяции волнового поля) чувствителен к погрешностям в задании априорной скоростной модели среды. Подобное исследование необходимо, но выходит за рамки настоящей работы, и мы ограничимся общими замечаниями.
Главным требованием к априорной модели среды является, прежде всего, то, чтобы она в целом отвечала кинематическим свойствам реальной среды. Выполнение этого условия можно обеспечить, так как наблюдения методом ВСП из нескольких пунктов возбуждения, а также привлечение дополнительной информации за счет использования совмещенных наземно-скважинных систем наблюдения [8] позволяют эффективно решать обратную кинематическую задачу и получать достаточно достоверную интервальную скоростную модель среды. Вместе с тем известно [37], что методы миграции, использующие при экстраполяции волнового поля интервальную скоростную модель (методы глубинной миграции, к которым относится и рассматриваемый в работе подход), более чувствительны к ошибкам задания скоростей, чем алгоритмы временной миграции, задействующие средние скорости. Средние скорости определяются по принципу наилучшей фокусировки рассеянной сейсмической энергии [39] и, вообще говоря, не описывают реального распределения параметров среды, в отличие от модели интервальных скоростей. Именно поэтому, хотя и являясь менее устойчивыми, в отсутствие резких неизвестных латеральных неоднородностей среды и серьезных ошибок в определении кинематических параметров методы глубинной миграции предпочтительнее, так как обеспечивают построение более корректного сейсмического изображения в смысле положения, структуры и динамических характеристик отображаемых горизонтов [59].
Представленный метод миграции включает экстраполяцию волнового поля на основе решения задачи распространения упругих колебаний в обратном времени и базируется на предположении о том, что используемая опорная скоростная модель правильно отражает основные кинематические характеристики разреза. Как уже отмечалось, рассматриваемый подход относится к группе методов глубинной миграции, и, следовательно, к нему относятся высказанные соображения относительно требований к точности задания опорной модели среды.
Важной задачей при построении продолженного поля является понимание того, куда именно возможно экстраполировать наблюдения с профиля регистрации. Формально продолжение поля может быть получено во всем исследуемом пространстве Q. Однако реально размеры области, в которой волновое поле восстанавливается корректно, ограничены и в первую очередь определяются апертурой наблюдений.
Это ограничение не самое существенное при миграции данных наземной сейсморазведки, когда апертурой наблюдений фактически служит земная поверхность, и продолжение волнового поля осуществляется в нижнее полупространство. В случае данных ВСП ситуация изменяется в силу того, что апертура наблюдений ВСП значительно меньше, чем у сейсморазведки на поверхности.
Эффекты, связанные с апертурными ограничениями, легко проиллюстрировать на простейшем примере. На рис. 13 изображена схема наблюдений ВСП. Среда предполагается однородной, освещаемая область выделена цветом.
Определение кинематических параметров среды и нормали к границе
Функционал (3.5) достигает минимума в точке, где его частные производные по параметрам обращаются в нуль: дФ - = 0, / = 1,2,3. (3.7) ш, Условие (3.7) задает систему нелинейных уравнений относительно вектора параметров d. Для ее решения воспользуемся методом Ньютона [11], в соответствии с которым вектор неизвестных системы уравнений F(x) = 0 находится как результат итерационного процесса x(«i) = хм _ [j(xW)] F(xU)), где J(x) = ( (x)) = L(x) fdF, Л ydxj матрица Якоби системы.
Во избежание необходимости обращения матрицы Якоби на каждой итерации метода алгоритм вычисления следующего приближения х(Л+) удобно записать в виде системы линейных уравнений: J(x(i))x(i+1) = J(x(s))x(i) -F(x(J)), так как решение системы линейных уравнений - процедура, обеспечить устойчивость которой легче, чем устойчивость вычисления обратной матрицы в силу доступности большего числа разнообразных алгоритмов ее выполнения. \дФ дФ дФ] Таким образом, для х = d = {dv d2, d3) и F = , , имеем: \ddx dd2 ddij J(d(5,)d(5+I) = J(d(,,)d ,) -F(d(i)), где (3.8) J(d)= (3.9) в предположении, что функционал (3.5), по крайней мере, дважды дифференцируем по параметрам J,, d2, d3 (справедливость этого предположения будет показана ниже).
Процесс минимизации продолжается до тех пор, пока относительная погрешность є между двумя последовательными приближениями не достигнет заданной достаточно малой величины, то есть пока не будет выполнено условие d( +1)-d ( ) lis) (3.10)
Метод Ньютона отличается высокой, квадратичной скоростью сходимости, то есть для последовательности приближений {d(5) выполнено условие ( +) с (S) где d — точное решение, С — положительная константа. Для рассматриваемой задачи при удачном начальном приближении приемлемая точность (относительная погрешность, не превышающая 0.1%) достигается за две-три итерации.
В качестве начального приближения — вектора d — могут быть заданы априорные отношения скоростей (если они известны) или значения, типичные для однородной среды: например, d = {і, 2, 2}.
Представление тензора градиента смещений в инвариантной форме
Сначала рассмотрим вектор смещений падающей продольной волны UDP = A0/(/-ar,(p0,r)), где А0 = {4,Л ,А30}, р0 = {р1р20,р30} и r = {x,y,z},u вычислим его градиент, исключая координаты точек пространства. Далее, опуская общий фазовый множитель / () запишем тензор Л Л Л градиента полного вектора смещений U = U, - U2. Для этого вспомним, что Ul = UDP + UUrr + UUPS u2 = uDPP + uDPS и, следовательно, и = -ог1(А0р0 + А,р1) + ог2А2р2-ДВ1 q, + j32B2q2. (3.11)
Окончательно, выразим тензор U через искомые параметры dlt d2, d3, разделив (3.11) на -or, (такая операция правомерна, так как не нарушает справедливости условия (3.3)): U = (i0(A0(8)p0 + A1pl)- 1A2(8 p2+J2B1(8)q1-J3B2q2. (3.12)
Как видно из (3.12), тензор градиента смещений линейно зависит от кинематических параметров. Это значит, что первое слагаемое минимизируемого функционала (3.5), содержащее квадраты миноров второго порядка матрицы U, относительно этих параметров представляет собой многочлен четвертой степени. Второе и третье слагаемое функционала является многочленом второй степени по переменным d(, d2, di. Таким образом, существование вторых производных функционала (3.5) заведомо обеспечено, что делает возможным применение метода Ньютона для поиска его минимума.
После определения искомых кинематических параметров среды, доставляющих минимум функционалу (3.5), вектор нормали п может быть вычислен либо по матрице градиента смещений U, либо по фазовой матрице 4у. Ранее было показано, что для этих матриц и для двух векторов т _1_ 1, касательных к границе раздела двух сред, выполнены равенства UT = UI = 0, Фт = 1 = 0, (3.13) из которых следует, что каждая строка тензоров U и Ч ортогональна границе и, следовательно, коллинеарна вектору нормали. Таким образом, для нахождения нормали достаточно взять любую ненулевую строку матрицы U или 4і, нормировать ее и задать направление, удовлетворяющее требованию п, О.
Алгоритм векторной инверсии, исследование устойчивости
Конечной целью процедуры миграции сейсмических данных является построение изображения среды. Это изображение должно отражать структурные особенности строения сейсмических границ, а динамически представлять упругие свойства горных пород, то есть амплитуда в каждой точке получаемого разреза количественно должна описывать тот или иной физический параметр материала.
Первое условие - адекватное отображение структурных особенностей -достигается путем использования при миграции как можно более полной и точной скоростной модели среды. Если модель среды, задаваемая при продолжении наблюденного волнового поля, кинематически отвечает реальной картине, все элементы исходного поля мигрируют к правильным положениям породивших их рассеивающих неоднородностей.
Обеспечение второго условия — восстановление распределения параметров упругой среды в исследуемой области - является задачей динамической инверсии. Отметим, что для одномерных сред обратная динамическая задача полностью решается в акустическом варианте методами оптимизационной динамической инверсии отраженных волн, в результате чего может быть получено распределение акустических импедансов [33, 72]. Для двух- и трехмерного варианта также разработаны методы упругой динамической инверсии. Одни из них базируются на определенных допущениях о решении прямой задачи распространения колебаний [45, 64], другие используют численное решение [65, 73]. Однако все такие методы, как правило, совмещают процесс инверсии с процедурой продолжения волнового поля. Это относят к недостаткам алгоритмов трех- и двухмерной обратной динамической задачи, так как приходится ограничиваться приближениями, снижающими точность результатов и устойчивость решений [22].
В настоящей работе предлагается метод построения сейсмического изображения околоскважинного пространства, основанный на разделении процедур экстраполяции волнового поля и инверсии и осуществляемый по следующей схеме (рис. 21).
Волновое поле, зарегистрированное на профиле наблюдения, продолжается во внутренние точки среды, описываемой опорной скоростной моделью, на основе решения задачи обратного распространения упругих колебаний методом конечно-разностных схем (см. гл. 2). В процессе экстраполяции в каждый момент времени поле смещений сохраняется для множества вертикальных профилей, располагаемых в исследуемой области с некоторым заданным шагом. По завершении процесса продолжения наблюденного сейсмического поля будет получен набор волновых полей ВСП на профилях, покрывающих исследуемую область.
Далее необходимо провести селекцию рассчитанных описанным способом волновых полей с целью получения многокомпонентных трасс однократных отражений для всех основных типов волн. Эти трассы подаются на вход процедуре векторной инверсии, результатом которой является распределение вдоль профиля коэффициентов отражения по нормалям к сейсмическим границам.
Таким образом, из всех заданных внутри исследуемого пространства вертикальных профилей могут быть составлены разрезы нормальных коэффициентов отражения продольных и поперечных волн. Отличительной чертой получаемых сейсмических изображений является то, что динамически они представляют реальный перепад акустических свойств среды в независимости от углов падения лучей на освещаемые границы раздела.
