CC BY
15
УДК 538.951 + 550.34.013
Характеристики и свойства динамической системы в диссипативной модели землетрясений Олами-Федера-Кристенсена
А.С. Черепанцев
Институт нанотехнологий, электроники и приборостроения, Южный федеральный университет, Таганрог, 347928, Россия
Приближение клеточных автоматов двумерной механической модели Burridge-Knopoff эволюции сейсмического процесса позволяет, наряду с традиционным рассмотрением дискретного потока сбросов, провести анализ временного ряда с эквидистантными отсчетами параметра энергии системы. Для динамической системы, формирующей наблюдаемые временные вариации энергии, проведены оценки ее инвариантных параметров. Полученные оценки значений параметров размерностей и энтропии динамической системы позволяют выделить диапазон значений параметра взаимодействия соседних элементов а, при котором поведение системы близко к критическому. Для данного диапазона а характерны инвариантность динамических параметров в широком диапазоне временных масштабов и связь степенных показателей распределений сбросов в пространственной и энергетической областях.
Ключевые слова: распределение сейсмичности в энергетической и пространственной областях, модель Олами-Федера-Крис-тенсена, самоорганизованное критическое состояние, динамическая система, корреляционная размерность, масштабное самоподобие
Characteristics and properties of dynamic system in the dissipative Olami-Feder-Christensen earthquake model
A.S. Cherepantsev
Institute of Nanotechnology, Electronics and Instrument Engineering, Southern Federal University, Taganrog, 347928, Russia
A cellular automaton approximation of the two-dimensional Burridge-Knopoff mechanical model of seismic process evolution can be applied, along with the conventional consideration of the discrete series of relaxations, to analyze a time series with equally spaced samples of energy system parameter. Invariant parameters for a dynamic system that forms observed time variations of energy are estimated. The obtained estimates for the dimension and entropy parameters of the dynamic system are used to determine the range of the interaction parameter of neighboring elements а at which the system behavior is close to critical. The given range of а is characterized by the invariance of dynamic parameters in a broad range of time scales and by the relation of relaxation distribution exponents in the spatial and energy domains.
Keywords: seismicity distribution in energy and spatial domain, Olami-Feder-Christensen earthquake model, self-organized critical state, dynamic system, correlation dimension, scalar self-similarity
1. Введение
В соответствии с концепцией самоорганизованного критического состояния сейсмического процесса [1, 2], система блоковых структур земной коры под действием активных тектонических процессов эволюционирует к статистически устойчивому состоянию со степенным поведением пространственной и временной корреляционных функций. При этом система приходит к самоорганизованному критическому состоянию вне зависимости от локальных особенностей блоковых структур и внешних воздействующих процессов. Наиболее
известным степенным распределением в сейсмичности является закон повторяемости землетрясений [3]. Он определяет связь числа наблюдаемых землетрясений с энергией большей Е0 как N(Е > Е0) = Щ0Е-В. Степенной показатель В в законе варьируется для широкого диапазона землетрясений в пределах 0.80-1.05.
В ряде модельных представлений сейсмичности получено распределение амплитуд сброса упругой энергии аналогичное наблюдаемому в реальном процессе. Модель блоков с упругими связями [4] предполагает, что поверхность разлома может быть представлена сис-
© Черепанцев А.С., 2015
темой блоков, упруго связанных друг с другом. Они расположены на основании с трением и медленно управляются внешней силой. В дальнейшем Z. Olami, H. Feder, K. Christensen [5] показали, что разработанная ими диссипативная модель (модель Олами-Федера-Крис-тенсена) клеточных автоматов эквивалентна двумерной модели блоков с упругими связями в квазистатическом приближении. Она сочетает в себе простоту построения и нетривиальный характер поведения. В отличие от начальных консервативных моделей возникновения самоорганизованного критического состояния, модель Олами-Федера-Кристенсена является диссипативной. Получаемая в модели оценка показателя степени в распределении амплитуд сбросов B оказывается близкой к наблюдаемой в сейсмичности.
Модель Олами-Федера-Кристенсена представляет собой двумерную клеточную решетку взаимодействующих элементов размером LxL. Для каждого элемента задан параметр, определяемый как «напряжение» или «энергия» Et (Et > 0). В начальном состоянии Et имеют случайные значения в интервале [0; 1]. На каждом шаге итерации, определяющем временной ход, Et возрастает на постоянную величину А до тех пор, пока в одной из ячеек величина напряжения не достигнет предельного значения Ek = Ec = 1. Тогда в k-й ячейке напряжение сбрасывается в нулевое значение, а часть сброшенного напряжения aEk (0<а<0.25) передается четырем ближайшим соседним ячейкам. Если для напряжения в соседней ячейке j за счет добавки aEk будет выполнено условие сброса (Ej > Ec = 1), то происходят также ее сброс и передача части сброшенного напряжения aEj ближайшим четырем соседям. Такая последовательность сбросов продолжается до тех пор, пока напряжения во всех ячейках не будут превышать предельного значения Ec = 1. После этого эволюция системы продолжается путем следующего приращения напряжений в ячейках (Et +А) до момента достижения в одной из ячеек предельного значения Ec. Сброшенные ячейки на одном шаге итерации, имеющие общую границу, формируют событие-сброс с амплитудой равной количеству таких ячеек. Модель является консервативной при a = 0.25 и неконсервативной при a < 0.25. В случае выбора открытых граничных условий система после конечного числа итераций переходит в устойчивое состояние, характеризуемое степенным распределением сбросов по размерам [5, 6].
Диссипация в системе определяется параметром a взаимодействия соседних ячеек. В ряде работ [7, 8] указывается на «ограниченную» критичность системы в неконсервативном режиме при a < 0.25. При этом выполнение соотношения конечно-размерного скейлинга оказывается справедливым лишь для сеток малых размеров L.
При анализе моделей сейсмичности, демонстрирующих критическое поведение, значительное внимание исследователей уделяется изучению потока дискретных событий сбросов-землетрясений, закономерностям получаемых распределений, сравнительному анализу структуры сбросов и особенностей реального сейсмического процесса, таких как афтершоковые и форшо-ковые последовательности. Вместе с тем практика изучения напряженного состояния земной коры включает в себя как контроль сейсмического режима, так и непрерывные наблюдения широкого круга геофизических полей, чувствительных к деформационным процессам. Модель Олами-Федера-Кристенсена позволяет получить непрерывные ряды параметра напряжения Е(/) в локальной области сетки и интегральной характеристики напряжения во всей системе. Выделение закономерностей поведения динамической системы, сформированной в модели Олами-Федера-Кристенсена, является в этом случае актуальным для понимания особенностей поведения временных рядов параметров различной природы, информативных в отношении напряженно-деформированного состояния геологической среды.
Разработанные к настоящему времени методы оценки инвариантных характеристик нелинейных динамических систем по временному ряду отдельной компоненты дают возможность оценить параметры размерности системы уравнений, сложности портрета системы в фазовом пространстве, степени предсказуемости ее поведения.
2. Особенности получения временных рядов в модели Олами-Федера-Кристенсена
Для оценки параметров динамической системы в качестве временного ряда выбраны вариации среднего значения напряжения на рассматриваемой решетке элементов [9] а = (Е) = 1/1'} £ Е\ j. Шкала времени
определяется величиной приращения напряжения А, предельным значением напряжения Ес и параметром диссипации а. При Ес = 1 шаг по времени составляет величину St = А/ (1 - 4 а). Выбор фрагмента временного ряда вариаций а(^) в режиме устойчивого предельного поведения динамической системы модели Олами-Фе-дера-Кристенсена предполагает проведение предварительного расчета переходного процесса достижения такого состояния. Как отмечено в [10], рост радиуса корреляции поведения соседних элементов растет со временем степенным образом: у(1)~ tв. Степенной показатель в при этом зависит от а и не зависит от L. С целью оценки времени достижения предельного состояния в системе размером L проведен расчет скорости роста величины максимального сброса в скользящем временном окне Аt ~ 3Т для различных значений а.
Рис. 1. Изменение величины максимального сброса с ростом времени эволюции системы L = 100 при различных значениях а. Окно сглаживания Ь = 0.07^
На основе расчетных зависимостей, представленных на рис. 1, рост £тах(^) может быть приближенно аппроксимирован как ^тах(а, г) = 3 г6 3а. Пунктирными линиями на рис. 1 показана данная аппроксимация для каждой расчетной зависимости с различным параметром связи.
Полученная оценка указывает на необходимость проведения длительного счета переходного состояния системы при малых а < 0.1. Для достижения предельных значений Smax ~ 4 • 103, полученных для а = 0.12-0.22 на сетке Ь = 100 (рис. 1), при а = 0.07 необходим расчет эволюции модели в течение t ~ 1.2 -10 Т. Данная оценка близка к оценке, полученной в [10] для числа расчетных сбросов N ~1010 при а = 0.07 и Ь = 128.
Другим важным параметром расчета при формировании временного ряда является выбор временного интервала дискретизации 5t и соответственно приращения энергии на отдельной итерации Д = (1 - 4а) 8г. Для информативности отсчетов временного ряда естественно предположить, что 5t не должно быть меньше среднего интервала времени между сбросами: 8г < T|nS, где nS — среднее число сбросов за период Т = 1. Зависимость nS (а) определяется различной величиной накло-
—В
на В(а) плотности распределения р(£)~ S [10]. Предполагая максимально возможный сброс Smax ~ L и условие однократного сброса каждой ячейки на сетке в течение периода Т = 1 - 4 а, для приближенной оценки 5t справедливо:
1
8г =
17
Е s
S=1
Г-В +1
|1LS -
S=1
(1)
Исходя из представленной в [9] зависимости В(а) изменение интервала дискретизации 8t в зависимости от параметра взаимодействия а представлено на рис. 2.
Данная зависимость является лишь приближенной оценкой, т.к. не учитывает временное группирование сбросов на старших масштабах, отмеченное в [9]. Улучшение временного разрешения временного ряда путем увеличения размера сетки Ь ограничено существенным увеличением временных затрат счета для достижения предельного устойчивого состояния системы.
Внимание к выбору интервала дискретизации при проведении расчета предельного состояния системы обусловлено и особенностью внутренней настройки уровней напряжений отдельных элементов Е при различных 8^ Для ускорения времени счета переходного режима можно рассмотреть алгоритм, в соответствии с которым необходимое время предварительного счета проводится с приращением напряжения Д1, существенно превышающим последующее приращение Д на этапе расчета временного ряда. На рис. 3 представлены фрагменты временных рядов, полученных по распределению напряжений после предварительного расчета с Д< Д1. Возникающая высокочастотная составляющая вариаций параметра среднего напряжения при расчете временного ряда с выбранным периодом дискретизации 8* = 0.18*1 указывает на существование дискретности распределений напряжений на сетке, кратных приращению напряжения на отдельном шаге временных итераций.
Подобная «настройка» значений напряжений отдельных элементов решетки для исходно случайно распределенной величины Е может быть показана в простейшем случае двух взаимодействующих ячеек. Пусть
Рис. 2. Зависимость В(а) (из работы [9]) и соответствующая зависимость 8^а) при размере сетки Ь = 150
Рис. 3. Фрагменты временного ряда среднего напряжения (Е) с различной величиной приращения напряжения на каждом шаге итераций Д = (1- 4а^ в модели Ь = 100, а = 0.07 после предварительного расчета в течение начального времени г = 2.5-105
Рис. 4. Распределение напряжений на сетке элементов при различных значениях а после начального расчета времени т: а = 0.07, L = 100, т =1.5 -105 (а); а = 0.12, L = 150, т =104 (б); а = 0.22, L = 500, т =1.3-103 (в)
Е1 (0) = 81 и Е2 (0) = 82 — два начальных случайных состояния элементов, таких что сброс одной ячейки не приводит к сбросу второй, т.е. [Е^к) = 81 + кА>1, [Е2(к) = 82 +k А< 1 [El(k) = 0,
[Е2 (к) = 82 + кА + а(81 + кА) < 1.
Если считать, что после этого на т-м шаге произошел сброс во второй ячейке, а далее после этого на ^м шаге — снова в первой, то
Е1(к + т + К) = (т + к)А + - + а[82 + (к + т)А + а(81 + кА)] >1, ^ Е2(к + т + К) = КА < 1
Е1(к + т + К) = 0,
^ 8 Е2 (к) = КА + а(т + К)А +
+ а2[82 + (к + т)А + а(81 + кА)] < 1.
Таким образом период сброса в первой ячейке составляет т + h итераций, т.е. (т + К)А = 1 - а. Тогда после г-го сброса первой ячейки, произошедшего на пЕ1 итерации,
Е2(г) = Ц -пЕ-!)А + а(1 - а) + а2(1 - а) +... + + аг-1 (1 - а) + аг[ 82 + (к + т) А + а( 81 + кА)] = = («Е - )А + а-аг + + аг [82 + (к + т) А + а(81 + кА)]. С учетом того что а < 1:
E2( i) = (nEi - n-1) A +
а.
(2)
Аналогичное соотношение справедливо и для напряжений в первом элементе: Е1(г) = (ПЕг - «Е-1)А + а.
Соотношение (2) указывает, что при малых а и отсутствии совместных сбросов на отдельной итерации напряжения в ячейках при достаточно большом числе итераций принимают значения кратные приращению внешнего напряжения на отдельной итерации А.
Для снятия выделенного выше ограничения по временному разрешению вариаций параметра (E) рассмотрена модель с приращением напряжения в каждой ячейке на отдельном шаге решетки А ^ 0. Для такой модификации модели основным механизмом формирования сброса является эффект взаимодействия соседних ячеек при передаче части напряжения а. При этом минимизируется эффект формирования макросброса за счет слияния границ отдельных сбросов с различными эпицентрами. Данный эффект может существенно искажать характер распределения сбросов по размеру при больших 8г и а. Выполнение условия А ^ 0 означает, что приращение на отдельной итерации не является постоянным значением, а определяется условием достижения элемента с максимальным значением Ei значения 1:
А = 1 -max Ei.
i
При этом напряжения в соседних ячейках после сброса соседнего элемента могут превышать единичное значение. Результат практического расчета модели в этом случае представляет собой значения отсчетов параметра (E) после сбросов, вызванных отдельным эпицентром и неэквидистантным шагом по времени btk = = Ак/(1 - 4 а), где k — номер сброса во временной серии. Такое представление может быть преобразовано к каталогу сбросов по амплитудам путем преобразования
Si =-
-«E(ti_i)> - <E(ti)> + (ti - ti_i)(1 - 4а)).
1 - 4а
В этом случае искажения полученных значений амплитуд возникают лишь для сбросов граничных элементов модели.
Распределение напряжений на сетке элементов после проведения предварительного расчета переходного процесса при различных значениях параметра диссипации а показано на рис. 4. Примеры полученных распределений, использованных в качестве карт начальных напряжений при расчете временных рядов вариаций (Е), представлены на рис. 5.
Возникновение сбросов групп соседних ячеек в рассматриваемой модели определяется возможностью за-
Рис. 5. Фрагменты временных рядов среднего напряжения при различных значениях параметра взаимодействия а
хвата соседних ячеек сформировавшимся кластером элементов, имеющих близкие значения напряжений. При этом скорость распространения таких захватов и роста кластера определяется величиной параметра связи а (рис. 1). Чтобы это показать, рассмотрим простейшую одномерную систему элементов. Пусть имеется группа п соседних элементов с величиной напряжений в них Е. В момент t она достигает предельного значения Ес = 1. Тогда после последовательных сбросов новые значения напряжений имеют вид Е1 =а(1 + а) = а + а2,
Е2 =а(1 + а(1 + а)) = а + а2 + а3,
(3)
Еп-1 =а + а2 +... + ап, Е = 0.
При этом соседняя (п + 1)-я ячейка, не входящая в кластер, получает приращение ДЕп+1 = а + а2 +... + ап.
Учтем, что а < 0.5 и при п >> 1 ДЕп+1 «а/(1 - а). В случае если
" (4)
Еп+1 >1 -
1-а
то (п + 1)-й элемент также сбрасывается, присоединяясь тем самым к кластеру: Еп+1 ^ 0, Еп ^ а. Ширина захвата возможных значений соседнего элемента определяется параметром связи а (4). В консервативном пределе а = 0.5, элемент с любым значением напряжения будет захвачен.
Рассмотрим устойчивость нового кластера размером п + 1 при следующем сбросе. Предполагая, что максимальное значение напряжения среди всех п + 1 элементов кластера имеет элемент с номером п - 1, последо-
вательность сброса, начиная с этого элемента, продолжится в сторону элементов с уменьшающимися индексами. Исходя из того что максимальная разность напряжений элементов в кластере (3) Еп-1 - Е1 = ап + +ап-1 +... + а3 и минимальное приращение ДЕк-1 при сбросе £-го элемента ^ = 2, ..., п - 1) ДЕк-1 > а, условие сброса всех элементов после сброса (п - 1)-го элемента: 1 -а3(1 + а + ... + ап-3) + а>1. Для а е [0; 0.5] неравенство справедливо, т.к. а 2/(1 -а) < 1. Для обоснования сброса элемента Еп, обусловленного сбросом Еп-1, учтем, что за интервал времени Т между сбросами кластера соседняя ячейка Еп+1 должна также предварительно сброситься и передать часть напряжения ячейке Еп: ДЕп >а.
Тогда условием сброса Еп в кластере является достижение значения
/ \
Е =1-
а
1-а
а
+ а>1,
которое справедливо для всего диапазона значений 0 <а<0.5. Аналогично рассуждая, можно показать, что и в случае захвата кластером при сбросе (п + 1)-го элемента при следующем сбросе кластера он также синхронно сбросится.
При малых значениях а, как показывают модельные расчеты, скорость формирования кластеров заданных размеров существенно ниже, чем наблюдаемая при а ^ 0.25 (рис. 1). Такое поведение может быть обусловлено как сужением диапазона возможных напряжений захвата (4), так и уменьшением скорости изменений напряжений в соседних элементах, обусловленным различной величиной диссипации энергии, т.е. если два кластера имеют различную среднюю величину приращения энергии, то через определенное время напряжение в кластере с более быстро изменяющимся циклическим процессом приблизится к напряжению в более медленном кластере и произойдет захват и объединение кластеров. Именно в областях максимальной изменчивости приращений напряжения в граничных областях решетки начинается захват элементов и формирование кластеров. Рассмотрим аналитическую оценку изменения приращения напряжений по мере удаления от границы одномерной модели размером N и при различных значениях параметра диссипации а.
Предполагая потери за счет превышения напряжений порогового значения Ес = 1 малыми, для числа сбросов в ячейке после прохождения времени т>> 1 справедливо: п1 = Е + п2 а, п2 = Е + п1а + п3а,
... (5)
пЩ-1 = Е + пЩ-2 а + пЩ а, пЩ = Е + пЩ -1а.
Рис. 6. Зависимость (7) скорости приращения напряжений в граничной области одномерной решетки при различных значениях параметра взаимодействия а
Системы уравнений (5) с трехдиагональной матрицей при N >> 1 представляет собой разностную форму неоднородного дифференциального уравнения:
d2«(х) 1 - 2^ Е
v - +-п (х) =--,
а а
ёх2 d«
ёп
х=0
1 -а
--п
а
1 -а
Ъ=0
х=Ь
а
« Т =~
1х=Ь
_ е
а
Е а
(6)
Решение этого уравнения представляется в виде
«(х)
= 1 --
е
-Хх - е_Х(ь -х)
-х (7)
Е/ (1 - 2а) Х + Р-(Х-Р)е-Хь'
где в = (1 - а)/а, Х = д/р-Г. Оно может служить оценкой средней скорости К = п/Т = п/[Е/(1 - 2а)], роста напряжения в элементе с координатой х = 0, 1, ..., Ь. На рис. 6 представлены зависимости изменения скорости приращения напряжений V(х, а), обусловленные внешним приращением Е и приращением, связанным с передачей напряжений соседними элементами при сбросе. Указанные зависимости демонстрируют, что с ростом параметра связи между состоянием соседних элементов, определяемого параметром а, увеличивается изменчивость АV в граничной области и расширяется зона такой изменчивости по х. Рост АV при увеличении а определяет в этом случае и большую скорость захвата кластерами с синхронным поведением ячеек соседних элементов и тем самым повышение скорости возникновения сбросов больших размеров. Полученное соотношение для скорости формирования сброса заданного размера в зависимости от параметра а сопоставлено с расчетными значениями. Оценки времени возникновения сброса заданного размера в одномерной модели проведены по аналогии с ранее рассмотренным двумерным случаем (рис. 1). Зависимость величины максимального сброса в скользящем временном окне в одномерной модели с начальным случайным распределе-
нием напряжений от времени аппроксимируется степенной зависимостью Smax(t, а) ~ гп. При этом показатель п имеет величину, меньшую, чем в двумерной модели, изменяясь от п = 0.32 при а = 0.07 до п = 0.75 при а = 0.45. Исходя из выбора степенной аппроксимации, время эволюции системы из исходного состояния в состояние с возникновением сброса величины £ определяется как
г (а) =
г ^ а)
С учетом того что при оценках п за единицу времени выбран период Т = 1/(1 - 2а), нормировка времени зависит от а. Поэтому для общего временного масштаба при сопоставлении параметров систем с различными значениями а можно определить время возникновения сброса размером £ как
„ \1/п(а)
г\а) = г (а)(1 - 2а) = (1-2а)
(8)
Вместе с тем оценку времени формирования сброса £ соседних элементов можно получить исходя из полученных выше соотношений для скорости приращения напряжений. Пусть путем захвата соседних элементов за счет разной скорости роста их напряжений сформирован сброс размера (£ - 1). Тогда время необходимое для захвата соседнего элемента и формирования сброса размером £ определяется разностью скоростей роста напряжений в кластере элементов размером (£ - 1) и захватываемого соседнего элемента. Если пг — число сбросов г-го элемента при величине приращения внешнего напряжения Е, то для Ес = 1 средняя скорость роста напряжений г-го элемента V = пгЕс/Е = щ/Е. В качестве приближенной оценки скорости роста напряжений синхронизованного кластера элементов рассмотрим среднее значение скоростей приращения напряжений отдельных элементов, входящих в кластер: = = п/Е. Тогда верхняя оценка времени формирования сброса размером £ определяется временем т, необходимым для достижения равенства напряжений кластера и соседнего элемента. Оно равно отношению максимально возможной величины расхождения АЕ = 1 и разности средних скоростей роста напряжений:
1Е
т =-= (9)
(пс1 - п VЕ пс1 - п
На рис. 7 представлена оценка времени достижения сброса величиной £ = 30, исходя из соотношений (8) и (9). Следует заметить, что для £ < 30 теоретическая оценка времени дает завышенные величины для а < 0.3, а для £ > 30 теоретическая оценка в том же диапазоне а оказывается несколько заниженной. Это может быть связано как с приближенностью выбранной оценки средней скорости роста напряжений в кластере, так и ограниченностью аппроксимации степенной зависимо-
0.0 0.1 0.2 0.3 0.4 а
Рис. 7. Зависимость времени достижения сброса величиной £ = 30 от значения параметра взаимодействия а между элементами линейной модели. Расчет проведен для модели Ь = 500, Е = 105
стью полученных модельных расчетов времени достижения сброса заданной величины.
3. Параметры динамической системы в модели Олами-Федера-Кристенсена
Оценки основных параметров динамической системы, формируемой в модели Олами-Федера-Кристенсена, проведены по данным значений среднего напряжения (Е) для режима предельного устойчивого поведения модели со степенным распределением сбросов по величине.
При отсутствии взаимодействия элементов динамической системы размерность системы определяется количеством таких элементов. Однако наличие взаимодействия между ними может приводить, как отмечалось выше, к синхронизации поведения соседних элементов, образованию кластеров и соответственно уменьшению размерности динамической системы. Как следует из рис. 5, характер вариаций (Е) зависит как от параметра взаимодействия а, так и от выбранного временного масштаба рассмотрения. В работе [9] на основе анализа коэффициента вариации временных интервалов между сбросами различной амплитуды показано существование группирования сбросов большой амплитуды при параметре взаимодействия а = 0.2 в модели Ола-ми-Федера-Кристенсена. Учет временного масштаба обусловлен и наличием периодической составляющей на масштабах процесса т ~ 1 - 4а при малых а.
Вид временных рядов указывает на ослабление периодической составляющей с периодом Т = 1 - 4а при росте а. Подобная закономерность ослабления амплитуды спектрального максимума при анализе временного хода сбросов указана в [11]. При приближении эволюции модели к режиму консервативной системы (а = = 0.25) растет хаотизация процесса. Для описания систем с хаотическим поведением широко используются подходы нелинейной динамики. Они позволяют по огра-
ниченным данным поведения динамической системы получить основные ее характеристики.
Под динамической системой понимается система п уравнений вида НХ
— = G (X), Х(г) = {*(г), Х2(г),..., Хп (г)}, (10) Нг
определяющая динамику во времени п независимых компонент, в общем случае нелинейно взаимодействующих, что определяется правой частью системы G(X). Основной интерес представляют закономерности поведения (10) в фазовом пространстве, включающем все возможные состояния системы. При этом разным состояниям соответствуют различные точки в нем. Путем непрерывных и однозначных преобразований параметров Х1(г), х2(г),..., хп(г) можно получить иное от (10) описание того же процесса. Фундаментальным свойством рассматриваемой системы является инвариантность ряда ее параметров при таких преобразованиях [12].
Теория динамических систем указывает на возможность построить образ неизвестной динамической системы (10) по временным вариациям единственного параметра хк ^ = 1, ..., п), имеющей те же инвариантные параметры, что и система (10). Как правило, для построения образа используется преобразование
X(г) = {Х1(г), Х2(г),..., Хп (г)} ^ ^ У (г) = {Хк (г), Хк (г + Д),..., Хк (г + (п - 1)Д)}.
Оцениваемые инвариантные параметры динамической системы включают размерность динамической системы п, корреляционную размерность d2, являющуюся нижней оценкой размерности Хаусдорфа в фазовом пространстве системы (d2 < п) и корреляционную энтропию К2, являющуюся нижней оценкой колмогоровс-кой энтропии и характеризующую степень предсказуемости поведения системы.
В соответствии с процедурой, предложенной в [13], оценку размерности фазового пространства, корреляционной размерности и корреляционной энтропии можно получить путем расчета корреляционного интеграла — функции распределения всевозможных межточечных расстояний амплитуд имеющегося временного ряда:
С2(г) = Х9(г-1 У(г,) - У(гj )|), (11)
N I
где 9() — единичная функция Хевисайда.
С целью анализа поведения динамической системы на различных временных масштабах исходные временные ряды отдельных сбросов преобразовывались к временным рядам средней энергии в системе с различными интервалами дискретизации Дt. Нижнее ограничение определялось средней величиной времени между последовательными сбросами. Верхнее ограничение определялось характерным периодом Т = 1 - 4а. Для
Рис. 8. Оценка размерностей d2, n по корреляционному интегралу C2(г), который рассчитан по выборке объемом N = 3 • 104 точек при a = 0.22, L = 500, компоненты ряда (E)(t) во временном диапазоне (0.0067-0.067)T. Фрагмент С2(г) для различных уровней вложения m = 1, ..., 14 (а); изменение наклона С2(г) в диапазоне амплитуд (0.7-1.5) • 102 г/ а (б)
получения выборки заданного временного масштаба исходные временные ряды подвергались фильтрации полосовым фильтром Баттерворта 16-го порядка и временному сжатию с интервалом дискретизации Аt = = 0.3//ир, где /ир — верхняя частота отсекания фильтра. Для компенсации фазовых искажений данная процедура представляла собой последовательное применение фильтров 8-го порядка к временному ряду при прямом и обратном временном ходе.
На рис. 8 представлен пример оценки размерностей d2, п по расчетным кривым корреляционного интеграла С2(г). Процедура оценки d2, п включает в себя последовательный расчет выборочной функции распределения межточечных расстояний (11) между отсчетами У(^) = {хк), хк(г,. + А),..., хк(г,. + (т -1)А)} при увеличении уровня вложения т. Величина т, начиная с которой линейный наклон у выборочной функции распределения, представленной в двойном логарифмическом масштабе, перестает изменяться, служит оценкой величины п (рис. 8, б). При этом показатель степени у зависимости С2(г) ~ гу является оценкой d2 (рис. 8, а).
Проведенный расчет размерностей позволяет оценить их зависимость от параметра взаимодействия а на различных временных масштабах т (рис. 9). В качестве т выбрано среднее геометрическое значение частотного окна полосового фильтра т = (/тп /тах)-^2, применявшегося при получении анализируемой выборки.
На основе временных вариаций среднего по элементам значения напряжения проведена оценка инвариантного параметра динамической системы корреляционной энтропии К 2. Обобщенный корреляционный интеграл С2 (г) определяет вероятность попадания в ячейку размером г групп из последовательных 5 точек фазовой траектории. При 5 = 1 он тождественен корреляционному интегралу (11), используемому при оценке d2 и п. Оценка К2 получается исходя из справедливости соотношения К2 = Нт1п[С5 (г)/С2+1( г)]/Аг, где Аt —
г ^0
интервал временной дискретизации временного ряда. Значение К2 позволяет получить оценку среднего времени предсказуемости поведения системы, которое определяется величиной скорости разбегания двух близ-
Рис. 9. Связь размерностей й2 (а) и п (б) с временным масштабом т вариаций напряжений в модели Олами-Федера-Кристенсена при различных значениях параметра взаимодействия а. п — ширина временного окна отдельной компоненты
Рис. 10. Зависимость параметра К2Дг от временного масштаба при различных параметрах взаимодействия а модели Олами-Федера-Кристенсена
ких траекторий (Iг(tt)-г(tj)1 < Дг) динамической системы £~ln(1/ Дг )/K2.
Значения K2 получены для различных временных масштабов вариаций напряжения и параметров взаимодействия а выбранной модели. Результат проведенных расчетов параметра K2Дt представлен на рис. 10.
На основе полученных оценок можно выделить особенности поведения динамической системы при различных значениях а. При слабом взаимодействии между соседними элементами (а < 0.17) поведение динамической системы на старших и младших временных масштабах имеет различный характер. Размерности d2, n уменьшаются с ростом масштаба т рассмотрения процесса во времени. Это сопровождается ростом де-терминизации процесса, что отражает уменьшение энтропии K2 динамической системы, приведенной к единичной временной шкале Дt = 1.
Поведение динамической системы с параметром связи а в диапазоне значений 0.175 < а < 0.220 демонстрирует близкое к постоянному значению размерностей d2, n при значениях т < (1 - 4а). При а = 0.22 динамическая система рассматриваемой модели в диапазоне масштабов т = 2 • 10-3 - 2 • 10-1 демонстрирует также и постоянное значение K2 Дt. Такое поведение указывает на самоподобный характер поведения динамической системы во временной области. При масштабном временном преобразовании характер поведения системы остается постоянным. Для указанного диапазона а оценка среднего времени предсказания поведения системы определяется лишь временным масштабом рассмотрения £ ~ const • Дt.
При приближении к консервативному пределу а = = 0.25 характерно увеличение d2 с ростом т. При этом размерность системы n и параметр K2Дt близки к постоянным значениям на рассматриваемых временных масштабах.
Полученное уменьшение размерностей d2, п с ростом временного масштаба при малых значениях а < 0.17 указывает на различный характер динамических процессов на рассматриваемых временных масштабах. Младшему временному масштабу соответствует хаотический процесс. С ростом же временного масштаба он стремится к периодическому детерминированному процессу, когда система переходит в исходное состояние через временной интервал Т = 1 - 4 а. В диапазоне значений параметра взаимодействия 0.175 < а < 0.220 оценки размерностей указывают на переход динамической системы в состояние, близкое к масштабному самоподобию, характеризуемому общим поведением системы на всех рассматриваемых временных масштабах. Параметры остаются близкими к постоянному значению d2 = 6.5 ± 0.5, п = 11 ± 1 в рассматриваемом диапазоне масштабов. Именно поведение модели Олами-Федера-Кристенсена в данном диапазоне параметра а рассматривается в большинстве работ по анализу самоорганизованного критического состояния [5, 6]. Это обусловлено прежде всего тем, что для данного диапазона параметра диссипации а значение наклона функции распределения В ~ 0.91 оказывается близким к наблюдаемым значениям в сейсмичности [6]. Как показано в [5], значение а = 0.2 соответствует выполнению условия изотропности упругих свойств взаимодействующих блоков в модели сейсмичности Burridge-Knopoff [4]. Рост размерности d2 с увеличением временного масштаба при значениях а, приближающихся к консервативному пределу а = 0.25, демонстрирует рост сложности динамической системы, включающей сбросы большой амплитуды.
Особенностью процесса сброса напряжений при а ^ 0.25 является рост числа повторных сбросов элементов в течение формирования отдельного макросброса. Это нарушает характер соотношения между площадью сброса и его величиной и, соответственно, может усложнить поведение динамической системы. Рассмотрим возможность возникновения множественных сбросов в одном и том же элементе на отдельном шаге итерации.
Элемент, удаленный от границ, после собственного сброса и вызванной этим лавины сбросов соседних ячеек изменяет напряжение на величину
ДЕ = 1 -Е (+а-1),
]
где Ej — значение энергии в соседних ячейках до сброса центральной ячейки, или, определив отклонение от Ес = 1 через 51 = 1 - Е,,
(12)
ДЕ = 1 -4 а(1 + а) + аУ5,.
В случае ДЕ < 0 ячейка осуществляет повторный сброс, которых может быть несколько, до тех пор пока не будет выполнено условие ДЕ > 0. Получим оценку значения
Р--
0.50.40.30.20.10.0 ^—9--1-1-1-1-10.20 0.22 0.24 а
Рис. 11. Доля повторных сбросов при различных значениях параметра взаимодействия а в модели Олами-Федера-Крис-тенсена на сетке размером Ь х Ь = 40 х 40
а, начиная с которого данный эффект будет значимо присутствовать. В качестве верхней оценки 8 выберем 8 тах = 1 - 4а, определяющей максимально возможное отклонение в ячейке после сброса четырех соседних ячеек. Тогда справедливо: а> 1Д/20 ~ 0.224. На рис. 11 показан модельный расчет доли повторных сбросов среди их общего числа при различных значениях а. Учет повторных сбросов можно определить путем задания дополнительной степени свободы в описании модели и перехода от двумерного распределения сбросов на сетке к трехмерному: S(х, у) ^ S(х, у, а ), где а — число сбросов. Такой рост размерности динамической системы может сопровождаться и более сложным ее поведением. Увеличение значения корреляционной размерности d2 с ростом временного масштаба при а = 0.245 демонстрирует повышение размерности фазового портрета динамической системы (рис. 9, а).
Рост хаотичности вариаций динамической системы при приближении а к консервативному пределу может определяться и изменением механизма диссипации энергии в системе. Наряду с механизмом диссипации во внутренних ячейках при сбросе присутствует и диссипация на границе решетки с иной величиной потерь при сбросе граничного элемента. Несложно оценить, что доля диссипации энергии в граничных элементах на решетке элементов размером Ь х Ь = 500 х 500 при а = 0.07 составляет лишь 0.8 %, увеличивается до 10 % при а = 0.245 и достигает 100 % при а = 0.25.
Условия самоподобия динамических процессов на различных амплитудно-временных масштабах и «плоская» модель сбросов на решетке элементов с малым числом повторных сбросов, наблюдаемые при значениях параметра взаимодействия а ~ 0.2, позволяют обосновать существование связи степенных показателей пространственного и амплитудного распределений сбросов в модели Олами-Федера-Кристенсена.
Возникновение самоорганизованного критического состояния в диссипативной модели Олами-Федера-
Кристенсена принято характеризовать степенным видом распределения сбросов по размерам. Предполагая отсутствие повторных сбросов, определим размер сброса как количество связанных (имеющих общую границу) ячеек, испытывающих сброс на отдельном временном шаге приращения энергии. Функция распределения амплитуд сбросов, имеет вид
Е (5) ~ S- в
Как показано в работах [11, 14], распределение эпицентров сбросов в пространственной области двух- и трехмерной неконсервативной модели Олами-Федера-Кристенсена имеет фрактальную организацию со степенным характером распределения расстояний между сбросами. Оценка фрактальной размерности для трехмерного случая составляет В2 -1.8 ± 0.1 [14]. Наблюдаемое в реальном сейсмическом процессе пространственное распределение эпицентров предложено [15] соотнести с фрактальной размерностью В длин разломов, приуроченных к областям сейсмичности N ~ I- п. В этом случае, исходя из представлений самоподобного характера сейсмического процесса в пространственной и энергетической областях и соотношения сейсмического момента и магнитуды землетрясения, для связи В и В справедливо соотношение Акк В = 2В [16]. Оно связывает организацию сейсмичности в пространственной и энергетической областях. Для рассматриваемой модели Олами-Федера-Кристенсена, предполагая, что величина сброса имеет геометрический смысл и определяется площадью сброшенной области, также можно определить соотношение параметров пространственного распределения сбросов и их размеров.
- в
Пусть N(5 > £0) = С150 — наблюдаемая зависимость числа сбросов больше £0 при однократном цикле сбросов в каждой ячейке. В двумерной модели при отсутствии повторных сбросов величина сброса численно равна площади области сброса. В общем случае зависимость площади сброса от линейного размера £(К) можно определить как £ ~ Я 8. Пусть в качестве линейного размера (размера очага) определено максимальное линейное расстояние между ячейками, принадлежащими данному сбросу. Тогда для распределения числа сбросов справедливо
N(£ > £0) = С2Я~в8. (13)
При рассмотрении пространственного распределения эпицентров сброса на заданной сетке эпицентром будем считать среднее значение координат сброшенных ячеек. Учтем, что наличие в одной ячейке размером I х I (в двумерном случае) двух и более эпицентров означает, что величины соответствующих сбросов не превышают £. < I . Тогда число эпицентров >2 в ячейках размером I определяет число сбросов, не превышающих £0 < 12: Й(£ < £0) = Й(£1: |г - г]| < I). В случае фрактального распределения эпицентров
Рис. 12. Степенные показатели амплитудного и пространственного распределения сбросов при различных значениях параметрах взаимодействия а модели Олами-Федера-Кристенсена. Показатель В функции распределения амплитуд сбросов, показатель d функции распределения межточечных расстояний между эпицентрами сбросов, показатель 5-распределения площади сброса по его линейному размеру (а); изменение параметра d - В5 с ростом параметра взаимодействия а модели Олами-Федера-Кристенсена (б)
N ^: | г - г] | < /)~ ^. (14)
В случае двумерной решетки d < 2.
Вид полученных функций распределения (13), (14) определяет связь степенных показателей:
d = В 5. (15)
В соответствии с (15), амплитудное распределение сбросов определяется фрактальной размерностью очага сброса и распределением областей сброса в пространстве.
Результаты оценки Ь, ^ 5 в модели Олами-Федера-Кристенсена для различных значений параметра взаимодействия а показаны на рис. 12, а. При оценке фрактальной размерности эпицентров d использована методика расчета корреляционного интеграла (11), где Уг- = = {х, у(} — координаты эпицентров и С2 ~ И [14].
Величина наклона графика повторяемости В в рассмотренной модифицированной неконсервативной модели Олами-Федера-Кристенсена уменьшается с ростом параметра взаимодействия а аналогично поведению данного параметра в исходной модели Олами-Фе-дера-Кристенсена с фиксированным ростом напряжений на каждом шаге итераций [9]. При этом характерен меньший диапазон изменения В, расширяющий диапазон возможных значений а для наблюдаемых в сейсмичности значений наклона графика повторяемости землетрясений [17]. Диапазон В = 0.80-1.05 соответствует диапазону а = 0.12-0.22.
Наличие связи В, d в фоновой сейсмичности отмечена в [18]. В рассматриваемой модели фрактальная размерность эпицентров сбросов d(а), подобно В(а), уменьшается с ростом параметра взаимодействия а. Это может указывать на переход системы от состояния слабо связанных групп сбросов к формированию синхронизованных сбросов большой величины. Предельным
максимальным значением размерности может служить значение d=2, соответствующее равномерному случайному распределению эпицентров на поверхности. Величина размерности d в диапазоне а, соответствующем B, наблюдаемому в сейсмичности, также имеет значения d = 1.5-1.7, соответствующие наблюдаемому распределению эпицентров землетрясений.
Оценка фрактальной размерности очага сброса по соотношению площади сброса и его линейного размера 5 = 1.8-1.9 демонстрирует отсутствие зависимости 5(а) при выбранной геометрии передачи части напряжений при сбросе соседним элементам. Данный результат может служить подтверждением предложенной K.Aki [15] прямой связи d = 3 В/c, где c = const.
Для анализа условий выполнения соотношения (15) на рис. 12, б представлена зависимость изменения значения d - В5 при различных значениях параметра взаимодействия а. Как следует из представленной кривой, поведение модели с наклоном графика повторяемости, соответствующего наблюдаемому в сейсмичности, отвечает приближенному (в пределах погрешности) выполнению условия (15). Как следует из рис. 12, б, при приближении модели к консервативному пределу (а = = 0.245) нарушение соотношения (15) обусловлено слишком малым значением B, а в случае слабого взаимодействия элементов (а = 0.07) нарушение (15) обусловлено слишком большим значением B. В рассматриваемой модели B является наиболее чувствительным параметром, характеризующим состояние взаимодействия элементов. Возможной причиной нарушения (15) при а ^ 0.25 может быть то, что рассчитываемая фрактальная размерность площади очага 5 не учитывает повторяющиеся сбросы в одной и той же ячейке, при этом они учитываются в величине сброса напряжения при оценке B.
4. Заключение
В работе рассмотрена возможность описания эволюции диссипативной модели Олами-Федера-Кристенсе-на как динамической системы со сложным хаотическим поведением на различных временных масштабах. Оценка инвариантных параметров динамической системы проведена по временному ряду средней энергии в системе. С целью повышения временного разрешения моделируемого временного ряда предложена модификация модели с различными значениями приращения энергии на отдельных итерациях. Для обоснования достижения предельного устойчивого состояния эволюции модели проведена оценка времени достижения сбросов заданной величины при различных значениях параметра взаимодействия а. Для простейшей двумерной модели рассмотрен механизм формирования кластера элементов с синхронным поведением. Показано, что различные величины градиента скорости приращения энергии при различных значениях а могут определять скорость формирования предельного критического состояния. Полученные оценки размерности динамической системы n, корреляционной размерности d2 и энтропии K2 позволяют выделить диапазон значений параметра связи: 0.175 <а<0.22, при котором динамическая система демонстрирует самоподобный характер на различных временных масштабах. Показано, что для данного диапазона справедливо соотношение связи между фрактальной размерностью эпицентров D, показателем B наклона графика повторяемости и фрактальной размерностью очага 8: d = B 8.
Работа выполнена при поддержке гранта РФФИ № 14-05-00521.
Литература
1. BakP., Tang C. Eartquakes as a self-organized critical phenomenon // J. Geophys. Res. B. - 1989. - V. 94. - No. 15. - P. 635-637.
2. Sornette A., Sornette D. Self-organized criticality and eartquakes // Europhys. Lett. - 1989. - V. 9. - P. 197-202.
3. Gutenberg B., Richter C.F. Magnitude and energy of eartquakes // Ann. Geophys. - 1956. - V. 9. - P. 1-15.
4. Burridge R., Knopoff L. Model and theoretical seismicity // Bull. Seism.
Soc. Am. - 1967. - V. 57. - P. 341-371.
5. Olami Z, Feder H.J.S., Christensen K. Self-organized criticality in a continuous, nonconservative cellular automaton modeling earthquakes // Phys. Rev. Lett. - 1992. - V. 68. - P. 1244-1247.
6. Lise S., Paczuski M. Self-organized criticality in a nonconservative earthquake model // Phys. Rev. E. - 2001. - V. 63. - P. 36111.
7. de Carvalho J.X., Prado C.P.C. Self-organized criticality in the Olami-Feder-Christensen model // Phys. Rev. Lett. - 2000. - V. 84. - P. 4006.
8. Miller G., Boulter C.J. Measurements of criticality in the Olami-Feder-
Christensen model // Phys. Rev. E. - 2002. - V. 66. - P. 016123.
9. Christensen K. Self-Organization in Models of Sandpiles, Earthquakes and Flashing Fireflies: PhD Thesis. - Oslo: Oslo University, 1992.
10. Middleton A.A., Tang C. Self-organized criticality in nonconserved systems // Phys. Rev. Lett. - 1995. - V. 74(5). - P. 742-745.
11. Peixoto T.P., Prado C. Statistics of epicenters in the Olami-Feder-Christensen model in two and three dimensions // Physica A. - 2004. -V. 342(1-2). - P. 171-177.
12. Schuster H.-G. Deterministic Chaos: An Introduction. - Weinheim: Physik Verlag, 1988.
13. Grassberger P., Procaccia I. On the characterization of strange attractors // Phys. Rev. Lett. - 1983. - V. 50. - P. 346-349.
14. Jansen F., Hergarten S. Basic properties of a three-dimensional springblock model with long-range stress transfer // Phys. Rev. E. - 2006. -V. 73. - P. 026124.
15. Aki K. A Probabilistic Synthesis of Precursory Phenomena // Earthquake Prediction / Ed. by D. Simpson, P. Richards. - Washington: Am. Geophys. Union, 1981. - P. 566-574.
16. LegrandD. Fractal dimensions of small, intermediate and large earthquakes // Bull. Seism. Soc. Am. - 2002. - V. 92. - P. 3318-3320.
17. Pacheco J.F., Scholz C.H., Sykes L.R. Changes in frequency-size relationship from small to large earthquakes // Nature. - 1992. -V. 355. - P. 71-73.
18. Caneva A., Smirnov V. Using the fractal dimension of earthquake distributions and the slope of the recurrence curve to forecast earthquakes in Colombia // Earth Sci. Res. J. - 2004. - V. 8. - No. 1. -P. 3-9.
Поступила в редакцию 10.06.2015 г.
Сведения об авторе
Черепанцев Александр Сергеевич, к.ф.-м.н., доц. ЮФУ, [email protected]
