Описание основных планируемых работ (поэтапно): Первый этап исследований:
Второй этап исследований: 1. Разработка программного кода на языках высокого уровня (С и/или Fortran) и на языке JCL системы SNDA, а также иными средствами, предоставляемыми системой UNIX для интеграции существующих алгоритмов в подсистему реального времени SNDA. 2. Адаптация системы SNDA. 3. Доработка программного кода. 4. Ре-инжиниринг встраиваемого ПО при наличии в нем коммерческого кода, замена коммерческого кода библиотечными функциями SNDA или разработка новых программных модулей (при обоснованной необходимости). Интеграция с базой данной, используемой в системе SNDA. Третий этап исследований: 1. Разработка тестовых наборов данных (на основе реальных данных сейсмических сетей) для интерактивного тестирования интегрированной, расширенной системы SNDA и дальнейшее тестирование SNDA 2. Доработка расширенной системы. 3. Тестирование SNDA как системы реального времени на основе симулирования среды реального времени, создания потоков данных и их ввода в подсистему реального времени. 4. Разработка элементов сейсмической группы. 5. Окончательная интеграция системы на основе одной из операционных систем – Linux или Solaris на платформах Intel или SPARC. Расширенное тестирование системы на реальных потоках данных от существующих сейсмических сетей (напр. GEOPHONE, или иных). Создание элементов сейсмической группы. ^ Как следует из главы 1, для эффективной оценки сейсмической опасности и последующего сейсмического мониторинга необходим достоверный расчет механизма очага, как для исторических землетрясений, так и для моделируемых, что самое важное – для происходящих в реальном времени. Расчет механизма очага – очень сложная процедура, требующего большого количества входных данных, такие расчеты еще совсем недавно казались, неосуществимы без вмешательства оператора и требовали значительных временных затрат. Однако вычислительные мощности достигли уровня, приемлемого не просто для автоматического расчета всех параметров, но и для их расчета в реальном времени с большим количеством входных данных. Для такого расчета в реальном времени требуется определение входных параметров также в реальном времени, а также соответствующим образом спроектированные сейсмические сети и знание скоростной модели земной коры в месте, где произошло событие. В данной работе мы предложим два программных продукта, предназначенных для расчета механизма очага. ^ Данная программа использовалась в университете Калифорнии, Сейсмологической лаборатории Беркли с 1993 года и автоматически исследовала все события с ML>3,5 в северной Калифорнии. В качестве входной информации программа использует файлы в формате SAC. Вместе с тем возможности SNDP позволяют создавать SAC-файлы с записью событий в реальном времени, таким образом, это не должно стать проблемой. Другим источником входных данных для данной программы является хорошо известный набор программ и алгоритмов – Numerical Recipes for C. SNDA, при соответствующей настройке, позволяет как использовать Numerical Recipes в качестве внешнего модуля, так и создавать необходимые входные файлы своими средствами. Остановимся более подробно на TDMT INVC. Основы методологии и основные предпосылки. В самом общем приближении сейсмический источник упрощается рассмотрением одновременно временного и пространственного точечного источника: ![]() Un, наблюденная n компонента смещения, Gni,j – n-ая компонента функции Грина для определенной ориентации пары сил, а Mij – скаляр сейсмического тензора момента, описывающий величину пары сил. Общие пары сил для девиаторного тензора момента могут быть представлены тремя фундаментными сбросами, а именно – вертикальным смещением по простиранию сброса, вертикальным смещением по погружению сброса и 450 смещением по погружению сброса. Индексы i и j определяют географические направления. Это уравнение решено с использованием линейных наименьших квадратов для данной глубины источника. Для этого распределения решен только девиаторный тензор сейсмического момента, а инверсия приводит к Mij, который раскладывается в скаляр сейсмического момента, тензор момента двойной пары и тензор момента диполя компенсированного линейного вектора. Разложение представляется в виде доли двойной пары (Pdc) и доли диполя компенсированного линейного вектора (PCLVD). Изотропная доля всегда равна нулю для данного применения девиатора. Двойная пара далее представляется в виде удара, наклона и погружения двух опорных плоскостей. Основы методологии и разложение тензора сейсмического момента описаны в Jost and Herrmann (1989). Глубина источника находится итеративно через нахождение решения с наименьшей дисперсией, ![]() где, data и sinth – данные и временные серии функции Грина соответственно, а суммирование применено для всех станций и компонент. Другое измерение, полезное для определения глубины источника в зонах, где взрывные события почти исключены, это RES/Pdc, дисперсия, поделенная на долю двойной пары, где: ![]() Отношение дисперсии к доле двойной пары должно стремиться к минимуму. Предполагается, что локация события хорошо представлена высокочастотной локацией гипоцентра, а низкочастотная локация центроида не определена. Также описанное упрощенное представление предполагает, что временная история источника синхронна для всех элементов тензора момента и что она может быть аппроксимирована в дельта-функцию. Эти предположения по большей части разумны для событий с Mw<7.5, поскольку в этом случае используются длиннопериодные волны (>10-20с.). Однако замечено, что для более крупных событий такие точечные представления источника разбивают диапазон задействованных периодов, и необходимы альтернативные конечные приближения сбросов (напр. Dreger and Kaverina, 2000 г.), либо волны большего периода или большие расстояния источник-станция (напр. Fukuyama and Dreger, 2000 г.). И, наконец, предполагается, что модель коры достаточно хорошо изучена для описания распространения низкочастотных волн. ^ Пример 1Цель этого примера – ознакомить с вычислением функции Грина, применением пост-процессинговой фильтрации и использованием программы инверсии тензора сейсмического момента, инвертировав синтетические волновые формы. Файлы, необходимые для Примера 1 расположены в MTPACKAGE/EXAMPLE_1. Необходимо использовать эту директорию как рабочую. В ней находятся следующие файлы: tdmt.config.linux, tdmt.config.sun, MODEL.socal, run_parallel, run_fkrsort, run_filter, b2s.par, mt_inv.in, and testdata[1-3]. Перед началом работы вам необходимо отредактировать файл tdmt.config.linux или tdmt.config.sun под соответствующую платформу. Список необходимых системных переменных представлен ниже: setenv REDI_MT_BINDIR исполняемая директория setenv REDI_MT_PROG_1 запускающая программа setenv REDI_MT_DATDIR директория в которой создаются рабочие директории по отдельным событиям setenv REDI_MT_LOGDIR директория с журналами результатов setenv REDI_MT_SYNTHDIR директория где хранятся функции Грина setenv REDI_MT_EXTRACT_OPTIONS рабочие флажки для программы извлечения qdata setenv REDI_MT_DEBUG_OPTION включение режима «debug» setenv REDI_MT_STATLIST файл координат станции setenv REDI_MT_RESP отклик инструмента (файл с нулями и полюсами) setenv REDI_MT_DATASTREAM используемый поток данных (в данном релизе поддерживаются BH и LH) setenv REDI_MT_PLOT флажок вывода печати (1=да; 0=нет) # Специфические программные переменные setenv REDI_MT_PROG1_PAGE флажок разбиения на страницы (1=да; 0=нет) setenv REDI_MT_PROG1_GFLOC директория где хранятся функции Грина (то же что REDI_MT_SYNTHDIR) setenv REDI_MT_PROG1_STATMAX число подтвержденных станций REDI setenv REDI_MT_PROG1_GETLIST список подтвержденных станций REDI (станции вне этого списка не будут обрабатываться, можно включать станции, не подтвержденные для расчета тензора момента) setenv REDI_MT_PROG1_NSTAT число станций, подтвержденных для расчета тензора момента setenv REDI_MT_PROG1_STATLIST список станций, подтвержденных для расчета тензора момента После внесения изменений необходимо выполнить программу source для этих файлов (а именно source tdmt.config.linux). Шаг 1. Вычисление функции Грина.Файл MODEL.socal содержит пример одномерной скоростной структуры, подходящей для Сьерра Невада и южной Калифорнии. Содержание файла описано ниже. Первая строчка в фале – флажок отладки - .F. Вторая строчка определяет диапазон частот для отладки, эта строчка не меняется. GREEN.X – имя выходного файла. Четвертая строчка определяет alpha (маленькое комплексное число для стабильности интеграции), глубину источника (км), индекс стартовой частоты (X), общее число частот (степень двойки), общее число временных отсчетов (в два раза больше числа частот), период дискретизации (секунды), общее число слоев, и число процессоров, на которых будет исполняться программа (Y). MODEL.socal: .F. 0 64 GREEN.X 6.0 8.00 X 512 1024 0.500 5 Y 1 1 1 1 1 1 1 1 1 1 0 0.5500E+01 0.5500E+01 0.3180E+01 0.2400E+01 600.00 300.00 0.2500E+01 0.6300E+01 0.3640E+01 0.2670E+01 600.00 300.00 0.8000E+01 0.6300E+01 0.3640E+01 0.2670E+01 600.00 300.00 0.1900E+02 0.6700E+01 0.3870E+01 0.2800E+01 600.00 300.00 0.4000E+03 0.7800E+01 0.4500E+01 0.3300E+01 600.00 300.00 3 0.4000000E+03 1.500000E+00 0 1 10000.0 30.0 2.9 2.5 100.00 0.0 10.0 Значения X и Y устанавливаются предварительным скриптом run_parallel, который описывается ниже. Строка 5 – серия флагов, запускающих различные вычисления для разломов фундамента. Строки 6-10 – список таких параметров модели как мощность пласта (км.), скорость p-волн (км/с), s-волн (км/с), плотность (г/см3), Q-альфа и Q-бета. Обращаем внимание, что источник должен находиться на искусственной границе, в которой скорости выше и ниже должны быть одинаковыми (пласты 2 и 3, строки 7 и 8 в примере). 11 строка дает число слоев ниже источника. Строку 12 менять не нужно. 13 строка дает число станций, включенных в вычисления (1 в нашем примере), и окно фазовой скорости для интегрирования. Не стоит менять значения фазовой скорости 10000 и 30 км/с., но необходимо убедиться, что другие два значения меньше чем скорость волн Релея в модели. Последующие строчки определяют расстояние до станции, время задержи (не учитывается в данной версии), и понижение скорости. Рекомендуется оставить значение понижения скорости 8 км/с. Функции Грина начнут вычисляться этой программой от t=расстояние/(понижение скорости). Файл, в том виде, в каком он здесь представлен, сгенерирует функции Грина для станций на расстоянии 100 км. Первой задачей данного примера будет запуск частотно-интегральной программы (FKRPROG) через скрипт run_parallel для вычисления функций Грина, используемых в инверсии синтетических волновых форм. Чтобы функции Грина сначала запускается скрипт run_parallel, который генерирует файлы GREEN: run_parallel <число процессоров> <входной файл>. К примеру, если выполнить команду: run_parallel 1 MODEL.socal, сначала сгенерируется новый входной файл MODEL1 со значениями X и Y равными 1. Затем запускается FKRPROG, производя выходной файл GREEN.1. Если число процессоров больше 1 (максимум 4), то создается несколько входных файлов MODEL[1-4] и FKRPROG выполняется несколько раз, производя соответствующие GREEN.[1-4] файлы. Шаг 2. Пост-обработка для создания функций Грина в применимом формате ASCII.Следующий шаг после FKRPROG – создание ASCII файлов с функциями Грина во временной области, для их последующего использования программой инверсии. Для этой цели используется скрипт run_fkrsort. Этот скрипт вначале вызывает программу wvint9, которая читает файлы GREEN. и осуществляет обратное преобразование Фурье, а затем записывает все функции смещений (см) в бинарный файл vec. Обращаем внимание, что все опции для wvint9 скрыты в скрипте. Существует возможность вывода функций скоростей (см/с) вместо функций смещений (см), если поменять d на v в стандартном выражении ввода для wvint9. Следом за wvint9 из vec файла извлекается определенная временная серия и создается 8-компонентный ASCII файл. Порядок компонент этого файла следующий (и такой же, как Jost и Herrmann (1989): тангенциальные компоненты вертикального смещения по простиранию сброса (tss) и вертикального смещения по погружению сброса (tds), радиальные компоненты вертикального смещения по простиранию сброса (xss), вертикального смещения по погружению сброса (xds)и 450 смещения по погружению сброса (xdd) и вертикальные компоненты вертикального смещения по простиранию сброса (zss), вертикального смещения по погружению сброса (zds)и 450 смещения по погружению сброса (zdd). Ниже пример использования run_fkrsort: run_fkrsort <префикс имени файла> <расстояние> <глубина> <число расстояний>. В нашем примере: run_fkrsort socal 100 8 1. Этот скрипт создаст файл с именем в формате {префикс_имени_файла}{расстояние}d{глубина}.disp. Если число расстояний не равно 1, то скрипт полагает что заданное расстояние – ближайшее, а последующие идут с интервалом 5 км. Обращаем внимание, что если изменить выход скрипта run_fkrsort на скорости, то необходимо также изменить разширения файлов с функциями Грина на vel. Шаг 3. Фильтрация функций Грина.Как только созданы файлы .disp , необходимо произвести полосовую фильтрацию, используя скрипт run_filter: run_filter <имя входного файла> <имя выходного файла> <верхняя граница фильтра> <нижняя граница>. Входной файл - disp. Имя выходного файла может быть любым, о обычно используется имя входного файла без disp верхний и нижний границы фильтра указываются в герцах. Для фильтрации используются команды SAC, встроенные в скрипт. Фильтрация – два прохождения фильтра Баттерворта 4го порядка. В данном примере необходимо отфильтровать функции Грина в полосе от 0.02 до 0.10 Гц: run_filter socal100d8.disp socal100d8 0.02 0.10. Шаг 4. Инверсия тензора момента.Как только вычислены отфильтрованные функции Грина, их можно использовать для инверсии тестовых данных для Примера 1. В рабочей директории находятся три файла данных. Эти файлы перечислены в прилагающемся входном файле для tdmt_invc, который называется mt_inv.in. Рисунок 5 иллюстрирует лучший доступный результат. ![]() Рисунок 5 - Пример графического вывода программы tdmt_invc с использованием функций Грина для источника на глубине 8 км. Данные показаны сплошной линией, синтетика – прерывистой. Для каждой станции сравниваются трёхкомпонентные данные и синтетика и выводится азимут (10, 40 и 50 градусов), максимальная амплитуда и уменьшение дисперсии. Решение включает величины strike, rake и dip для двух возможных плоскостей двойной пары, скаляр сейсмического момента и Mw. Также приводится информация о разложении тензора момента – Pdc, PCLVD и PISO. Также указаны подходящие величины дисперсии (Variance), понижения дисперсии (Var. Red) и дисперсии относительно Pdc (RES/Pdc.). Выходные параметры VAR, VR, Var/Pdc и Quality используются для оценки качества инверсии. VAR – общая оценка дисперсии, VR –уменьшение дисперсии (оценка с учетом и без учета весов расстояний), Var/Pdc – отношение дисперсии к доле двойной пары, quality – субъективная оценка, 4 – лучшая, 1 – худшая. Чем выше значение VR, тем лучше решение. Var/Pdc – также очень полезная характеристика для зон, где решения без двойной пары маловероятны. Для нахождения глубины источника, инверсия проводится по нескольким глубинам и берется лучшее – с наибольшим общим уменьшением дисперсии либо с наименьшим отношением Var/Pdc. Рисунок 6 иллюстрирует уменьшение дисперсии для разных глубин источника для данного примера. ![]() Рисунок 6 - Уменьшение дисперсии для разных глубин. Для каждой инверсии параметр выравнивания функции Грина был оптимизирован. Пример 2В данном примере будут проиллюстрированы возможности вычисления функции Грина и инверсии для широкополосных сейсмических данных. В поддиректории MTPACKAGE/EXAMPLE_2 вы найдете следующие файлы: tdmt.config.linux, tdmt.config.sun, mt_inv.in, run_parallel, run_fkrsort, run_filter, MODEL.gil7, b2s.par. Дополнительно там находятся две поддиректории: DATA_PC и DATA_SUN, содержащие бинарные SAC-файлы для PC и SUN платформ соответственно. Файлы имеют то же назначение, что и в Примере 1. MODEL.gil7 – входной файл для FKRPROG для вычисления функций Грина, подходящих для береговой линии центральной и северной Калифорнии. Перед работой необходимо отредактировать один из файлов tdmt.config.linux либо tdmt.config.sun как было описано в примере 1. Бинарные файлы с данными для станций BKS, CMB, KCC и PKD необходимо поместить в рабочую директорию. Эти станции проиллюстрируют использование программы. Необходимо скопировать следующие файлы: 19980812141000.BKS.BHE 19980812141000.BKS.BHN 19980812141000.BKS.BHZ 19980812141000.CMB.BHE 19980812141000.CMB.BHN19980812141000.CMB.BHZ 19980812141000.KCC.BHE 19980812141000.KCC.BHN 19980812141000.KCC.BHZ 19980812141000.PKD.BHE 19980812141000.PKD.BHN 19980812141000.PKD.BHZ Для проверки файлов на отсутствие повреждений необходимо использовать SAC. Шаг 1. Коррекция инструмента, фильтрация данных и запись данных в ASCII файл.Первым шагом является создание ASCII файлов с трёхкомпонентными данными, которые будут использоваться для tdmt_invc. Этот включает нахождение нулей и полюсов инструмента, использование SAC для вычитания среднего, деконволюции отклика инструмента, интеграции к смещению (см), поворота к тангенциальной и радиальной компонентам, полосовая фильтрация, передискретизация до 1 Гц, и, наконец, запись в ASCII файлы. Всё это делается одним скриптом - tdmt_redi_prepdata – расположенным в директории MTPACKAGE/MTCODE/BIN. Детального описания обрабатывающего скрипта не приводится, но пользователю предоставлено самому опробовать его для понимания обрабатывающих шагов SAC. Обращаем внимание, что файл отклика инструмента – instr.resp – читается программой get_resp, пишущей на выходе SAC файлы нулей и полюсов инструмента, которые впоследствии используются для деконволюции отклика инструмента. Путь к instr.resp устанавливается системной переменной REDI_MT_RESP, описанной в Примере 1. Программа предполагает, что отклики нулей и полюсов в файле instr.resp конвертируют из цифровых единиц в скорость в м/с. Использование tdmt_redi_prepdata для станции BKS (REDI_MT_BINDIR – системная переменная, указывающая путь к исполняемым файлам) показано ниже: $REDI_MT_BINDIR/tdmt_redi_prepdata 19980812141000 BKS 1998 224 36.755 -121.464 0.02 0.05. Первый аргумент командной строки – префикс имени файла с датой начала записи (ггггммддччммсс). Затем имя станции, год и день года (224), широта, долгота события, верхняя (0,02 Гц) и нижняя (0,05 Гц) границы фильтра. В данном случае используется волны от 0,02 Гц до 0,05 Гц, которые считаются эффективными для региональных событий (напр. Pasyanos et al., 1996 г.; Fukuyama и Dreger, 2000 г.). Возможно использование и других фильтров при условии, что функции Грина и данные профильтрованы одинаково и скоростная структура, использованная для вычисления функций Грина дает хорошо представлена теми длинами волн, каковые используются при моделировании. На практике используются полосы фильтрации в зависимости от магнитуды – 0,02 - 0.1 Гц для M< 4,0 Гц; 0,02 – 0,05 Гц для 4,0 Гц <= M < 5,0 Гц; 0,01 – 0,05 Гц для M>= 5,0 Гц. Для очень крупных событий (М> 7,5 Гц) желателен фильтр от 0,005 до 0,02 Гц (Fukuyama и Dreger, 2000 г.). После выполнения tdmt_redi_prepdata для станции BKS отображаются следующие азимут и расстояние: az = 3,314721e+02 dist = 1,420443e+02 Скрипт необходимо выполнить для всех трёх станций, после чего создаются файлы: BKS_f0.05.data, CMB_f0.05.data, KCC_f0.05.data, PKD_f0.05.data. Шаг 2. Вычисление функций Грина подходящих для данных расстояний и глубин источников.Следующий шаг – создание функций Грина в некотором диапазоне глубин источника для каждого расстояния источник-станция. Файл MODEL.gil7 был настроен для вычисления функций Грина на 16 расстояниях. Обычно расстояния округляются до ближайших 5 км. Возможно сгенерировать синтетику gil7, использую следующую последовательность команд: run_parallel 1 MODEL.gil7 run_fkrsort gil7_ 125 8 16 run_filter gil7_*d*.disp gil7_*d* 0.02 0.05. Обращаем внимание, что необходимо четко указать расстояние и глубину вместо звездочек при фильтрации функций Грина и что сохранение каждой вычисленной функции Грина может пригодиться для создания каталога отфильтрованных функций Грина. Использование вычисленных функций Грина, обработанных волновых форм и процедура использования tdmt_invc, описанная в примере 1 определяют наиболее подходящее решение тензора момента для данного события. На рисунке 7 показаны результаты автоматических взаимно корреляционных функций для трёх станций (BKS, CMB и PKD), а на рисунке 8 показан лучший результат при добавлении станции KCC и при инверсии, проведенном в некотором диапазоне глубин источника. ![]() Рисунок 7 - Результат, получаемый при использовании функций Грина для источника на глубине 8 км и возможности автоматического нахождения взаимно корреляционной функции программой для трёх станций: BKS, CMB и PKD. ![]() Рисунок 8 - Лучшее решение, полученное при переборе глубин источника – 4,5; 8; 11 и 14 км для станций BKS, CMB, PKD and KCC. Подобранная глубина – 11 км (хотя решение плохо сошлось для одной станции).
|