Содержание к диссертации
Введение
Глава 1 Локальная сейсмическая томография 11
1.1 Обзор исследований в области обратных кинематических задач сейсмики 11
2.1 Определение глубинного строения но данным от близких землетря сений 12
2.1.1 Отбор и предварительная обработка входных данных . 14
2.1.2 Построение начальной референтной скоростной модели . 15
2.1.3 Параметризация модели 15
2.1.4 Решение прямой кинематической задачи 16
2.1.5 Разделение системы на дискретную и непрерывную части . 17
2.1.6 Обращение систем линейных уравнений, регуляризация . 18
2.1.7 Оценка разрешающей способности метода и ошибок в решении 20
Глава 2 Определение скоростной структуры среды с известными и неизвестными источниками 22
1.2 Численное решение задачи линейной сейсмической томографии . 23
1.2.1 Исследование прямого томографического оператора в двухмерном случае 23
1.2.2 Компактность линейного томографического оператора . 26
1.2.3 Выбор численного алгоритма для решения задачи линейной сейсмической томографии 27
1.2.4 Определение трехмерной скоростной структуры среды методом сейсмической томографии без трассировки лучей . 29
1.2.5 Метод LSQRA, получение оценок матриц разрешающей способности, информационной плотности и ковариации 39
2.2 Смешанная задача: определение скоростной модели среды с одно временным уточнением параметров гипоцентров землетрясений . 41
2.2.1 Постановка задачи 41
2.2.2 Линеаризация системы 42
2.2.3 Разделение системы 43
Глава 3 Обращение синтетических и реальных данных 45
1.3 Обработка данных физического моделирования 45
1.3.1 Физическое моделирование 45
1.3.2 Аппроксимация линейного томографического оператора. Его сингулярный спектр 50
1.3.3 Результаты обработки данных физического моделирования . 51
2.3 Сейсмическая томография эпицентралыюй зоны Чуйского землетрясения 57
2.3.1 Численные эксперименты 57
2.3.2 Обращение данных афтершоков Чуйского землетрясения . 68
Заключение 87
Литература
- Определение глубинного строения но данным от близких землетря сений
- Решение прямой кинематической задачи
- Исследование прямого томографического оператора в двухмерном случае
- Аппроксимация линейного томографического оператора. Его сингулярный спектр
Введение к работе
Объект исследования — решение обратных кинематических задач сейсмики для случаев известных и неизвестных источников сейсмических волн, на предмет определения скоростного строения среды.
Актуальность темы
Задачи определения скоростного строения среды по кинематическим данным имеют широкую область приложения в науках о Земле: начиная от применения томографических методов в инженерной сейсморазведке до построения глобальных томографических моделей земной коры и мантии по тслесейсмическим данным. Несмотря на бурное развитие этой области, произошедшее за последние несколько десятилетий, вопросы об оценке надежности получаемых решений и об уменьшении объемов вычислений при обращении данных остаются острыми и по сей день. Особенно актуальными эти вопросы являются в задачах локальной сейсмо-томографии, связанной с одновременным определением скоростной структуры и параметров источников сейсмических волн. Неизвестные положения гипоцентров землетрясений ведут к большей неопределенности в скоростных моделях по сравнению с задачами, где источники известны, а обращение больших объемов данных требует применения мощных вычислительных систем (затрат оперативной памяти и процессорного времени). Оценка точности параметров получаемых скоростных моделей традиционно сопряжена с многократной инверсией данных или зачастую с непомерно дорогими процедурами полного обращения возникающей при этом матрицы, поэтому остается слабым местом большинства подходов.
В работе предлагается один из возможных способов инверсии больших объемов данных, основанный на применении регуляризующих итерационных процедур решения больших систем линейных уравнений с несимметричными матрицами. Точность скоростных моделей оценивается с помощью матриц ковариации и разрешающей способности, определяемых попутно с обращением данных, без привлечения дополнительных вычислительных ресурсов.
Цель исследований — разработка нового томографического алгоритма определения скоростного строения среды по кинематическим данным, не зависящего от используемых скоростных моделей, а также повышение достоверности решений путем увеличения детальности за счет обращения больших объемов данных с помощью итерационных алгоритмов решения систем линейных уравнений, и мак-
симально полного использования дополнительной априорной информации о стартовой модели среды и входных данных.
Научная задача исследований — определение скоростного строения среды но временам прихода волн для случаев известных и неизвестных источников. Поставленная задача предполагает выполнение следующих этапов:
Изучить область определения оператора линейной сейсмической томографии и выяснить возможность его матричного представления.
Разработать и программно реализовать универсальный метод определения скоростного строения среды и параметров источников по кинематическим данным, основанный на применении итерационных алгоритмов обращения систем линейных уравнений.
Адаптировать итерационную методику оценки точности скоростных моделей LSQRA к разработанному методу сейсмической томографии без трассировки лучей.
Произнести обращение реальных данных, зарегистрированных как от известных, так и от неизвестных источников сейсмических волн с последующей оценкой точности полученных решений.
Фактический материал и методы исследования
Определение скоростного строения среды и параметров источников, а также оценка точности полученных решений основывались на теории обратных задач и применении регрессионного анализа. На каждом этане разработки метода производилось тестирование его численной реализации. Точность получаемых решений оценивалась с применением синтетических моделей. Кроме того, для верификации разработанного алгоритма использовались данные физического моделирования, позволившие дать дополнительную оценку качества получаемых решений.
С помощью разработанного метода была проведена обработка
данных, полученых в результате физического моделирования;
синтетических данных, позволивших оценить разрешающую способность данных, зарегистрированных па редкой сети;
времен вступлений Р- и S-волн, зарегистрированных временной сетью цифровых станций Актагаского полигона от афтершоков Чуйского землетрясения;
Основной метод исследования — математическое моделирование. Исследование включало и себя: построение разностных схем решения систем дифференциальных уравнений, решение плохо обусловленных систем линейных уравнений с
разреженными несимметричными матрицами с помощью LSQR-алгоритма, определение параметров гипоцентров землетрясений методом Гейгера, применение пакета программ VELEST для определения одномерной оптимальной скоростной модели, нелинейная минимизация методом наименьших квадратов с изменяющимися весами, решение прямой кинематической задачи конечно-разностным способом, вычисление сингулярного разложения матриц и расчет матриц коварпацин, проектирование баз данных сейсмологической информации и разработку пользовательского интерфейса к ним.
Защищаемые научные результаты
Доказана ограниченность и компактность линейного томографического оператора в двухмерном случае, в предположении о регулярности поля лучей.
Разработан и программно реализован способ определения скоростного строения среды и параметров источников, основанный на решениях как систем линейных уравнения итерационным методом, так и дифференциальных уравнений для прямого и сопряженного томографических операторов конечными разностями, позволяющий производить инверсию кинематических данных без применения громоздких процедур полного обращения матриц, не требующий трассировки лучей.
Адаптирована к разработанному алгоритму и программно реализована итерационная методика оценки точности скоростных моделей LSQRA (Zhang, McMcchan, 1995), основанная на вычислении матриц ковариацин и разрешающей способности.
Алгоритм опробован на данных от афтершоков Чуйского землетрясения: выделены аномалии скорости продольных (до 0.15-0.25 км/с) и поперечных (до 0.1 км/с) волн, со средним горизонтальным размером около 30 км при разрешающей способности обусловленной неравномерной плотностью расположения станций не более 0.16.
Новизна работы. Личный вклад
1. Предложен оригинальный алгоритм обращения времен вступления волн и получению оценок точности рассчитываемых скоростных моделей:
- с использованием интегрального представления линейного томографического оператора в непрерывной постановке, показана его ограниченность и компактность, и как следствие, возможность его конечномерной аппроксимации;
в результате применения конечно-разностного метода решения уравнения эйконала, удалось избежать трассировки лучей при решении обратной кинематической задачи в линеаризованном виде и тем самым обеспечить надежность работы алгоритма в сложных трехмерных средах;
используя тот факт, что на каждой итерации по методу LSQR достаточно знания действия прямого и сопряженного операторов на конкретный вектор, развит подход, не требующий знания матричного представления томографического оператора;
с использованием представления граничных узлов сетки в качестве вторичных источников сейсмических волн для конечно-разностного метода решения уравнения эйконала разработан блочный способ решения прямой кинематической задачи;
на основе метода разделения задачи на дискретную и непрерывную части метод сейсмической томографии без трассировки лучей адаптирован к задаче определения скоростного строения с неизвестными источниками;
на основе решения на каждом шаге линеаризованной обратной кинематической задачи и введения процедуры взвешивания данных по определенному критерию реализован нелинейный итерационный метод наименьших квадратов с изменяющимися весами;
с применением гауссовского сглаживания разработатапы регуляризую-щис процедуры, максимально использующие априорную информацию;
с использованием модификации метода LSQR реализован алгоритм расчета матриц ковариации и разрешающей способности.
Разработанный подход реализован в виде комплекса прикладных программ и библиотек.
С применением разработанного программного обеспечения получены трехмерные скоростные модели Р- и S-волн эпицентралыюй зоны Чуйского землетрясения, а также впервые проведена оценка разрешающей способности данных, зарегистрированных на редкой сети.
Теоретическая и практическая значимость результатов
Предложенный томографический алгоритм является инструментом, позволяющим определять глубинное строение по кинематическим данным в сложных трехмерных средах применительно к широкому диапазону задач геофизики: от построения инженерных сейсмических разрезов до получения трехмерных скоростных моделей по сейсмологическим данным. Применение итерационного метода LSQR
позволяет производить обработку больших объемов данных без значительных затрат оперативной памяти и ресурсов процессора. Адаптированная итерационная методика вычисления матриц ковариацин и разрешающей способности делает возможной оценку точности трехмерных скоростных моделей для задач сейсмической томографии с большим количеством наблюдений и неизвестных параметров.
Разработанное программное обеспечение представляет из себя набор классов объектно-ориентированного языка программирования C++, реализующих основные понятия томографического подхода: распределение скорости и среде, поля времен, годографы, системы наблюдений, процедуры сглаживания модели и взвешивания данных, процедуры инверсии, и может быть собрано практически на любой платформе как в виде динамических библиотек, так и в виде исполняемых программ, предоставляя конечному пользователю широкий набор инструментов.
Найденные трехмерные скоростные модели и параметры очагов землетрясений являются первым шагом к построению детальной скоростной модели очаговой зоны Чуйского землетрясения, которая могла бы пролить свет на развитие сейсмичности в афтершоковой зоне, а проведенный анализ разрешающей способности данных, зарегистрированных на редкой сети временных станций Акташского полигона АСФ ГС СО РАН, определяет доверительные интервалы для получаемых скоростных моделей и параметров гипоцентров землетрясений.
Апробация работы
Основные результаты работы докладывались на различных российских и международных конференциях, наиболее значимыми из которых являются: Генеральная ассамблея IUGG (Саппоро, Япония, 2003 г.), Международная конференция "Проблемы сейсмологии Ш-ого тысячелетня" (Новосибирск, 2003 г.), Международная конференция "Математические методы в геофизике" (Новосибирск, 2003 г.), Международная научная конференция, посвященная 90-летию академика Н.Н. Пузырена (Новосибирск, 2004 г.), а также на геофизическом семинаре ИГФ ОИГГМ СО РАН, 29 сентября 2004 г.
Объем и структура работы
Данная диссертационная работа состоит из трех глав, заключения и приложения.
В первой главе проводится анализ работ, посвященных задачам одновременного определения скоростной структуры и параметров гипоцентров землетрясений -задачам локальной сейсмической томографии. Материал излагается в виде последовательности основных этапов, выполняемых при решении задачи, с подведением кратких итогов в конце.
Определение глубинного строения но данным от близких землетря сений
Решая задачу об определении глубинного строения земной коры по данным близких землетрясений мы попадаем в заколдованный круг: с одной стороны, для того чтобы определить скоростное строение некоторой области, нам необходимы достаточно точные координаты гипоцентров землетрясений; с другой стороны, для точного определения гипоцентров землетрясений, нам необходимо хорошо знать скоростное глубинное строение Земли в данной области. Здравый смысл подсказывает, что нам следует зафиксировать" источники сейсмических волн или скоростную модель. Однако на практике, зафиксировать, например, источники (используя промышленные взрывы для определения скоростной структуры) не удается, в связи с тем, что их слишком мало и к тому же сеть взрывов как правило, не покрывает исследуемую область. При этом у нас может не быть априорной информации о скоростной модели среды.
Тем не менее, времена пробега волн от значительного числа землетрясений содержат достаточно информации как для определения трехмерного скоростного строения Земли, так и для определения параметров гипоцентров. Эта задача является нелинейной, и обычно решается итерационным методом Ньютона с линеаризацией на каждом шаге. Впервые, решение линейной задачи было предложено Павлисом и Букером в 1980 году, в работе [63], где задачи подобного типа были названы смешанными, т.к. они состоят из двух типов неизвестных параметров: параметры гипоцентров - дискретная часть модели; распределение скоростей - непрерывная часть модели.
Дискретная часть модели характеризуется конечным числом неизвестных параметров. В нашем случае, у каждого гипоцентра их четыре: (xo,yo,zo,to); это пространственные координаты и время в очаге. Непрерывная часть модели содержит бесконечное число неизвестных параметров, поэтому задача об их нахождении по конечному числу наблюдений недоопределена в принципе. Павлис и Букер впервые показали, что для смешанной линейной обратной задачи, когда число неизвестных дискретных параметров модели р меньше числа наблюдений d, можно подобрать такое преобразование исходной системы линейных уравнений, что по крайней мере d — р полученных уравнений не будут зависеть от дискретных параметров и могут быть использованы для определения скоростной модели среды. Иначе говоря, мы мооїсем рассматривать задачи об определении гипоцентров и скоростной модели среды независимо друг от друга. Таким образом, существующие подходы к решению смешанной задачи можно классифицировать по методам определения скоростной структуры и координат гипоцентров, а также по способам оценки разрешающей способности и погрешности решений.
Все существующие па данный момент методы определения скоростной структуры по данным от близких землетрясений, имеют несколько обязательных этапов, которые необходимо выполнить для получения решения: Отбор и предварительная обработка входных данных Построение начальной референтной скоростной модели Параметризация модели Решение прямой кинематической задачи Разделение системы па дискретную и непрерывную части Обращение систем линейных уравнений Оценка разрешающей способности метода и ошибок в решении Рассмотрим каждый из этих этапов более подробно.
Отбор и предварительная обработка входных данных
Для построения достоверной скоростной модели необходимо использовать только качественные данные. Иными словами, из каталога должны выбираться землетрясения, удовлетворяющие следующим условиям: Глубина очага землетрясения не должна превышать hmax Невязки времен пробега не должны превышать сгтах Число зарегистрированных па одной станции вступлений должно быть больше і тіп Эпицентр землетрясения должен находиться в пределах сети станций. Т.е. если соединить эпицентр с каждой станцией сети отрезком, максимальный угол фдар среди всех углов между соседними отрезками не должен превышать Фтах (РИС 1.1) Ближайшая к эпицентру станция должна быть на расстоянии, не большем чем Атах Магпитуда землетрясения должна быть больше Mmt„ Ошибка определения горизонтальных координат эпицентра землетрясения не должна превышать аг Ошибка определения глубины очага землетрясения не должна превышать ст/,
Заданные в этих условиях величины являются управляющими параметрами, от правильного выбора которых зависит качество решения.
Помимо землетрясений, использование промышленных взрывов, с известными координатами и временами в очагах позволяет откалибровать скоростную модель.
Большинство моделей, расчитанных по временам пробега от близких землетрясений, используют времена, зарегистрированные на плотных сетях сейсмологических станций, со средними удалениями Д = 50 — 100 км. При этом наблюдения ведутся в сейсмоактивных зонах, что позволяет получать значительное число наблюдений (порядка 5000-50000 времен пробега) при относительно небольшом количестве станций (порядка 50). Для получения устойчивого гипоцентрального решения, необходимо чтобы землетрясение находилось в пределах сети станций, т.е. чтобы угол фдаР не превышал 180 [58], или даже 135 [62].
Решение прямой кинематической задачи
Использование двух типов неизвестных параметров в системе линейных уравнений существенно отличает задачу определения скоростной структуры по временам пробега волн близких землетрясений от других задач сейсмической томографии, где известны положения источников воли. Впервые времена пробега от близких землетрясений были использованы для определения скоростной структуры в [36], затем, в [63] метод разделения неизвестных скоростных параметров и параметров гипоцентров был формализован, и используется до сих пор практически в первоначальном виде. Для смешанной системы линейных уравнений
Aa; + Gv = d, где A - матрица, связанная с параметрами гипоцентров, подбирается такое преобразование U, чтобы выполнялось равенство UA = 0. Долученая в итоге система будет содержать только неизвестные скоростные параметры v. Определив поправки к скоростным параметрам, мы можем уточнить положения гипоцентров, на следующей итерации в свою очередь уточнить скоростные параметры. Итерационный процесс продолжается до тех пор, пока дисперсия невязок времен пробега не станет меньше некоторой пороговой величины.
Обращение систем линейных уравнений, регуляризация Система линейных уравнений Ах = Ь, (1.1) свазанная с определением поправок к скоростной модели среды по невязкам времен пробега волн обладает следующими особенностями: [31]: она несимметричная; разреженная, так как только сравнительно небольшая часть элементов матрицы А отличны от нуля; прямоугольная, число строк значительно превышает число столбцов; при этом ранг матрицы А меньше числа столбцов; несовместная, т.е. не существует вектора х, точно удовлетворяющего системе (1.1).
Поскольку (1.1) не имеет точного решения, прибегают к решению в смысле наименьших квадратов, минимизируя \\Ах — Ь2. Если гапк(А) п, существует бесконечно много векторов, минимизирующих \\Ах — Ь2. В этом случае, среди всех решений (1.1) ищут минимальное по норме решение (обобщенное обратное решение).
В классическом учебнике К. Аки и П. Ричардса "Количественая сейсмология" [1], авторами было дано определение понятий обобщенной обратной матрицы и обобщенного нормального решения с использованием понятий SVD-разложения матриц (SVD-разложение и обобщенное обратное решение уравнения (1.1) подробно описаны в Приложении А).
Несмотря на то, что SVD-разложснис даст минимальное по норме решение и позволяет оценить его погрешность, оно практически никогда не используется для обращения реальных данных. Причина в том, что SVD-разложение большой, разреженной матрицы является очень дорогой процедурой в смысле компьютерного времени и памяти.
Для получения обобщенного обратного решения, необходим метод способный решать системы линейных алгебраических уравнений с нсссимст-ричной матрицей не требующий храпения всей матрицы А в оперативной памяти компьютера по возможности являющийся регуляризующим
Всем этим условиям удовлетворяет метод LSQR Пейджа-Сауидерса [74], являющийся модификацией метода сопряженных градиентов для несимметричных матриц, который получил широкое распространение в последнее время. На каждой итерации LSQR требуется вычислить действие матриц на векторы: у = Ах и х = Ату. Таким образом, нет необходимости хранить всю матрицу А в оперативной памяти компьютера. Более того, не обязательно вычислять саму матрицу А, достаточно вычислить действие прямого и сопряженного линейных операторов на соответствующие векторы. Регуляризующим параметром для LSQR является число итераций: но мере роста числа итераций увеличивается разрешенпость, но падает устойчивость решения (краткое описание и последовательность шагов LSQR приведены в 2.1.5 и Приложении В).
Описанные методы получения решений в смысле наименьших квадратов работают без дополнительных ограничений только для "идеальных" задач: обращение результатов синтетических экспериментов, обращение данных, неотягощенных ошибками. На практике же, регуляризующих свойств LSQR и SVD оказывается недостаточно. Кроме того, решение в смысле наименьших квадратов очень чувствительно к так называемым "отскокам" - небольшим порциям данных, значительно отличающимся от данных, предсказанных моделью [72], [26]. Для получения устойчивого решения и уменьшения влияния отскоков, необходимо привлечение дополнительной априорной информации. В связи с этим, вместо минимизации Лх — Ь2, ищется минимум следующей функции [72]:
Априорная матрица ковариации модели Сд/ содержит информацию о модельных параметрах, которая используется для регуляризации. Самым простым способом регуляризации будет задание одинаковых дисперсий для всех модельных параметров: Сд/ = 0д/1, что ведет к уменьшению нормы решения х в зависимости от значения дисперсии o\t. Чуть более сложный метод регуляризации, предложенный Нолетом [26] состоит во взвешивании модельных параметров суммами длин лучей, проходящих через блоки. Однако, самое широкое распространение получило использование сглаживающих операторов в качестве Сд/. Например, если принять ограничение, что модель характеризуется изотропным радиусом корреляции иеоднородностей L, то можно положить [26]:СД/(ХІ,:Г,-) = a2exp(—l/2\xi — Xj\2/L). В работах [41], [62], [58], [70] для сглаживания к модели применяется оператор, обратный лапласиану. Причем, частные производные вычисляются только по горизонтальным пространственным координатам.
Ковариационная матрица данных Сд(г), с элементами, зависящими от величин невязок соответствующих наблюдений г, призвана уменьшить влияние больших невязок (в идеале - обеспечить одинаковую дисперсию ошибок в данных), связанных с недостаточно точно заданной стартовой скоростной моделью среды или со случайными ошибками в данных. Обычно матрица С/ (г) считается диагональной, что отражает предположение об отсутствии корреляции в данных [45], [54]. Таким образом, на каждой итерации процесса (2.34), с помощью матрицы Сд(г) мы можем уменьшить (в том числе до нуля) влияние данных, отличающихся от предсказанных моделью. При этом, на следующей итерации эти данные могут быть приняты обратно. Этот метод получил название итерационного метода наименьших квадратов с изменяющимися весами [26].
Исследование прямого томографического оператора в двухмерном случае
Глава посвящена задаче определения скоростного строения среды по кинематическим данным для случаев известных и неизвестных источников сейсмических волн. В первой части рассматривается задача сейсмической томографии в линеаризованном виде. Описываются свойства линейного томографического оператора для случая двух пространственных переменных; показывается возможность его конечномерной аппроксимации. Однако, основное внимание в данной главе уделено изложению метода сейсмической томографии без трассировки лучей, основанного на применении итерационного метода LSQR: описывается способ вычисления прямого и сопряженного томографических операторов, а затем способ получения их конечномерной аппроксимации. Приводится адаптированный к разработанному подходу способ оценки точности скоростных моделей, основанный на итерационном методе расчета матриц разрешающей способности и ковариации [81]. Во второй части главы изложенный подход к определению скоростного строения среды обобщается на случай неизвестных источников: рассматривается задача одновременного определения скоростного строения среды и координат гипоцентров землетрясений. Следуя [63], эта задача называется смешанной, и решается путем разделения системы линейных уравнений на дискретную и непрерывную части.
Постановка задачи и основные свойства оператора Будем считать, что полоса 0 х Н заполнена средой, скорость распространения волн в которой представляется в виде: с(х, z) = co(x,z) + ci(x,z) где Со(х, z) априорно заданная функция, а с\(х, z) - добавка, которая должна быть определена по известным временам пробега между источниками, расположенными на отрезке [О, Н\ прямой х = 0 и приемниками, расположенными на таком же отрезке прямой х = L (Рис. 2.1). В дальнейшем вместо времен пробега мы будем рассматривать невязки между полным временем пробега и временем пробега для среды с априорно заданной скоростью. В предположении малости добавки, то есть в предположении CI(X,Z)/CQ(X,Z) « 1, получим следующее соотношение, связывающее невязку времени пробега между источником с координатами (0, s) и приемником с координатами (L,r):
Таким образом, обращение невязок времен пробега в рамках линейной сейсмической томографии сводится к необходимости решения линейного операторного уравнения первого рода с компактным оператором. Самым первым шагом для получения численного решения этой задачи является ее дискретизация, в самом общем виде сводящаяся к представлению искомого решения в виде конечной линейной комбинации некоторых базисных функций. Компактность оператора обеспечивает существование его матричного представления и позволяет использовать проекционный метод для отыскания численного решения уравнения (2.2). Напомним, что проекционный метод состоит в представлении решения в виде линейной комбинации базисных векторов {фі,..., фм} из областип определения оператора, в то время как правая часть заменяется своей проекцией па подпространство, являющееся линейной обололчкой базисных векторов { і,... ,фы} из образа опреато-ра. Выбирая эти базисы ортонормированными, получим матричное представление операторного уравнения:
Компактность оператора обеспечивает сходимость решения матричного уравнения к решению исходного операторного уравнения при М, N —+ со если последнее существует. Основная трудность па этом пути состоит в том, что чем точнее приближается оператор, то есть чем больше размерности подпространств iV, М, тем хуже обусловленность матрицы (ноль принадлежит сингулярному спектру компактного оператора!) и, следовательно, тем более неустойчивым, чувствительным к погрешностям, является получаемое решение. Одним из способов преодоления этого затруднения является описанный в работах [22], [55] метод построения г-решения, сводящийся к использованию усеченного SVD-разложения матрицы. Однако, при работе с реалистичными размерностями процедура отыскания SVD-разложения становится чрезмерно дорогой и не может применяться для построения численного решения. Здесь необходимо отметить, что дискретизация уже сама по себе является регуляризирующей процедурой и, вообще говоря, идеальный выбор базиса - правые сингулярные векторы оператора - решает все проблемы, связанные с построением численного решения. Конечно, такой выбора базиса за редчайшим исключением, невозможен, так как задача отыскания сингулярного спектра оператора гораздо сложнее задачи построения конкретного решения. Поэтому переход к дискретной постановке обязательно должен сопровождаться и некоторой регуляризующей процедурой. Кроме того, как уже было отмечено выше, спецификой рассматриваемой задачи является потенциально чрезвычайно большая ее размерность, поэтому численный алгорити решения должен быть итерационным и не требовать хранения в оперативной памяти всей матрицы. Для симметричных матриц идеально отвечающим поставленным требованиям является метод сопряженных градиентов ([48]) - он является итерационным, требует только вычисления действия матрицы на вектор и не требует храпения всей матрицы и, что имеет принципиальное значение, является регуляризующим алгоритмом с параметром регуляризации "число итераций". Для систем уравнений с несимметричной матрицей в работе [74] предложена и обоснована модификация метода сопряженных градиентов, получившая название LSQR. Как показано в [81] этот метод обладает всеми теми же свойствами, что и метод сопряженных градиентов, и, прежде все го, является регуляризущим по числу итераций. К сожалению, нам не известны результаты по монотонности этого параметра регуляризации для LSQR, однако все проведенные нами численные эксперименты се подтверждают - если изображение после некоторой итерации начало разрушаться, то это разрушение шаг от шага только усиливается. Поэтому критерием остановки итерационного процесса является визуальное сравнение результатов двух последовательных итераций - если отмечается разрушение изображения, за окончательное решение берется результат, полученный па предыдущей итерации.
Согласно выводам, сделаным в предыдущем разделе, итерационный метод LSQR Пейджа-Саундерса [74] наиболее удачно подходит для численного решения линейной задачи сейсмической томографии. LSQR позволяет находить решение систем линейных уравнений н смысле наименьших квадратов, при этом па каждой итерации алгоритма требуется вычислить действие прямого и сопряженного операторов с последующей нормировкой (последовательность шагов LSQR изложена в Приложении В). Таким образом, для получения решения методом LSQR нет необходимости даже определять матричное представление оператора А. Нужно научиться вычислять действие прямого томографического оператора А на элемент пространства моделей П\ и действие сопряженного к нему оператора А на элемент пространства данных Т\. Данный раздел посвящен получению выражений для прямого и сопряженного операторов в случае трех пространственных переменных.
Вычисление прямого и сопряженного операторов в случае непрерывных функций
Вычисление прямого оператора Пусть D - некоторая область, имеющая форму параллелепипеда с размерами Нх, Ну и Hz. Введем в нашу область систему декартовых координат, таким образом, чтобы оси совпадали с ребрами параллелепипеда D (Рис. 2.2).
Аппроксимация линейного томографического оператора. Его сингулярный спектр
Аппроксимация линейного томографического оператора производилась стандартным образом - целевая область покрывалась регулярной прямоугольной сеткой, в каждой ячейке которой искомая скорость предполагалась постоянной. Пространство входных данных, полученных і! ходе проведения эксперимента, само по себе обязательно является дискретным, так как состоит из времен пробега, измеренных для конечного числа источников в конечном числе приемников. Количество источников, умноженное на количество приемников дает размерность пространства данных. Размерность пространства моделей определяется количеством базисных функций, используемых для аппроксимации искомой скорости распространения волн в целевой области. Для аппроксимации томографического оператора для второго эксперимента, целевая область была покрыта прямоугольной сеткой разбита на 50x95 равных прямоугольников, расположенных таким образом, чтобы источники и приемники располагались бы в серединах соответствующих сторон. Интегральный оператор (2.2) аппроксимировался по методу прямоугольников. В итоге исходное интегральное уравнение свелось к системе линейных алгебраических уравнений с прямоугольной матрицей, число строк которой (число уравнений) равно 45x95=4275 (число источников умножить на число приемников), а число неизвестных (размерность пространства моделей) равно 4750 (50x95). Первые 500 старших сингулярных чисел приведены на Рис. 3.3. О поведении правых сингулярных векторов позволяет судить Рис. 3.4. Как можно заметить, самый старший - первый - правый сингулярный вектор отображает так называемую "освещенность" изучаемой области, то есть плотность пересекающих ее лучей. Этот вектор является знакопостоянным и монотонным на всей целевой области. Приведенный па следующем рисунке сотый правый сингулярный вектор является весьма сильно осциллирующим, но, тем не менее сосредоточенным в весьма "яркой" области, то есть в области, достаточно плотно покрытой лучами, соединяющими источники и приемники.Заметим, что правые сингулярные векторы, принадлежащие "почти" ядру опретора (то есть соответствующие отношениям а(1)/а(п) -г-104 и выше) сосредоточены в основном в слабо освещенной области и являются быстро осциллирующими.
Основное заключение, которое можно сделать из анализа поведения сингулярных чисел и правых сингулярных векторов оператора формулируется следующим образом: при решении задачи линейной сейсмической томографии наиболее устойчиво определяются плавные неоднородности, расположенные в наиболее ярко освещенной части целевой области.
Результаты обработки данных физического моделирования
Па Рис. 3.5 приведены результаты применения метода LSQR для обращения невязок времен пробега для первой модели, имитирующей задачу определения скоростного строения земной коры и верхней мантии по данным о невязках времен пробега от локальных землетрясений. Как можно видеть, уже после двух итераций достаточно уверенно выделяются контуры объекта, размещенные на фоне неоднородной яркости освещения целевой области. Характерной особенностью полученного изображения является его "размытость" вплоть до границ целевой области, что связано с отсутствием лучей, проходящих между верхним (нижним) краем объекта и системой приемников (источников), то есть отсутствием лучей, "отсекающих" объект сверху и снизу. Естетственно, что при выполнении последующих итераций объект остается "размытым" сверху и снизу. Неоднородность яркости постепенно исчезает, но начинает разрушаться изображение объекта. Разрушение накапливается сверху и снизу, постепенно распространяясь к центру. Весьма сажным является и тот факт, что кроме восстановления геометрии объекта, восстанавливается и значение скорости распространения волн внутри него.
Результат применения метода LSQR для обработки данных, полученных для второй модели, представлен на Рис. 3.6. Здесь мы сразу приводим результат после семи итераций, как наиболее информативный. Последующее выполнение итераций приводит к разрушению изображений. Сравнивая полученное изображение с приведенным на Рис. 3.2 строением среды можно отметить, что:
отчетливо прослеживается слой пониженной скорости в интервале глубин от 1140 м до примерно 1200 м, подпираемый клиновидным слоем с повышенной скоростью распространения волн (соответствие вмещающей среде в соответствии с приведенной шкалой оттенков серого достаточно хорошее);
в интервале глубин от 1200 м до 1250 м угадывается клиновидный слой с пониженной скоростью распространения волн, сразу за которым прослеживается слой с повышенной скоростью;
куполообразная структура располагается, как п должно, на глубине порядка 1300 м, а в интрвале от 1500 м до примерно 1350 м можно заметить достаточно четко выраженное поднятие слоев слева направо.
Таким образом, можно констатировать, что в общих чертах структура скоростного строения среды восстановлена правильно. Кроме того, сравнивая распределение оттенков серого на изображении, можно заключить, что и абсолютные значения скоростей находятся верно.
В заключение данного раздела необходимо подчеркнуть, что рассматриваемая модель весьма хорошо соответствует понятию "малое возмущение" однородной вмещающей среды, так как отклонение скорости внутри объекта от скорости во вмещающей среде составляет менее 3%. Тем самым как входные данные, так и сам оператор достаточно хорошо описываются линейным приближением, что позволяет выполнить достаточно много итераций по методу LSQR до начала разрушения изображения объекта.
Похожие диссертации на Определение скоростной структуры среды с известными и неизвестными источниками методом сейсмической томографии без трассировки лучей
