Сейсмическая томография земли (стр. 2 )

Задача томографии по первым вступлениям проходящих волн решается, как правило, в рамках модели геометрической сейсмики. Процесс распространения волн рассматривается как распространение слабых разрывов в решении уравнений теории упругости. При этом слабые разрывы могут существовать только на характеристиках исходных уравнений. Характеристическое уравнение для уравнений теории упругости, записанное для времён пробега, называется уравнением эйконала. Время прихода волны из источника в приёмник зависит лишь от свойств среды вдоль фиксированной криволинейной траектории, называемой лучом. Лучи являются характеристиками уравнения эйконала и бихарактеристиками уравнений теории упругости. Сейсмические волны, возникающие в реальных физических процессах, вообще говоря, не удовлетворяют указанному выше математическому определению волны. Законы геометрической сейсмики подходят для описания реальных физических волн в случае высоких частот. Чем меньше длина волны по сравнению с характерным размером неоднородности среды - тем точнее процесс распространения волн может быть описан с помощью лучей. Глубокое, математически обоснованное понимание связи геометрической сейсмики с реальными волновыми процессами пришло вместе с развитием лучевого метода , , [Алексеев и др., 1958; Алексеев, Гельчинский, 1959] в середине 50-х годов прошлого столетия. Линейные соотношения, связывающие возмущения времён пробега с возмущениями параметров среды, могут быть получены с помощью линеаризации уравнения эйконала. При этом оказывается, что возмущения времён связаны с возмущениями скорости вдоль лучевых траекторий. Для изотропной среды требуемые соотношения сводятся к интегрированию приращения величины обратной скорости вдоль лучевой траектории. Для анизотропной среды данные соотношения выведены для уравнений акустики [Романов, 1972]. Для уравнений теории упругости в анизотропной среде нам не удалось найти требуемые соотношения в литературе. В приложениях чаще всего рассматривается трансверсально-изотропная среда.

При решении обратной кинематической задачи возникает необходимость решения прямой кинематической задачи расчёта времён пробега сейсмических волн в среде с заданными свойствами. При решении задачи сейсмической томографии в рамках геометрической сейсмики достаточно построить все возможные лучи, соединяющие источники с приёмниками, и вычислить время пробега вдоль каждого луча. Традиционно прямая кинематическая задача решается методами, использующими трассировку лучей. Возможны два варианта постановки: 1) задача Коши, заключающаяся в построении луча, выходящего из заданной точки в заданном направлении; 2) двухточечная краевая задача о построении луча, соединяющего две фиксированные точки. При решении задачи Коши возникает система уравнений луча, которая быстро и точно интегрируется методом Рунге-Кутта.

Для задачи сейсмической томографии необходимо строить лучи, соединяющие источники с приёмниками, то есть много раз решать задачу двухточечного трассирования. Решение двухточечной задачи является более сложной задачей. Простейший способ её решения - метод пристрелки, суть которого заключается в многократном решении задачи Коши для источника с целью подбора направления, для которого луч попадёт в приёмник [Peryra et al., 1980; Sun, 1993]. Более широкое распространение получили методы изгиба лучей, позволяющие подбирать геометрию луча для фиксированной пары источник-приёмник с помощью минимизации функционала Ферма [Оболенцева, 1988; Um, Thurber, 1987; Moser et al., 1992; Thurber, Kissling, 2000]. Плохая устойчивость в средах с большими контрастами скорости и большие вычислительные затраты не позволяют применять традиционные алгоритмы двухточечного трассирования в сложных трёхмерных средах [Peryra et al., 1980; Lions, 1982]. Конечно-разностные методы [Trier, Symes, 1991; Osher, Shu, 1991; Sethian, 1999] решения уравнения эйконала лишены недостатков, присущих методам двухточечного трассирования. Чтобы решить прямую кинематическую задачу, то есть вычислить времена прихода волны, необходимо решить уравнение эйконала, и для этого не обязательно строить лучи. После того как в некоторой области уже известны времена прихода, построение лучей, также необходимое для решения задачи сейсмической томографии, представляет собой техническую задачу.

Математическое обоснование численного решения уравнения эйконала с помощью конечно-разностных схем дано в работе французского математика [Crandall, Lions, 1983], в которой рассматривается обобщённое решение уравнения эйконала. Использование математической концепции обобщенного решения является существенным, поскольку уравнение эйконала может не иметь всюду дифференцируемого классического решения. Обобщённое решение не единственно, например, при наличии сильной рефракции. Лионс и другие авторы формально выделяют особое обобщённое решение уравнения эйконала, которое соответствует временам первых вступлений. Поскольку данное обобщённое решение обладает устойчивостью по отношению к параметрам среды, положению источника и т. д., оно может быть найдено с помощью конечно-разносных схем. Это особое обобщённое решение называется вязким решением (viscosity solution) [Lions, 1982; Crandall, Lions, 1983]. Результаты Лионса [1982] относятся к изотропному случаю, однако логично предполагать [Qian, Symes, 2002], что для уравнения эйконала, описывающего распространение квазипродольных волн в случае трансверсально-изотропной среды, первые времена также соответствуют устойчивому обобщённому решению уравнения эйконала и поэтому могут быть найдены с помощью конечно-разностной схемы.

Конечно-разностные методы решения уравнения эйконала можно разделить на два типа: методы, распространяющие решение из заданной точки во всех направлениях, в иностранной литературе называющиеся "point expanding methods" [Kim, 2002], и методы, продолжающие в глубину заданное поле времён на поверхности, соответственно "box expanding methods". Формально начальное условия для уравнения эйконала в двумерном и трёхмерном случаях ставить в точке нельзя: задача имеет множество решений, поэтому в "point expanding"-методах начальные значения решения задаются на некоторой бесконечно малой сфере, окружающей источник.

Наиболее популярным из "point expanding"-методов является метод, получивший в иностранной литературе название fast marching [Sethian, 1996]. Суть алгоритма решения этим методом заключается в специальном алгоритме пересчёта значений поля времён в узлах расчётной сетки. Область с уже вычисленными временами расширяется за счёт добаления к ней той точки, прилегающей к её границе, значение времени прихода в которой, вычисленное на основе специальной разностной аппроксимации уравнения эйконала, минимально. Хотя в большинстве узлов возникает необходимость пересчитывать значения по нескольку раз, а также постоянно находить точку с минимальным временем из всех соседних к рассчитанной области, алгоритм требует O(Nln(N))-действий для вычисления времени в N узлах сетки. Для поиска минимума при этом используется эффективный алгоритм двоичной сортировки. Алгоритм, лежащий в основе метода fast marching, был первоначально предложен в 1959 г. нидерландским математиком Э. Дейкстрой для решения задачи о нахождении минимальных путей в графах [Dijkstra, 1959; Cormen et al., 2001]. Конечно-разностный метод fast marching решения уравнения эйконала был предложен американским математиком Джеймсом Сезианом [Sethian, 1996; Sethian,1999а]. В ряде работ рассматривается множество приложений этого метода в задачах теории оптимального управления, при обработке изображений, а также при решении задач геофизики [Falcone et al., 1994; Sethian, 1999b]. Существует множество модификаций метода fast marching, например, модификации метода для цилиндрической и сферической систем координат[Alkhalifah, Fomel, 1995]. Достоинствами метода являются возможность распространять временное поле во все стороны от источника при относительно небольших вычислительных затратах и универсальность метода: fast marching работает практически в любой среде, например в случае присутствия высококонтрастных включений типа соляных тел. На основе fast marching создан алгоритм быстрой сейсмической миграции [Cameron et al., 2006]. У метода fast marching два существенных недостатка. Метод вычисляет времена с невысоким первым порядком точности. При вычислении времён используется совпадение направления градиента временного поля с нормалью волнового фронта, а это верно только для изотропной среды. Fast marching не удаётся применить для анизотропной среды. В работе E. Cristiani [2009] предложено обобщение метода fast marching, получившее название buffered fast marching, которое может работать для анизотропных сред, но алгоритм получается весьма громоздким и более медленным. В работе E. Cristiani [2009] показаны примеры расчётов только для простейшего случая эллиптической анизотропии, применить данную методику для расчёта времён в более сложных случаях анизотропии затруднительно.

В большинстве задач сейсморазведки источники расположены на поверхности. Как правило, достаточно рассматривать только лучи, уходящие в глубину z. Поэтому применяют "box expanding"-методы, в которых временное поле, заданное на поверхности при z=0, пересчитывается в глубь среды. Если рассматривать только лучи, уходящие в полупространство, уравнение эйконала можно однозначно алгебраически разрешить относительно пространственной производной времени прихода по глубине. В результате получается уравнение Гамильтона-Якоби. Записывается разностная схема, в которой переменная играет роль времени, времена пробега пересчитываются вниз по слоям.

Широко распространённым в настоящее время является конечно-разностный метод решения [Vidale, 1988; Trier, Symes, 1991; Kim, Cook, 1999], основанный на ENO-аппроксимации (essentially nonoscillatory - существенно неосциллирующие) пространственных производных по горизонтальным переменным. Идея ENO-апроксимации состоит в использовании адаптивного шаблона, обеспечивающего высокий порядок в областях с гладким изменением скорости и избегающего по возможности больших градиентов скорости. При ENO-аппроксимации вычисляются значения производной для различных возможных шаблонов, а затем из них выбирается один наиболее подходящий. Выбор шаблона требует дополнительных вычислительных затрат и для схем высокого порядка может существенно замедлить вычисления. WENO-аппроксимация (weighted essentially nonoscillatory - взвешенная существенно неосциллирующая) производных по свойствам аналогична ENO-, но при WENO-аппроксимации значения производной на разных шаблонах складываются с определёнными, вычисляемыми по некоторому правилу весами, следовательно, нет необходимости сравнивать значения на различных шаблонах между собой. В результате этого WENO-схемы оказываются быстрее. После того как значения горизонтальных производных вычислены, для перехода на следующий слой по глубине используется метод интегрирования Рунге-Кутта, причём на каждом шаге приближённо решается задача Римана [Bardi, Osher, 1991]. Делать это можно разными способами [Osher, Shu, 1991; Bardi, Osher, 1991]. Наиболее простым способом решения задачи Римана является использование функции потока Ошера [Osher, Shu,1991; Bardi, Osher, 1991], как это и делается в работах [Kim, Cook, 1999; Osher, Shu, 1991]. Изложенная схема решения принадлежит классу сеточно-характеристических конечно-разностных методов, в основе которых лежит идея, использованная в схеме [Годунов, 1959; Годунов и др., 1976]. При данном подходе можно вычислить только времена, соответствующие лучам, уходящим в глубину и лежащим в некотором растворе углов. Эта особенность является следствием условия устойчивости подобных разностных схем, которое заключается в том, что характеристики должны лежать внутри расчётных ячеек. Выбирая шаг сетки по вертикальной оси значительно меньшим чем по горизонтальным направлениям, можно расширять диапазон углов, но это требует и дополнительных вычислительных затрат. В работе J. Qian, W. W.Symes [2002] предложен алгоритм, при котором шаг сетки на каждом слое выбирается адаптивно с использованием оценки невязки путём сравнения результатов вычислений для схем второго и четвёртого порядков точности. В работе S. Kim, R. Cook [1999] рассматривается модификация рассматриваемого конечно-разностного метода, позволяющая рассчитывать поле времён во все стороны от источника, последовательно решая уравнение в разных направлениях. К сожалению, для сред со сложными высококонтрастными включениями (соляные тела) этот метод всё же не применим: лучи слишком много раз меняют направление.