Книги, научные публикации Pages:     | 1 | 2 | -- [ Страница 1 ] --

ИНСТИТУТ ДИНАМИКИ ГЕОСФЕР РОССИЙСКОЙ АКАДЕМИИ НАУК

На правах рукописи

Усольцева Ольга Алексеевна ТРЕХМЕРНЫЕ СКОРОСТНЫЕ МОДЕЛИ ЗЕМНОЙ КОРЫ ТЯНЬ-ШАНЯ НА ОСНОВЕ БИ-СПЛАЙН ПАРАМЕТРИЗАЦИИ И ТРИАНГУЛЯЦИИ

ДЕЛОНЕ Специальность 25.00.10 - геофизика, геофизические методы поисков полезных ископаемых ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

Научный руководитель - доктор физико-математических наук Ирина Альфатовна Санина Москва - 2004 2 Оглавление.

СОКРАЩЕНИЯ, ТЕРМИНЫ И ОБОЗНАЧЕНИЯ.........................................................5 ВВЕДЕНИЕ....................................................................................................................6 ГЛАВА 1 ИСТОРИЯ РАЗВИТИЯ СЕЙСМИЧЕСКОЙ ТОМОГРАФИИ И СУЩЕСТВУЮЩИЕ ТОМОГРАФИЧЕСКИЕ МОДЕЛИ..............................................12 1.1 Предыстория создания первых трехмерных моделей скоростей сейсмических волн....................................................................................................................................... 1.2 Некоторые существующие трехмерные модели скоростей Р и S волн.................. Выводы к главе........................................................................................................................ ГЛАВА 2 МЕТОДИКА ПОСТРОЕНИЯ ТРЕХМЕРНЫХ РЕГИОНАЛЬНЫХ ТОМОГРАФИЧЕСКИХ МОДЕЛЕЙ ПО ДАННЫМ ВРЕМЕН ПРОБЕГА ОТ ЛОКАЛЬНЫХ ЗЕМЛЕТРЯСЕНИЙ............................................................................. 2.1 Томография по временам пробега объемных сейсмических волн от землетрясений. Постановка задачи. Введение понятия параметризации функции.......... 2.2 Формулы для вычисления частных производных. Различные способы параметризации модели.......................................................................................................... 2.2.1 Блоковая параметризация................................................................................... 2.2.2 Параметризация с использованием Би-сплайнов............................................. 2.2.3 Параметризация модели с помощью тетраэдров, внутри которых скорость меняется линейно непрерывно........................................................................................... 2.2.4 Несколько слов о других способах параметризации модели.......................... 2.3 Используемые стандартные статистические характеристики для оценки качества полученных результатов......................................................................................................... 2.4 Алгоритмы построения сейсмических лучей........................................................... 2.4.1 Построение лучей в сферически-симметричной среде (одномерный луч)... 2.4.2 Построение лучей в горизонтально-неоднородных средах............................ 2.5 Методы обращения матриц. Вычисление матрицы разрешения. Решение обратной задачи....................................................................................................................... 2.5.1 Нахождение решения уравнения (34) с помощью построения системы нормальных уравнений....................................................................................................... 2.5.2 Нахождение решения уравнения (34) с помощью сингулярного разложения матрицы B............................................................................................................................ 2.5.3 Нахождение решения уравнения (34) методом LSQR..................................... 2.6 Методы локации землетрясений в одномерных скоростных моделях................... 2.7 Различные способы введения весовых коэффициентов.......................................... 2.8 Основные этапы любого сейсмотомографического исследования........................ 2.8.1 Сортировка данных............................................................................................. 2.8.2 Построение оптимальной одномерной модели и перелокация событий в полученной одномерной модели. Вычисление станционных поправок........................ 2.8.3 Выбор алгоритма и его тестирование................................................................ 2.8.4 Построение трехмерной модели и перелокация событий, проверка правильности проведенных расчетов................................................................................ Выводы к главе........................................................................................................................ ГЛАВА 3 СРАВНИТЕЛЬНЫЙ И ОЦЕНОЧНЫЙ АНАЛИЗ.................................. 3.1 Сравнение различных численных алгоритмов для расчета времени пробега луча и вычисления траекторий сейсмических лучей.................................................................... 3.2 Сравнение различных способов параметризации модели....................................... 3.3 Сравнение различных способов обращения матриц................................................ 3.4 Оценка влияния процедуры пересчета лучей и перелокации событий после каждой итерации на значение RMSw..................................................................................... 3.5 Тестирование различных алгоритмов построения сейсмотомографических моделей..................................................................................................................................... Выводы к главе........................................................................................................................ ГЛАВА 4 ОБЩАЯ ГЕОФИЗИЧЕСКАЯ ХАРАКТЕРИСТИКА ТЯНЬ-ШАНЯ....... 4.1 Расположение исследуемого региона на топографической и тектонической картах........................................................................................................................................ 4.2 История развития региона.......................................................................................... 4.3 Магнитные аномалии.................................................................................................. 4.4 Аномалии силы тяжести, изостатическое состояние земной коры, плотность пород........................................................................................................................................ 4.5 Тепловое поле.............................................................................................................. 4.6 Активные разломы, скорости современных деформаций....................................... 4.6.1 Скорости горизонтальных движений................................................................ 4.6.2 Скорости вертикальных движений.................................................................... 4.7 Исследования Тянь-Шаня по сейсмологическим данным...................................... 4.7.1 Общая информация............................................................................................. 4.7.2 Сейсмичность..................................................................................................... 4.7.3 Исследования методами глубинного сейсмического зондирования (ГСЗ). 4.7.4 Исследования земной коры всего Тянь-Шаня методом сейсмической томографии по данным объемных волн.......................................................................... 4.7.5 Исследования земной коры Cеверного Тянь-Шаня методом сейсмической томографии по данным объемных волн.......................................................................... 4.8 Детальное изучение зоны сочленения Чуйской впадины и Киргизского хребта.....

..................................................................................................................................... Выводы к главе...................................................................................................................... ГЛАВА 5 ОПИСАНИЕ ИСПОЛЬЗУЕМЫХ В РАБОТЕ ДАННЫХ И АНАЛИЗ СЕЙСМИЧНОСТИ ТЯНЬ-ШАНЯ.............................................................................. 5.1 Общие сведения об используемых данных............................................................. 5.2 Анализ сейсмичности и месторасположения станций на территории Тянь-Шаня по имеющимся данным......................................................................................................... 5.3 Отбор данных для построения одномерных и трехмерных томографических моделей................................................................................................................................... 5.3.1 Отбор данных для территории Cеверного Тянь-Шаня (41.9-43.4 с.ш., 73.5 76.5 в.д.).............................................................................................................................. 5.3.2 Отбор данных для всей территории Тянь-Шаня............................................ Выводы к главе...................................................................................................................... ГЛАВА 6 ПОЛУЧЕННЫЕ ОДНОМЕРНЫЕ И ТРЕХМЕРНЫЕ СКОРОСТНЫЕ МОДЕЛИ ТЯНЬ-ШАНЯ............................................................................................. 6.1 Скоростные модели Северного Тянь-Шаня............................................................ 6.1.1 Построение одномерных моделей, расчет станционных поправок и перелокация событий в одномерной модели для территории Северного Тянь-Шаня.....

............................................................................................................................. 6.1.2 Трехмерные модели Р и S волн для территории Северного Тянь-Шаня, перелокация событий в трехмерной модели................................................................... 6.1.3 Геологотектоническая интерпретация полученных результатов и сравнение с существующими моделями для территории Северного Тянь-Шаня......................... 6.2 Исследование всей территории Тянь-Шаня с использованием локальных данных.

..................................................................................................................................... 6.2.1 Построение одномерных моделей, расчет станционных поправок и перелокация событий в одномерной модели для территории всего Тянь-Шаня........ 6.2.2 Трехмерные модели Р и S волн для территории всего Тянь-Шаня, перелокация событий в трехмерной модели................................................................... 6.2.3 Геологотектоническая интерпретация полученных результатов и сравнение с существующими моделями территории всего Тянь-Шаня........................................ Выводы к главе...................................................................................................................... ВЫВОДЫ И ЗАКЛЮЧЕНИЕ..................................................................................... ПРИЛОЖЕНИЕ 1 СУЩЕСТВУЮЩЕЕ ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ ДЛЯ РЕШЕНИЯ СЕЙСМОТОМОГРАФИЧЕСКИХ ЗАДАЧ.............................................. 1. Программа Velest [21]................................................................................................... 2. Программа Sphypit90 [43]............................................................................................. 3 Программа Simulps [55]................................................................................................ 4 Программа Fatomo [19]................................................................................................. ЛИТЕРАТУРА............................................................................................................ Сокращения, термины и обозначения ТФР - Таласо-Ферганский разлом ГСЗ - глубинное сейсмическое зондирование МОСМ - минимальная одномерная скоростная модель DSW - derivative weighted sum (сумма взвешенных значений частных производных по всем лучам для данного параметра модели) RDE - resolution diagonal elements (диагональные элементы матрицы разрешения) Vp - стандартная ошибка для скоростей Р волн (измеряется в км/с) dVp/Vp - стандартная ошибка для возмущений скоростей Р волн (измеряется в %) Vp/Vs - стандартная ошибка для трехмерной модели распределения отношения Vp/Vs Р волна - продольная сейсмическая волна S волна - поперечная сейсмическая волна T=tреальн-tрассчетн - невязка времен пробега (разница между реальным и расчетным временем пробега луча);

- эпицентральное расстояние между двумя точками по дуге большого круга;

RMS - среднеквадратичная невязка по всем лучам RMSw - взвешенная среднеквадратичная невязка по всем лучам RMSсоб - среднеквадратичная невязка для данного события RMSwсоб - взвешенная среднеквадратичная невязка для данного события Введение Актуальность работы.

Несколько последних десятилетий сейсмическая томография является одним из наиболее распространенных и эффективных методов получения информации о скоростных свойствах пород внутри Земли. Существует большое количество различных сейсмотомографических алгоритмов.

Наличие этих алгоритмов с одной стороны существенно упрощает и убыстряет решение задачи восстановления трехмерной скоростной структуры в Земле. С другой стороны, не зная как устроены эти алгоритмы и какова их область применимости, трудно получить удовлетворительное решение. В этой связи особенно актуальным на сегодняшний день является сравнительное изучение физических основ различных алгоритмов, освоение этих алгоритмов и их усовершенствование для конкретного набора экспериментальных данных и особенностей геолого-тектонического строения изучаемого региона.

На данном этапе развития сейсмической томографии наиболее актуальным является исследование сложно-построенных регионов, например, такого горного массива, как Тянь-Шань. К настоящему моменту для территории Тянь-Шаня очевидно, что кора и мантия существенно неоднородны. Характер этих неоднородностей более сложный, с различной степенью контрастности, чем, например, в зонах субдукции и в местах расположения плюмов (нет ярко выраженной высокоскоростной области, связанной с погружающейся пластиной океанской литосферы, как в районах субдукции, или низкоскоростного канала, связанного с восходящей струей разогретой мантии, как над плюмами). Поэтому нужны более совершенные методы для восстановления этих неоднородностей. До сих пор при построении сейсмотомографических моделей территории Тянь-Шаня изучаемая область разбивалась на прямоугольные блоки с постоянной скоростью внутри или скоростная функция раскладывалась по полиномам Лежандра (гармонически-полиномиальный способ разложения). Оба эти способа имеют ряд существенных недостатков. При использовании гармонически-полиномиального способа часто сталкиваются с проблемой ложной экстраполяции искомой функции в слабо изученных районах. При разбивке исследуемой территории на блоки с постоянной скоростью не всегда удается правильно установить границы разноскоростных областей.

Сравнительно недавно в Тянь-Шаньском регионе наряду с аналоговой формой записи сейсмического сигнала, стало возможным производить запись сейсмического сигнала в цифровом виде. Первые цифровые станции на территории Тянь-Шаня появились в 1991 году. Количество цифровых станций на территории Тянь-Шаня с каждым годом увеличивается. В 1997 2000 гг. их было около 40. Точность определения времени вступления различных волн по цифровой волновой форме значительно выше, чем с использованием аналоговой сейсмограммы. Очевидно, что из двух скоростных моделей для одного и того же региона, та которая получена с использованием данных на цифровых станциях является более точной и детальной. Построение томографических моделей с использованием цифровых данных описано в работах [11] и [67]. В данной работе представлен анализ цифровых данных за более длительный период времени и при построении моделей используются эти данные в совокупности с наиболее совершенными сейсмотомографическими методиками.

Разработка и усовершенствование имеющихся сейсмотомографических алгоритмов для сложно-построенных сред, а также расчет трехмерных скоростных моделей с помощью этих алгоритмов и новых массивов данных для территории Тянь-Шаня актуальны при проведении сейсмологических исследований горного сооружения Тянь-Шаня. Новая информация, полученная на основе анализа рассчитанных вновь скоростных моделей совместно с новыми уточненными тектоническими картами [86;

109], в значительной степени дополняет существующие представления об особенностях геологического строения и динамических процессах, происходящих в этом регионе.

Основная цель.

Целью работы является изучение пространственных скоростных неоднородностей строения земной коры Тянь-Шаня на основе непрерывного способа параметризации модели по данным цифровых и аналоговых сейсмических станций.

Основные задачи

исследования.

- провести сравнительный анализ наиболее часто используемых различных сейсмотомографических методов, а также алгоритмов, разработанных автором на тестовых примерах;

- для конкретных условий изучаемого региона (объем экспериментального материала, относительное расположение источников и приемников, размеры территории и др.) разработать усовершенствованные алгоритмы томографической инверсии;

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

- построить наиболее адекватную скоростную модель земной коры по Р волнам для территории всего Тянь-Шаня.

Научная новизна.

Предложены два новых алгоритма TomoTetraFD и TomoCubeFD для построения непрерывных и квазинепрерывных (функция скорости претерпевает разрыв на конкретных глубинах) моделей скоростей продольных (Р) и поперечных (S) сейсмических волн. Алгоритм TomoTetraFD особенно эффективен в районах с существенно неравномерным расположением источников сейсмических волн. В отличие от многих существующих алгоритмов применение данных возможно не только для территорий локального масштаба (200*200 км), но и регионального (1000*1000 км).

Впервые проведен сравнительный анализ пяти различных сейсмотомографических программ: Sphypit90, Simulps14, Fatomo, TomoCubeFD, TomoTetraFD, в которых в качестве исходных данных использованы времена пробега продольных (Р) и поперечных (S) сейсмических волн от локальных землетрясений. Три первых алгоритма (Sphypit90, Simulps14, Fatomo) активно применяются для построения трехмерных моделей скоростей объемных волн различных регионов. Два других (TomoCubeFD и TomoTetraFD) разработаны автором на базе существующего сейсмотомографического программного обеспечения. Даны рекомендации по поводу того, в каких случаях удобнее всего использовать тот или иной алгоритм.

По данным времен пробега от местных землетрясений, зарегистрированных на аналоговых и цифровых станциях, впервые построена квазинепрерывная трехмерная скоростная модель Р волн всей территории Тянь-Шаня. Предложен и реализован на практике ряд тестов, подтверждающих устойчивость полученного результата. Проведена геолого геофизическая интерпретация полученного результата. Впервые проведена локация большого количества землетрясений и взрывов, произошедших на всей территории Тянь-Шаня методом сеточного поиска [34] в трехмерной скоростной модели.

Применен новый подход к построению трехмерных скоростных моделей P и S волн для верхней части коры под Северным Тянь-Шанем, а также для анализа отношения Vp/Vs в этом регионе. Этот подход позволяет выявить наиболее достоверные скоростные неоднородности, т.к. включает в себя построение целого множества трехмерных скоростных моделей с использованием различных наборов данных, различных сейсмотомографических алгоритмов для Р волн и совместно для Р и S волн.

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

Защищаемые научные положения.

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

Обнаруженные скоростные неоднородности в верхней коре под Тянь Шанем хорошо согласуются с геологической и тектонической картой данного региона.

Практическая значимость.

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

Установленные скоростные модели могут быть также использованы в интересах Международной системы сейсмического мониторинга (IMS) для контроля за соблюдением Договора о всеобъемлющем запрещении ядерных испытаний. Полученные скоростные разрезы могут быть использованы для расчета очагово-станционных сейсмических поправок SSSC для станций IMS Тянь-Шаньского региона. В работе [42] показано, как на основе высокоточной трехмерной скоростной модели для Индо-Пакистанского региона произведен расчет поправок SSSC и с их помощью проведена более уточненная локация событий по станциям IMS.

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

Глава 1 История развития сейсмической томографии и существующие томографические модели В первой главе описаны некоторые исторические факты, предшествующие построению первых трехмерных моделей скоростей Р и S сейсмических волн. Проводится краткий обзор работ, в которых построены сейсмотомографические модели, сильно повлиявшие на существующие сейчас представления о строении Земли и процессах, происходящих в ее недрах.

1.1 Предыстория создания первых трехмерных моделей скоростей сейсмических волн К концу шестидесятых годов прошлого века развитие наук о Земле вышло на ту стадию, когда стало возможным и необходимым построение трехмерных горизонтально неоднородных моделей скоростей Р и S сейсмических волн.

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

Возможность построения таких моделей объясняется рядом причин.

Во-первых, был накоплен обширный экспериментальный материал, который позволил выявить основные черты внутреннего строения Земли в рамках одномерных скоростных моделей. Была создана довольно точная сферически симметричная модель внутреннего строения Земли и на ее основе рассчитаны таблицы времен прихода сейсмических волн в зависимости от эпицентрального расстояния - хорошо известные таблицы Джеффриса Буллена и Херрина, EASPEI-91. Во-вторых, была создана плотная мировая сеть сейсмологических наблюдений, большие сейсмические антенны (LASA в США в 1965г., NORSAR в Норвегии в 1972г.), развиты национальные системы сейсмологических наблюдений. В-третьих, была разработана методика построения одномерных скоростных моделей, а также методы построения сейсмических лучей. В-четвертых, развитие вычислительной техники позволило реализовать предложенные алгоритмы построения сейсмических лучей и расчета трехмерных моделей скоростного строения Земли по объемным волнам.

1.2 Некоторые существующие трехмерные модели скоростей Р и S волн Первые сейсмотомографические модели появились около 30 лет назад.

На сегодняшний день существует большое количество глобальных и региональных скоростных моделей Земли. Для построения скоростных моделей отдельных регионов используют данные времен пробега объемных (Р и S) и высокочастотных поверхностных волн. В глобальных моделях в качестве начальных данных часто используются: короткопериодные продольные объемные волны, прошедшие без отражений, отраженные от земного ядра, а также прошедшие через внешнее земное ядро;

длиннопериодные объемные волны с периодами от 45 до 200 сек;

поверхностные волны с периодами от 35 до 250 сек.;

собственные колебания всей планеты (3-55 мин). Региональные скоростные модели строятся по данным времен пробега сейсмических волн как от локальных землетрясений или взрывов (т.е. источников, расположенных на исследуемой территории), так и по данным времен пробега от телесейсмических землетрясений или взрывов (удаленных на большие расстояния от региона).

Одна из первых работ, посвященных построению региональной скоростной модели по локальным данным, была проведена K.Aki и W.H.K.Lee в 1976 году [3]. Они построили трехмерную модель скоростей Р волн для верхних 15 км коры в районе разломов Сан Андреас и Калаверас в Калифорнии и произвели перелокацию источников в этом районе.

Впоследствии глубинные разломы в Калифорнии методом сейсмической томографии изучались в работах [7;

8;

57;

64]. После перелокации событий в трехмерной скоростной модели большинство гипоцентров землетрясений располагаются ближе к тектоническим разломным структурам. Границы высоко- или низкоскоростных аномалий часто совпадают с линиями разломов, которые были обнаружены с помощью геологических методов. В перечисленных работах для расчетов использовались данные от землетрясений, произошедших в верхней коре до глубин 10-15 км. В работе [46] использовались данные от землетрясений, произошедших на значительно больших глубинах до 250 км. С помощью этих данных С.Рокеру [46] удалось выявить низкоскоростную область, простирающуюся под Памиро-Гиндукушской зоной в Средней Азии до глубин 250 км.

По данным удаленных землетрясений в нашей стране первая трехмерная скоростная модель была построена А.С. Алексеевым и М.М.

Лаврентьевым на основе наблюдений времен пробега на сейсмологическом профиле Памир-Байкал [69]. Одной из основополагающих работ телесейсмической томографии, в которой подробно описывается методика построения скоростных моделей верхней мантии с использованием данных от удаленных землетрясений, также является работа [2].

По количеству работ в области глобальной сейсмической томографии особое место занимают работы американской школы Адама Дзивонского [5].

На сегодняшний день для большинства глобальных сейсмотомографических моделей в верхней мантии на глубинах 78-148 км с сейсмически активными районами связаны низкоскоростные аномалии, а с сейсмически спокойными регионами высокоскоростные аномалии.

Исследованию пространственного строения литосферы по данным поверхностных волн посвящены работы Т.Б. Яновской [122].

Помимо трехмерных моделей скоростей Р и S волн с помощью сейсмической томографии по локальным данным также строятся карты отношения Vp/Vs для различных глубин [64], которые несут в себе информацию о химическом составе, слагающих земные недра пород.

В последнее время для регионального изучения крупных тектонических структур Европейским Геофизическим Союзом проводился ряд исследований с привлечением различных стран. На территории Германии, Дании, Швеции проводился глобальный эксперимент TOR (Transeuropaen Suture Zone) в 1996/97 годах, на территории Финляндии и России эксперимент Svekalapko в 1999 году. Основная цель этих экспериментов - более детальное изучение внутреннего строения в этих регионах за счет создания плотной сети переносных и стационарных сейсмических станций (при проведении эксперимента Svekapko станции располагались через каждые 50 км). Полученные данные были использованы для изучения горизонтальных неоднородностей не только методом сейсмической томографии, а также методом обменных волн и другими.

Наиболее яркими сейсмотомографическими результатами являются трехмерные скоростные модели для зон субдукции, разломных зон и областей плюмов. Ниже на Рис. 1а представлено вертикальное сечение трехмерной скоростной модели, полученной А.Горбатовым для зоны субдукции под Камчатским полуостровом [13]. Скоростные сейсмотомографические модели для зон субдукции под Новой Зеландией построены в работе [64], для Идзу-Бонинской зоны (Тихий океан) в [63], для района Памир-Гиндукуш в [46]. С помощью сейсмической томографии во всех работах (аналогично Рис. 1а) удается проследить, как субдуцирующая холодная высокоскоростная литосферная плита пересекает астеносферный слой и следует наклонно до поверхности нижней мантии (670 км). На Рис. 1а построена скоростная томографическая модель в районе субдукции до глубин 200 км. Трехмерные скоростные неоднородности над плюмами обсуждаются в работах [54] (Гавайский плюм) и [10] (Исландский плюм).

Рис. 1 Вертикальные сечения трехмерных скоростных моделей: а) под Камчатским полуостровом из [13], б) сейсмотомографическая модель Исландского плюма [10].

В данных работах удается проследить, как огромная капля перегретого низкоскоростного вещества поднимается из недр к коре (см. Рис. 1б).

В последнее время с помощью сейсмотомографического метода также производят расчет трехмерных моделей параметра добротности (характеристика затухания сейсмических волн) [51;

27], анизотропных скоростных моделей [18] и трехмерных моделей температур вещества, слагающего земные недра.

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

Глава 2 Методика построения трехмерных региональных томографических моделей по данным времен пробега от локальных землетрясений В первой главе приводились примеры расчетов различных томографических моделей. Вторая глава посвящена теоретическим основам сейсмотомографического алгоритма. В данной главе с позиции автора подробно описана суть предложенного K.Aki и W.H.K.Lee в 1976 году подхода и дальнейшие модернизации этого подхода, проведенные большим количеством ученых. При описании каждой компоненты метода отмечается, проводилось ли автором только изучение этой компоненты, или же дополнительно тестирование и совершенствование.

2.1 Томография по временам пробега объемных сейсмических волн от землетрясений. Постановка задачи. Введение понятия параметризации функции.

Как уже упоминалось во введении, сейсмическая томография позволяет, используя данные времен пробега объемных сейсмических волн, рассчитывать трехмерные статистически обоснованные модели распределения скоростей сейсмических волн в коре и мантии.

Время пробега сейсмической волны определяется функционалом Ферма, в котором интегрирование выполняется вдоль экстремали (луча).

dl T =,.

(1) v (r ) L v(r) - скорость сейсмической волны в точке r, L - траектория сейсмического луча, также зависящая от скорости сейсмических волн в среде, координат источника. Координаты землетрясений и время, когда оно произошло, точно не известны, поэтому эти характеристики вносят ошибку в определение траектории луча. Координаты приемников обычно известны с высокой точностью. В более общем виде время пробега сейсмической волны T, распространяющейся от некоторого землетрясения с координатами xист yист zист и временем в источнике tист на станцию, можно представить как некую функцию от значений скоростей v(r) в различных точках пространства, координат землетрясений xист yист zист и времени в источнике tист.

T = f (v (r ), xист, yист, zист, tист ),.

(2) Основная задача заключается в определении скорости v(r), координат землетрясений и времени в источнике по множеству измерений времен пробега T на поверхности для различных лучей.

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

M v ( r ) = h ( r ) (3) k k k = hk(r) - М базисных функций, - различные М коэффициентов при этих k функциях. Базисные функции представляют собой набор стандартных функций, значения которых обычно изменяются от нуля до единицы. В сейсмической томографии часто осуществляется параметризация не функции скорости, а функции медленности (величины, обратной скорости) [44], или квадрата медленности [11]. Все дальнейшие выкладки этого пункта верны в не зависимости от того, какая из функций (скорость, медленность или квадрат медленности) раскладывается по набору базисных функций.

Учитывая (3) равенство (2) можно переписать следующим образом:

T = f (,K,, xист, yист, z, tист ),.

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

В основу сейсмотомографического метода заложены принципы теории возмущений.

Предположим, что известна начальная приближенная скоростная модель v0(r), приближенные координаты события xист0 yист0 zист0 и время в источнике tист0. При параметризации начальной скоростной модели возможно определить начальные коэффициенты, соответствующие заданным k базисным функциям. Время пробега сейсмической волны от землетрясений с координатами xист0 yист0 zист0 и временем в источнике tист0 в начальной скоростной модели обозначим T0.

0 0 0 0 0 T0 = f (,K,, xист, yист, zист, tист ),.

(5) 1 k Если искомая скоростная модель, координаты источников и время в источнике незначительно отличаются от нулевого приближения, то можно воспользоваться разложением Тейлора по малому параметру для функции многих переменных. Здесь необходимо отметить, что используются слагаемые только первого порядка малости и пренебрегается слагаемыми более высокого порядка малости. Получаем:

0 0 0 0 0 T = f ( +,K, +, xист + x, yист + y, zист + z, tист + t ) 1 1 k k 0 0 0 0 0 f (,K,, xист, yист, zист, tист ) +.

1 k (6) M T T T T + + x + y + z + t k x y z k = k Разница между временем пробега T и T0 равна:

M T T T T T - T = T = + x + y + z + t.

(7) k x y z k = k Основная сейсмотомографическая задача свелась к определению поправок к коэффициентам, соответствующим нулевой модели, поправок к координатам источника и поправки ко времени в очаге по разнице времен пробега T для различных лучей. Другими словами, необходимо решить систему линейных уравнений, количество уравнений в которой равно количеству сейсмических лучей, а неизвестными являются величины k, x, y, z и t.

Как было сказано выше, задачу нахождения поправок k, x, y, z и t, зная невязки T, возможно свести к линейной системе уравнений только в том случае, если при разложении функции T в (6) пренебречь членами с k2, x2, y2, z2, t2, k3, x3 и.т.д. Т.к. слагаемые, в которые искомые поправки входят во второй, третьей и более высоких степенях не учитывается, решение системы линейных уравнений типа (7) не является решением, поставленной задачи. Чтобы получить истинные значения k, x, y, z и t процедуру решения линейной системы уравнений типа (7) необходимо повторить итеративно несколько раз согласно формуле (8) до тех пор, пока значения искомых поправок на очередной итерации не будут пренебрежимо малы.

i i i i M T T T T i i +1 i +1 i +1 i +1 i + T = + x + y + z + t, (8) k x y z k = k где i и i+1 - номера итераций. В этом случае искомые поправки k, x, y, z и t равны сумме поправок, полученных после каждой итерации (9).

i i i i i = x = x y = y z = z t = t, k k (9) i i i i i Далее рассматриваются различные способы решения системы уравнений типа (7) только на первой итерации, т.к. методы, применяемые при решении систем на каждой итерации, одинаковы.

В матричном виде система уравнений типа (7) имеет вид K T = A + H hn, где n n = T T = T1,K, TN (10) T T =,K,, hn = tn, xn, yn, zn, 1 M Ti Ti Ti Ti Ti A = Aik =, H =, n tn xn yn zn k N - количество лучей, М - количество коэффициентов при базисных функциях, К - количество землетрясений. Матрица А, элементами которой являются частные производные времени пробега по коэффициентам (производные Фреше), в дальнейшем будет называться матрицей из частных производных.

Ti Ti Ti Ti Ti,,,, Вычисление частных производных tn xn yn zn для всех k лучей i=1,Е,N есть результат решения прямой задачи. Обратная задача, hn включает в себя нахождение поправок, удовлетворяющих k уравнению (10), зная матрицы A, Hn и вектор T.

2.2 Формулы для вычисления частных производных.

Различные способы параметризации модели.

Обозначим функцию скорости - v(r), функцию медленности - s(r) и функцию квадрата медленности - s2(r). Выражение для времени пробега через функцию скорости приведено в (1), ниже приведены выражения для времени пробега через функцию медленности (11) и через функцию квадрата медленности (11):

T = s (r )dl (a ) L.

(11) T = s (r )dl (б ) L При параметризации модели искомая функция раскладывается по некоторому набору базисных функций с определенными коэффициентами. В (12) приведен общий вид разложения функций s(r) и s2(r) по набору базисных функций.

M ` s ( r ) = h ( r ) k k k = (12), M 2 `` s ( r ) = h ( r ) k k k = Эти разложения аналогичны разложению (3) для функции v(r). Исходя из (3) и (12) с учетом (1) и (11) общий вид формул для вычисления производных Фреше для трех видов функций следующий:

v (r ) Ti hk (r ) = dl = - dl i i v (r ) k k Li Li Ti s (r ).

(13) = dl = hk (r )dl i i ` ` k k Li Li Ti s (r ) 1 hk (r ) = dl = dl i i `` `` s (r ) k k Li Li Из формул (13) видно, что проще всего вычислить производные, если параметризовать функцию медленности. Производные, вычисляемые для функций s(r) и s2(r) всегда положительны, следовательно, матрица A в (10), которую они составляют, является положительно определенной. Если Ti используется функция v(r), то все производные - отрицательные и k положительно определенной является матрица ЦA.

Параметризацию модели возможно проводить в классе кусочно постоянных, в классе непрерывных, в классе один раз непрерывно дифференцируемых, дважды-непрерывно дифференцируемых (и.т.д. до бесконечности) функций. Ниже более подробно описаны те виды параметризаций, которые использовались автором в ходе работы.

2.2.1 Блоковая параметризация.

Блоковая параметризация применяется и в декартовой и в сферической системе координат и является наиболее распространенной, как в глобальной, так и в локальной томографии. Исследуемый объем разбивается на конечное число блоков. В декартовой системе координат блоки представляют собой параллелепипеды, в сферической системе координат сферические сегменты.

Считается, что искомая функция (v(r), s(r) или s2(r)) постоянна внутри каждого блока. Количество базисных функций равно количеству блоков.

Каждая базисная функция имеет вид:

1, если r k - ому блоку hk (r ) =.

(14) 0, если r k - ому блоку вокруг любого построенного тетраэдра, не попадало ни одной из заданных точек [105]. Считается, что разбиение на тетраэдры, удовлетворяющее условию Делоне, является наиболее равномерным из всех возможных.

Расчет базисных функций производится по формуле [49]:

(r - r1) {(r2 - r1) (r3 - r1)} hk (r) = (18) (rk - r1) {(r2 - r1) (r3 - r1)}, где rk-положение точки k, r1, r2, r3 - положение трех других точек тетраэдра, содержащего точки r и rk. Параметризация модели с помощью тетраэдров применена в работах И.Ю.Кулакова [93;

92], а также в [64].

2.2.4 Несколько слов о других способах параметризации модели.

Возможно производить параметризацию с использованием ортогональных многочленов Лежандра или Чебышева [73]. Такой способ параметризации применен в [108]. При построении глобальных моделей часто используется разложение функции в конечный ряд полностью нормированных сферических гармоник [95]. Глобальные модели с использованием этого типа параметризации построены А.Морелли и А.Дзевонски [95].

Автором работы подробно изучены 4 различных типа параметризации, а также написано несколько программ для того, чтобы использовать тот или иной вид параметризации в процессе сейсмотомографической инверсии.

Краткая информация о степени участия автора в разработке того или иного способа параметризации представлена в Табл. 1.

Табл. 1 Степень участия автора в разработке четырех различных способов параметризации модели (И - изучен;

П - запрограммирован автором;

С - используется при расчетах реальных скоростных моделей (может быть использован и не в программной реализации автора);

У - усовершенствован автором;

Т - использовался при тестировании.) Способ параметризации модели Степень участия автора Блоковая (2.2.1) ИС Линейными Би-сплайнами (2.2.2) ИС Кубическими Би-сплайнами (2.2.2) ИПСУТ С помощью тетраэдеров (2.2.3) ИПС 2.3 Используемые стандартные статистические характеристики для оценки качества полученных результатов Обычная и взвешенная среднеквадратичные невязки по всем лучам RMS и RMSw рассчитываются по формулам:

N Ti i= RMS = N (19), N WiTi i= RMS =, Wi = W1i W2i W3i W4i, W N Wi i= где Ti - разница между реальным и расчетным временем пробега (невязка) для i-того луча, N - количество лучей, Wi - итоговый весовой коэффициент, W1i - весовой коэффициент, зависящий от качества вступления волны, W2i - весовой коэффициент, зависящий от величины невязки Ti, W3i - весовой коэффициент, зависящий от эпицентрального расстояния, W4i - весовой коэффициент, зависящий от качества регистрации на станции.

Обычная и взвешенная среднеквадратичные невязки для одного события RMSсоб и RMSwсоб вычисляются по формулам:

N k Ti соб i = RMS = N k (20) N л W T i i соб i = RMS =, W N л W i i = где Nk - количество лучей, зарегистрированных от k-го события, остальные обозначения такие же как и для (19).

Количество степеней свободы системы из N уравнений с M неизвестными равно:

, N = N - M (21) св Вариация данных вычисляется по формуле (22) и имеет размерность сек2:

N W Ti i i = var =, (22) N св Вариация параметров модели вычисляется по формуле (23) и имеет размерность (км/сек)2:

K V i i = var =, (23) Д N св где Vi - найденная i-ая поправка к первоначальному значению скорости, К - количество неизвестных скоростных параметров.

2.4 Алгоритмы построения сейсмических лучей 2.4.1 Построение лучей в сферически-симметричной среде (одномерный луч) Для однородной слоистой среды траекторией луча является ломаная линия, состоящая из отрезков прямых линий, соответствующих каждому слою. На границах слоев луч преломляется по закону Снеллиуса r sin = const (24), c где - угол между лучом и радиусом, r- радиус Земли, с - скорость сейсмической волны в слое. Чтобы упростить расчеты используют плоскопараллельное преобразование для перехода из сферической системы координат в декартову [68]. Rz и vсфvпл переводятся по формулам:

z r RЗ RЗ (25) e = vпл (z) = vсф (r) RЗ r Если задана непрерывная зависимость скорости от глубины, тогда для расчета траектории луча и времени пробега в сферической среде используют формулы R 2dr T ( p) = r(2 - p2)1/ rmin (26) R, pdr ( p) = r(2 - p2)1/ rmin где =r/c, p=rsin/c - лучевой параметр, r - радиус Земли, (p)-угол в радианах. Зависимости (26) выводятся из закона Снеллиуса в предположении бесконечно тонких слоев.

Вышеприведенные выражения широко применяются в телесейсмической томографии при построении лучей в стандартных одномерных скоростных моделях - IASPEI91, PREM, Herrin model.

Например, в программе С.Рокера Sphypit90 (подпрограмма delcal2.f). Также построение лучей в горизонтально-неоднородных средах в локальной томографии иногда сводится к формулам (24) и (25). Например, метод, описанный в [56].

2.4.2 Построение лучей в горизонтально-неоднородных средах В неоднородных средах фронт (x,y,z) сейсмической волны в лучевом приближении в некоторый фиксированный момент времени t описывается уравнением эйконала:

2 (27) + + =, x y z v2 (x, y, z) 2.4.2.1 Расчет лучевых трасс аналитически 1.Для некоторых функций v [72;

114] уравнение (27) решается аналитически. В [72] в аналитической форме приведена запись решения уравнения (27) для скоростных функций следующего вида:

v(x, y, z) = Px + Qy + Rz (P2 + Q2 + R2 > 0) (28), v(x, y, z) = k(x2 + y2 + z2 ) ( > 0) Если скоростная функция представима в одном из видов (28), тогда фронт волны, распространяющейся от точечного источника с координатами (x0,y0,z0) является сферическим, а траектория сейсмического луча представляет собой дугу окружности. В [114] представлено аналитическое решение уравнения (27) для функции v вида (29).

v(x, y, z) =, (29) 1/ (S + Px + Qy + Rz) 2.4.2.2 Расчет лучевых трасс численно. Методы изгиба траектории 1. Численный метод, предложенный С.H. Thurber и W.L.Ellsworth в 1980 году, сводится к построению лучей в одномерной слоистой среде. Суть метода изложена в [56], а также на русском языке в [90]. В трехмерной горизонтально-неоднородной скоростной модели для каждого луча вычисляется своя специальная усредненная слоистая скоростная модель и затем рассчитывается луч в слоистой скоростной модели. Данный метод используется в программе С. Рокера SPHYPIT90. Ограниченность использования метода состоит в том, что его возможно использовать только в случае блоковой параметризации скоростной функции или функции медленности. Недостатком является также тот факт, что даже в горизонтально неоднородном скоростном поле траектория сейсмического луча, рассчитанная эти методом, является плоской, т.е. существует плоскость, которой принадлежат все точки рассчитанной траектории. С другой стороны скорость работы алгоритма очень высокая и даже если эпицентральные расстояния достигают 1000 км, расчет траекторий лучей на ПК Athlon 1,33ГГц и ОП 256 МБ занимает меньше минуты.

2. Метод окружностей (approximate ray tracing ART). Суть метода описана в работе [57]. Использовался в [108] для вычисления трехмерной скоростной модели коры Тянь-Шаня. Этот метод также запрограммирован в Simulps (подпрограмма rayweb.f), используется для построения и приближенной лучевой траектории. В работе не тестировался. Все точки результирующей траектории также как и в предыдущем методе принадлежат одной плоскости.

3. Метод изгиба траектории был предложен в работе [59], он широко распространен благодаря своей простоте и быстроте вычислений. Суть метода на русском языке также подробно описывается в [90]. Этот алгоритм используется при построении лучей в программе Simulps (подпрограмма minima.f). Также автором написана собственная программа для расчета сейсмических лучей эти методом. Большое количество сейсмотомографических моделей получено с использованием именно обсуждаемого подхода для вычисления траекторий сейсмических лучей, например в [7]. В отличие от методов 1 и 2, рассчитанная траектория сейсмического луча является трехмерной. Необходимым условием метода, описанного в [59] является существование первой производной функции скорости в любой точке пространства, что ограничивает множество скоростных полей, в которых можно использовать этот метод.

При упоминании этого метода в следующей главе будут использоваться его две различные модификации: М3-1 и М3-2. В М3-1 в качестве нулевого приближения при построении лучей используется прямая линия. В М3-2 в качестве нулевого приближения используется траектория, рассчитанная с использованием метода, описанного в параграфе 2.4.2. (пункте 2).

В 1998 году K.Koketsu и S.Sekino [23] метод был усовершенствован.

Появилась возможность использовать данный алгоритм для построения лучей не в декартовых, а в сферических координатах, а также в моделях, где функция скорости в каких-то точках претерпевает разрыв, т.е. в этих точках первая производная функции скорости не определена. В этом обновленном варианте метод использовался в работе [14].

Общий недостаток, свойственный группе методов изгиба траектории, заключается в недостаточной чувствительности этих методов к резким изменениям скорости. Если в модели имеется скоростное включение, существенно изменяющее ход луча, при использовании метода изгиба существует вероятность того, что это резкое изменение скорости не отразится в траектории. Частичная демонстрация этого отрицательного свойства методов изгиба приведена в следующей главе.

2.4.2.3 Расчет лучевых трасс численно. Методы пристрелки (shooting) 4. Суть метода Перейры для расчета траекторий сейсмических лучей [39] состоит в численном решении системы дифференциальных уравнений первого порядка с нелинейными граничными условиями. Теория метода описана в [90]. Расчеты с использованием этого метода в работе не проводились. Существенным недостатком по сравнению с 1-3 и 5 является большая громоздкость вычислений.

5. Суть метода Рунге-Кутта описывается в [68]. В этом методе для системы дифференциальных уравнений вида (30) dxi dmi = vmi, = (1/ v) ds ds xi (30) где m1 =, m2 =, m3 =, x y z в которой (m1,m2,m3) Цвектор медленности, (x1(s),x2(s),x3(s))- траектория луча, заданная параметрически, s - длина дуги луча, остальные обозначения такие же как и в (27), составляется некоторая конечно-разностная схема. Система дифференциальных уравнений (30) выводится из (27) [114]. Вообще метод Рунге-Кутта n-го порядка предполагает ошибку ограничения порядка hn+1, в работе поиск решения производился методом Рунге-Кутта второго порядка.

Сначала с помощью метода Рунге-Кутта рассчитывались траектории трех сейсмических лучей, вышедших под разными углами из источника. Зная точки выхода этих трех лучей на поверхность, можно вычислить частные производные функции изменения местоположения конечных точек луча на поверхности от углов выхода и найти поправки к углам выхода.

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

6. Параксиальный метод построения лучей был изучен и применялся в расчетах в ходе работы. Основы метода изложены в [619]. Программная реализация этого метода выполнена J.Viriux. F.Haslinger модифицировал существующую программу Simulps, таким образом, чтобы в ней параксиальный метод построения лучей также мог быть применен [17].

Наиболее удобно производить вычисления по этому методу, если осуществлять параметризацию функции квадрата медленности с помощью набора кубических Би-сплайнов, т.к. для использования этого метода необходимо существование первой и второй производной квадрата медленности в любой точке пространства. В процессе работы алгоритма сначала производится вычисление отсчетного луча, выходящего из источника и приходящего в какую-то точку на поверхности Земли. Затем корректируются углы выхода луча по отношению к отсчетным углам и рассчитываются поправки для всех точек траектории отсчетного луча. Далее производится проверка насколько точка выхода скорректированной траектории на поверхность ближе к приемнику, чем точка выхода, которая была до коррекции. Недостаток метода 6 в том, что в любой точке пространства необходимо знание не только первой, но и второй производной заданной функции скорости или квадрата медленности.

В методах 4 и 5 искомая траектория сейсмического луча состоит из конечного числа прямолинейных отрезков и для каждого отрезка составляется своя система дифференциальных уравнений, набор этих систем дифференциальных уравнений в 4 решается одновременно, а в последовательно. При использовании параксиального метода построения луча нет необходимости решать большое количество систем дифференциальных уравнений, т.к. выведены формулы, по которым, зная поправки к рассчитанным углам выхода, можно вычислить пространственную поправку к любой точке траектории отсчетного луча.

Благодаря этим формулам существенно сокращается время вычислений.

2.4.2.4 Расчет лучевых трасс численно. Методы конечных разностей 7. В вышеописанных методах 1-6 сначала производится построение сейсмических лучей от источника к приемнику, а затем расчет времени пробега вдоль этих лучей. В методах 7 и 8 используется другая тактика.

Сначала определяется время пробега от источника до любой точки в заданном объеме. Зная время пробега в каждой точке пространства, появляется возможность построить поверхность с одинаковым временем пробега (волновой фронт). Затем можно вычислить траекторию сейсмического луча, руководствуясь правилом, что сейсмический луч в каждый момент времени перпендикулярен волновому фронту.

Для того, чтобы вычислить время пробега от заданного источника до любой точки в пространстве, сначала весь исследуемый объем разбивается на достаточно мелкие блоки-кубики. Потом вычисляется значение медленности, характерное для каждого блока. Затем, представляя частные производные в (27) в виде конечных разностей, вычисляется время пробега в одной из вершин куба, зная времена пробега в других вершинах.

Для большинства представленных в данной работе скоростных моделей для расчета невязки и траектории лучей используется именно метод конечных разностей, а точнее программа time3d.c, написанная P.Podvin и I.Lecomte в 1991 году [41]. Согласно [64] также существует другая программная реализация метода конечных разностей, авторами которой являются Hole и Zelt (1995) (модифицирована Vidale (1990)).

Ошибка вычисления времени пробега, которая возникает в процессе применения метода конечных разностей, главным образом связана с размером блоков, на которые разбивается объем. Чем блоки мельче, тем ошибка меньше. С другой стороны при увеличении количества блоков требуется больше машинного времени для расчета поля времен пробега для данного источника и больше дискового пространства для хранения массивов.

8. Метод расчета сейсмического луча ненулевой толщины (fat ray).

Вычисляется множество F тех точек, которые при распространении сейсмической волны от источника, находятся внутри первой зоны Френеля.

Первой зоной Френеля называется область, в которой разность хода волн меньше чем половина длины волны. Согласно [19;

19] точка х принадлежит множеству F, если она удовлетворяет неравенству (31).

T tsx + trx - tsr (31), где Т - период сейсмической волны, tsr - время пробега вдоль сейсмического луча, tsx - время пробега от источника к точке x, trx - время пробега от приемника к точке x. Легко проверить какие точки принадлежат, а какие не принадлежат множеству F, если знать для всего исследуемого объема поле времен пробега, соответствующих данному источнику и поле времен пробега, соответствующих данному приемнику. Чтобы вычислить поля времен пробега, в данной работе применен метод [41], изложенной пунктом ранее.

Количество вычислений, производимых с использованием методов 7- значительно больше, чем с использованием методов 1-6. В этом заключается основная трудность работы с этими методами. Заметим дополнительно, что для построения толстого луча (метод 8) требуется приблизительно в два раза больше машинного времени, чем при построении луча методом 7. Т.о.

для проведения быстрых оценочных вычислений методы 7-8 не подходят.

Они нужны при проведении детального исследования высокой точности.

Большим плюсом методов 7 и 8 является также возможность использовать их в скоростных полях, которые имеют разрывы первого рода. Возможность вычислить достаточно быстро с помощью метода 7 поле времен пробега от данного источника для всех точек сетки очень важна при определении координат и времени в источнике, т.к. позволяет максимальным образом уменьшить значения RMSсоб или RMSwсоб (20). Ниже в Табл. 2 представлена степень участия автора в развитии восьми разных методов построения лучей.

Табл. 2 Степень участия автора в развитии восьми разных методов построения лучей (И - изучен;

П - запрограммирован автором;

С - используется при расчетах реальных скоростных моделей (может быть использован и не в программной реализации автора);

Т - использовался при тестировании.) Название метода построения лучей Степень участия автора М1 Метод Ellswoth ИС М2 Метод окружностей ИСТ M3 Метод изгиба траектории ИПСТ М4 Параксиальный метод ИТ М5 Метод Pereyra И М6 Метод Runge-Kutta ИП М7 Метод конечных разностей ИСТ М8 Построение толстого луча ИСТ 2.5 Методы обращения матриц. Вычисление матрицы разрешения. Решение обратной задачи Матрицу (10) иначе можно переписать в виде (32). Смысл обозначений такой же, как и в (10).

A`x = T T K* 64748 4M T T = T1,K,TN, x = 1,K,,t1,K,zK M K* M 6 6 (32), L 0 Ti A` = 0 0 k 0 0 L Первые М строк матрицы А` представляют собой разреженную матрицу, остальные 4*K строк имеют вид блочно диагональной матрицы.

Матрица А` является переопределенной, т.е. N>M+4K. Т.о. точного решения системы (32) не существуют. Чаще всего решение системы (32) ищется по норме L2, т.е. находится такой вектор x из всех возможных, который удовлетворяет условию:

N M +K* ` min Aik xk - Ti.

(33) i=1 k = Из-за наличия ошибок в данных и ошибок, возникших при решении прямой задачи, удовлетворяющие условию (33) поправки к значениям скоростей могут во много раз превосходить начальные значения скоростей, а поправки к параметрам гипоцентров могут быть значительно больше реально допустимых. Поэтому были разработаны методы для учета имеющейся информации об ошибках в данных и сведений о скоростных свойствах в исследуемом регионе, а также методы для нахождения сглаженного решения, т.е. решения в котором все искомые поправки значительно меньше самих параметров, к которым они относятся.

Учет ошибок во входных данных становится возможен с введением ковариационной матрицы данных Сd, которая для (32) имеет размер (M+4K)(M+4K) и умножается слева на матрицу A` и вектор невязок T.

Процедура умножения на матрицу Сd иначе называют взвешиванием каждого уравнения в системе. Если матрица Сd является диагональной, то из правил матричного умножения следует, что произведение Сd на A` эквивалентно умножению всей строки матрицы A` на соответствующий диагональный элемент. Т.о. если известно, что точность данного измерения низкая, то все члены уравнения, соответствующего этому измерению домножаются на некоторое число и это уравнение, влияет меньше на общий ход решения, чем остальные. Более подробно различные способы взвешивания, применяемые в сейсмической томографии, описаны в разделе 2.7. В данной работе используется ковариационная матрица данных только диагонального вида.

Также существует ковариационная матрица параметров модели Сm, которая имеет размер NN. Эта матрица умножается справа на вектор данных. Если матрица Сm имеет диагональный вид, то умножение справа эквивалентно умножению целого столбца матрицы A` на соответствующий диагональный элемент. В работе используется только диагональная ковариационная матрица параметров модели.

Для получения сглаженного решения вводится система дополнительных уравнений. После проведения вышеописанных модификаций системы (32) для получения более устойчивого и гладкого решения имеем систему (34):

B b Dz = 0, где B = Cd A`Cm, (34) b = Cd T, z = Cm1x.

Условие для нахождения решения z по норме L2 теперь с использованием обозначений из (34) выглядит следующим образом:

2 N M +K*4 S M +K*, min zk - b + zk где B = Bik, D = Dik (35).

Bik Dik i=1 k=1 i=1 k= Матрица D называется матрицей демфирования. В работе используются демфирующие матрицы двух разных типов. Первая демфирующая матрица D1 является диагональной (36) и имеет размер (M+4K)(M+4K).

1 0 D1 = 0 i (36), 0 0 M +4K где 1,, М+4K - некоторые заданные числа, которые называются демферами, они могут быть равны или не равны между собой. Введение этой матрицы позволяет найти решение незначительно отличающееся от нуля.

Смысл введения второй демфирующей матрицы D2 состоит в том, чтобы сгладить решение за счет уменьшения разницы между искомыми возмущениями, соответствующими двум соседним блокам или двум соседним точкам сетки. Эта матрица является переопределенной и имеет следующий вид:

1 j k l 1 0 -1 1 0 0 - D2 = (37) 0 L L L, 0 1 0 - 0 0 1 - где - демфирующий коэффициент, 1 и k, 1 и l, j и l, k и l - номера соседних блоков или соседних точек сетки.

После введения демфирующих матриц имеет смысл вводить параметр, оценивающий насколько полученное решение близко к идеальному по норме L2, т.е. к решению, получаемому при выполнении условия (33). Такая оценка обеспечивается матрицей разрешения. Если решение уравнения (32) обозначить xист, а решение уравнения (34) xвыч, тогда для матрицы разрешения R будет верно равенство:

Rxист = xвыч, (38) При нахождении решения (34) методом нормальных уравнений (см. 2.5.1) или с помощью сингулярного разложения (см.2.5.2) xвыч представимо в виде xвыч = Fb, (39) где F - некоторая матрица. С использованием (39) для матрицы R будет верно равенство:

(40) R = FA, Матрица разрешения R имеет размер (M+4K)(M+4K), значение каждого элемента матрицы R всегда лежит на отрезке [0,1]. Чтобы информацию, содержащуюся в матрице разрешения представить наглядно, строят карты диагональных элементов матрицы разрешения. Другой способ представления информации, скрытой в матрице R, - это расчет спрединговой функции [58;

33].

Также часто производится для оценки качества полученных результатов производится вычисление матрицы ковариации. С использованием (39) выражение для матрицы ковариации cov имеет вид T T (41) cov = F b b F, где b - вектор ошибок, связанный с неточностями во входных данных.

Размер матрицы ковариации такой же, как и размер матрицы разрешения.

Ниже в 6.1.2 и в 6.2.2 построены карты стандартной ошибки полученного решения dVp/Vp, Vp/Vs. Под стандартной ошибкой в данном случае понимается квадратный корень диагональных элементов матрицы ковариации. При расчете стандартной ошибки точность определения времени вступления для всех лучей считается равной 0.1 сек. Т.о., в общем виде вычисляется по формуле:

T (42) = 0.1 ( FF ), Из формул (32) и (10) видно, что первые M столбцов матрицы A` образуют разреженную матрицу A, а последние 4*K столбца образуют блочно-диагональную матрицу, каждый блок которой представляет собой матрицу Hn.

K A' x = A + hn,, Hn (43) n= Hn Каждая матрица состоит из 4 столбцов и Nk строк (Nk - количество лучей Hn от k-го землетрясения). Согласно [94] каждая матрица может быть Hn сведена к верхней треугольной путем домножения на некоторую Hn ортогональную матрицу Qn. Т.к. размер матрицы 4Nk, то, следовательно, после приведения ее к треугольному виду только 4 строки будут ненулевыми и Nk-4 будут нулевыми. Домножим первые N1 строк матрицы B на соответствующую матрицу Q1, следующие N2 строк на Q2 и.т.д. В результате получится некоторая матрица, у которой в N-4*K строках элементы, соответствующие переменным hn, равны нулю. С помощью остальных 4*K строк матрицы переменные hn определяются однозначно. Т.о. чтобы найти решение x, сначала необходимо обратить матрицу размером M(N 4*K) и найти вектор, а затем последовательно, уже зная, вычислить hn для всех землетрясений (для всех n). Процедура нахождения ортогональных матриц Qn и расчет матрицы далее в работе будет называться разделением переменных [38]. При использовании разделения переменных система (32) преобразуется в x = r, где T 4M (44), x = 1,K, M размер M(N-4*K), размер r N-4*K. Ковариационные матрицы данных и параметров модели, а также демфирующие матрицы для уравнения (44) вводятся аналогичным образом, также как и для (32). Все приведенные ниже формулы соответствуют случаю, когда разделение переменных не производилось, т.е. решалась система уравнений (32). Аналогичные формулы легко могут быть выведены и для системы уравнений (44), т.е. случая, когда разделение переменных имело место. Метод разделения переменных [38] запрограммирован в программах Sphypit90 [43] (подпрограмма buildg2.f и u2.f), Simulps (автор С. Thurber, подпрограмма parsep.for). Этот метод применялся при расчете скоростных моделей в [44;

7] и в других работах.

2.5.1 Нахождение решения уравнения (34) с помощью построения системы нормальных уравнений Как было замечено выше, чтобы найти решение системы (32) по норме L2, необходимо найти такой вектор x, который удовлетворяет условию (33).

Дифференцируя условие (33) по xk получаем систему уравнений, которая называется системой нормальных уравнений.

A`T A`x = A`T T.

(45) В системе (45) количество неизвестных равно количеству переменных и равно M+4K. Эта система имеет единственное решение:

x = A`+T (46).

A`+ = (A`T A`)-1 A`T Матрица A`+ имеет M+4K строк и N столбцов. Эта матрица называется псевдообратной матрицей или матрицей Мура-Пенроуза.

Для модифицированного условия (35), если использовать обозначения, введенные в (34), система нормальных уравнений имеет вид:

(BT B + DT D)z = BTb,.

(47) Псевдообратная матрица A`+ с использованием обозначений (34) равна - (48) A`+ = Cm(BT B + DT D) BT, а вектор решения x с использованием обозначений (34) и (48) (49) x = A`+ b, Из (38), (45), (34) и (48) следует, что матрица разрешения равна:

R = A`+ BCm1 (50).

Из (48)и (49) видно, что решение x можно вычислить, только зная обратную матрицу к BT B + DT D. Для обращения этой матрицы возможно использовать различные методы. В рамках данной работы применялось разложение Холесского. Суть этого разложения в следующем.

Какой бы ни была матрица B (34) в любом случае матрица BT B + DT D является квадратной симметричной и неотрицательно определенной. Для любой квадратной симметричной и неотрицательно определенной матрицы P согласно [94;

87] существует верхняя треугольная матрица U, такая что верно равенство T U U = P.

(51) Если, например, представить матрицу BT B + DT D из (47) в виде (51):

T BT B + DT D = U U, (52) то вектор z из (47) можно получить, решая две треугольные системы:

T U y = BTb.

(53) Uz = y Матрицу разрешения (50) можно посчитать с помощью решения M+4K уравнений таким двухэтапным методом.

Если обозначить элементы матрицы U uij, а элементы матрицы P из (51) pij, тогда практически uij вычисляются по формуле из [94]:

i- pij - ukj uki k= (54) uij =, j = i +1,K,n, i = 1,K,n.

uii В формуле (54) сумма считается равной нулю при i=1.

При использовании этого метода точное решение, удовлетворяющее условию (35), находится за конечное число шагов.

2.5.2 Нахождение решения уравнения (34) с помощью сингулярного разложения матрицы B В данной работе этот метод обращения матрицы использовался только для случая, когда D = D1 = I,.

(55) где I - единичная матрица. Из курса матричного анализа известна следующая теорема [76]:

Для любой RT (R>T) матрицы S существует ортогональная RR матрица U, ортогональная TT матрица V и RT-диагональная матрица с диагональными элементами 12ЕЕ0, такие что выполняется равенство (56).

T S = UV, (56) Ортогональной называется такая матрица, у которой транспонированная и обратная матрицы равны между собой. Элементы матрицы 1, 2 и.т.д.

называются сингулярными числами. Если известно сингулярное разложение матрицы S и S - матрица полного ранга, найти обратную просто по формуле:

-1 T S = V-1U, (57) Т.о. представляя матрицу B в виде (55) и используя (56) уравнение (47) можно переписать в следующем виде:

T T (2 +2I )V z = U b.

(58) Из (58) легко вывести формулу для вычисления z, а также с учетом (38) формулу для вычисления матрицы разрешения.

T z = V U b +.

(59) T R = V V + В работе для сингулярного разложения матрицы B использовался программный код R.C. Singleton, созданный в 1973 г (Simulps, подпрограмма fksvd.for). В указанной подпрограмме матрица B сначала преобразуется к верхней двудиагональной посредством последовательных умножений B справа на специально вычисленные матрицы. А затем производится сингулярное разложение полученной двудиагональной матрицы. [94] Также как и для предыдущего способа обращения в данном случае точное решение, удовлетворяющее условию (35) находится за конечное число шагов.

2.5.3 Нахождение решения уравнения (34) методом LSQR Этот способ решения системы (34) также используется только для случая, когда демфирующая матрица D представима в виде (55): Данный метод решения уравнения (34) является итерационным, т.е. решение ищется с помощью серии последовательных итераций. Итерации продолжаются до тех пор, пока норма L2 не будет меньше некоторого числа. Следовательно, решение, найденное с помощью этого метода является приближенным в отличие от двух предыдущих методов. В описываемом алгоритме в явном виде псевдообратная матрица A`+ не вычисляется, а сразу находится вектор решения. Т.о. невозможно вычислить матрицу разрешения R и матрицу ковариации cov, как это можно было сделать при использовании двух предыдущих методов. Этот недостаток описываемого метода хорошо известен и существуют работы [35], в которых предлагаются различные алгоритмы для приближенного вычисления матрицы разрешения. С другой стороны очень эффективным является использование метода LSQR для разреженных переопределенных матриц больших размеров (более неизвестных), т.к. скорость работы этого алгоритма при обращении таких матриц существенно выше, чем скорость работы двух предыдущих алгоритмов, а необходимый для хранения массивов объем памяти значительно меньше. Для получения решения этим методом используется программа, написанная C.C. Paige и M.A. Saunders в 1982 году. В этой программе согласно [53;

36] для нахождения решения на каждой итерации сначала необходимо вычислить некоторую систему базисных векторов.

Затем на p-ой итерации найти решение по норме L2 простейшей трехдиагональной системы из p+1 уравнений с p неизвестными, а затем провести еще одну процедуру умножения матрицы на вектор.

Табл. 3 Степень участия автора в развитии трех различных методов обращения матриц (И - изучен;

С - используется при расчетах реальных скоростных моделей (может быть использован и не в программной реализации автора);

Т - использовался при тестировании.) Название метода обращения матрицы Степень участия автора Нормальных уравнений ИС Сингулярного разложения ИСТ LSQR ИСТ Теоретические основы решения обратных задач в геофизике вообще изложены в работе [121].

2.6 Методы локации землетрясений в одномерных скоростных моделях При определении координат вводится упрощающее предположение, что землетрясение - это точечный источник возбуждения поперечных и продольных сейсмических волн. Выше в 2.1 и далее в работе землетрясение рассматривается именно как точечный источник возбуждения сейсмических волн.

Задачи определения параметров гипоцентра и одномерной или трехмерной модели скоростей сейсмических волн в Земле, как видно из 2.1, являются взаимосвязанными и взаимозависимыми. Если неизвестны скорости сейсмических волн в исследуемом регионе, то нельзя с высокой точностью определить координаты землетрясения. Наоборот, если координаты землетрясений определены достаточно грубо, то траектории сейсмических лучей являются неверными и полученная скоростная модель абсолютно не соответствует реально существующей. Далее речь пойдет о методах, применяемых для локации событий в одномерной и трехмерной скоростных моделях. Более детально будут освещены методы, используемые автором в процессе работы по локации событий.

1. Самый простой и грубый способ определения эпицентрального расстояния до землетрясения по данным времен пробега волны Р и S, зарегистрированным на одной станции, изложен в [99]. Он заключается в вычислении эпицентрального расстояния по формуле:

= 1.73*(tP - tS ), (60) где tP и tS - это время пробега волны Pи S соответственно. Если станция является трехкомпонентной, то азимут на событие возможно определить по отношению смещений продольной волны, зарегистрированной двумя перпендикулярными горизонтальными сейсмографами. Эпицентральное расстояние и азимут однозначно определяют точку, которая является гипоцентром данного землетрясения.

2. Большинство программ локации событий (Hypo71, Hypocenter и др.) основаны на использовании метода Гейгера, предложенного в 1910 году.

Например, в программе Hypoellipse сначала приближенные координаты событий и время в источнике определяются в однородной скоростной модели (скорость во всем пространстве предполагается равной V0) с помощью системы из К нелинейных уравнений вида:

(x - x1)2 + (y - y1)2 + (z - z1)2 = (t1 - t)2 V L, (61) (x - xK )2 + (y - yK )2 + (z - zK )2 = (tK - t)2 V где К - количество измерений одного события на разных станциях, x,y,z,t - неизвестные координаты события и время в источнике, xi,yi,zi,ti - координаты станции, регистрирующие i-ую волну и время вступления i-ой волны на станции (могут быть использованы различные типы волн). Эта система нелинейных уравнений может быть сведена к линейной путем последовательного вычитания всех нижних уравнений из первого. Далее поправки к приближенно определенным координатам находятся с помощью системы уравнений (62).

T = Hh, где T T = T1,K, TN k T h =, x, y, z, T1 T1 T, (62) x y z H = 1..., TN TN TN k k k x y z где Nk - количество лучей от одного события на разные станции, Ti -.разница между реальным и расчетным временем пробега для i-го луча,, x, y, z - поправки к приближенно определенным координатам события и времени в источнике. Системы уравнений (61) и (62) чаще всего являются переопределенным, поэтому решение ищется по норме L2 (33) или какой нибудь другой норме. Если число измерений меньше 4-х, то неизвестные Ti Ti Ti параметры найти невозможно. Для нахождения производных,, x y z необходимо знать траектории сейсмических лучей в заданной слоистой скоростной модели. Эти производные связаны с углом выхода луча и скоростью в источнике следующим соотношением:

Ti cosi Ti cos i Ti cos i =, =, =,, (63) x Vист y Vист z Vист где i, i, i - углы между касательной к i-ому лучу в источнике и осью x, y, z соответственно, Vист - скорость сейсмической волны в источнике.

Программы, вычисляющие координаты и время в источнике по методу Гейгера, отличаются друг от друга введением дополнительных различных весовых матриц и весовых коэффициентов, различными способами обращения матриц. В ходе работы автор использовал программы Hypoellipse [26], Sphrel3d [43]. В Simulps возможно проведение локации с использованием метода Гейгера в трехмерной скоростной модели.

3. Также для локации событий применялся метод сеточного поиска. С развитием компьютерной техники, а также с появлением конечно разностных методов для вычисления времени пробега от источника до различных точек исследуемой области стало возможным не обращать матрицу H из (62), а искать абсолютный минимум обычным методом перебора. Метод перебора является более простым и надежным по сравнению с численным решением переопределенной системы линейных уравнений. Автором работы был запрограммирован алгоритм, предложенный в [34]. Локация событий производится как в одномерной, так и в трехмерной скоростной модели. Вычисленные функции RMSсоб или RMSwсоб после проведения локации этим методом действительно имеют меньшие значения, чем после проведения локации методом 1.

4. Также заслуживает внимания метод двойных разностей для определения локации событий, т.к. в некоторых случаях он является более эффективным, чем описанные выше методы. Он был предложен в 2000 году F. Waldhauser и W.Ellsworth [62]. Метод основан на решении линейной системы, каждое уравнение в которой имеет вид:

Tki Tki Tki Tkj Tkj Tkj i j j j j xi + yi + zi + - x - y - z - = drkij x y z x y z, (64) i i i j где drkij = (tk - tkj ) - (tk - tkj ) = T - T реальн рассчетн i - номер одного события, j - номер другого события, k - порядковый номер уравнения, drijk - двойная разность, остальные обозначения аналогичны (62).

Т.о. в предыдущих методах производился поиск таких значений x,y,z и t, которые уменьшают максимальным образом функцию RMSсоб или RMSwсоб, в этом методе предлагается минимизировать функцию всех возможных комбинаций двойных разностей. Этот метод позволяет производить более точную локацию событий, локализованных в одной области (кластере). Чаще всего наибольшее скопление землетрясений наблюдается в зонах разломов. В [62] методом двойных разностей произведена более точная локация событий, произошедших в районе Северного Хауварского (Hayward) разлома в Калифорнии. Афтершоки и форшоки [99] какого-либо крупного землетрясения также образуют кластер. В работе [47] произведена перелокация событий методом двойных разностей для землетрясения, произошедшего 26 июля 2001 года на острове Скирос (Греция). При использовании метода двойных разностей контуры разломных зон или областей, в которых произошла разрядка накопленного за некоторый период геологического времени напряжения, удается выявить более точно. Заметим, что абсолютная точность локации всей этой области или группы событий по прежнему зависит от количества станций, регистрирующих его.

Преимущества этого метода перед более ранними в случаях применения в областях наибольшего скопления событий заключаются также в том, что 1. исчезает ошибка, связанная с неоднородной скоростной структурой непосредственно под станцией 2. локация не зависит от скоростных неоднородностей на пути от станции к паре событий, а зависит только от скоростных неоднородностей непосредственно вблизи перелоцируемых событий.

Авторами статьи [62] также была написана программа HypoDD для локации событий с использованием предложенного ими метода. В ходе данного исследования программа HypoDD была усовершенствована для случая локации событий в трехмерной среде с помощью построения лучей в скоростной модели с горизонтальными неоднородностями. Анализ сейсмичности территории Тянь-Шаня показал, что в этом районе существуют зоны, в которых плотность событий в несколько раз больше, чем во всем остальном пространстве. Для некоторых из этих зон в настоящей работе была произведена попытка локации событий методом двойных разностей. Ниже в Табл. 4 представлена степень участия автора в развитии трех различных методов локации событий.

Табл. 4 Степень участия автора в развитии трех различных методов локации событий (И - изучен;

П - запрограммирован автором;

С - используется при расчетах реальных скоростных моделей (может быть использован и не в программной реализации автора);

У - усовершенствован автором;

Т - использовался при тестировании.) Название метода локации событий Степень участия автора Метод Гейгера ИС Метод сеточного поиска ИПС Метод двойных разностей ИУТ 2.7 Различные способы введения весовых коэффициентов Введение весовых коэффициентов является неотъемлемой частью любой томографической программы. Каждой станции может быть присвоен свой вес в зависимости от качества регистрации вступлений на этой станции, каждому событию может быть присвоен также свой весовой коэффициент в зависимости от качества вступлений сейсмических волн, пришедших только от этого события. Более распространено присвоение собственного веса каждому лучу в зависимости от того, насколько точно определено время вступления волны. Обозначим параметр, характеризующий точность определения времени вступления волны t. Наиболее часто используемым способом взвешивания является разделение всех лучей на весовые группы.

Например, если t0.07сек, то целочисленный вес этого луча wt=0, если 0.07

число, fw =0 =1, fw =1 = 5, fw =2 =10,. и.т.д.;

в программе Simulps:

t t t t 2w Введение целочисленных весов удобно, если не известно конкретное значение точности определения вступления волны, а известны лишь качество вступления (I или E). В случае данной работы дана точность определения времени вступления волны в секундах. Автором предлагается своя формула для вычисления W1i, зная t:

W1i =. (65) t Значение W1i непрерывным образом зависит от t. Для случая, когда известна точность определения времени вступления сейсмических волн, дискретный способ взвешивания (деление лучей на различные весовые группы) менее эффективен, чем непрерывный способ взвешивания. Использование дискретного способа взвешивания в данной ситуации может привести к искажению конечного результата, т.к. два времени вступления, имеющие почти одинаковую точность вступления, например 0.34 и 0.36 сек, могут оказаться в разных весовых группах. Таких проблем в случае применения непрерывного способа взвешивания не возникает. В дальнейшем в работе при использовании стандартных программ для расчетов (Hypoellipse, Velest, Fatomo) деление всех данных на весовые группы производилось, а при расчете трехмерной модели по собственным программам весовой коэффициент №1 определялся по формуле (65).

Далее при расчетах согласно (19) в программах используется итоговый весовой коэффициент, который зависит не только от точности определения вступления, но и от эпицентрального расстояния (W3i) и невязки (W2i). В Sphypit90 и Simulps W2i и W3i вычисляются по-разному. Общая закономерность в такая, чем больше величина невязки и больше эпицентральное расстояние, тем W2i и W3i соответственно меньше.

Зависимость W3i от расстояния и W2i от невязки легко представить графически (Рис. 4).

Рис. 4 а) Зависимость W2i от значения невязки, б) Зависимость W3i от эпицентрального расстояния для программ Sphypit90 и Simulps.

2.8 Основные этапы любого сейсмотомографического исследования В этом разделе выделено четыре этапа работы, необходимых для получения достоверных трехмерных скоростных моделей.

2.8.1 Сортировка данных Времена вступления сейсмических волн всегда определяются с некоторой погрешностью. Для аналоговых станций считается, что максимальная точность определения сейсмической волны 0.1 сек, а для цифровых 0.01 сек. Далеко не для всех сейсмических волн удается определить времена вступления с указанной точностью. От одних землетрясений сейсмическая запись является более четкой, от других - менее. Точность определения времени вступления волны также зависит от местоположения сейсмической станции. На станциях, которые расположены вдали от населенных пунктов, уровень сейсмического шума меньше и полезный сигнал более четкий. Частой ошибкой является неправильное определение фазы волны. Интерпретатор пытается зарегистрировать время вступления первой волны, а регистрирует время вступления вторичной фазы.

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

Анализ скоростных свойств в исследуемом регионе, построение одномерных и трехмерных скоростных моделей следует проводить только, опираясь на данные времен вступлений, определенные с хорошей точностью и только от тех землетрясений, для которых ошибка вычисления координат и времени в источнике невелики. Использование всех имеющихся данных в томографическом исследовании может привести к получению неверных одномерных и трехмерных скоростных моделей. Наиболее грубые ошибки в определении времени вступления удается зафиксировать при построении индивидуальных и сводных годографов. Чтобы из общей массы землетрясений выявить те события, координаты и время в источнике которых определены достаточно надежно, Э.Кислингом предлагается использовать следующие критерии.

1. Выбирать только те события, которые зарегистрированы достаточно большим количеством станций.

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

По мнению С.Рокера [1;

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

1 - число обусловленности матрицы из частных производных времени пробега по координатам, т.е. отношение максимального собственного числа матрицы к минимальному, не должно быть слишком большим;

2 - процесс локации должен быть сходящимся, т.е., например, на последней итерации координаты должны сдвигаться не более, чем на 3 км;

3 - координаты событий должны быть слабо чувствительны к локации в различных скоростных моделях. Т.е., например, координаты событий не должны отличаться более, чем на 10 км при локации в скоростных моделях, различающихся на 10%;

4 - стандартная ошибка определения координат также не должна быть слишком большой.

Конкретные числа в критериях для каждого набора данных необходимо подбирать свои в зависимости от специфики данных и региона, также с использованием метода проб и ошибок.

И Э.Кисслингом, и С.Рокером также отмечается, что при сортировке данных необходимо обратить внимание на значение RMSсоб и на расстояние до ближайшей станции.

Критерии Э.Кисслинга и С.Рокера применены для сортировки событий во многих томографических исследованиях [45;

40]. Результаты этих исследований подтверждают действенность описанных критериев при сортировки событий. По этой причине автором в данной работе также использовались эти критерии для отбора качественных данных.

2.8.2 Построение оптимальной одномерной модели и перелокация событий в полученной одномерной модели. Вычисление станционных поправок.

Важность нахождения минимальной одномерной скоростной модели (МОСМ) подчеркивалась во многих работах, например в [21]. Под МОСМ понимается такая скоростная модель, относительно которой среднеквадратичная невязка по всем лучам минимальна, среди всех возможных скоростных моделей, т.е. модель обеспечивающая минимум функционала RMSw (19). Использование для трехмерных расчетов в качестве нулевого приближения неоптимальной одномерной модели может привести к появлению артефактов (несуществующих неоднородностей) в трехмерной модели. Подбор оптимальной одномерной модели должен производиться для каждого набора данных отдельно с учетом существующей априорной информации о скоростных границах и значениях скоростей в интересующем регионе. Поэтому поиск оптимальной одномерной скоростной модели является неотъемлемой частью любого сейсмотомографического исследования. Процесс нахождения наилучшей одномерной скоростной модели, как и трехмерной, состоит из решения прямой и обратной задачи.

Следовательно, почти любая сейсмотомографическая программа, предназначенная для расчета оптимальной трехмерной модели, пригодна для поиска оптимальной одномерной модели, если иначе произвести параметризацию среды. Процесс поиска оптимальной одномерной модели является более неустойчивым, чем трехмерной. Наилучшей сходимостью при построении оптимальной одномерной модели по мнению автора обладает программа Velest (автор E.Kissling). Ранее для построения скоростного разреза, например по данным ГСЗ или просто по данным площадных измерений, использовалась формула Герглотца-Вихерта, а также формулы геометрической оптики.

Построение оптимальных скоростных моделей в работе производилось с учетом рекомендаций, данных в [21]. Достаточно детально процесс построения и тестирования оптимальной одномерной модели описан также в диссертации S.Husen [19]. Часть тестов, использованные S.Husen, повторены автором для своего набора данных.

Основными критериями оптимальности одномерной модели считается близость (взвешенного) среднего арифметического по всем невязкам к нулю и минимальность (взвешенной) среднеквадратичной невязки для всех лучей RMSw.

Одновременно со скоростями обычно рассчитываются станционные временные поправки для каждой станции. Станционные поправки в центре сети обычно несут в себе информацию о геологическом строении приповерхностных слоев. Временные поправки для станций, расположенных на окраинах, отражают не только особенности структуры осадочного чехла в зоне приемника, но и другие эффекты, связанные с неравномерным расположением событий в зависимости от азимута для данной станции.

2.8.3 Выбор алгоритма и его тестирование.

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

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

Эта модель называется нулевым приближением. Затем по рассчитанным временам пробега пытаются восстановить скоростные неоднородности с использованием выбранного алгоритма.

2.8.4 Построение трехмерной модели и перелокация событий, проверка правильности проведенных расчетов.

Когда алгоритм достаточно оттестирован, данные качественно отсортированы и найдена минимальная скоростная модель и станционные поправки, а также проведена перелокация событий в этой одномерной модели, переходят к построению трехмерных скоростных моделей по реальному набору данных.

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

2.8.4.1 Критерии:

1. Уменьшение взвешенной среднеквадратичной невязки RMSw (19). Рассчитанная в полученной минимальной одномерной или трехмерной скоростной модели среднеквадратичная невязка RMSw (19) должна быть меньше невязки, рассчитанной в начальной одномерной скоростной модели.

2. Диагональные элементы матрицы разрешения (38) должны быть больше 0.5 для большинства искомых переменных, а рассчитанная стандартная ошибка в несколько раз меньше вычисленных поправок.

2.8.4.2 Оценки:

1. Оценка устойчивости полученного результата. Из имеющегося списка землетрясений выбирают только четные события и проводят расчеты только с этой группой событий. Затем выбирают только нечетные события и снова проводят расчеты. Два расчета, проведенные по разным наборам данных, должны быть схожи. При построении трехмерных томографических моделей эта оценка была выполнена, например, в [93]. В данной работе с помощью описанной процедуры проверялась устойчивость полученной одномерной скоростной модели Р волн для территории Северного Тянь Шаня.

2. Оценка разрешающей способности имеющегося набора данных.

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

Доказательство возможности восстановления полученной трехмерной скоростной структуры (restoring resolution test). Предположим, что после проведения расчетов с использованием реальных данных, получена некоторая трехмерная скоростная модель Мреал. Чтобы проверить, возможно ли восстановление именно этой скоростной модели с помощью данного набора данных, рассчитывают времена пробега для всех возможных пар источник-приемник в имеющейся трехмерной скоростной модели (trrt) Мреал.

Затем, возвращаясь обратно к отсчетной скоростной модели, проводится новая серия вычислений, в которых реальные значения времен пробега заменены временами trrt. Если, полученная в результате этих расчетов новая скоростная модель сходна с моделью Мреал, значит, имеющийся набор данных способен восстановить модель Мреал.

3. Если на исследуемой территории производились взрывы, то оценить качество полученной трехмерной модели можно с помощью проведения локации этих событий в рассчитанных скоростных моделях.

Координаты места проведения взрыва и время взрыва обычно точно известно. В то же время определить координаты этого события и времени в источнике можно сейсмологическими методами. В трехмерных скоростных моделях, которые действительно отражают истинный характер скоростных неоднородностей в Земле точность локации должна быть выше, чем в различных одномерных моделях.

Выводы к главе Все сейсмотомографические алгоритмы могут отличаться друг от друга:

1. Методами параметризации среды;

2. Методами построения лучей;

3. Методами обращения больших разреженных матриц;

4. Методами локации событий;

5. Способом взвешивания данных.

Автором проанализировано 4 различных метода параметризации модели, 6 различных методов построения лучей (М1-М4, М6-М8), 3 различных способа обращения больших разреженных матриц и 3 различных метода локации событий. Предложено использование не дискретного, а непрерывного способа взвешивания данных.

Глава 3 Сравнительный и оценочный анализ В предыдущей главе показано, насколько широк спектр вычислительных алгоритмов, применяемых в рамках метода региональной сейсмичской томографии, по данным от локальных землетрясений. В большинстве работ, посвященных сейсмотомографии, обычно приводятся результаты, полученные с помощью одного какого-либо алгоритма. Между тем, выявление случаев, когда один метод построения лучей является более точным, чем другой, или один способ параметризации среды более адекватно описывает среду, чем другой имеет большое значение при построении реальных томографических моделей. В данной главе приводятся результаты сравнения различных методов построения лучей, различных способов параметризации среды, различных способов обращения матриц на тестовых примерах. Автору не известны работы, в которых проводились большинство из описываемых тестов. Также важным и необходимым по мнению автора является оценка влияния сферичности Земли на различных и оценка влияния пересчета лучей на значение RMSw, представленные в этой главе.

3.1 Сравнение различных численных алгоритмов для расчета времени пробега луча и вычисления траекторий сейсмических лучей.

В процессе работы проводилось сравнение описанных в 2.4 методов построения лучей на тестовых примерах. Сравнение было выполнено для того, чтобы выбрать наиболее точный и несложный в использовании метод построения лучей.

В статьях [37;

28;

17] сравнивались различные методы построения лучей. Тесты, проведенные в [37], демонстрируют преимущество вычисления трехмерных траекторий сейсмических лучей (методы 3-8) по сравнению с плоскими траекториями (методы 1-2). Если на каждом шаге итеративного процесса (8) производить заново пересчет трехмерной траектории сейсмического луча, восстановленные контуры скоростных неоднородностей будут ближе к истинным контурам, чем в случае расчета плоских траекторий лучей. Сравнение М3_2 (см. 2.4.2.2) с М6 (см. 2.4.2.3) проводилось в работе [17]. Э.Кислингом и Хаслингером был сделан вывод, что метод изгиба траектории (М3_2) обеспечивает высокую точность только для небольших расстояний, до 80 км. Для больших расстояний лучше использовать параксиальный метод построения лучей (М6) вместо метода изгиба траектории (М3_2).

Первый тест, проведенный автором, выявляет недостатки методов 2 и по сравнению с 7. Для простоты представления тест рассматривается не в 3-х мерном, а в 2-х мерном пространстве. Задана некоторой скоростная модель, содержащая низкоскоростную неоднородность (Рис. 5). Заданы координаты источника и приемника. Вычисляется траектория сейсмического луча и время пробега вдоль этого луча. Анализ полученных результатов (Рис. 5) показывает, что наименьшее время пробега соответствует траектории, рассчитанной с помощью метода M7 (см. 2.4.2.4). Следовательно, эту траекторию стоит считать истинной.

Рис. 5 Сравнение четырех методов построения лучей в скоростной модели, содержащей низкоскоростное включение.

Наихудший результат получается при использовании метода M3_1, когда в качестве нулевого приближения используется прямая линия между источником и приемником. Если же в качестве нулевого приближения использовать сейсмический луч, рассчитанный с помощью M2, а затем производить уточнение траектории с помощью методики, описанной в [59] (метода М3_2), результирующая кривая ближе к истинной, но не совпадает с ней полностью. Разница между временами пробега, вычисленными вдоль траекторий М7 и М3_2 1.3%. При использовании метода М7 вся исследуемая площадь была разбита на квадраты со стороной 0.5 км. Лучи под различными углами выходили из источника до тех пор, пока луч не попал в ячейку, содержащую приемник. Когда производился расчет с использованием М3_2, вся траектория состояла из 129 точек, каждая из которых (за исключением конечных) последовательно возмущалась. Итеративный процесс подбора луча продолжался до тех пор, пока разница времен пробега, рассчитанных на двух последовательных итерациях, не стала меньше некоторого заданного tlim, в данном случае равном 0.001 сек. Этот тест показывает, что для данной скоростной модели метод М7 более точно вычисляет траекторию сейсмического луча и время пробега, чем методы М2 и М3. Другой итоговый результат заключается в том, что использование в качестве нулевого приближения для метода М3 прямой линии в данном случае является значительно менее точным, чем использование в качестве нулевого приближения траектории, рассчитанной с помощью M2.

Второй тест выявляет особенности, свойственные различным алгоритмам расчета лучей, на различных эпицентральных расстояниях.

Предположим, что горизонтально-неоднородная тестовая скоростная модель задана формулой x y v = v0 (1 + 0,06 sin( ) sin( )), (66) 35 где v0=5,8. Также предположим, что координаты источников и приемников распределены в объеме так, как представлено на Рис. 6.

Рис. 6 Карта гипоцентров землетрясений и координаты станций (слева), зависимость глубины землетрясения от координаты y (справа).

Произведем расчет времен пробега лучей для всех пар источник приемник с использованием метода изгиба траектории М3_2, параксиального метода построения луча (М6) и с помощью конечно-разностного метода с шагом 2 км (М7(2км)) и с шагом 0.5 км (М7(0.5км)). Скорость вычислений максимальна при использовании метода М3_2, немного медленнее расчеты производятся методами М6 и М7(2км), наибольшее количество вычислений требуется при использовании метода М7(0.5км). Результаты расчетов представлены на Рис. 7.

Рис. 7 Зависимость разности времен пробега, рассчитанных с помощью различных методов, от эпицентрального расстояния для тестовой модели, заданной формулой (66) и координат источников и приемников, представленных на Рис. 6.

Из Рис. 7 видно, что наименьшие времена пробега рассчитываются с использованием метода М6. Важным является тот факт, что при расчетах с использованием метода М3_2 после 100 км образуется группа лучей, для которой времена пробега более, чем на 0.1 сек, отличаются от времен пробега, рассчитанных с помощью М6, М7 (2 км) или М7 (0.5 км). Одно из возможных объяснений этого факта состоит в том, что с помощью метода М3_2 для этой группы лучей найдена неоптимальная траектория сейсмических лучей. Следовательно, данный тест показывает, что метод М3_2 менее точен при эпицентральных расстояниях более 100 км, чем методы М6, М7 (2км) и М7 (0.5км). Разницы времен пробега tМ6- tМ7(0.5км) и tМ7(2км)- tМ7(0.5км) для большинства лучей с эпицентральными расстояниями 0 200 км находятся в пределах 0.1 сек. Т.к. скорость расчета времени пробега с использованием метода конечных разностей М7 существенно зависит от размера ячейки, на которые разбивается исследуемый объем, представляет интерес оценить насколько более точно производятся вычисления времен пробега при использовании более мелких ячеек. Из Рис. 7в видно, что при разбиении на более мелкие ячейки вычисленные времена пробега всегда имеют меньшее время, чем при разбиении на более крупные. При уменьшении стороны каждой ячейки-куба в 4 раза максимальное отличие расчетных времен пробега равно 0.13 сек. Следовательно, когда в процессе томографической инверсии разница между реальным и расчетным временем пробега уменьшилась до значений 0.10-0.15 сек, вычисление времени пробега методом конечных разностей с помощью разбиения на ячейки-кубы со стороной 2 км производить бессмысленно, т.к точность вычисления времени пробега величина того же порядка, что и невязка. Необходимо продолжить расчеты с использованием более мелкого разбиения на блоки.

В третьем тесте расчитывалось время пробега в одномерной слоистой скоростной модели с использованием плоскопараллельного преобразования (25), учитывающего кривизну Земли, и без использования. Для проведения теста было взято некоторое реальное расположение источников и приемников. На Рис. 8 представлены результаты теста. По оси ординат откладывается значение разности времен пробега, рассчитанной для одной и той же пары источник-приемник с использованием плоскопараллельного преобразования и без. По оси абсцисс откладывается эпицентральное расстояние между источником и приемником. Все лучи разбиты на группы в зависимости от того, является ли луч прямым или преломленным на какой либо границе.

Рис. 8 Зависимость разности времен пробега Р волн, рассчитанных с использованием плоскопараллельного преобразования и без в слоистой скоростной модели, представленной на врезке, от эпицентрального расстояния.

Из Рис. 8 видно, что время пробега, рассчитанное с использованием плоскопараллельного преобразования, всегда меньше времени пробега, рассчитанного без использования преобразования. На эпицентральных расстояниях 200 км разброс поправок, связанных со сферичностью Земли, для лучей, принадлежащих различным группам, наибольший. Так поправка для лучей, преломленных на границе 50 км, на 0.1 сек больше поправки для лучей, преломленных на границе 10 км. При эпицентральных расстояниях ~600 км разность tсф-tпл по модулю составляет около 1% от значения времени пробега и равна 0.6 сек. Величины поправок на эпицентральных расстояниях более 200 км уже сравнимы со значениями невязок, рассчитанными для сейсмических лучей с данным эпицентральным расстоянием в некоторой одномерной скоростной модели. Отсюда следует, что отсутствие корректировки, связанной со сферичностью Земли, в исследованиях, в которых расстояния между источником и приемником для большинства лучей больше 200 км, может привести к существенным искажениям при построении трехмерных скоростных моделей. В программах Simulps и Fatomo учет кривизны не производится, поэтому эти программы не желательно использовать для расчета трехмерных скоростных моделей областей более 200 км в поперечнике. В алгоритме Sphypit90, напротив, осуществляется коррекция глубины и скорости согласно (25), поэтому Sphypit90 подходит для получения трехмерных скоростных моделей регионального масштаба (более 200 км в поперечнике).

Аналогичный тест проводился также с использованием одномерной непрерывной скоростной модели. Полученные результаты имеют похожий вид.

3.2 Сравнение различных способов параметризации модели В данном разделе описан тест, демонстрирующий преимущества алгоритмов, использующих непрерывные способы параметризации модели, по сравнению с алгоритмами, использующими блоковый способ параметризации.

В качестве тестовой скоростной модели использовалась скоростная модель, заданная формулой (66) при v0=5.8 км/с. Горизонтальное сечение тестовой скоростной модели для всех глубин представлено на Рис. 9. В тесте используется расположение источников и приемников, изображенное на Рис.

6. Для всех пар источник-приемник в тестовой скоростной модели проведен расчет времен пробега. В качестве нулевого приближения была задана модель с постоянной скоростью во всем объеме, равной 5.8 км/с (Рис. 9).

Далее с использованием различных способов параметризации среды проведено восстановление истинной скоростной структуру. Решение обратной задачи проводилось с помощью разделения переменных [38].

Сейсмический шум (ошибка, связанная с неправильным определением времени вступления), а также ошибка, связанные с неправильным определение координат источника, в данном случае не вводились.

Как видно из Рис. 10 границы высоко- и низкоскоростных неоднородностей в тестовой модели расположены на расстоянии 35 км.

Восстановление блочной скоростной модели проведено для случаев, когда весь исследуемый объем разбивался на блоки с горизонтальными размерами 2020 км и 1010 км. При использовании Би-сплайн параметризации узлы сетки не совпадали с центрами высоко- и низкоскоростных клеток тестовой модели, а располагались равномерно по осям x и y через каждые 20 км. При восстановлении скоростной структуры с помощью триангуляции Делоне использовалась неравномерная сетка, наибольшее количество узлов было сосредоточено в зонах максимальной плотности сейсмических лучей.

Для того, чтобы произвести сравнение влияния именно типа параметризации на Рис. 10 и Рис. 11 приведены результаты расчетов после первой итерации. В этом случае траектории лучей для всех четырех вариантов одинаковы и совпадают с прямыми линиями, а следовательно, решение не зависит от способа построения луча. Процесс восстановления трехмерной скоростной структуры является итеративным, но характер основных неоднородностей обычно проявляется уже после первой итерации.

Рис. 9 Горизонтальные сечения модели скоростей Р волн для всех глубин для тестовой скоростной модели (слева) и для нулевого приближения (справа).

Рис. 10 Горизонтальные сечения модели скоростей Р волн на глубине 3 км восстановленных моделей при триангуляции Делоне, параметризации линейными и кубическими Би-сплайнами. На модели нанесены контуры истинных скоростных неоднородностей, представленных на Рис. 9 (слева).

Рис. 11 Горизонтальные сечения модели скоростей Р волн на глубине 3 км при использовании блоковой параметризации среды:

а1) модель, рассчитанная при использовании разбиения на блоки 20 на 20 км, до сглаживания, а2) модель, рассчитанная при использовании разбиения на блоки 20 на 20 км после сглаживания с помощью линейных Би-сплайнов, б1) модель, рассчитанная при использовании разбиения на блоки 10 на 10 км, до сглаживания, б2) модель, рассчитанная при использовании разбиения на блоки 10 на 10 км после сглаживания с помощью линейных Би-сплайнов.

На сглаженные модели нанесены контуры истинных скоростных неоднородностей, представленных на Рис. 9 (слева).

Для восстановления скоростной структуры с использование блоковой параметризации среды использовалась программа Sphypit90, с использованием линейной Би-сплайн параметризации - программа Simulps, с использованием кубической Би-сплайн параметризации и тетраэдеров Делоне - программы, написанные автором работы. На Рис. 11 результаты расчетов с использованием блоковой параметризации представлены в двух разных вариантах. Оба эти варианта используются в сейсмической томографии для представления конечного результата. В некоторых работах результаты расчетов изображаются в виде а1, б1, а в других, например [44;

67], в виде а2, б2. На Рис. 11 а1 и б1 представлены скоростные модели, которые действительно были посчитаны. Скорость сейсмической волны внутри каждого блока постоянна. На Рис. 11 а2 и б2 представлены скоростные модели, полученные после сглаживания результатов, изображенных на Рис. 11 а1 и б1. Та скорость, которая соответствует данному блоку, присваивалась центральной точке этого блока. Затем, зная скорости в центрах всех блоков, производилась аппроксимация скоростной функции линейными Би-сплайнами по значениям скоростей в конкретных точках.

Сравнение результатов, представленных на Рис. 10 и Рис. 11, показывает, что для данного расположения источников и приемников и для описанного выше разбиения среды на блоки и задания точек сетки, параметризация линейными или кубическими Би-сплайнами, а также триангуляция Делоне лучше восстанавливают клеточную скоростную структуру, чем блоковая параметризация. На Рис. 11 а1 и б1 клеточная структура неоднородностей еле уловима. На Рис. 11 а2 и б2 она проявляется более четко, но все равно хуже, чем на Рис. 10. Для случая б на Рис. 11, когда горизонтальные размеры блоков 1010 км, рассчитанная скоростная модель является более детальной. Если сопоставить контуры истинных скоростных неоднородностей с рассчитанными, то обнаружится, что второй снизу ряд клеток абсолютно неверно восстанавливается при использовании блоковой параметризации (Рис. 10), в то время как при использовании Би-сплайн параметризации и триангуляции Делоне этот ряд клеток восстанавливается достаточно точно (Рис. 11).

Проверка разрешающей способности имеющегося набора данных с помощью клеточного теста при сейсмотомографическом исследовании также проводится при использовании блокового способа параметризации, например в [15]. Из сказанного выше, также как и из работы [15] следует, что клеточная структура в целом восстанавливается при блоковом способе параметризации модели. Однако, при использовании непрерывных способов параметризации в большинстве случаев удается добиться лучшего восстановления клеточной скоростной структуры. Особенно отчетливо разница в восстановлении клеточной структуры при использовании блоковых и непрерывных способов параметризации среды заметна для случаев систем наблюдения с малым количеством регистрирующих станций и неравномерного расположения источников и приемников.

3.3 Сравнение различных способов обращения матриц.

В работе проведено сравнение 2-х способов обращения переопределенной разреженной матрицы: методом сингулярного разложения и методом LSQR. Рассмотрена одна из матриц, которая образовалась в ходе проведения расчетов. Эта матрица состоит из 2163 строк. Количество столбцов в матрице равно 1128, из них 904 (226*4) соответствуют неизвестным параметрам гипоцентров, а 224 неизвестным параметрам модели. Около 90% элементов в исследуемой матрице нулевые, но при этом в каждом столбце, который соответствует неизвестным параметрам модели, содержится как минимум 500 ненулевых элементов. Для обращения такой матрицы с помощью LSQR потребовалось 250 итераций. При использовании метода сингулярного разложения было проведено вычисление собственных чисел, соответствующих данной матрице. Отношение максимального собственного числа к минимальному (число обусловленности) равно 52015. В обоих случаях использовался демфирующий параметр 0.5. Результаты расчетов, с использованием обоих методов, абсолютно идентичны. Среднеквадратичная невязка RMS понизилась со значения 0.19 сек до значения 0.11 сек и в том и в другом случае. Среднее арифметическое найденных параметров модели равно 0. км/с, среднее арифметрическое поправок к горизонтальным координатам гипоцентров 0,5 км, к вертикальной координате 0.3 км, ко времени в очаге 0.06 сек.

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

3.4 Оценка влияния процедуры пересчета лучей и перелокации событий после каждой итерации на значение RMSw.

Как отмечается в [91] подбор минимальной одномерном скоростной модели по данным времен пробега объемных сейсмических волн является задачей существенно неустойчивой, т.е. большие вариации скоростной модели могут приводить к малым изменениям значения RMSw. При подборе минимальных одномерных моделей в ходе данного исследования такого рода неустойчивость также неоднократно отмечались.

При построении минимальной трехмерной скоростной модели в одних областях (там где плотность лучей низкая) введение больших вариаций скорости не приводит к сильным изменениям значения RMSw после пересчета лучей. Напротив, в других областях (там где плотность лучей высокая) небольшие изменения скорости сильно изменяют значение RMSw, которое рассчитывается после пересчета лучей в новой трехмерной скоростной среде. Ниже на примере реального набора данных для территории Северного Тянь-Шаня вычислено значение RMSw в различных скоростных моделях тремя различными способами. Для анализа используются 5 скоростных моделей, полученных после первой итерации с использованием демфирующих параметров (см. 2.5) D1 0.5;

5;

10;

50 и 100 и одного и того же набора данных. Как известно, с помощью демфирующего параметра можно регулировать уровень рассчитанных скоростных вариаций, поэтому при введении большого демфера контрастность результирующей скоростной модели значительно меньше, чем при введении маленького. В Табл. 5 в строке 1 значение RMSw рассчитано на основе старых траекторий лучей, координат событий и времени в источнике. В строке 2 RMSw рассчитано после пересчета траекторий лучей в новой скоростной модели с использованием старых координат источников и приемников. В строке RMSw рассчитано после пересчета траекторий лучей и перелокации событий, проведенной методом сеточного поиска. Начальное значение RMSw равно 0.152 сек.

Табл. 5 Значения RMSw в сек, рассчитанные в трехмерных скоростных моделях, полученных в результате обращения одной и той же матрицы с различными демфирующими параметрами. Тр0, t0, x0, y0, z0 Цтраектории сейсмических лучей, время в источнике и координаты землетрясений в начальной скоростной модели. Тр1, t1, x1, y1, z1 Цтраектории сейсмических лучей, время в источнике и координаты землетрясений в рассчитанной скоростной модели.

d0.5 d5 d10 d50 d RMSw нач 0.152 0.152 0.152 0.152 0. тр0, t0, x0, y0, z0 0.111 0.131 0.138 0.149 0. тр1, t0, x0, y0, z0 0.204 0.139 0.141 0.145 0. тр1, t1, x1, y1, z1 0.145 0.133 0.137 0.147 0. dVpсред.ариф (км/с) 0.44 0.11 0.06 0.008 0. Из Табл. 5 видно, что при демфере 0.5 после обращения матрицы для траекторий Тр0 и параметров эпицентров t0, x0, y0, z0 рассчитанная RMSw минимальна по сравнению с другими рассматриваемыми моделями. Однако, w кон RMS после пересчета лучей в трехмерной скоростной модели RMSw становится максимальной по сравнению с другими четырьмя моделями и больше, чем RMSw нач. После проведения перелокации событий в новой трехмерной скоростной модели рассчитанное значение RMSw уменьшается и становится приблизительно равной взвешенной среднеквадратичной невязке для случая использования демфера 50, когда вычисленные поправки к значениям скоростей в несколько раз меньше. Поэтому на данной стадии выбора оптимального демфирующего параметра возникает задача - найти такие неизвестные параметры модели, которые несильно отличаются от нуля, но при этом обеспечивают минимальное значение невязки, рассчитанной обязательно после пересчета траекторий сейсмических лучей в новой скоростной модели с использованием новых параметров гипоцентров. Из табл.4 видно, что в данном представленном расчете оптимальным демфером является число 5.

Тема выбора оптимального демфирующего параметра затрагивается во многих сейсмотомографических работах, например в [8]. В отличие от выше представленного исследования в них не приводятся численные оценки для различных случаев, когда пересчет лучей не производится или пересчет производится, но перелокация не осуществляется. Eberhart-Phillips [8] предложен выбор оптимального демфирующего параметра с помощью построения trade off кривой. При построении этой кривой по оси х откладывается вариация данных (23), а по оси y вариация параметров модели (22). Эмпирическим путем установлено, что такая кривая должна иметь абсолютный минимум. Оптимальным считают то значение демфера, которое соответствует этому абсолютному минимуму.

3.5 Тестирование различных алгоритмов построения сейсмотомографических моделей.

Во всех предыдущих разделах этой главы были описаны результаты тестирования различных компонент, из которых составляется любой сейсмотомографический алгоритм. В данной разделе описаны результаты тестирования четырех различных алгоритмов восстановления трехмерной скоростной модели в целом. Эти тесты помогают разобраться в том, как в совокупности использование различных методов построения лучей, различных способов параметризации модели и различных методов обращения матрицы, сказывается на точности восстановления трехмерной скоростной модели. Сравнение алгоритмов в целом интересно также и потому, что два из четырех сравниваемых алгоритмов: Simulps14 и Fatomo, - достаточно активно применяются в мире для построения сейсмотомографических моделей различных регионов, два других:

TomoCubeFD и TomoTetraFD созданы автором на основе существующего сейсмотомографического программного обеспечения. Информация о четырех используемых алгоритмах и об их составных элементах приведена в Табл. 6.

Табл. 6 Информация о четырех томографических алгоритмах, принимающих участие в тестировании.

A1 A2 A3 A Назв.

SIMULPS14 FATOMO TOMOCUBEFD TOMOTETRAFD программы Thurber Husen Усольцева Усольцева Автор Метод М построения М3_2 (см.2.4.2.2) М7 (см.2.4.2.4) М7 (см.2.4.2.4) (см.2.4.2.4) лучей Способ Линейные Би- Кубические Би- Триангуляция параметр. Блоки сплайны сплайны Делоне Модели Функция, параметр.

Относительная Относительна Относительная Относительная которой скорость я медленность скорость медленность производитс я 1. Обращение матрицы, Способ содержащей 1.разделение решения частные переменных [38] Так же как в обратной производные Так же как в А 2. факторизация А задачи параметров модели Холесского и всех гипоцентров с помощью методов SVD или LSQR Методы Метод сеточного Метод сеточного локации Метод Гейгера поиска поиска событий Из Табл. 6 видно, что программы автора отличаются от существующих методом построения лучей, способом параметризации модели, методом решения обратной задачи и методом перелокации событий. Отличительными особенностями авторских алгоритмов, которые не были включены в таблицу, также является возможность построения с их помощью не только непрерывных, но и квазинепрерывных (функция скорости или медленности претерпевает разрыв на конкретных глубинах) скоростных моделей.

Алгоритмы автора возможно применять для исследования территорий в несколько раз большего размера, чем алгоритмы Simulps14 и Fatomo, т.к. в них вводится коррекция, связанная со сферичностью Земли, которая отсутствует в программах Simulps14 и Fatomo. Как уже отмечалось 2.2.3, при использовании триангуляции Делоне, узлами сетки может быть любой случайным образом заданный набор точек. Данная особенность является существенным преимуществом этого вида параметризации по сравнению с другими способами параметризации при работе с чрезвычайно неравномерно расположенными источниками сейсмических волн. Удается добиться более детального восстановления скоростной структуры в зонах наибольшего скопления сейсмических лучей. Поэтому алгоритм TomoTetraFD, в котором параметризация производится с помощью триангуляции Делоне, особенно эффективен при неравномерном покрытии исследуемой области лучами.

Все описываемые алгоритмы также могут быть использованы в начале работы для нахождения минимальной одномерной скоростной модели.

Для выявления недостатков различных алгоритмов и условий, при которых эти недостатки проявляются, тестирование проводилось по следующей схеме:

1. Расчеты при достаточно равномерном покрытии исследуемой территории лучами.

2. Расчеты при неравномерном покрытии исследуемой территории лучами.

Каждый пункт включал в себя четыре разных теста:

а. шахматный тест для случая, когда все источники являются взрывами;

б. теста пятна для случая, когда все источники являются взрывами;

в. шахматный тест для случая, когда все источники являются землетрясениями;

г. теста пятна для случая, когда все источники являются землетрясениями.

В шахматном тесте тестовая скоростная модель имеет вид типа Рис.

9, т.е. по всему объему равномерно расположены высоко- и низкоскоростные неоднородности в шахматном порядке. В тесте пятна задается тестовая скоростная модель типа Рис. 15, в этой модели присутствуют одна или две обособленные скоростные неоднородности некоторой произвольной формы.

Уровень вводимого при тестировании сейсмического шума примерно соответствовал точности определения времени вступления для используемых в дальнейшем в работе реальных данных (0.1 сек).

Для обеспечения достаточно равномерного покрытия исследуемой области лучами использовалась схема эксперимента, предложенная в [19].

Координаты станций, а также горизонтальная и две перпендикулярные вертикальные проекции расположения источников представлены на Рис. 12.

На Рис. 13 представлена тестовая скоростная модель и результаты расчетов с помощью программ Simulps, Fatomo и TomoCubeFD.

Рис. 12 Карта расположения источников и приемников из [19]. Случай равномерного покрытия исследуемой территории лучами.

Рис. 13 Горизонтальные сечения на различных глубинах тестовой скоростной модели (а) и восстановленных с использованием различных алгоритмов: б) Simulps, в) Fatomo, г) TomoСubeFd.

Рис. 14 Dws на различных глубинах при расчетах с использованием различных алгоритмов: а) Simulps, б) Fatomo, в) TomoCubeFD.

В качестве нулевого приближения используется некоторая одномерная скоростная модель. На Рис. 14 представлены карты сумм взвешенных производных (DWS), которые были рассчитаны при восстановлении скоростной структуры с помощью различных алгоритмов. Карты DWS позволяют определить области наибольшего и наименьшего скопления сейсмических лучей, они необходимы для того, чтобы правильным образом расположить узлы сетки или правильно разбить исследуемый объем на блоки. Из Рис. 13 видно, что при расположении источников и приемников, представленном на Рис. 12, скоростная структура восстанавливается достаточно отчетливо особенно на глубине 15 км. Анализ карт DWS (Рис. 14) показывает, что именно на глубинах 10-15 км исследуемая область достаточно равномерно принизывается лучами. В то же время, возвращаясь к опять к Рис. 12, хотелось бы отметить, что на глубинах 10 и 20 км во всех моделях присутствуют артефактные (несуществующие) неоднородности. На глубине 5 км, скоростная модель, полученная с помощью TomoCubeFD является более контрастной и следовательно более схожей с тестовой моделью, чем скоростные модели, полученные с помощью Simulps и Fatomo.

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

Также в работе проведено сравнение алгоритмов Simulps и Fatomo с использованием теста пятна для случая, когда все источники являются землетрясениями и сейсмические лучи равномерным образом пронизывают весь исследуемый объем. С помощью тестов, проведенных автором для случая равномерного покрытия исследуемой территории лучами, трудно выявить лучший из четырех рассматриваемых методов. Все алгоритмы достаточно хорошо восстанавливают заданные скоростные неоднородности.

При проведении тестирования для случая неравномерно расположенных источников и приемников использовалась схема эксперимента, представленная на Рис. 6. Частично, результаты шахматного теста для случая неравномерного покрытия исследуемой территории лучами описаны в 3.2. Клеточная структура наиболее отчетливо проявляется при использовании параметризации с помощью кубических Би-сплайнов (TomoCubeFD), достаточно хорошо восстанавливается при использовании параметризации с помощью линейных Би-сплайнов (Simulps) и немного хуже при использовании триангуляции Делоне (TomoTetraFD).

Тестовая модель пятна для случая неравномерного покрытия исследуемой территории лучами представлена на Рис. 15. Результаты восстановления этой модели с использованием четырех различных алгоритмов приведены на Рис. 16 и Рис. 17. Видно, что существующее в тестовой модели высокоскоростное включение на глубина 6-9 км присутствует и в восстановленных с помощью различных программ моделях.

Рис. 15 Горизонтальные сечения на различных глубинах тестовой скоростной модели пятна.

На восстановленных скоростных моделях это включение начинается раньше продолжается на большую глубину. Низкоскоростная аномалия, появившаяся в тестовой модели на 9 км, восстанавливается значительно хуже с использованием всех алгоритмов.

Наилучшего восстановления этой аномалии удается добиться с помощью программы TomoCubeFD. Вообще, модели, построенные с помощью TomoCubeFD, в данном случае содержат меньше артефактных скоростных неоднородностей.

пятн Рис. Результаты восстановления тестовой скоростной модели л а ( Рис. 15) с помощью различных сейсмотомографических алгоритмов :

а ) Simulps, б ) Fatomo, в ) TomoCubeFD, г ) TomoTetraFD.

Представлены горизонтальные сечения скоростных моделей на глубинах 0, и км.

Выводы к главе.

1. Результаты тестирования различных методов построения лучей показали, что в большинстве случае удобнее использовать метод конечных разностей для построения лучей:

пятна Рис. Результаты восстановления тестовой скоростной модели л ( Рис. 15) с помощью различных сейсмотомографических алгоритмов :

а ) Simulps, б ) Fatomo, в ) TomoCubeFD, г ) TomoTetraFD.

Представлены горизонтальные сечения скоростных моделей на глубинах 9, и км.

a. он дает минимальное время пробега по сравнению с большинством тестируемых методов;

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

Pages:     | 1 | 2 |    Книги, научные публикации