CC BY
42
В. М. Урицкий, Н. А. Смирнова, В. Н. Троян
КРИТИЧЕСКАЯ ДИНАМИКА ФРАКТАЛЬНЫХ СИСТЕМ И ЕЕ ВОЗМОЖНАЯ РОЛЬ В ФОРМИРОВАНИИ ЗОН РАЗЛОМОВ И ГЕНЕРАЦИИ ЭЛЕКТРОМАГНИТНЫХ ПРЕДВЕСТНИКОВ ЗЕМЛЕТРЯСЕНИЙ*)
Введение. Теория самоорганизованной критичности (СОК)[1, 2] - один из наиболее универсальных подходов к описанию масштабно-инвариантных процессов и структур в природных системах [3]. К настоящему времени разработан ряд моделей СОК, описывающих степенную статистику энергий землетрясений в различных по строению сейсмоактивных зонах [4]. Общий недостаток большинства таких моделей заключается в их неспособности отразить эффект пространственно-временной фрактальной кластеризации сейсмических событий [5]. Первой моделью СОК, позволившей исследовать связь степенной формы плотности вероятности энергии землетрясений с фрактальной динамикой матрицы тектонических разломов, стала модель Д. Хьюза и М/Пацзус-ки [6]. Ключевым ее компонентом является нарушенная абелева'симметрия локальных взаимодействий [7], что приводит к появлению долговременной памяти и крупномасштабным корреляциям флуктуаций решетки. Возникающая под действием СОК-неустойчивостей стохастическая структура, состоящая из субкритических элементов модели, определяет наиболее вероятные направления распространения возмущений. Она подвержена медленной самосогласованной эволюции, обеспечивающей сохранение критического поведения модели в широком диапазоне значений управляющих параметров.
Авторы модели [6] высказали предположение о возможности ее использования для изучения связи статистики землетрясений с крупномасштабной динамикой тектонических систем. В настоящей работе проведен количественный анализ корреляционной структуры модели Хьюза-Пацзуски, позволяющий оценить перспективы ее применения к описанию сейсмических явлений. Исследуются энтропия решетки модели, фрактальная структура временных рядов параметра порядка, а также зависимость энерговыделения от скорости измененияч перколяционной проводимости. Рассматривается также возможность моделирования сейсмоэлектромагнитных предвестников землетрясений, связанных с процессами фрактального трещинообразования.
Описание модели. Исследованная нами версия модели Хьюза-Пацзуски определена на 2-мерной ромбической решетке, каждый узел которой описывается целочисленными координатами ж € 0, 1,— 1, у = 0, 1,..., — 1 и энергией г(х, у) е Д. Если в момент времени
4 локальное значение энергии превышает порог развития неустойчивости г — гс, некоторая ее
часть, обозначенная (1г, передается двум соседним элементам (рис. 1, а):
г4+1 (х,у) = 2* (х,у) - <1г, '
г±+1(х + 1,у + 1) = гг(х + 1,у + 1) +рйг, (1)
24+1 {х,у + 1) =* гг (х,у + 1) + (1 -р) йг. * '
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-05-64938) и программы «Интергеофизика».
© В. М. Урицкий, Н. А. Смирнова, В. Н. Троян, 2004
40СМ
Рис. 1. Фрагмент решетки 2-мерного клеточного автомата Хъюза-ПаЦзуски, поясняющий механизм перемещения лавины неустойчивостей с энергией г > 2С (большие черные кружки) по сети субкритических элементов с 2: € [гс — <1г, гс] (маленькие черные кружки) (а) и характерный пример пространственного распределения энергии в неабелевой версии модели (б) (черные области решетки соответствуют активным1 элементам (лавинам), серые - субкрити-ческим элементам).
Здесь р - случайная величина с однородным распределением в интервале [0,1]. На следующем шаге каждый из получивших энергию элементов может также перейти в активное состояние.
. что приводит к возникновению лавины неустойчивостей, распространяющихся сверху вниз в направлении оси у. Лавина останавливается, достигнув нижней поглощающей границы у — Ny — 1 или попав в область решетки с низкими значениями z. На левую и правую границы решетки наложено периодическое условие z (0, у) = z{Nx — 1, у). В описанных далее численных экспериментах использовался следующий набор параметров: Nx = Ny = 400, h = 1, Zc — 1.
Неравновесное состояние модели поддерживается добавлением малого по сравнению с порогом 2С количества энергии к случайно выбранным элементам Верхнего ряда (у = 0), которое продолжается до тех пор, пока один из них не потеряет устойчивость (г > zc), что служит стартовой точкой лавины. Критическая динамика модели выражается в степенной форме распределения лавин по размеру .s и времени жизни Т: .
: ' p{s) = а_т* f{s/sc),
■ ’ - (2)
р(Т) = Т ~Tt д (Т/Тс),
где / и д - скейлинговые функции, описывающие закон спадания распределений в окрестности пороговых значений s = sc и Т = Тс; т3 и rt - критические индексы лавин [8]. Степенная статистика лавин в моделях СОК, как правило, не предполагает наличия корреляций в поведении элементов модели на масштабах времени, превышающих длительность отдельных лавин. . Однако лавинная динамика, описываемая соотношениями (1), при определенных условиях приводит к возникновению дальнего порядка между элементами модели и способствует ее крупномасштабной самоорганизации на периодах времени, значительно превышающих время жизни неустойчивостей. .
Ключевым параметром модели, определяющим наличие или отсутствие крупномасштабных корреляций, является количество энергии dz, перераспределяемой при взаимодействии элементов [6]. В случае dz = const клеточный алгоритм (1) становится абелевым, т.е. не зависящим от порядка обновления переменной решетки [7], что вызывает пространственную хаотизацию субкритических состояний модели, определяемых из условия 2 ^ 2С. В случае, когда dz - переменная величина, равная или пропорциональная энергии возбужденного элемента, происходит нарушение абелевой симметрии и поведение системы становится сложно коррелированным. При этом модель приобретает масштабно-инвариантную структуру, образованную фрактальной сетью взаимосвязанных субкритических элементов (рис. 1, б).
Результаты и их обсуждение. Для количественного изучения фрактальных корреляций в неабелевом режиме модели была исследована масштабная зависимость информационной энтропии 5 распределения состояний решетки (в битах)
/
00 ' .
' 8(1}'=- £&{3)\о&2(3)
о ,
энергии решетки, оцененного в ячейках с длиной стороны I. Принимая гауссову модель распределения, что оправдано при достаточно большом числе ячеек, приходим к приближенному соотношению
5(1) *» ^1о§2(27ге сгр), (4)
в котором в качестве сгг используется стандартное отклонение значений полученное при усреднении по всем ячейкам размера I. . .
В абелевой версии модели величина г(х,у) как функция пространственных координат имеет вид некоррелированного шума, описываемого соотношением <Т{ ос I ~ а/2, где ё — 2 - размерность решетки, что после подстановки в (4) дает
вое-^1о?2!= -^^108,01. (5)
Такое поведение энтропии означает, что при больших I флуктуации энергии усредняются и становятся пренебрежимо малы. В неабелевой версии модели влияние крупномасштабных корреляций должно приводить к значительно более медленному спаданию 5(/), поскольку в этом случае усреднение флуктуаций на все большем масштабе не гарантирует снижение <т( и сходимость средних оценок энергии г. Действительно, было обнаружено, что режим с сохраненной абелевой симметрией отличается быстрым снижением 5 при увеличении размера ячеек. Угловой коэффициент наклона этой зависимости, построенной в полулогарифмическом масштабе, хорошо согласуется с теоретическим значением —3,332, полученным по формуле (5). Режим с нарушенной абелевой симметрией характеризуется значительно более плавным уменьшением энтропии, значение которой на несколько порядков выше, чем в абелевом режиме (рис. 2, а). Образование крупномасштабных структур в неабелевой версии модели, вызывающее аномальный скейлинг информационной энтропии, является одним из проявлений самоорганизации критического состояния модели. Для иллюстрации этого эффекта на рис. 2, 5 представлен пример начальной релаксации фрактальной сложности модели
^гпах .
Е{1) = I 5(/, «)сй. ' (6)
Исходное состояние характеризовалось равномерным случайным распределением энергии, среднее значение которой выбиралось близким к величине средней энергии конечного стационарного состояния. При этом условии наблюдаемая динамика определяется главным образом структурными изменениями в системе. Как следует из рис. 2, б, рост фрактальной сложности в неабелевом режиме продолжается существенно дольше, что связано с более продолжительным временем, требуемым для установления крупномасштабных корреляций. При этом зависимость Е(4) имеет существенно более гладкий
20 40 60 80 ИЮ 120 МО Ш> 1 $0 ^
О
>хт ч [ см кю
Рис. 2. Распределение информационной энтропии 5 по масштабам в двух режимах модели (а) и начальная релаксация фрактальной сложности Е(Ь) при формировании стационарного состояния СОК (б). :
1 - неабелева версия модели, 2 - абелева.
вид, свидетельствуя о наличии долговременных корреляций в системе. И, наконец, абсолютная величина фрактальной сложности по окончании переходного процесса оказывается в таком режиме существенно выше.
На качественном уровне неабелев режим эволюции модели позволяет воспроизвести важную для сейсмологии связь между степенной статистикой землетрясений и масштабно-инвариантной структурой разломов земной коры. Действительно, крупномасштабные пространственные корреляции энергии элементов модели являются следствием ее лавинной динамики, что имеет определенное сходство с поведением тектонических систем, структура которых модифицируется под влиянием землетрясений. В обоих случаях скорость изменения крупномасштабного ландшафта оказывается значительно ниже скорости развития отдельных неустойчивостей. На временном масштабе отдельной лавины модели или отдельного землетрясения может создаться впёчатле-ние, что не затронутая неустойчивостью геометрия распределения свободной энергии не подвержена изменениям. Однако при более длительном наблюдении становится.очевидно, что она непрерывно эволюционирует, определяя наиболее медленный временной масштаб системы, отражающий кооперативные свойства исследуемых систем. Механизмом сопряжения медленного и быстрого временных масштабов в модели является стационарная критическая динамика в окрестности точки СОК, что позволяет предположить действие такого же механизма в сейсмически активных регионах, описываемых рассматриваемым классом моделей.
Связь медленного и быстрого временных масштабов в состоянии СОК интересна тем, что содержит возможность прогноза катастрофических землетрясений по долговременным изменениям фрактальной геометрии разломов коры. В частности, можно ожидать, что средний размер сейсмических событий и вероятность катастрофических ' землетрясений уменьшатся, если масштабно-инвариантная структура разломов будет подвержена действию дополнительных возмущений, выводящих систему из окрестности точки СОК в область более спокойной субкритической динамики. В противоположном случае, при переходе системы в суперкритическое состояние, средний размер событий должен возрасти, а вероятность сейсмических катастроф увеличиться. Описанная
3
' 100
80 60
х 40
Рис. 3. Зависимость среднего размера лавин Ъ от фрактальной сложности Е в стационарном состоянии неабелевой версии моде-
закономерность подтверждается результатами фрактального анализа пространственного распределения региональной сейсмичности, обнаружившего снижение размерности кластеров эпицентров средних и малых землетрясений перед крупными землетрясениями [9]. Здесь мы кратко рассмотрим возможности численного моделирования этого эффекта. _
Пусть з (Е) характеризует средний размер лавин, порождаемых при разных уровнях возмущенности стационарного состояния СОК, соответствующих различным значениям Е.~ Для определения этой зависимости исследуемая конфигурация решетки сохранялась в памяти компьютера, после чего определялось начальное значение в как результат усреднения по всем лавинам, инициированным последовательно в каждом из элементов верхнего ряда решетки. Далее энергии определенной части случайно вы- , бранных пар элементов обменивались между собой, что приводило к частичному разрушению крупномасштабных корреляций и снижению фрактальной сложности. Процедура последовательной рандомизации решетки и оценки 5 продолжалась до тех пор, пока значение Е не достигало нижнего порога насыщения, соответствующего полностью неупорядоченному распределению энергии в модели. ,
На рис. 3 приведена зависимость 1(Е), установленная в результате усреднения по 1100 начальным конфигурациям. Кривая на рисунке показывает, что наибольший размер лавин наблюдается в исходном состоянии модели, соответствующем максимальному значению Е. С падением фрактальной сложности система переходит в субкритичес-кий режим, и размер лавин закономерно уменьшается. Параметр I может рассматриваться в качестве оценки критической восприимчивости модели к слабым возмущениям. Согласно анализу, проведенному в рамках динамической теории среднего поля [10], восприимчивость СОК-систем ограничена внутренними диссипативными процессами, которые в исследуемой модели представлены поглощением энергии на открытой нижней границе. При стремлении скорости диссипации к нулю восприимчивость, а вместе с ней и средний размер лавин стремятся к бесконечности. При конечных уровнях диссипации поведение модели с максимальными значениями «и Е является катастрофичным в том смысле, что обеспечивает максимальную вероятность неустойчивостей, соизмеримых по масштабу с размерами решетки. При понижении фрактальной сложности, связанном с разрушением крупномасштабных корреляций в пространственном распределении энергии элементов системы, эта вероятность также быстро снижается. Поведение модели указывает на то, что размер землетрясений в сейсмических систе-
мах может предсказываться на основе анализа фрактальной структуры матрицы разломов земной коры в рассматриваемом регионе. Увеличение уровня фрактальной слож- .
ности матрицы разломов, оцененного с использованием формул (3)-(6), может служить индикатором приближающегося характеристического землетрясения. Для прогноза его магнитуды следует построить эмпирическую зависимость, аналогичную показанной на рис. 3. /
Фрактальная структура модели способна генерировать временные сигналы, несущие информацию об изменении внутреннего состояния сейсмической системы, в частности о накоплении тектонических напряжений, обусловливающих уменьшение информационной энтропии. Эта возможность основана на том, что развитие неустойчивостей в неабелевой модели Хьюза-Пацзуски определяется характером распределения энергии, на различных пространственных масштабах. В стационарном состоянии модели суб-критические элементы, энергия которых близка к порогу возбуждения 2С, образуют' системы ветвящихся «каналов» (см. рис. 1, б), которые определяют пути прохождения последующих лавин и меняют свою форму под их воздействием. Если лавина начинает свой рост на длинной цепочке таких элементов, она имеёт существенно более высокую вероятность дорасти до размеров катастрофического события по сравнению с лавиной, возбужденной на короткой цепочке.
В том случае, если энергия взаимодействия элементов определяет их электропроводность, существует возможность выработать схему сейсмоэлектромагнитных измерений, позволяющую получать информацию об условиях прохождения крупномасштабных лавин, связанных с риском катастрофических землетрясений. Электропроводность рассчитывалась исходя из предположения, что устойчивые [г — 0) и субкрити-ческие (0 < 2 ^ гс) элементы соответствуют двум фазам земной коры - твердой горной породе и жидкости, наполняющей имеющиеся в ней пустоты. Для простоты считалось, что проводимость устойчивых элементов равна 0, в то время как проводимость субкри-тических элементов равна 1. Концентрация С проводящей фазы в пространственной области П как функция дискретного времени 2 может быть оценена по формуле ' ■ • '
С(«) = <0(1-2, (®,»)))Я>1>6П, . . (7)
где символ ()х у еП означает усреднение по всем узлам решетки в пределах в -функция Хевисайда. Известно, что эффективная электропроводность неупорядочен- . ной смеси проводящего материала с диэлектриком меняется в зависимости от концентрации как (С — Со)*1 [11]. Порог перколяции Со и критический индекс ц чувствительны к пространственным корреляциям, связанным с самоаффинным распределением энер- / гии в неабелевом режиме модели, и на данном этапе исследования не могут быть получены аналитически. Вместо этого была проведена прямая оценка перколяционной проводимости ар. При этом предполагалось, что используемое для измерения сгр электрическое поле создается постоянной разностью потенциалов, приложенной к области решетки $7 в направлении у. Электрический ток через эту область обусловливается существованием связного перколяционного кластера П С П субкритических элементов, соединяющего левую и правую границы П. На периодах времени, когда такой кластер существует, величина ар рассчитывалась по формуле
"»(*) = (£,€Л£,еп[0(*'!')><',(1~2‘(х’!'))]) ') ■
(8)
где
Q [х v) = I 1 ^x' ^ e
^ . 1 0 {x,y) iU. , .
Для остальных периодов полагалось стр = 0. Соотношения (8) ‘дают верхнюю оценку проводимости, при которой учитываются все топологические структуры перколяцион-ного кластера П, в том числе такие, которые не влияют на его способность к пропусканию тока. По нашим оценкам, такое приближение оправдано, когда соотношение вертикального и горизонтального размеров области Г? мало. В обсуждаемых ниже результатах это условие было соблюдено.
Фрактальная топология перколяционных кластеров проводящей фазы является одним из главных факторов, определяющих электрические свойства неупорядоченных систем [12]. Было высказано предположение [13], что изменение ее свойств перед наступлением землетрясений способно вызывать генерацию предвестниковых сейсмоэлек-тромагнитных сигналов. Для моделирования этого явления были рассмотрены временные ряды флуктуаций концентрации (7) и электропроводности (8). Было обнаружено, что значения производных по времени от этих параметров достоверно коррелируют с размером лавин и образуют область регрессии с четко выраженной верхней границей. В основе этой зависимости лежит тот факт, что чем больше лавина, тем более существенные изменения она вызывает в распределении энергии между элементами. Таким образом, прогнозируя временною эволюцию параметров электропроводности, можно оценивать и максимальный размер лавин, ожидаемых в следующие моменты времени.
Горизонты прогноза временных рядов С и av определяются их корреляционной структурой, которая; была исследована методом Пенга (detrended fluctuation analysis - метод DFA), состоящим в вычислении среднеквадратичного отклонения F исследуемого сигнала относительно его локальных трендов. . Для фрактального временного ряда оно зависит от временного масштаба т как F ~ та. В случае, когда исследуемый процесс представим в виде броуновского шума с некоррелированными приращениями, степенной индекс а равен 1.5. Другие значения а сигнализируют о положительных или отрицательных корреляциях приращений. Преимуществом метода DFA является . его устойчивость к наличию во временном ряде нестационарных помех. '
Как показал анализ временных рядов модели, динамика концентраций проводящей фазы характеризуется значением а = 1,5 (рис. 4), практически не допускающим возможность прогноза динамики модели на основе C(t). Однако динамика проводимости ap(t) характеризуется другими значениями.а, указывающими на наличие долговременных корреляций приращений этого процесса. На масштабах т < 102 индекс DFA достоверно ниже 1,5, что соответствует отрицательным корреляциям приращений проводимости и ее производной по времени. На более длинных масштабах он достоверно выше этого уровня, что указывает на положительные корреляции приращений. Таким образом, в обоих интервалах т существует возможность прогноза временной эволюции процесса daPldt и определения на этой основе вероятности характеристических лавин, соответствующих крупным землетрясениям в сейсмических системах.
Заключение. Полученные данные указывают на перспективность использования модели Хьюза-Пацзуски с нарушенной абелевой симметрией для изучения эффекта взаимодействия между статистикой Гуттенберга-Рихтера и масштабно-инвариантной геометрией матриц тектонических разломов в зоне землетрясений. Преимущество использованной модели состоит в крупномасштабной самосогласованной корреляционной структуре распределения энергии локальных взаимодействий. Предложена количественная мера фрактальной сложности Е, позволяющая оценивать пространственную
. 1 , ' 2 3 1 2 - Ъ \%Г
Рис. 4• Примеры флуктуаций концентрации проводящей фазы С(4) и электропроводности ор(1) в области х = 0-200, у = 20-40 (а); соответствующие им функции БРА (слева - для С(*), справа - для ар{{)) (б).
самоорганизацию модели в окрестности стационарной критической точки. Значение Е определяет чувствительность модели к слабым возмущениям, а также средний размер лавин, что указывает на возможность разработки алгоритмов прогноза сильных землетрясений, основанных на анализе самоаффинной структуры разломов.
Было также показано, что крупномасштабные корреляции в пространственном распределении энергии приводят к фрактальной временной эволюции характеристик элек-
тропроводности модели. Временные ряды электропроводности, зависящей от переменных перколяционных свойств кластеров субкритических элементов, демонстрируют долговременные корреляции, которые могут использоваться в целях диагностики текущего состояния системы и прогноза ее будущих состояний. На качественном уровне динамика величины dap(t)/dt сопоставима с электромагнитными сигналами диапазона УНЧ, регистрируемыми перед крупными землетрясениями [14]. Для построения их количественной модели необходимы включение более реалистичного механизма ло-. кальных взаимодействий, а также учет связи критической динамики с формированием сейсмических циклов и другими кооперативными явлениями, наблюдаемыми в сейсмических системах [15, 16]. ,
Авторы выражают признательность М. Пацзуски (Имперский колледж, Лондон) за обсуждение алгоритма модели, А. Тзанису (Афинский университет) и Ф. Валлианатосу (Технологический образовательный институт Крита) за полезные замечания при работе над рукописью.
Summary
Uritsky V. М., Smirnova N. A., Troyan V. N. Criticality evolution of fractal systems as a possible mechanism of the formation of the fault zones and the generation of seismo-electromagnetic precursors.
Regional seismicity is known to demonstrate scale-invariant properties in different ways. Some typical examples are fractal spatial distributions of hypocent^rs, Gutenberg-Richter magnitude statistics, fractal clustering of earthquake onset times, power-law decay of aftershock sequences, as well as scale-invariant geometry of fault systems. In some regions, the observed scale-free Affects are likely to be connected to a cooperative behavior of interacting tectonic plates and can be described in terms of the self-organized criticality (SOC) concept. In this work, we investigate a new SOC model incorporating short-term fractal dynamics of seismic instabilities and slowly evolving matrix of cracks (faults) reflecting long-term history of preceding events. The model is based on a non-Abelian directed sandpile algorithm proposed recently by Hughes and Paczuski (2002), and displays a self-organizing fractal network of occupied grid sites similar to the structure of stress fields in seismic active regions. Depending on the geometry of local stress distribution, some places on the model grid have higher probability of major events compared to the others. This dependence makes it possible to consider a time-dependent structure of the background earth crust geometry as a sensitive seismic risk indicator. We also propose a simple framework for modeling ultra-low frequency (ULF) electromagnetic emissions associated with abrupt changes in the large-scale geometry of the stress distribution before characteristic seismic events. ■
Литература / ^
1. Bah P., Tang C-, Wiesenfeld K. j / Phys. Rev. Letters. 1987. Vol. 59, N 4. P. 381-384. 2. Bak P., Tang C., Wiesenfeld K. // Phys. Rev. A. 1988. Vol. 38, N 1. P. 364-374. 3. Takayasu H. Fractals in the physical science (Nonlinear Science: Theory and application series). Manchester, 1990. 4. Turcotte D. L. // Physics of the Earth and Planetary Interiors. 1999. Vol. Ill, N 3-4. P. 275-293. 5. Pavlides S. B:, Zhang P., Pantosti D. jj Tectonophysics. 1999. Vol. 308: N 1-2. P. vii-x. 6. Hughes D., Paczuski M. // Phys. Rev. Letters. 2002. Vol. 88, N 5. P. 054302-1054302-4. 7. Dhar D. j j Physica A. 1999. Vol. 263, N 1-4. P. 4-25. 8. Bak P. How nature works: The science of self-organized criticality. Oxford, 1997. 9. Урицкий В. М., Троян В. Н. // Вопросы геофизики. 1998. Т. 35, № 1. С. 39-42. 10. Vespignani A., Zapperi S. // Phys. Rev. Е. 1998. Vol. 57, N 6. P. 6345-6362. 11. Fractals and disordered systems / Eds. A. Bund, S. Havin. Berlin,
1996. 12. Bahr K. jj Geophys. J. Intern. 1997. Vol. 130, N 2. P. 649-660. 13. Smirnova N. // Atmospheric and ionospheric electromagnetic phenomena associated with earthquakes / Ed. by , M. Hayakawa. Tokyo, 1999. P. 215-232. 14. Uritsky V., Smirnova N., Troyan V. et al. // Physics and Chemistry of Earth. 2004. Vol. 29. P. 473-480. 15. Tzanis A., Vallianat6s F. // Natural Hazards and Earth System Sciences. 2003- Vol. 3, N 1. P. 1-17. 16. Rundle J. B., Klein W-> Turcotte. D. L. et al. // Pure and Applied Geophysics. 2000. Vol. 157, N 5. P. 2165-2182. '
Статья поступила в редакцию 29 января 2004 г.
