Содержание к диссертации
Введение
Глава 1. Математические модели распространения волн в трещиноватых средах 10
1.1 Моделирование рассеяния акустических волн на неоднородностях типа трещин 12
1.2 Моделирование рассеянного поля в длинноволновой асимптотике ... 14
1.3 Затухание падающего поля в трещиноватой среде 17
1.4 Комплексная модель рассеяния 20
Глава 2. Численное моделирование процессов рассеяния волн на неоднородностях типа трещин 28
2.1 Численные методы решения скалярных волновых уравнений 29
2.2 Поглощающие граничные условия для разностных схем 33
2.3 Моделирование рассеяния волн от областей со случайно распределенными трещинами 35
2.4 Определение областей повышенной трещиноватости методом сейсмической локации бокового обзора (СЛБО) 40
2.5 Влияние параметров трещин на результаты исследований методом СЛБО 45
Глава 3. Поиск локальных неоднородностей на основе решения задачи определения коэффициента поглощения среды 52
3.1 Задача распространения акустических волн в неоднородной поглощающей среде в обобщенной постановке 54
3.2 Обратная задача об определении коэффициента поглощения 66
3.3 Решение обратной задачи вариационным методом 73
3.4 Результаты численного моделирования 77
Глава 4. Томографический подход к решению задачи определения локальных неоднородностей 83
4.1 Постановка прямой и обратной задачи лучевой томографии 91
4.2 Решение прямых задач для уравнения эйконала и переноса 94
4.3 Решение задачи определения коэффициента преломления среды 98
4.4 Результаты численного моделирования задач сейсмики в лучевом приближении 103
Заключение 109
Список литературы
- Моделирование рассеянного поля в длинноволновой асимптотике
- Поглощающие граничные условия для разностных схем
- Обратная задача об определении коэффициента поглощения
- Решение задачи определения коэффициента преломления среды
Введение к работе
Одной из основных задач разведочной геофизики является поиск мест скопления углеводородного сырья в геосреде [20]. Такие области обычно называют углеводородными коллекторами. Они представляют интерес главным образом для определения мест бурения промысловых скважин с целью извлечения повышенных для данного района объемов углеводородного сырья. Поэтому при разведке нового месторождения основной задачей реализуемых на нем методов является поиск геологической структуры, которая с большой вероятностью может быть коллектором [20]. На рис.1, приведен результат работы одного из геофизических методов (метод общей глубинной точки [74]), применяемого для изучения структуры геосреды. Здесь хорошо видно, что часть одной из сейсмических границ образует небольшой купол, который может служить местом скопления углеводородов. Как правило, при разработке нового месторождения первые скважины закладывают именно над такими куполами.
Рис. 1. Изображение среды, полученное по методу общей глубинной точки.
предполагаемое место скопления
Размер и форма углеводородного коллектора сильно зависит от геологических условий района. Обычно коллектор с повышенным содержанием углеводородов является локальным образованием [16], т. е. основной объем коллектора сосредоточен в сравнительно ограниченном пространстве. В случае, если же коллектор имеет протяженную структуру, основные запасы углеводородов концентрируются в локальных областях его объема.
Многие современные методы прямо или косвенно направлены на поиск мест залегания углеводородного сырья. Подходы реализуемые для решения этой задачи преимущественно основываются на использовании свойств (физических, химических, и. т. д.) коллекторов, отличающих их от других объектов геосреды. Как видно из рис. 1, одним из наиболее простых способов является выделение коллектора по его форме. Другим отличительным свойством коллектора {например, карбонатного) является повышенная концентрация трещин [11] по сравнению с близлежащими участками геосреды. На основе этого факта разработан ряд методов направленных на поиск областей повышенной трещиноватости [47], рис. 2.
Рис. 2. Вертикальный разрез поля трещиноватости по методу СЛБО. 4
Следует так же отметить, что в пространстве коллектора происходит повышенное затухание амплитуды проходящей через него волны. Прежде всего это связано с рассеянием энергии волны на неоднород-ностях [48], но также существует ряд иных механизмов ее поглощения [81]. Наиболее характерным свойством коллектора является значительное изменение скорости распространения волн (как правило понижение) в его объеме по сравнению соседними участками геосреды, рис. 3. По этой причине результаты томографических методов [19, 59] представляют повышенный интерес с точки зрения определения местоположения и формы коллектора.
Свойства коллектора зависят в основном от геологических характеристик и условий района, но описанные выше особенности коллектора выделяются в большинстве случаев. Таким образом, можно ввести формализованное понятие углеводородного коллектора как локальной неоднородности, характеризующейся следующими свойствами: резкое изменение скорости распространения волн (как правило уменьшение), высокая концентрация неоднородностей типа трещин, повышенное затухание волн по сравнению с окружающими областями геосреды. Такая формализация не является строгой (так как учитывает только некоторые из свойств коллектора) и носит качественный характер. Однако она позволяет сформулировать ряд математических задач по определению локальной неоднородности геосреды. Более того, рассмотрение и решение таких задач позволяет как смоделировать некоторые из существующих методов по поиску коллекторов, так и оценить потенциальные возможности новых подходов.
О 20 40 SO вО
Рис, 3. Пример скоростного разреза, полученного по методу межскважинного прозвучивания.
Цель данной работы заключается в анализе возможностей определения локальных неоднородностей среды путем использования каждого из введенного в формализации свойств в отдельности.
В первой главе работы рассмотрены наиболее актуальные на сегодняшний день подходы к моделированию трещин, размеры которых много меньше длины падающей волны. В общем случае именно такие трещины заполняют объем карбонатного коллектора. Рассмотренные методы позволяют провести моделирование рассеянного поля, образующегося на ансамблях трещин, при прохождении волны через коллектор. Предложенные подходы рассматриваются применительно к двумерной среде, что обусловлено невозможностью провести численные эксперименты в трехмерном случае по причине больших вычислительных затрат. Трещины представляются линейными отрезками и характеризуются своей длиной и ориентацией на плоскости.
Во второй части первой главы предложен подход к моделированию распространения волн в трещиноватых средах, позволяющий перейти от конечного числа рассеивателей к их континуальному распределению. Для этого в рассмотрение вводится функция являющаяся удельным объемом рассеивателей и входящая в волновое уравнение как ко-
эффициент. При этом сами трещины представляются в виде бесконечно тонких дисков, характеризующихся своей ориентацией в пространстве и радиусом. В основу модели положен механизм рассеяния на одиночных трещинах, описанный в [50]. Следует отметить, что в рассматриваемой модели можно учесть диссипативные потери в падающем поле, вызванные рассеянием на трещинах, которыми пренебрегали в предыдущих работах [57, 75]. Новизной предложенного метода является механизм взаимодействия падающего и рассеянного полей, который приводит к закономерностям, наблюдаемым при физическом моделировании.
Вторая глава исследования посвящена математическому моделированию метода сейсмической локации бокового обзора (СЛБО) [47]. Этот метод используется в геофизических исследованиях для определения областей повышенной трещиноватости в геосреде путем накопления рассеянной от скоплений мелких (10-100 раз меньше длины падающей волны) трещин энергии, зарегистрированной на земной поверхности.
Разработчики метода СЛБО напрямую связывают энергию рассеянных волн со степенью трещиноватости, что, как будет показано в работе, справедливо на качественном уровне. Помимо этого существует ряд нерешенных вопросов относительно использования метода, которые препятствуют его применению в промышленных масштабах. Среди основных можно выделить: возможность выделения трещиноватых областей по методу СЛБО, соответствие полученного результата реальному распределению трещиноватости, влияние ориентации, размера и концентрации трещин на энергию рассеянных волн. Ответы на эти вопросы помогут понять как научную, так и промышленную ценность метода СЛБО, а с точки зрения данного исследования оценить возможность выделения локальных неоднородностей по повышенной в них концентрации трещин. Для этого в первой главе проводится как моделирование процессов рассеяния на ансамблях мелких трещин (что соответствует решению прямой задачи в методе СЛБО), так и моде-
лирование алгоритма восстановления трещиноватости среды по методу СЛБО (что соответствует решению обратной задачи в методе СЛБО).
В третьей главе работы проводится исследование задачи распространения волн в однородной среде с поглощением. Механизм поглощения энергии в пространстве коллектора объясняется различными гипотезами. Однако экспериментально установлено [48], что одной из причин такого поглощения является рассеяние волн на скоплениях трещин. Во введении ко второй главе показано, как можно свести задачу о рассеянии волн на трещинах к задаче распространения волн в среде с поглощением. Установлено, что если пренебречь рассеянием на жестких трещинах, то коэффициент в волновом уравнении характеризующий интенсивность рассеяния совпадает с коэффициентом поглощения. Таким образом, повышенные значения последнего свидетельствуют о возможном наличии локальных неоднородностей в области зондирования. Так как поглощение энергии волны в пространстве коллектора происходит неравномерно (дискретно), то естественно предположить, что коэффициент поглощения принадлежит пространству интегрируемых (с квадратом) функций. По этой причине прямая и обратная задачи формулируются в обобщенной постановке. Обратная задача заключается в определении коэффициента поглощения по рассеянному полю, зарегистрированному на земной поверхности. Для ее решения применяется вариационный метод минимизации функционала невязки между зарегистрированным и полученным в результате решения прямой задачи волновыми полями.
В четвертой главе работы сформулирована задача сейсмической томографии и предложен один из подходов к ее решению. Как было отмечено ранее, наиболее характерным свойством коллектора является резкое изменение скорости распространения волн в его объеме по сравнению с соседними участками геосреды. На сегодняшний день широко распространены томографические подходы основанные на межсква-
жинном прозвучивании [19], то есть когда излучение и прием сигнала производятся в скважинах. Однако такие методы очень затратны и дают информацию о скоростном разрезе лишь между двумя скважинами. В данной главе рассмотрен подход, в котором излучение можно проводить с земной поверхности, при этом регистрация сигнала производится в скважине. Используемый подход, с одной стороны, основан на регулярности поля лучей, с другой, позволяет отказаться от системы прямых, рассматриваемых в методе радоновской [30] томографии. Для нахождения скоростной характеристики среды применяется метод линеаризации нелинейных операторных уравнений. По схеме решения задачи Радона на основе проекционной теоремы построен численный алгоритм, который, в случае вырождения среды в однородную, обладает всеми свойствами проекционной схемы. Поскольку уже в случае радоновской томографии обратная задача становится некорректной, то не приходится ожидать получения результатов решения обратной задачи для неоднородных сред, характеризуемых удовлетворительной точностью. Однако для решения задач поиска локальных неоднородностей геосреды достаточно лишь качественного описания скоростной модели. Все поставленные в работе задачи решаются способами математического моделирования: для их решения на компьютере моделируется физический процесс, который наблюдается в реальных условиях, и регистрируются параметры, использующиеся для решения обратной задачи. В качестве регистрируемого параметра в нашем случае выступает волновое поле, а в качестве восстанавливаемого — характеристики среды. Решение обратной задачи реализуется на компьютере в виде отдельной программы, на вход которой подаются данные из прямой задачи. Таким образом, математическое моделирование представляет собой процесс повторения физического моделирования на компьютере, являясь при этом более гибким и дешевым способом.
Моделирование рассеянного поля в длинноволновой асимптотике
Более адекватной природным условиям с точки зрения частотных зависимостей является модель закрепленной (жесткой) трещины. В этом случае задачу рассеяния на трещинах можно сформулировать следующим образом: V2psc(r, и) + psc{v, и) = 0, vD (1.9) nmVpsc(v, и) = -nmVpinc(v, w), г Є D, (1.10) где индекс m соответствует m-ой трещине. Как было показано в п.1.1.2. рассеянное поле от закрепленной трещины ищется в виде psc{r)= J nVG(r, i (r ) dr . (1.11) r eD
В силу того, что мы рассматриваем трещины, размер которых много меньше длины падающей волны можно пренебречь изменением функции Грина по поверхности трещины и вынести ее за знак интеграла [17].
В результате получим r(r) = nmVG(r,rm) / (r )dr . (1.12) VeDm Так как функция Грина G(r,r ) является решением уравнения Гель-мгольца (1.1) с правой частью s(r, г ) = 5(г — г ), то из представления (1.12) следует, что рассеянное поле psc удовлетворяет следующему уравнению V2psc(r,o;) + sc(r,a;) = n-V(5(r-r-) / ) - (1-13) m " m г єД
В работе [75] показано, что при выполнении следующих условий: 1) ка «С 1, где к = UJ/C, а — размер трещины, т. е. размер трещины мал по сравнению с длиной падающей на нее волны, 2) рассеяние второго порядка и выше не учитывается, можно вычислить интеграл, входящий в правую часть уравнения (1.13): / p(r ) dr й шІЇт, (1.14) r eDm где qm(u) = nmWnc(rm,w). (1.15) После перехода из частотной области во временную уравнение (1.13) преобразуется в уравнение: VVC(r, ) - Со2- гМ) = X m fa( )nmV (r _ rJ} (U6) m где qm(t) = nmVpinc(rm,t). (1.17) Отсюда следует, что задача рассеяния волн сводится к решению волнового уравнения (1.16) с соответствующим источником, который зависит от ргпс. Алгоритм для расчета рассеянного поля выглядит следующим образом: а) вычисляется падающее поле ргпс, б) по полю ргпс для каждой m-ой трещины рассчитывается функция 9m(t), в) вычисляется рассеянное поле psc(r, t) путем решения волнового урав нения (1.16).
Таким образом, предложенный алгоритм позволяет явно не учитывать граничное условие на трещине, а сама трещина выражается в виде точечного источника в волновом уравнении. Однако источник, входящий в уравнение (1.16), невозможно использовать для численного моделирования, так как в него входит градиент дельта-функции. Поэтому требуется найти представление, позволяющее аппроксимировать его на сетке.
В работе [75] показано, что при выполнении условий 1) и 2) функции р(т,ш) и ф(г,и), входящие в уравнения (1.5) и (1.8), можно представить в виде произведения двух функций а (г) и b(w), где функция а(г) известна и определяется в каждом отдельном случае. Функция b{uS) находится из соответствующего интегрального уравнения. Поэтому можно выписать явное представление для рассеянного поля, которое в случае закрепленной трещины выглядит следующим образом: й(г, и) = ш1П f-jf H?\ka\v - г„)ь,М, (1.18) где &о — w/c, Щ - функция Ханкеля 1-го рода. Рассмотрим в частотной области разностный волновой оператор, имеющий порядок аппроксимации 0(Аш\ Аг10) [63]: eh = _к2 + ( _ [дЛ + Mf!!(Afc)2 12 - - , v , , (1.19) где As = Ах = Az, fiij = CijAt/As. В работе [57] показано, как путем применения этого оператора к равенству (1.18) и последующим переходом во временную область можно построить разностную схему для расчета psc с источником следующего вида: SM) = 2ш2гп[ Mr m)qm{t) + A2(r,rm) - f() + J44(r,rm)i (i)), (1.20) который и используется при проведении численных экспериментов. Коэффициенты AQ, А2, АІ не зависят от времени и быстро убывают от центра трещины гт. Их вид представлен в работе [57].
Для получения рассеянного поля прямую задачу приходится решать дважды: первый раз - для расчета волнового поля в отсутствии трещин ргпс, второй - для расчета поля рассеянного на трещинах psc. При расчете ргпс источник в волновом уравнении задается произвольно, при расчете psc - определяется формулой (1.20).
Поглощающие граничные условия для разностных схем
При моделировании волновых процессов возникает проблема появления отраженных от границ вычислительной области волн. Так как эти отражения являются следствием ограниченности вычислительной сетки и не образуются в реальных условиях, требуется применять методы их подавления. Наиболее распространенным методом подавления отраженной энергии является задание на границах вычислительной области условия следующего вида [67]: dp др (2.7) - + с = 0. от ох Условие (2.7) описывает распространение волн только в одном направлении. Замена волнового уравнения на это условие в приграничных точках вычислительной области гарантирует отсутствие отраженной энергии. Следует заметить, что условие (2.7) обеспечивает полное поглощение волн только лишь по направлению нормали к границе, и чем больше угол падения волны составляет с нормалью, тем слабее поглощается энергия.
В случае, когда волна падает на границу под известным углом можно использовать условие ! + c(cosa)- g = 0; (2,) которое поглощает энергию в направлении (cos a, sin а), где а — угол между нормалью и направлением падения волны.
Так как в вычислительных экспериментах, как правило, моделируется распространение сферической волны, которая падает на границу области под всевозможными углами, то условия (2.7) и (2.8) не обеспечивают удовлетворительного поглощения энергии по всем направлениям. Поэтому следует дополнительно применять другие способы поглощения энергии на границах вычислительной области. Эффективным является условие поглощения, предложенное в работе [71], которое в случае рассмотрения правой границы сетки выглядит следующим образом: М=Лм + №)-Є{х){ії+Сд-Х)- (2-9)
Коэффициент поглощения є(х) должен принимать наименьшее значение внутри области вычисления и постепенно увеличиваться до границы области, где он достигает своего максимума. Обычно для его задания используют 15 — 20 точек сетки. Максимальное значение коэффициента є(х) выбирается в зависимости от числа Куранта и, как правило, меняется от 0.05 до 0.1.
Одним из наиболее простых в смысле реализации условий поглощения является применение на границах области вычисления так называ емого "слоя-губки"[73]. В каждой точке этого слоя рассчитанное волновое поле умножается на убывающую функцию, максимальное значение которой достигается на внутренней границе слоя, а минимальное — на границе вычислительной области. Толщина слоя (обычно 20 — 40 точек) выбирается в зависимости от вида убывающей функции. Хотя такой способ и является наиболее удобным в смысле реализации, он требует задания большого числа точек поглощающего слоя, что сильно сказывается на времени вычислений.
При моделировании волновых процессов приходиться применять несколько способов подавления отражений от границ сетки одновременно, так как каждый из них в отдельности имеет свои недостатки.
Моделирование проводилось на сетке с шагом по пространственным координатам Ах = Az = 5 м и по временной координате At = 0,001 с, размер сетки равнялся 300 X 300 точек по координатам х и z соответственно, скорость распространения волн была равна 2500 м/с. В области с трещинами каждый узел сетки содержал трещину, длина которой менялась случайным образом от 1 см до 200 см, также случайно менялась ориентировка трещин (на плоскости). В первом приближении такая область соответствовала модели среды со случайно распределенными трещинами. Размеры областей рассеяния (здесь и дальше так будут обозначаться области с трещинами) приведены в таблице 2.2.
Схема проведения экспериментов приведена на рис. 2.2. Центры областей рассеяния на всех моделях находятся в центре вычислительной сетки. Положение излучателя относительно верхней границы области менялось для каждой модели так, чтобы расстояние от излучателя до верхней границы области рассеяния оставалось постоянным. В качестве импульса падающей волны (в дальнейшем будет называться зондирующим импульсом) использовался импульс Риккера с доминирующей частотой 25 Гц, определяющийся по следующей формуле:
Обратная задача об определении коэффициента поглощения
Подобное рассмотрение обратной задачи возможно и в том случае, когда источник в краевом условии (3.11) имеет вид г;(0, t) = s(t), s(t) Є І 2[0,Т]. В сейсморазведке, как правило, используются импульсные источники сравнительно малой длительности, удовлетворяющие последнему условию. Пусть также известен не точный след решения f(t) задачи с источником s(t), а его приближение f(t) такое, что ИЛ /L2[O,T] єі гДе величина є известна. Задача минимизации функционала «/(//) в 1 2[0, Т/2] в этом случае является некорректной, поскольку точка минимума может определяться неустойчиво.
Рассмотрим вариант метода регуляризации А. Н. Тихонова [40] для решения операторных уравнений с нелинейным оператором, переводящим слабо сходящиеся последовательности в сходящиеся. Для задачи (3.9 - 3.12) при v(0,t) = s(t) таким оператором является оператор А : 1/2[0,Т/2] — Ьг[0, Т], определяемый равенством А\ц{х)] = «"( ). При этом функционал Тихонова имеет вид Ма( ) = ) + а/2МІ2[от (3.27) Из Теорем 2.2.3-2.2.4 [15] вытекает Следствие 3.9. Функционал Ma(ii) достигает в г[0, Г/2] своей точной нижней грани.
Если { } — множество точек минимума функционала Ма, то каждая из них является слабым приближенным решением обратной задачи (3.9 - 3.12) при соответствующем согласовании параметров а я є. Последняя проблема подробно изучена в литературе, см. [41, 45] и др.
Для решения обратной задачи может быть использована априорная информация об искомом решении вида L2[O,T/2] R или 1И1с[о,г/2] R. Соответствующие множества являются слабыми компактами в пространстве 1/2[0,Г/2], что позволяет построить слабые квазирешения /І как решения вариационной задачи /z = arg miii J(fi), m\\ R при этом существование точки минимума следует из приведенных при обосновании Следствия 3.9 рассуждений. Заметим, что в силу природы рассматриваемой задачи (3.9 - 3.12) не представляется возможным делать более сильные предположения о функции (i(x) такие как монотонность, выпуклость, дифференцируемость или ограниченность в соответствующих множествах.
Для практического решения обратной задачи (3.9 - 3.12), т. е. в условно - корректной постановке, был выбран метод проекции градиента [7] для вариационной задачи /і = argmin J([i), MR = {її Є L2[0,T/2] : 0 ц R). MR
Решение этой задачи строится итерационно по следующему правилу: Д +і = Рмп&к - РкУць + Qtfc/ jfe)), где PMR — проектор на множество MR, а последовательности а , р& таковы, что ak,Pk,pk/ctk - 0 и \ak+i - 0 14 - 0 при к -» оо. Например, а ; — fc-1/3, pk = fc-1/2, см. [7]. Проектор PMR реализуется в виде срезки, т. е. Рмпу {х) = min{max{0,/i}, R}.
Для вычисления J x) используется приближенная формула, следующая из (3.26): х Т - J x) « д(2х) ехр{-2 / fi(y) dy} + / [2v(x,t - х) О 2х t-x + / u(x,t — T)u(x,r)dT]g(i)dt, х учитывающая члены до второго порядка по степеням //, где, соответственно, х й{х, t) = -fi((t + )/2)/2, v(x, t) = \J KtMt + ( - )/2) de На практике при решении задач с распределенным источником рассматриваются фурье-образы приведенных выражений, которые вычисляются с помощью алгоритма быстрого преобразования Фурье, а параметры итерационного процесса согласуются между собой и погрешностью входных данных в соответствии с методом итеративной регуляризации [7].
Для оценки предложенного метода был проведен ряд вычислительных экспериментов по восстановлению модельных коэффициентов поглощения. При этом рассматривались как типичные для геофизики коэффициенты, так и те, которые представляют интерес в математическом смысле. В силу того, что на практике выбор дельта функции в качестве источника невозможен, в некоторых экспериментах использовался финитный импульс сравнительно малой длительности.
Решение задачи определения коэффициента преломления среды
При проведении томографического эксперимента на апертуре приема регистрируются времена и амплитуды первых вступлений волн. Например, при мо делировании метода межскважинного прозвучивания известными можно считать функции т(х,0) и А(х,0). Покажем, как по этим функциям можно найти функцию t(l,6). Так как в кинематической сейсмике функция т(х,0) называется годографом первых вступлений, то задача нахождения функции t(l,9) является задачей параметризации годографа, которая, в свою очередь, является весьма важной и самостоятельной задачей интерпретации данных.
При излучении плоской волны известен угол ее выхода в, поэтому задача параметризации годографа сводится к нахождению точки выхода I. Ранее было установлено, что Цх,у;в) = А\х,у-в)ту = ±А2{х,у;9) п2(х,у) - т2(х,у;9), где обозначения / = 1(х, у;9),т = т(х, у;9),А = А(х, у; 9) подчеркивают параметрическую зависимость от 9. Так как 1Х знакопостоянна в силу регулярности поля лучей, то остается лишь одно значение корня, т. е. 1х(х,у; 9) = А2{х, у; 9)ту = А2(х,у; 9) /п2(х,у)-т2(х,у;9). Пусть для определенности /(0,0; 9) = 0, тогда 1х(х, 0; 0) = А2(х, 0; 0)л/п2(х, 0) - т2{х, 0; 0), откуда следует, что х l(x, 0;9) = J А% 0; 0) ,0)- ,0 . (4.8) о Так как функция 1(х, 0; 9) монотонна по х, то существует обратная к ней х = х(1,9). Искомое параметрическое представление годографа имеет вид , 0) = ф(/, 0),0).
Заметим, что при излучении сферической волны известен вынос /, а угол 0 находится аналогично нахождению I для плоской волны.
Расчет поля времен. Рассмотрим задачу нахождения функций т(х,у) и А(х,у) по заданной функции п(х,у). Один из способов заключается в численном решении волнового уравнения [63]. Практика вычислений показывает, что такой подход оказывается неэффективным с точки зрения затрат вычислительного времени, поэтому для нахождения т(х, у) использовался численный метод решения уравнения эйконала [83, 82, 72]. Заметим, что в случае плоских волн уравнение эйконала решалось в декартовых координатах [72], а для сферических волн в полярной системе координат; рассчитанное при этом поле времен затем переводилось в декартову систему координат.
Перейдем к рассмотрению алгоритма численного решения уравнения эйконала. В полярной системе координат уравнение выглядит следующим образом: rr+ = ( ). Выражая тг (с учетом роста г вдоль луча), находим Tr = \jn2{r,ip)- T$.
Пусть г = гАг, ip = jA(p, где Аг — шаг по переменной г, Аїр — шаг по переменной ip , г и j — целочисленные индексы, изменяющиеся от 0 до N. Заменив частные производные по г и р центральными разностными соотношениями, получаем / Аг 1 Ц+и = П-ij + W4n?,Ar2 - д Кя-і - -i)2p 2 + (Af2 + А 2) Для использования этой явной схемы необходимо задать значения функций ТфТ{ -1,то ,Т1 . На рис.4.2а изображена скоростная характеристика среды v(x,y) и полученное для нее поле времен, рис. 4.26.
Расчет поля лучей. Для решения обратной задачи и нахождения амплитуд в прямой задаче требуется знание поля лучей. Из того, что изолинии поля времен ортогональны изолиниям поля лучей, следует, что VrV0 = О, где в(х,у) — поле лучей. В полярной системе координат это условие выглядит следующим образом: тгвг + - Тфвц, = 0. Выразим из этого уравнения $г как вг = — —вф. г1 тг Заменим частные производные по г и ip центральными разностными соотношениями. В результате получаем
Для использования этой явной схемы необходимо задать значения функций 9oj,6ij,9ito, Qi,N-\- На рис. 4.2а изображена модель среды и полученное для нее поле лучей, рис. 4.2в.
Расчет поля амплитуд. Поле амплитуд можно рассчитать несколькими способами. Во-первых, его можно найти из уравнения переноса, зная решение уравнения эйконала, [78]. Однако такой способ на практике используется редко, так как в уравнение переноса входят вторые частные производные поля времен и их вычисление приводит к значительным ошибкам. Гораздо более эффективным оказывается определение поля амплитуд из уравнения (4.7). При этом сначала численно рассчитывается поле времен, затем поле лучей и находится поле амплитуд из уравнения (4.7). Пример поля амплитуд, полученного для модели, изображенной на рис. 4.2а, приведен на рис. 4.2г.
