CC BY
56
Асминг В.Э. и др. Моделирование сейсмической локации в трёхмерных средах
УДК 550.34.01 : 519.245 : 004.436
Моделирование сейсмической локации в трёхмерных средах
В.Э. Асминг, Ю.А. Виноградов, А.В. Прокудина
Кольский филиал Геофизической службы РАН
Аннотация. Дано описание и представлен принцип работы программы "Сейсмоконфигуратор", разработанной в Кольском филиале Геофизической службы РАН (КФ ГС РАН). Показано, что данная программа позволяет задать трёхмерную скоростную модель сейсмической среды, лоцировать в ней сейсмические события и моделировать ошибки локации, вызванные неточностью скоростной модели и погрешностями измерений времён приходов сейсмических волн на станции. Применён новый подход к моделированию, в котором проводящая среда заменяется стохастическим графом из большого числа вершин и рёбер. Для описания моделей предложен простой язык задания параметров.
Abstract. The description and principle of operation of the Seismic configurator programme developed in the Kola branch of Geophysical Survey of RAS have been presented. The programme allows a user to specify 3D seismic medium velocity models, to locate seismic events and to model location errors caused by uncertainties of travel time models and onset picking errors. A new approach to modeling (the seismic medium is replaced by a random graph consisting of large number of vertices and edges) has been applied. Models have been specified by a simple parameters description language.
Ключевые слова: сейсмическая локация, сеть сейсмостанций, скоростная модель, моделирование, стохастический граф Кеу words: seismic location, seismic network, velocity model, modeling, random graph
1. Введение
В ОАО "Апатит" проводится непрерывный сейсмический мониторинг полей объединённого Кировского и Расвумчоррского рудников (Хибинский горный массив) при помощи собственных локальных сейсмосетей. Кольский филиал Геофизической службы РАН (КФ ГС РАН) выполняет мониторинг сейсмичности всей Европейской Арктики. В 2010 г. обе системы мониторинга были объединены для наблюдения за сейсмичностью Хибинского массива (Аккуратов и др., 2011). Точность локации сейсмических событий в пятидесятикилометровой зоне производственной деятельности ОАО "Апатит" существенно ниже, чем в зоне подземных рудников, контролируемых системой сейсмических станций ОАО "Апатит". Так, например, в зоне рудников Кировский и Расвумчоррский ошибка локации эпицентров сейсмических событий составляет десятки метров, в зоне карьера Центральный она увеличивается до первых сотен метров, а в зоне карьеров рудника Восточный, хвостохранилищ АНОФ-2 и АНОФ-1 может достигать километра и более. Такая точность недостаточна для целей мониторинга. Повышение точности локации возможно двумя способами -увеличением количества станций наблюдения и уточнением используемой скоростной модели. В то же время, обеспечение максимальной точности на всей контролируемой территории экономически нецелесообразно и трудноосуществимо по техническим причинам. Для понимания, насколько установка новых станций и схема их расположения обеспечит достаточное качество наблюдений, необходимо математическое моделирование сейсмической локации в трёхмерных средах, в том числе с учётом изменяющейся в ходе проходки геометрии горных выработок.
В КФ ГС РАН разработана программа "Сейсмоконфигуратор". Программа позволяет задать трёхмерную скоростную модель сейсмической среды - горизонтально-слоистую среду с вкраплёнными в неё телами произвольной формы, в которых могут быть заданы сейсмические скорости, отличающиеся от скоростей в окружающей среде, в том числе и нулевые (пустоты); лоцировать в этой среде сейсмические события и моделировать ошибки локации, вызванные неточностью скоростной модели и погрешностями измерений времён приходов сейсмических волн на станции.
2. Принцип локации сейсмических событий
Стохастический граф. В данной модели трёхмерная проводящая сейсмические волны среда заменяется стохастическим графом (рис. 1). Это набор вершин - точек, случайным образом "высыпанных" в интересующую нас область Земли, и рёбер - связей между ними. Связи задаются между вершинами, находящимися не дальше заданного расстояния друг от друга. Для каждой пары вершин по скоростной модели рассчитываются времена пробегов волн между ними. Это - веса рёбер в нашем стохастическом графе.
644
Вестник МГТУ, том 16, № 4, 2013 г. стр.644-649
Рис. 1. Фрагмент стохастического графа Рис. 2. Пример стохастического графа с неравномерной
плотностью вершин. Показан фрагмент модели, предназначенной для локации событий в Хибинском и Ловозерском массивах. В них плотность вершин графа максимальна. Вне массивов плотность вершин выше в верхнем слое
Вершины стохастического графа могут располагаться в пространстве неравномерно (рис. 2). Пользователь задаёт плотности распределения вершин графа (единицы измерения - число вершин на кубический километр) в различных участках моделируемой области.
Связи (рёбра графа) задаются другим параметром - критическим расстоянием. Это максимальное расстояние между вершинами, при котором они ещё связываются. Оно может быть не постоянным, а являться функцией координат, задаваемой пользователем.
Таким образом, для задания стохастического графа пользователю достаточно задать две пространственно-распределённые функции - плотность (далее density) и критическое расстояние для установки связи (далее radius).
Алгоритм локации на графе. Существует быстрый алгоритм поиска кратчайших путей в графе, различные модификации которого известны как волновой алгоритм, алгоритм Ли (Lee, 1961), алгоритм Дейкстры (Dijkstra, 1959). В нашем случае этот алгоритм действует так. В некой начальной вершине графа задаётся время выхода оттуда сейсмической волны или время, когда там произошло сейсмическое событие. В результате работы алгоритма для всех остальных вершин графа определяются времена приходов в них сейсмической волны, вышедшей из исходной вершины.
Сейсмостанции включены в стохастический граф, т.е. если мы моделируем N сейсмостанций, то они соответствуют первым N вершинам графа.
Рис. 3. Иллюстрация локации методом инверсии. Если волны от события достигли станций в моменты времени t1, t2, t3, то вторичные волны, выпущенные из станций в моменты -t1; -t2, -t3 придут в точку -источник в один и тот же момент времени 0. Звёздочкой обозначено сейсмическое событие, а
треугольниками - станции
645
Асминг В.Э. и др. Моделирование сейсмической локации в трёхмерных средах
Локация производится методом инверсии времён приходов (рис. 3). Идея этого метода родилась в обсуждении с ныне покойным Александром Фёдоровичем Буяновым1. Она заключается в следующем. Пусть мы имеем времена приходов волны на N станций - 4 Тогда, если T - времена пробега волны от источника до i-й станции, для всех станций величины t-T совпадают и равны времени, когда произошло событие.
Переходя из среды в стохастический граф, можем сказать, что для вершины графа, наиболее близкой к результату локации, величины T-ti, где теперь - время пробега от i-й вершины (приёмника) до проверяемой вершины, для всех i должны быть близки между собой, т.е. их дисперсия минимальна.
Таким образом для локации события N раз применяется алгоритм Ли, за начальную вершину берётся каждый приёмник (i), а за начальное время - инвертированное время прихода -4 В каждой вершине графа по результатам алгоритма Ли накапливаются дисперсии времён. По окончании процесса ищется вершина с минимальной дисперсией.
3. Описание модели локации
Для моделирования и локации в программу должна быть введена исходная информация: скоростные модели, структура стохастического графа, сеть используемых сейсмостанций, рельеф и управляющие параметры. Совокупность этих данных мы будем называть моделью локации. Такая модель задаётся в виде текста на специализированном языке задания параметров (ЯЗАП).
Язык задания параметров ЯЗАП. Язык предназначен для унифицированного задания произвольных параметров. Язык декларативный, состоит исключительно из операторов присваивания. Синтаксис его крайне прост. Имена параметров специально не оговариваются, программа, в зависимости от контекста, сама извлекает нужные значения из текста. В этом ЯЗАП очень близок языку XML. Оператор присваивания значения имеет вид <имя>=<значение>. Операторы разделяются пробелами. Имя - идентификатор, состоящий из латинских букв и цифр. Значение может быть идентификатором, числом, строкой текста или списком (набором вложенных операторов присваивания, заключённых в круглые скобки). Использование списков позволяет реализовать иерархические структуры параметров. В текст в фигурных скобках могут быть включены комментарии.
Примеры операторов с комментариями приводятся ниже:
x=123 {Присвоение числового значения параметру х}
text=^TO текстовая строка' {Присвоение строкового значения параметру text}
color=(red=255 green=0 blue=100) {Присвоение значения - списка}
Надо отметить, что декларативный язык ЯЗАП, при всей его внешней простоте, оказался мощным средством задания произвольных параметров. Авторы используют его не только в данной системе, но и во множестве других приложений, в частности, для создания конфигурационных файлов различных программ.
Задание рабочей области моделирования. Описание модели начинается с задания рабочей области моделирования - трёхмерного фрагмента Земли, для которого моделирование, собственно, и проводится. Она задаётся в виде параллелепипеда по широтам и долготам (в градусах) и глубинам (вниз в километрах).
Задание произвольных объёмных тел. В данной программе можно описывать объёмные тела произвольной формы. Они используются для задания плотностей вершин стохастического графа (например, можно задать области, в которой вершины будут располагаться гуще, чем в окружающем пространстве и т.д.), а также для задания скоростных моделей - внутри объёмных тел скорости могут отличаться от скоростей, заданных в окружающей среде.
Объёмные тела задаются в виде систем замкнутых контуров, расположенных на горизонтальных плоскостях. Контуры увязываются друг с другом и образуют нечто вроде вертикально стоящего цилиндра (рис. 4).
Однако хранение объёмных тел в памяти компьютера в виде системы связанных контуров не оптимально. Используемые в данной программе алгоритмы требуют многократных проверок того, попала ли какая-нибудь точка внутрь объёмного тела. Такие проверки для тел, заданных контурами, заняли бы слишком много времени. Поэтому перед их использованием тела "дискретизируются". Каждое тело заключается в параллелепипед с осями, параллельными осям координат. Он разбивается на кубики. Дискретизируемое тело огрубляется - мы считаем, что если центр кубика попал в тело, то и весь кубик принадлежит телу. Параллелепипед с кубиками представляется в виде трёхмерного битового массива, в котором 0 обозначает, что кубик не принадлежит исходному телу, а 1 - что принадлежит (рис. 5).
1 Александр Фёдорович Буянов (1959-1997) — талантливый геофизик, работавший в Геологическом институте КНЦ РАН, кандидат физико-математических наук. Человек с обширной сферой интересов, которая включала, наряду со многим другим, сейсмологию и сейсмическую томографию.
646
Вестник МГТУ, том 16, № 4, 2013 г.
стр.644-649
Рис. 4. Пример задания объёмного тела в виде трёх связанных контуров, лежащих на горизонтальных параллельных плоскостях
Такое представление несколько огрубляет объёмные тела, но позволяет компактно хранить их в памяти компьютера и быстро проверять, попала ли в них точка.
Задание стохастического графа. Стохастический граф в данной программе задаётся двумя пространственно распределёнными функциями - мгновенной плотностью вершин (далее density), измеряемой в числе точек на кубический километр и максимальным расстоянием от данной вершины до другой, при которой между ними устанавливается связь. Далее эта функция обозначается как radius, измеряется в километрах.
Эти функции задаются в виде иерархической структуры. На нижнем уровне иерархии задаются базовые постоянные значения во всём пространстве. На следующем задаётся горизонтально-слоистая модель, в каждом слое значения функций постоянны. Затем задаются объёмные тела, в каждом из которых функции принимают другие (постоянные в пределах тел) значения.
Объёмные тела в описании графа задаются в порядке возрастания приоритета. Это значит, что если точка, в которой мы вычисляем значение параметров, попала в какое-либо тело, тела, которые были заданы до этого, не проверяются, берутся значения, связанные с данным телом. Если точка не попала ни в одно из тел, проверяется, не попала ли она в слоистую модель. Если да - берутся значения из соответствующего слоя. Если нет - берутся базовые значения.
Задание скоростных моделей. В данной программе под скоростной моделью понимается функция координат, описывающая значения скорости определённого типа сейсмической волны (P или S) в каждой точке рабочей области. В программе реализовано задание модели в виде горизонтально -слоистой среды с вкраплёнными в неё телами с изменёнными значениями скорости (структура и синтаксис задания полностью аналогичны заданию стохастического графа). Как и в случае задания стохастического графа, геометрические тела в моделях располагаются в порядке возрастания их приоритета. Последнее тело имеет наивысший приоритет.
Всего в программе используется четыре модели скорости (по две для P и S-волны). Так называемые прямые или "forward", обозначаются как P1 и S1. Их смысл при моделировании в том, что, когда мы выясняем возможную ошибку локации некоего сейсмического события, мы задаём его координаты и считаем времена приходов сейсмических волн на станции согласно этим моделям. Две других модели называются обратными или "backward" и обозначаются как P2 и S2. Их смысл в моделировании в том, что, когда рассчитаны времена приходов по моделям P1, S1, локация производится по моделям P2, S2. Расхождение между исходной точкой сейсмического события и результатом локации и будет оценкой погрешности за счёт разницы скоростных моделей.
4. Моделирование и локация
Задания на моделирование локации составляются на том же самом языке задания параметров ЯЗАП, что и модель среды. Реализованы два принципиально разных варианта моделирования - локация по временам и обратная локация, а также трассировка лучей из точки в точку.
Локация по временам. В этом варианте задаются времена приходов сейсмических волн на станции и, возможно, диапазон их предполагаемых ошибок. Задаётся также количество испытаний. Программа генерирует времена приходов, добавляя к заданным временам случайные ошибки, и лоцирует сейсмическое событие заданное количество раз. Результат отображается на карте и разрезах в виде облака точек (результатов локации) и эллипса ошибок (рис. 6).
По умолчанию локация производится по обеим парам скоростным моделей (P1, S1) и (P2, S2). На карту наносятся две точки локации и два эллипса ошибок. Оценивается расстояние между локациями. Можно лоцировать по одной модели 1 или 2, указав это явно в тексте задания.
Данный режим может быть использован как для моделирования ошибок локации, так и просто для локации сейсмических событий по трёхмерной скоростной модели.
647
Асминг В.Э. и др. Моделирование сейсмической локации в трёхмерных средах
Рис. 6. Отображение на карте результата исполнения задания на локацию. Кружком отмечен средний по всем вариантам результат локации, звёздочки - отдельные локации по временам вступлений с внесёнными случайными погрешностями. Число (3.1575) соответствует длине большой оси эллипса ошибок в километрах
Пример. Задание на локацию события в Хибинском массиве и результат его исполнения.
station=(name='APA' tp=(1999 8 13 11 48 7.6)
ts=(1999 8 13 11 48 10.3) errp=0.1 errs=0.2) station=(name='AP0' tp=(1999 8 13 11 48 10.2)
ts=(1999 8 13 11 48 14.4) errp=0.1 errs=0.2) station=(name='LVZ' tp=(1999 8 13 11 48 10.6)
ts=(1999 8 13 11 48 15.7) errp=0.1 errs=0.2) tests=20 {Число испытаний} model=1 {Номер скоростной модели 1/2}
Обратная локация. Данный режим предназначен для моделирования ошибок локации, вызванных как неточностью скоростной модели, так и погрешностями измерений времён приходов волн (рис. 7). В нём пользователь задаёт координаты точки, в которой как будто бы произошло сейсмическое событие. Программа генерирует времена приходов волн P и S на все станции, используя скоростные модели P1 и S1 (как будто бы это истинные скорости). Производит несколько испытаний, внося в эти рассчитанные времена случайные погрешности, и по этим, возмущённым, временам лоцирует сейсмическое событие, но используя скоростные модели P2 и S2 (которые могут отличаться от P1, S1, имитируя тем самым ошибку задания модели). В результате оценивается систематическое отклонение получаемых координат от заданной точки - известных координат события и параметры эллипса ошибок (рис. 8).
Рис. 7. Иллюстрация обратной локации. Для некой исходной точки согласно модели 1 рассчитаны времена приходов волн на станции t\, t2, t3. К ним добавлены случайные ошибки и проведены случайные испытания. Рассчитан эллипс ошибок и систематическая погрешность
648
Вестник МГТУ, том 16, № 4, 2013 г.
стр. 644-649
Рис. 8. Результат моделирования обратной локацией. Показано истинное положение лоцируемого события (звёздочка слева) и результаты локации по временам вступлений с внесёнными в них случайными погрешностями (звёздочки справа). По ним построен эллипс ошибок, большая ось которого имеет размер 6.85 км. Систематическая ошибка, вызванная различием исходной скоростной модели и модели локации, соответствует расстоянию между средней локацией и истинным положением события
и в данном случае составляет 2.9 км
Пример задания на обратную локацию и результат счёта приведены ниже.
{Исходная точка}
Fi=67.66000 Ld=33.93014 d=0
{Какие волны каких станций использовать и с какими ошибками} station=(name='APA' p=yes errp=0.1 s=yes errs=0.2) station=(name='AP0' p=yes errp=0.1 s=yes errs=0.2) tests=10 {Число случайных испытаний}
Трассировка лучей. В данном режиме пользователь задаёт две точки. Моделируется кратчайший путь сейсмической волны из одной точки в другую. Для этого исходные точки "сносятся" в стохастический граф, т.е. заменяются ближайшими вершинами графа. Путь сейсмической волны моделируется кратчайшим путём в графе. Он рисуется на карте и разрезах. Считается эффективная скорость пробега, т.е. отношение расстояния между точками ко времени пробега.
Такая возможность может оказаться полезной при отладке скоростных моделей, подборе моделей под наблюдения, проверки адекватности стохастического графа.
5. Заключение
Применение программного продукта "Сейсмоконфигуратор", разработанного в Кольском филиале Геофизической службы РАН, поможет выработать оптимальную конфигурацию сети сейсмостанций для мониторинга зоны производственной деятельности ОАО "Апатит", а также скоростную модель, достаточную для локации сейсмических событий в этой зоне.
Литература
Dijkstra E.W. A note on two problems in connexion with graphs. Numerische Mathematik, v. 1, p. 269-271, 1959.
Lee C.Y. An algorithm for path connections and its applications. IRE Transactions on Electronic Computers, v. EC-10, N 2, p. 364-365, 1961.
Аккуратов М.В., Асминг В.Э., Виноградов Ю.А., Корчак П.А. Объединённая система контроля состояния Хибинского горного массива на базе сетей сейсмических станций Кольского филиала
ГС РАН и ОАО "Апатит". Мат. 6-й Междун. сейсмологической школы "Современные методы обработки и интерпретации сейсмологических данных", Апатиты, 15-19 августа 2011 года. Обнинск, ГС РАН, с. 7-10, 2011.
649
