CC BY
30
УДК 550.34.(013.4+042.4)
ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ХАРАКТЕРИСТИК САМОПОДОБИЯ СЕЙСМИЧНОСТИ
1 9
© Ф.Л. Зуев1, А.В. Ключевский2
Институт земной коры СО РАН,
664033, Россия, г. Иркутск, ул. Лермонтова, 128.
Описаны методические особенности вычисления параметров самоподобия сейсмичности, примененные в практике сейсмологических исследований Байкальского региона и Монголии. Усовершенствованные алгоритмы для клеточной, информационной, корреляционной размерности показателя Херста позволяют улучшить результат при небольшом количестве исходных данных в выборках землетрясений этих регионов. Фрактальные оценки показывают точность до второго десятичного знака. Ошибка растет с увеличением фрактальной размерности, что открывает пути усовершенствования метода.
Ключевые слова: сейсмичность; самоподобие; фрактальные размерности; информационная размерность; корреляционная размерность; показатель Херста.
COMPUTATION ASPECTS OF SEISMICITY SELF-SIMILARITY CHARACTERISTICS F.L. Zuev, A.V. Klyuchevsky
Institute of the Earth Crust SB RAS, 128 Lermontov St., Irkutsk, 664033, Russia.
The paper describes the methodological aspects of seismicity self-similarity parameters computation as applied to seis-mological studies in the Baikal region and Mongolia. Improved algorithms for cellular, information and correlation dimensions of the Hurst exponent produce better results for decreased amounts of input data on earthquakes in these regions. Estimated fractal values show the precision up to second decimal digit. A systematic error increases proportionally to fractal dimension growth that opens possibilities for further improvements of the method.
Keywords: seismicity; self-similarity; fractal dimensions; information dimension; correlation dimension; Hurst exponent.
Проблемы описания структуры дискретного пространственно-временного распределения землетрясений привели к использованию методов теории самоподобия для обработки и интерпретации материалов каталогов сейсмических событий. При исследовании фрактальных характеристик сейсмичности по данным «Каталога землетрясений Прибайкалья» и «Каталога землетрясений Монголии» у авторов часто появлялась необходимость в уточнении и усовершенствовании некоторых общепринятых методических приемов и вычислительных процедур. В каталогах собрана информация примерно о ста тысячах сейсмических событий, и необходимость в уточнении возникала обычно при использовании небольших пространственно-временных выборок землетрясений и разделении толчков по энергетическим классам. В настоящей работе показано, как и почему авторами внесены некоторые изменения в методику оценки параметров самоподобия сейсмичности, не изменяя основных принципов фрактального анализа.
Известно [7], что фрактальная размерность, или показатель самоподобия, - специфическая характеристика геометрических свойств природных или модельных объектов. Любая геометрическая характеристика упрощает и идеализирует измеряемый объект, но ес-
ли классические евклидовы характеристики, такие как длина или площадь, подразумевают непрерывность и гладкость исследуемого объекта, то различные фрактальные характеристики позволяют оценить и описать его сложности - неоднородность, разрывность и неравномерность. В «прямой» задаче моделирования фрактального объекта эти сложности формируются достаточно просто - последовательными итерациями некоторого масштабированного «генератора», наложенного на «инициатор». Инвариантность относительно мультипликативных изменений масштаба природных объектов обусловлена самоподобием пространственно-временных процессов, но о самоподобии и параметрах исследуемой структуры можно говорить только в заданном масштабном диапазоне.
Периодические повторения в геологической истории Земли некоторых процессов, таких как суточное и годовое вращение планеты, аналогичны итерациям масштабированного «генератора» и вполне могут создать фрактальные структуры. Региональные и локальные фрактальные структуры формируются происходящими на этих территориях периодическими и квазипериодическими повторениями одних и тех же природных процессов. В литосфере Байкальского региона формирование фрактальных объектов происходит
1Зуев Федор Леонидович, ведущий инженер лаборатории современной геодинамики, тел.: (3952) 471417, 89501376016, e-mail: fedor@crust.irk.ru
Zuev Fedor, Leading Engineer of the Laboratory of Modern Geodynamics, tel.: (3952) 471417, 89501376016, e-mail: fedor@crust.irk.ru
2Ключевский Анатолий Васильевич, ведущий научный сотрудник лаборатории инженерной сейсмологии и сейсмогеологии, тел.: (3952) 427462, 89500766435, e-mail: akluchev@crust.irk.ru
Klyuchevsky Anatoly, Leading Researcher of the Laboratory of Engineering Seismology and Seismogeology, tel.: (3952) 427462, 89500766435, e-mail: akluchev@crust.irk.ru
преимущественно под влиянием структур-аттракторов рифтогенеза, ритмичное импульсное воздействие которых приводит к перестройкам напряженного состояния среды и изменению скорости потока землетрясений [9]. Поскольку это не идеальные, а реальные процессы с различного рода флуктуациями, то самоподобие и фрактальные оценки в приложении к случайным множествам, таким как сейсмичность, - понятия не столь строгие. В этом случае части не обязательно должны быть в точности подобны целому и достаточно того, что части и уменьшенное в масштабе целое имеют одинаковые распределения. Такие оценки характеризуют структуру сейсмичности в первом приближении и используются обычно для сравнительного анализа при условии относительной однородности и сопоставимости используемых пространственно-временных выборок данных.
Основные понятия фрактальной геометрии существовали в математике с конца XIX века и были обобщены Б. Мандельбротом [7]. Мы применяем математический аппарат фрактальной математики и развитые в его рамках метрики и алгоритмы для изучения самоподобия сейсмичности, т.е. законов сохранения ее пропорций при изменении пространственно-временных масштабов [2-6]. Определения фрактальной размерности сами по себе алгоритмичны (рекурсивны) и поэтому часто возникает соблазн буквально запрограммировать их и этим ограничиться. Однако такой алгоритм может создавать значительную ошибку (что легко проверяется на тестовых данных с заранее известной размерностью) и может быть неустойчив, когда небольшие случайные изменения исходных данных приводят к скачку результата. Поэтому в данной работе описаны методические особенности и протестированы результаты вычисления параметров самоподобия сейсмичности, примененные в практике сейсмологических исследований Байкальского региона и Монголии [2-6].
Теоретические основы и общие вычислительные приемы
Мы используем известные метрики фрактальной геометрии - клеточную (Хаусдорфа), информационную и корреляционную размерности и показатель Херста применительно к геофизическим явлениям, описываемым наборами точек, прежде всего эпицентрами землетрясений [4, 5]. С инженерно-измерительной точки зрения нас интересует достижимая точность оценки размерностей и способы ее улучшения. Можно выделить следующие источники погрешности:
1. Неопределенность измеряемой величины. В общем случае под фрактальной размерностью эпи-центрального поля землетрясений подразумевается размерность множества всех сейсмических событий, которые могли произойти на данной территории за большой период времени. Зарегистрированные землетрясения, представленные в каталогах, являются ограниченной пространственно-временной выборкой и параметры сейсмичности, определенные по такой выборке, отражают объективную неопределенность между полной и ограниченной выборкой.
2. Погрешность, обусловленная аппроксимацией природного объекта набором точек. Под
эпицентром землетрясения подразумевается проекция центра очага землетрясения на земную поверхность. Размеры проекций сильных и слабых землетрясений могут отличаться на четыре порядка. Область очага землетрясения обычно представляется эллипсоидом, а проекция ее на поверхность есть эллипс.
3. Погрешности и ошибки в исходных данных
зависят от расположения и плотности расстановки региональных сейсмических станций, точности калибровки сейсмографов и квалификации обслуживающего персонала. Ошибки в исходных данных также зависят от принятых скоростных моделей среды и законов затухания сейсмических волн с увеличением расстояния.
Эти три источника неопределенностей и погрешностей присутствуют всегда.
4. Погрешности, связанные с особенностями применяемого алгоритма, можно оценить при сопоставлении с тестовыми расчетами известной размерности и внести необходимые исправления.
Среднеквадратичная аппроксимация. В определениях характеристик самоподобия используются выражения вида:
/ (^)
D = lim
гlog r '
(1)
где условие г^0 не следует понимать буквально. Так, спектр размерностей Реньи, частными случаями которого являются клеточная О0, информационная 01 и корреляционная 02 размерности, определяется как [1]:
D
1
lim--—
q-1 r^o log r
(2)
где РI = Ы/Ы- доля точек в /'-ой ячейке от общего количества точек во множестве, или вероятность найти точку в /-ой ячейке сетки. Но при его вычислении не изучается значение 0Я в ближайших окрестностях нуля, а ищется такое О, при котором ^р в интересующем нас диапазоне наилучшим образом квадратично аппроксимируется степенной функцией г0. Геометрически это обычно интерпретируют как коэффициент при г или тангенс угла наклона линейного графика в двойном логарифмическом масштабе.
Разбиение на элементарные ячейки. Большинство определений масштабных характеристик начинается с покрытия изучаемого множества сеткой квадратных ячеек с последующей манипуляцией над «заполненными» клетками. При изучении фракталов сетку на плоскости обычно строят, описывая вокруг изучаемого множества квадрат, который затем двумя ортогональными линиями делят на четыре части, на следующем шаге каждую из них делят еще на четыре части вдвое меньшего размера и т.д. Такое построение удобно для математических доказательств, поскольку позволяет использовать рекурсию, однако не является единственным возможным и с вычислитель-
ной точки зрения обладает рядом недостатков. Так, изучаемая область должна быть квадратом, первые разбиения могут дать слишком грубые значения, ухудшающие результат. Кроме того, подобное разбиение приводит к расширению диапазона масштабов до нескольких десятичных порядков, и при фиксированных границах ячеек повышается влияние случайных флуктуаций. О.И. Шелухин и его соавторы [8] рекомендуют начинать с разбиения каждой стороны не меньше чем на пять частей. Из-за ограниченности имеющихся у нас выборок данных (иногда менее сотни точек) и ошибки в определении координат эпицентров землетрясений (около 10 км), мы начинали с разбиения на три части, при котором стандартное отклонение оценивается величиной до 5%.
Расположение ячеек сетки. Определения параметров самоподобия нигде не требуют ни того, чтобы разбиение начиналось с минимально возможного деления надвое, ни того, чтобы размер ячейки следующей сетки был кратен размеру ячейки предыдущего, ни того, чтобы покрытие множества сеткой было наиболее компактно. Поэтому для вычисления различных фрактальных размерностей нами разработан следующий алгоритм покрытия множества точек на плоскости последовательностью сеток с уменьшающимися ячейками:
1. Находится минимальный прямоугольник со сторонами параллельными осям координат, описанный вокруг изучаемого множества. Крайние точки множества оказываются на сторонах этого прямоугольника.
2. Размер стороны ячейки г1 первой сетки устанавливается в 1/3 от длины меньшей стороны описанного прямоугольника.
3. Сетка с ячейками размера rk накладывается на плоскость таким образом, что левый нижний угол первой ячейки сетки устанавливается на г/2 ниже и на г/2 правее левого нижнего угла описанного прямоугольника, другими словами, левый нижний угол описанного прямоугольника попадает в центр первой ячейки.
4. Над сеткой производятся необходимые вычисления, и создается некоторая (своя для каждого метода) функция Я(гк).
5. Размер ячейки следующей сетки устанавливается в 80% от размера предыдущей - = 0,8-гк. Число 80% было подобрано эмпирически, большие значения существенно удлиняли счет, ничего не добавляя к точности результата.
6. Начиная с п. 3, повторять до тех пор, пока размер ячейки гк не окажется меньше заданной точности входных данных в 10 км.
Выделение отрезка последовательности. Не все члены последовательности Я(гк), / = 1, 2, ... k полезны для аппроксимации. Крайние ее члены иногда ухудшают качество оценки, в частности, потому что:
- первые значения Я(г) часто характеризуют изучаемое множество как нерасчлененное целое, не давая информации о соразмерности ее частей;
- по мере того, как число непустых ячеек сетки становится сопоставимо с числом точек множества, Я(г) начинает вырождаться, отражая вместо размерностей связного объекта, представленного набором то-
чек, размерность (нулевую) самого набора точек;
- точность входных данных ограничена, и выходя за ее пределы, мы вносим ошибки в расчеты.
Авторы использовали следующий алгоритм отсечения крайних членов.
1. В начале последовательности пропускаются все члены, у которых все ячейки сетки оказываются непустыми. Отсчет начинается с того момента, когда у нас появляется хотя бы одна ячейка сетки, лежащая вне измеряемого множества.
2. Построение последовательности прекращается, когда число непустых клеток становится больше чем N/d, т.е. в каждой ячейке в среднем оказывается меньше чем d точек [10] (d - топологическая размерность).
3. Построение прекращается, если сторона ячейки сетки становится меньше точности входных данных.
Алгоритмы вычислений
Клеточная (Хаусдорфа) размерность D0. Самая простая и наиболее часто используемая из фрактальных метрик. Она выражается формулой
-log n(r)
(3)
D0 = lim
r^o log r
где n(r) - количество непустых, содержащих хотя бы одну точку, клеток сетки с длиной стороны r, накрывающей изучаемое множество. Геометрический смысл метрики - показатель изменения покрываемой множеством области в зависимости от единицы измерения и характерного масштаба. Алгоритм вычисления достаточно прост. На множество накладываются сетки с последовательно уменьшающимся размером ячейки rk до тех пор, пока rk не станет меньше точности измерения. Затем у последовательности пар (1/n(rk), rk), k = 1, 2, 3, ... отбрасываются неинформативные крайние члены. Оставшиеся пары аппроксимируются методом наименьших квадратов, а среднеквадратичное отклонение берется как оценка погрешности метода (что не совсем верно методически, но дает представление о порядке величины). Алгоритмы использовались для топологических размерностей один (временной сейсмический процесс) и два (поле эпицентров землетрясений) [2-6].
Информационная размерность D1. Вторая из спектра размерностей Реньи. Она выражается формулой
А
lim
r
-m
log r
IP log P
lim
r ^ 0
log r
(4)
где S(r) —^p l°g p - информационная энтро-
i
пия по Шеннону. Геометрический смысл метрики -показатель изменения информации, необходимой для определения положения точек в зависимости от единицы измерения [1]. Алгоритм вычисления аналогичен вычислению клеточной размерности с тем различием, что вместо log n(r) вычисляется S(r) [4, 5].
Корреляционная размерность D2. Определяется формулой
log с ( r )
(5)
А
lim
log r
где С(г) - корреляционный интеграл, то есть вероятность того, что расстояние между двумя произвольными точками множества меньше г. Интеграл может быть записан через функцию Хевисайда:
C(г) = Тд('-|(X -X
(6)
где суммирование проводится по всем парам точек x,, Xj изучаемого множества. Нами использовался алгоритм Грассбергера-Прокаччиа [1]. Для геометрической прогрессии масштабов rk рассчитываются значения C(rk), затем выполняется удаление крайних точек и аппроксимация аналогично предыдущему [4, 5].
Показатель Херста H. Является характеристикой случайного процесса и не применяется непосредственно для оценки множества точек. Построим случайную функцию Am(t) как число точек, попадающих в заданный отрезок времени (t, t+ût) (ее можно выразить аналогично (6) через функцию Хевисайда). Согласно [7] случайная функция будет иметь H = 2-D0, где D0 - размерность Хаусдорфа образующего множества.
Мы использовали алгоритм вычисления H через индекс дисперсии для отсчетов (Index of dispersion for counts, IDC) [8]. Индекс дисперсии IDC - отношение дисперсии к математическому ожиданию
A(r ))
IDC (г) = -
(7)
М(Л(г))'
Имеющийся временной отрезок длинной Т делится последовательно на к = 5, 6, 7, ... ячеек длины гк = Т/к до достижения пороговой точности либо до
того, как количество ячеек достигнет 1/5 от общего количества точек. Подсчитывается число точек в каждой ячейке, и вычисляются выборочное среднее и дисперсия. Из них получаем ЮС(гк). После аппроксимации экспонентой получаем показатель Херста как Н = (1+а)/2, где а - тангенс угла наклона прямой на графике в двойном логарифмическом масштабе [2, 3, 5, 6].
Результаты на искусственных множествах
Для тестирования и оценки качества описанных выше алгоритмов использовались искусственные множества топологических размерностей 1 и 2 с заданной фрактальной размерностью. Для каждого множества в формате сейсмических каталогов подготовлены наборы данных с 64, 256, 1024, 4096, 16384, 65536 точками.
В пространстве топологической размерности 1 расчеты выполнены для набора множеств с размерностью 0,2-0,95, синтезированных на основе «канто-ровой пыли». На рис.1 представлены графики зависимости клеточной размерности О0 от количества точек N в выборках данных. Графики зависимости информационной й1 и корреляционной й2 размерности от N имеют аналогичный вид и поэтому не рассматриваются. На рис.1 видно, что максимальные отклонения от заданной размерности происходят при минимальном числе точек N = 64 и не превышают нескольких процентов. С увеличением N различие между заданной и вычисленной размерностью уменьшается, и при больших N размерности фактически совпадают. Можно отметить, что стандартное отклонение максимально для минимального количества точек N = 64 и повышается с ростом О0, а различие между заданной и вычисленной размерностью обычно не превышает одного стандарта.
Рис. 1. Графики зависимости клеточной размерности Do от количества точек N в выборках данных
В пространстве топологической размерности 2
использовались следующие синтезированные множества:
1. «Канторова пыль», расположенная на прямой диагональной осям координат, й0 = 0,63.
2. Диагональная прямая линия с равномерно расположенными вдоль нее точками, й0 = 1,0.
3. Диагональная прямая линия со случайно расположенными на ней точками, й0 = 1,0.
4. Две пересекающиеся диагональные линии (крест), проведенные аналогично п. 3, й0 = 1,0.
5. Триадная кривая Коха, й0 = 1,26.
6. Ковер Серпинского, й0 = 1,89.
7. Квадрат, заполненный случайно расположенными точками, и две точки вне квадрата по диагонали,
ния от заданной размерности происходят при минимальном числе точек N = 64. С увеличением N различие между заданной и вычисленной размерностью уменьшается, и при больших N размерности фактически совпадают. Наиболее значительно изменяется размерность множества из двух пересекающихся диагональных линий (креста) - от й2 = 1,79 при N = 64 до й2 = 1,04 при N = 65536. Как и в предыдущем примере, стандартное отклонение максимально для минимального количества точек и повышается с ростом 02, однако различие между заданной и вычисленной размерностью может превышать стандартное отклонение. Такие вариации наблюдаются на графике, построенном для двух пересекающихся диагональных линий (рис. 2, график 4).
Рис. 2. Графики зависимости корреляционной размерности D2 от количества точек N в выборках данных:
1-7 - синтезированные множества
чтобы стороны квадрата были строго внутри множества, й0 = 2,0.
На рис. 2 представлены графики зависимости корреляционной размерности й2 от количества точек N в выборках данных. Размерность й2 выбрана для представления из-за своей существенной изменчивости; графики клеточной О0 и информационной й1 размерности имеют в целом аналогичный вид, но их вариации менее значительны и поэтому не рассматриваются. На рис. 2 видно, что максимальные отклоне-
Выводы
1. Фрактальные оценки позволяют охарактеризовать структуру сейсмичности с точностью до второго десятичного знака.
2. Ошибка в определении растет с увеличением фрактальной размерности и усложнением топологии исследуемого объекта.
3. Примененные алгоритмы имеют преимущество при небольшом количестве исходных данных.
Статья поступила 23.12.2014 г.
Библиографический список
1. Божокин С.В., Паршин Д.А. Фракталы и мультифракталы. Ижевск: НИЦ «Регулярная и хаотическая динамика», 2001. 128 с.
2. Ключевский А.В., Зуев Ф.Л. Исследование динамики сейсмичности в Байкальском регионе // Доклады Академии наук. 2006. Т. 409. № 2. С. 248-253.
3. Ключевский А.В., Зуев Ф.Л., Демьянович В.М., Баяр Г. Оценки самоподобия энергетической структуры и динамики сейсмичности Монголии: тр. VI Российско-Монгольской конференции по астрономии и геофизике. Иркутск: Изд-во ИЗК СО РАН, 2006. С. 57-62.
4. Ключевский А.В., Зуев Ф.Л. Структура поля эпицентров землетрясений Байкальского региона // Доклады Академии наук. 2007. Т. 415. № 5. С. 682-687.
5. Ключевский А.В., Зуев Ф.Л. [и др.]. Подобие, самоподобие и синхронизация сейсмичности Монголо-Байкальского реги-
она // Современная геодинамика и опасные природные процессы в Центральной Азии. Иркутск: Изд-во ИЗК СО РАН, 2010. Вып. 6. С. 29-40.
6. Ключевский А.В., Зуев Ф.Л. Фрактальные оценки сейсмического процесса в литосфере Байкальского региона // Литосфера. 2011. № 1. С. 143-149.
7. Мандельброт Б. Фрактальная геометрия природы. М.: Изд-во ИКИ, 2002. 656 с.
8. Шелухин О.И., Тенякшев А.М., Осин А.В. Фрактальные процессы в телекоммуникациях. М.: Радиотехника, 2003. 480 с.
9. Klyuchevskii A.V. Rifting Attractor Structures in the Baikal Rift System: Location and Effects // Journal of Asian Earth Sciences. 2014. V. 88. P. 246-256.
10. Theiler J. Estimating fractal dimension // Journal of Optical Society of America, (A). 1990. V. 7. № 6. P. 1055-1073.
