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



Скачать 440,31 Kb.
страница1/2
Дата29.03.2018
Размер440,31 Kb.
  1   2
СЕЙСМИЧЕСКАЯ ТОМОГРАФИЯ ЗЕМЛИ
Научно-популярное описание сейсмотомографического метода изучения внутреннего строения Земли. Примеры расшифровки внутренних структур Земли сейсмотомографическим методом. Материалы подготовлены при финансовой поддержке Минобрнауки РФ в рамках Соглашения о предоставлении гранта в форме субсидии №8668.

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

Кинематический подход исторически является более ранним. Исходными данными для обратной кинематической задачи являются времена пробега упругих волн из источников в приёмники, которые могут быть определены с высокой точностью. Первые результаты восстановления скоростного строения по временам пробега упругих волн от землетрясений получены в начале двадцатого века Герглотцем и Вихертом [Herglotz, 1907; Wiechert, Zoeppritz, 1907].Обратная кинематическая задача решалась аналитическими методами. Так, в работе Шлихтера [Schlichter, 1932] рассмотрена задача восстановления скоростного строения вертикально- неоднородной среды.

Дальнейшие усилия долгое время были направлены на развитие методов решения одномерной задачи. В частности, характер неоднозначности решения обратной кинематической задачи изучен М.Л. Гервером и В.М. Марушкевичем [Гервер, Марушкевич, 1965]. Одними из первых являются результаты для двумерной постановки обратной кинематической задачи, изложенные М.М.Лаврентьевым и В.Г.Романовым [Лаврентьев, Романов, 1966]. Основополагающими в области обратных кинематических задач являются работы А.С.Алексеева, М.М. Лаврентьева, В.Г.Романова и др., посвящённые исследованию кинематической задачи в линеаризованной постановке [Алексеев и др., 1969; Алексеев и др., 1979].

В современной научной литературе обратную кинематическую задачу часто называют задачей сейсмической томографии. Компьютерной томографией называется численное восстановление функций по их криволинейным или поверхностным интегралам, которое находит применение в различных областях науки и техники [Натеррер, 1990]. Впервые о компьютерной томографии заговорили в начале 70-х годов прошлого столетия в связи с рентгенодиагностикой в медицине [Натеррер, 1990; Ganzkörper – Computer - Tomographie, 1977].

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

В современной геофизике обратная задача кинематики решается численно. Применяется стандартная техника поиска минимума целевого функционала, представляющего собой норму разности измеренных данных и данных, рассчитанных с помощью решения прямой задачи для текущей модели среды. Первые работы в этом направлении появились в конце 70-х - начале 80-х годов [Dines, Lytle, 1979; Pederson et al., 1985].

Известно, что многие горные породы обладают анизотропными свойствами, и скорость распространения сейсмических волн зависит от направления [Thomsen, 1986]. Закон Гука, уравнения, описывающие распространение волн, законы преломления, трассировка лучей в анизотропной среде существенно отличаются от изотропного случая. В частности, широко распространённые глинистые сланцы обладают анизотропными свойствами вследствие своей микроструктуры, состоящей из тонких слоёв. Распространение сейсмических волн в сланцах может быть описано в рамках модели трансверсально-симметричной среды с вертикальной осью симметрии. Методика построения эффективных трансверсально-изотропных моделей, описывающих тонкослоистые среды, представлена в работах Backus [1962] и Berryman [1979]. Модели для глин и сланцев рассматриваются в [Sayers, 2005] и [Xu, White, 1995]. Методы введения анизотропии для трещиноватых сред разработаны Hudson [1981] и Thomsen [1995].

Традиционно обработка сейсмических данных проводится в рамках модели изотропной среды. Для небольших выносов (меньше 30°) и незначительной анизотропии "изотропные" методы работают достаточно хорошо, но в ряде случаев неучёт анизотропии приводит к значительным ошибкам. В работе Levin [1990], Гурвич [1940] проведены модельные исследования ошибки, возникающей при неучёте анизотропии в методе отражённых волн. Задача томографии в анизотропной среде изучена слабо и недостаточно представлена в литературе.

Известны три основных способа численного решения задачи томографии: методы матричного обращения, методы преобразования Фурье, алгебраические процедуры восстановления [Хаттон, 1989].

Методы матричного представления основаны на линеаризации задачи, в результате которой получается линеаризованный томографический оператор. Возмущения времён пробега при этом зависят от возмущений среды вдоль лучевых траекторий в невозмущённой среде. Целевая область дискретизуется, как правило, с помощью прямоугольной сетки с постоянным значением скорости в ячейке, берётся конечное число источников и приёмников. Исходный линейный томографический оператор, действующий в бесконечномерных пространствах, аппроксимируется конечномерным линейным оператором, действие которого представляется умножением соответствующей томографической матрицы на вектор из пространства моделей. Решением системы алгебраических уравнений находятся неизвестные возмущения параметров среды. После этого можно модифицировать начальную модель, заново построить лучи и невязки времён пробега, сформировать новую матрицу и вектор данных, решить систему линейных уравнений и получить следующее приближение. Процедура повторяется до тех пор, пока мера различия между наблюдёнными и рассчитанными временами пробега не станет меньше некоторой заданной величины. Достоинством матричного обращения является отсутствие ограничений на геометрию расстановки источников и приёмников и на форму лучевых траекторий. Недостатком метода, не позволявшим широко применять его в 80-х и начале 90-х годов прошлого века являются высокие затраты машинной памяти на хранение томографических матриц и время вычислений для трассировки лучей и решения соответствующих систем линейных уравнений. На современных компьютерах, даже на персональных, можно составлять и решать томографические системы достаточно больших размерностей.

Методы преобразования Фурье основаны на теореме о сечении спектра, которая утверждает, что одномерное преобразование Фурье проекции на направление, задаваемое углом относительно фиксированной оси некоторой функции двух переменных, совпадает с сечением двумерного преобразования Фурье этой функции, направленной под тем же углом плоскостью [Mersereau, Oppenheim, 1974]. Проекцией функции двух переменных на направление называется функция одной переменной, получающаяся в результате интегрирования исходной функции по прямым с одинаковым наклоном относительно фиксированной оси. В случае прямолинейных лучевых траекторий и полного набора данных, когда лучи пересекают область под всевозможными углами, на основе изложенной теоремы можно создать очень быстрый метод обращения. Вначале вычисляются преобразования Фурье всех проекций на различные направления и по ним строится двумерный спектр. Затем с помощью обратного двумерного преобразования Фурье находится исходное поле скоростей. На данном этапе приходится применять интерполяцию двумерного спектра при переходе от полярных к декартовым координатам. На практике более эффективным оказывается эквивалентный свёрточный подход, или восстановление по отфильтрованным проекциям, основанный на обратном преобразовании Радона [Хаттон, 1989; Наттерер, 1990]. Достоинствами данной группы методов являются быстрота и строгая математическая обоснованность. Недостатками методов, основанных на преобразовании Фурье, являются предположение о прямолинейности лучевых траекторий и требование хорошей освещённости области лучами. Эти ограничения делают методы мало применимыми в задачах сейсморазведки на проходящих волнах.

Алгебраические процедуры восстановления имеют дело с таким же матричным уравнением, как и метод матричного представления. Единственное отличие состоит в том, что решение системы уравнений выполняется специальным методом итераций. В пределах целевой области, изображение которой требуется построить, рассчитывается ход лучей и определяются времена пробега. Невязки между наблюдёнными и расчётными значениями времён пробега перераспределяются вдоль каждой лучевой траектории, отсюда получается новое приближение. Достоинствами метода являются быстрота решения линеаризованной задачи на каждой итерации и гибкость в отношении геометрии наблюдений и допустимой формы лучевых траекторий. Известно большое число модификаций изложенного алгоритма [Herman et al., 1976; Gordon, 1974; Scudder, 1978; Dines, Lytle, 1979]. Данный метод широко применялся в 80-х годах прошлого века в связи со слабыми мощностями вычислительной техники того времени. В некоторых современных работах этот метод также рассматривается. Процедура решения линеаризованной задачи с помощью размазывания возмущений вдоль лучей даёт грубое приближение решения. Восстановить достаточно сложные структуры этим методом не удаётся.

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

Задача томографии по первым вступлениям проходящих волн решается, как правило, в рамках модели геометрической сейсмики. Процесс распространения волн рассматривается как распространение слабых разрывов в решении уравнений теории упругости. При этом слабые разрывы могут существовать только на характеристиках исходных уравнений. Характеристическое уравнение для уравнений теории упругости, записанное для времён пробега, называется уравнением эйконала. Время прихода волны из источника в приёмник зависит лишь от свойств среды вдоль фиксированной криволинейной траектории, называемой лучом. Лучи являются характеристиками уравнения эйконала и бихарактеристиками уравнений теории упругости. Сейсмические волны, возникающие в реальных физических процессах, вообще говоря, не удовлетворяют указанному выше математическому определению волны. Законы геометрической сейсмики подходят для описания реальных физических волн в случае высоких частот. Чем меньше длина волны по сравнению с характерным размером неоднородности среды - тем точнее процесс распространения волн может быть описан с помощью лучей. Глубокое, математически обоснованное понимание связи геометрической сейсмики с реальными волновыми процессами пришло вместе с развитием лучевого метода В.М.Бабичем, А.С.Алексеевым, Б.Я.Гельчинским [Алексеев и др., 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] рассматривается модификация рассматриваемого конечно-разностного метода, позволяющая рассчитывать поле времён во все стороны от источника, последовательно решая уравнение в разных направлениях. К сожалению, для сред со сложными высококонтрастными включениями (соляные тела) этот метод всё же не применим: лучи слишком много раз меняют направление.

Рассмотрим методы дискретизации пространства моделей и построения томографической системы линейных алгебраических уравнений. Исходный томографический оператор действует в бесконечномерном пространстве скоростей. При решении задачи томографии матричным методом рассматриваются конечномерные пространства моделей и данных. Существуют различные методы дискретизации целевой области. Традиционным и наиболее распространённым методом является использование прямоугольной сетки: среда считается постоянной внутри ячейки, интегралы вдоль лучей аппроксимируются соответствующими суммами.

Возможны и другие способы дискретизации. В работе Michelena, Harris [1991] предложено разбиение области с помощью ячеек, растянутых вдоль лучевых траекторий. Такой способ дискретизации был назван natural pixels. В работе Vasco et al. [1995] этот подход развивается с использованием объёмов Френеля, при этом учитывается физическая природа волны: зоны Френеля зависят от длины волны. Преимуществом подхода является его быстрота. В работе Michelena, Harris [1991] на модельном примере показано, что для восстановления модельного объекта с одинаковым качеством подход, использующий natural pixels, требует существенно меньше времени, чем стандартный, основанный на использовании прямоугольной сетки. Дискретизация с помощью natural pixels, или объёмов Френеля, не подходит для анизотропных сред: это связано с тем, что при её использовании необходимо учитывать влияние ячеек друг на друга в зоне пересечений лучей. В изотропной среде это влияние легко учитывается. Использование прямоугольных сеток является более универсальным и простым при программной реализации алгоритмов.

В работе T.Bodin, M. Sambridge [2009] предложен алгоритм решения задачи томографии, на каждой итерации которого модифицируется не только скоростная модель, но и сетка дискретизации. При этом используется нерегулярная сетка на основе многоугольников Вороного. Такое разбиение определяется только количеством ячеек и их «центральными» точками. Этот подход является оригинальным и интересным, но очень сложен в реализации и требует большого количества вычислений на каждом шаге итераций.

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

В настоящей работе будет рассматриваться стандартная дискретизация целевой области на основе прямоугольной сетки.

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

Можно ли применять подобную аппроксимацию, насколько она точна и как связаны между собой томографическая матрица и исходный оператор, какой должна быть дискретизация? Эти вопросы в большинстве работ исследуются методом обращения синтетических данных для различных моделей. Понятно, что такой подход не является достаточно строгим и универсальным. Существует ряд теоретических результатов, связывающих дискретизацию и разрешение [Наттерер, 1990], вывод которых основан на преобразовании Фурье. Эти теоремы и оценки носят общий характер, они не учитывают компактности томографического оператора, не предполагают наличие ошибок в данных и ошибок аппроксимации оператора.

Оператор томографии является компактным, поэтому обратный оператор не является ограниченным. Получающаяся в результате дискретизации томографическая система наследует свойства исходного бесконечномерного оператора. Томографическая матрица плохо обусловлена. Известно, что оптимальным выбором, позволяющим получать устойчивые решения плохо обусловленной задачи, является базис в пространстве решений, образованный правыми сингулярными векторами исходного линейного оператора [Чеверда, Костин, 2010]. Устойчиво определяются проекции на подпространства, соответствующие старшим правым сингулярным векторам. В научной литературе эти проекции называются r-решениями [Годунов и др., 1992; Чеверда, Костин, 2010].

Существует множество программных пакетов, позволяющих проводить различные операции линейной алгебры, в частности находить сингулярное разложение матриц. Общепринятым стандартом является библиотека LAPACK. Эта библиотека находится в свободном доступе (netlib.org) и написана на фортране. Она входит в широко распространённые математические пакеты, например, в MATLAB и математическую библиотеку компании INTEL MKL (Math Kernel Library). Существуют различные модификации LAPACK: SCALAPACK-LAPACK, переписанный на языке C, PROPACK - параллельная версия и другие. LAPACK позволяет вычислять сингулярное разложение для комплексных матриц двойной точности с использованием процедуры ZGESVD, комплексных матриц одинарной точности CGESVD, вещественных матриц двойной точности DGESVD и вещественных матриц одинарной точности с помощью процедуры SGESVD. В MATLAB этим процедурам соответствует команда svd.

Во все перечисленные процедуры исходные матрицы необходимо передавать в "полном" виде (двумерные массивы соответствующих типов). Процедуры вычисляют сразу весь спектр. Указанные особенности алгоритмов накладывают существенные ограничение на размеры матриц и, следовательно, на размеры сетки дискретизации и количество пар источников-приёмников. Вместе с исходной матрицей необходимо также хранить вспомогательные массивы, которые служат «рабочим пространством» для процедуры SGESVD, а также получаемые сингулярные векторы и числа. Кроме этого, вычисление сингулярного разложения с применением процедур LAPACK требует большого количества времени, так как вычисляется сразу весь спектр.

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

LAPACK не позволяется вычислять сингулярное разложение для спарсовых матриц. Нет такой возможности и в MKL. Для вычисления сингулярного спектра спарсовых матриц применяется пакет ARPACK. Пакет ARPACK (http://www.caam.rice.edu/software/ARPACK) предназначен для вычисления части собственных чисел и векторов квадратных матриц. В пакете есть различные процедуры для вычисления собственных чисел и векторов соответственно для симметричных и несимметричных комплексных и вещественных матриц разной точности. Доступна и параллельная версия ARPAСK. Задача поиска части сингулярного разложения матрицы B может быть сведена к поиску собственных чисел симметричной квадратной матрицы BB*. В пакете ARPACK есть готовые подпрограммы ssvd, dsvd для вычисления первых сингулярных чисел и векторов соответствующих матриц ("s","d" вещественные одинарной и двойной точности). Пакет ARPACK используется в MATLAB для вычисления части сингулярного разложения спарсовых матриц (команда ssvd). Существенным недостатком MATLAB является невозможность использования спарсовых матриц одинарной точности. Если работать с библиотекой ARPACK, непосредственно требуется задавать лишь действие матрицы на вектор и действие сопряжённой матрицы в виде процедур. Таким образом, хранить матрицу вообще не требуется. Тем не менее в ходе выполнения сингулярного разложения действие матрицы и сопряжённой матрицы считается большое число раз, поэтому на практике всё же приходится вычислить элементы матрицы один раз, а затем хранить их в спарсовом формате для действия прямой и сопряжённых матриц.

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

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

При решении обратных задач геофизики на практике используются различные процедуры сглаживания решения за счёт введения специальных сглаживающих матриц, которые выступают в роли предобуславливателей для исходной системы линейных уравнений [Fomel, 2007]. Достоинством подхода является возможность учитывать априорную информацию о восстанавливаемой скоростной модели. Недостатком является неопределённость в выборе сглаживания. В общем случае, при отсутствии априорной информации неясно, какими должны быть параметры сглаживания.

Использование многомасштабных базисов [Woodward et al., 2008] также является одним из вариантов регуляризации. Решение задачи отыскивается итерационно с использованием последовательности измельчающихся прямоугольных сеток дискретизации целевой области. Решение, полученное на крупной сетке, используется в качестве начального приближения для решения задачи на более мелкой. При переходе с крупной сетки на мелкую используется сглаживание. Этот метод регуляризации по свойствам схож с предыдущим.

Универсальным и математически обоснованным методом является регуляризация на основе усечения сингулярного разложения. При решении линейных алгебраических уравнений устойчиво могут быть найдены r-решения, являющиеся проекциями точного решения на пространства правых сингулярных векторов. Существуют оценки, связывающие точность данных и задания линейного оператора с точностью получаемого решения [Чеверда, Костин 2010]. R-решения для старших сингулярных векторов содержат информацию о пространственном спектре решения для низких частот. Другими словами, регуляризация на основе сингулярного разложения эквивалентна сглаживанию. Преимущество в том, что единственным параметром является количество сингулярных векторов, которое определяется ошибкой и соответствующим числом обусловленности.

Существует ряд способов описания анизотропии упругих сред, наиболее распространённые из которых параметризации Томсена [Thomsen, 1986] и Шонберга [Schoenberg, de Hoop, 2000]. Предварительным этапом разработки надёжного численного метода восстановления строения среды с учётом анизотропных свойств является определение оптимальной параметризации. Надёжным способом выбора оптимальной параметризации является изучение проекций на пространства старших правых сингулярных векторов [Протасов, Сердюков, Чеверда 2010].

Изложенные подходы регуляризации решения систем линейных уравнений предполагают поиск решения методом наименьших квадратов. Решение определяется как точка минимума функционала среднеквадратичного уклонения действия линейного оператора на вектор моделей и данных. Такой подход устойчиво работает в условиях Гауссовского распределения ошибок в данных. Недостатком метода наименьших квадратов является его повышенная чувствительность к отдельным большим ошибкам в данных. Известно, что решение оказывается устойчивым к таким ошибкам, если в качестве критерия использовать не наименьшие квадраты, а минимизировать абсолютные отклонения [Claerbout, Muir, 1973].

Примеры расшифровки внутренних структур Земли сейсмотомографическим методом приведены на рисунках 1-3.




  1. Численные исследования на модели Северного моря:

Рисунок 1 - Слева – скоростная модель среды, справа – результат обращения синтетических данных (сейсмотомографическое изображение среды).



  1. Примеры сейсмотомографической расшифровки структуры верхней мантии Земли:


Рисунок 2 – Сейсмотомографическое изображение верхней мантии Земли



Поделитесь с Вашими друзьями:
  1   2


База данных защищена авторским правом ©grazit.ru 2017
обратиться к администрации

    Главная страница