Моделирование ультранизкочастотных электромагнитных эмиссий, возникающих перед и во время сильных землетрясений Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

МОДЕЛИРОВАНИЕ УЛЬТРАНИЗКОЧАСТОТНЫХ ЭЛЕКТРОМАГНИТНЫХ ЭМИССИЙ, ВОЗНИКАЮЩИХ ПЕРЕД И ВО ВРЕМЯ СИЛЬНЫХ ЗЕМЛЕТРЯСЕНИЙ

Н.А. Семенов

Научный руководитель - доктор технических наук, профессор А.Г. Коробейников

Перед и в течение сейсмически активного периода на земной поверхности наблюдается ультранизкочастотное (УНЧ) электромагнитное излучение (Р=0.001-10 Гц). Одним из источников этого излучения служит процесс образования микротрещин, возникающих в области очага будущего землетрясения при усилении тектонических движений. Моделированию этого процесса и посвящена данная работа.

Введение

Землетрясения - это явления, при которых земная кора и верхняя мантия постепенно подвергаются воздействию обширного поля сил [1]. В ослабленных участках внезапно образуются разрывы, земная кора как бы вспарывается на какой-то глубине, при этом единым импульсом высвобождаются напряжения и деформации и возбуждаются сейсмические волны. Как правило, в течение сейсмоактивного периода происходит много сейсмических толчков различной магнитуды. Собственно землетрясением обычно называется самый сильный сейсмический толчок в сейсмическом рое. В сейсмоактивных зонах ведутся многочисленные наблюдения и съемки с целью слежения за движениями коры, сейсмичностью, геомагнетизмом, грунтовыми водами, геохимическими изменениями, выделяющимися газами и т.д., однако проблема предсказания разрушительных землетрясений еще далека от своего решения. Как указывается в литературе, очень много сильных землетрясений (~50%) происходит без форшоковой активности, поэтому исследования предвестников землетрясений несейсмического характера, имеют большое значение для решения проблемы предсказания землетрясений [1].

Среди выделенных многочисленных среднесрочных и краткосрочных прогностических характеристик особое место занимают электромагнитные предвестники землетрясений, особенно в ультранизкочастотном (УНЧ) диапазоне (Р =10-0.001 Гц) [2]. Одним из источников этого УНЧ излучения служит процесс образования микротрещин размером ~0.5-10 мм, возникающих в области очага будущего землетрясения при усилении тектонических движений в земной коре [3]. В момент раскрытия микротрещин в лабораторных экспериментах наблюдалось возникновение коротких электромагнитных импульсов длительностью ~1 мкс [4, 5]. Данная работа посвящена моделированию этого процесса.

Описание математической модели

Для решения поставленной задачи была выбрана математическая модель, основанная на том, что участок земной среды, электромагнитное излучение которого моделируется, берется в форме параллелепипеда. Задаются координаты расположенной на земной поверхности точки наблюдения (той точки, в которой мы хотим наблюдать излучение), обозначаемые на рис. 1 как (Х,У,2).

Сам параллелепипед разбивается на параллелепипеды меньшего размера, каждый из которых считается в дальнейшем местом расположения источника электромагнитного излучения. Задается также некоторый период времени, в течение которого происходит тектонический процесс, вызывающий образование микротрещин размером 0.5-10 мм в породах земной коры. В момент своего образования микротрещины являются источниками, генерирующими электромагнитные импульсы очень малой длительности (~1 мкс) [1]. Этот период разбивается на две части:

первая часть (1), в которой происходит нарастание плотности излучения перед землетрясением;

вторая часть (2), в которой происходит снижение плотности излучения после землетрясения.

(Х,у,г)

Рис. 1. Геометрическая схема модели

Плотность временных отметок (моментов возникновения микротрещин) в первой части (Р1) рассчитывается по формуле

Р1 = Щв*+\

где и Ь1 - некоторые положительные коэффициенты Плотность временных отметок во второй части (р2) рассчитываются по формуле

Р2 = «2в +?2,

где «2 и ¿2 - некоторые положительные коэффициенты, t - время, а tl и ?2 - начальные отступы. Коэффициенты а1, Ъ1, а2 и Ъ2 выбираются таким образом, что плотности временных отметок в точке пересечения первой и второй частей совпадают. При этом спад плотности излучения происходит значительно быстрее, чем возрастание [5].

Сам процесс излучения представляется таким образом, что в каждый момент времени из тех, которые были рассчитаны заранее, одна из точек в изначальном параллелепипеде может стать источником направленного излучения (излучения, максимум амплитуды которого приходится на направление, перпендикулярное большей оси микротрещины). Значение амплитуды излучаемого сигнала в точке излучения представляется формулой

Б(0 = ав

= ав~к ^ ^о)

(1)

где а и к - некоторые положительные коэффициенты, а 10 - момент начала излучения.

При этом коэффициенты берутся таким образом, что через время, равное некоторому заранее заданному At, значение амплитуды сигнала можно принять равным нулю.

Поскольку амплитуда сигнала с расстоянием меняется в зависимости от частоты и проводимости земной коры, необходимо разложить сигнал из формулы (1) на частотные составляющие. Для разложения используется быстрое преобразование Фурье [6].

Рассмотрим значение амплитуды магнитного поля в точке, отдаленной от источника излучения на некоторое расстояние г в момент времени t:

Б^) = с(г)Б0 (^)в~вг ^(М - кг) * ^(а), где w - частота сигнала, к - волновой вектор, Б0 - амплитуда сигнала, с(г) - некоторый положительный коэффициент, зависящий от расстояния, в - некоторый положи-

тельный коэффициент, а - угол между направлением распространения сигнала и направлением на точку наблюдения. Применяя приближение плоской волны, окончательно получим:

B(t) = B0 (w)e~kr cos(wt - kr) * cos(a). (2)

Для плоской волны волновой вектор k рассчитывается по формуле

, 2п

k =—, к

где к = 107 Pk^ - длина волны, Т = - период сигнала, Pk - удельное электриче-

л1 р^Т - длина волны, Т = 2п

w

ское сопротивление среды [2, 7].

Амплитуда сигнала B0 (f) зависит от частоты и частотной плотности в соответствии со следующей формулой:

ВоИ = ,

где р - частотная плотность, а f - частота, причем w = 2я[ .

Следует также заметить, что, поскольку скорость распространения волны конечна, то сигнал достигнет точки наблюдения не мгновенно, а через некоторое время г

*р = V, (3)

где г - расстояние до точки наблюдения, а V - фазовая скорость. При этом

107Pk (4)

V = 4

Т V Т

Подставляя значение V из формулы (4) в формулу (3), окончательно получаем:

= г4Т

р л/Ю^РТ'

Каждый источник излучения (рис. 1) имеет свои собственные координаты (х ,, у,, ). Соответственно, координаты точки наблюдения относительно источника с

этими координатами будут (X - х,, У - у {, 2 -1 .,), а расстояние до источника будет вычисляться по формуле

Г =д/(X - х )2 + (У - у )2 + (2 -1, )2 . (5)

Поскольку каждый из источников имеет свое собственное направление распространения, необходимо отдельно рассчитывать горизонтальные и вертикальную составляющую магнитного поля. Вертикальная составляющая В2 характеризует составляющую магнитного поля, перпендикулярную земной поверхности и направленную вниз, а Ва - составляющую магнитного поля, лежащую в плоскости земной поверхности. Горизонтальная составляющая Ва, в свою очередь, раскладывается на компоненты по широте и по долготе Вн и Вв, как показано на рис. 2. Здесь BV - проекция вектора магнитного поля на плоскость, перпендикулярную земной поверхности, О - проекция направления распространения волны на земную поверхность, Н - ось координат, соответствующая широте, а1 - угол между В и В^,, а2 - угол между BV и В2, Д- угол

между Ва и О, в2 - угол между О и Н.

Рис. 2. Схема проекции вектора магнитного поля на оси координат

Как следует из рис. 2: В2 = В сов(а ) сов(а2),

Ва = В^ 1 - соб2 (а) СОБ2 (а2).

где В - величина магнитного поля в точке наблюдения. Соответственно, Вн = Ва СОБ(Д + Д), Вв = Ва мп( рх + Д2), СОБ(р1 + Р2) = СОБ(р1) СОБ(р2) - бЗД ) Мпф2) =

(6)

(7)

(8) (9)

А

1 -

б1п (а1)

2 2 1 - соб (а1)соБ (а2)

соб(Р2) -

Б1п(а1)

1 - соБ(а1)соБ(а 2)

81п(р2):

в1П(Р1 + Р2) = ) СОБ(р2 ) + СОБ(р1) БШ^ ) =

Б1п(а1)

■СО8(р2)+

1-

б1п (а1)

22 1 - соб (а1)соБ (а2)

81п(р2)

1 - соБ(а1)соБ(а 2)

Для получения графика компонент магнитного поля в точке наблюдения весь период времени наблюдения разбивается на равные части с шагом, величина которого такова, что значение компонент магнитного поля для каждого из отрезков можно считать постоянным. Тогда

Вг = X В2 , Вэ] = X Во , ВН = X Вн ,

¿ЕЬ. ¿ЕЬ. ¿ЕЬ.

где Ь. - отрезок времени с индексом .. При этом компоненты магнитного поля берутся из формул (6)-(9), а величина вектора магнитного поля - из формулы

В, ) = В0(ч)е-* соб(^ - ) - кг) * с08(аг). (10)

Реализация математической модели

Реализация математической модели строилась исходя из минимальных системных требований и требований высокой производительности. Соответственно, исключались все компоненты, использование которых могло привести к снижению скорости расчетов. Также проводился предварительный расчет некоторых величин, постоянных для всех частот. Это несколько усложняло структуру основного цикла, но давало ощутимые преимущества в скорости.

Структуру основного расчета можно представить следующим образом.

1. Разбиение исходного экспоненциального излучения на некоторое заранее заданное число частот. Поскольку для получения спектра использовалось быстрое преобразование Фурье, то число это должно было равняться степени двойки.

2. Расчет некоторого числа точек во времени, в которое происходило излучение. Расчет этот делался по схеме, описанной в предыдущей части. При этом учитывалось, что излучение происходит в случайные моменты времени, и, соответственно, были внесены случайные сдвиги от экспоненциальной схемы.

3. Период времени, за который велся расчет, разбивался на отрезки фиксированной длины. Число этих отрезков было достаточно велико (обычно бралось 103-105), но при этом на несколько порядков меньше, чем число источников излучения.

4. Основной цикл стартовал по заранее рассчитанным точкам времени, при этом имелся внутренний цикл по частотам. Происходило постепенное заполнение выходного потока. Следует отметить, что такая схема полностью аналогична схеме, непосредственно вытекающей из формул

Вг =£ В2 , Во] =£ Вп', ВН =£ ВН ,

где - отрезок времени с индексом / Но такая схема требует индексирования в

массиве по данным (косвенного индексирования). Выбранная же схема позволяет использовать прямое индексирование, хотя она и менее очевидна.

5. При вычислении компонент магнитного поля следовало учитывать, что для каждого источника излучения имеются свои параметры:

/2 2 2

• г, =у (X - х,) + (У - у,) + (2 -) - расстояние до точки наблюдения;

• - время, за которое сигнал достигает точки наблюдения;

• а - угол между направлением на точку наблюдения и основным направлением

распространения сигнала. Перед расчетом задаются удельное электрическое сопротивление среды, начальная амплитуда магнитного поля в микрофарадах, а также координаты точки, для которой производится расчет.

Результаты

На верхней панели рис. 3 показан результат проведенного моделирования (представлена вертикальная компонента магнитного поля). Момент землетрясения совпадает с максимумом амплитуды вариаций магнитного поля. Точка наблюдения располагалась на расстоянии 50 км от эпицентра

Очаг землетрясения (параллелепипед из раздела 2) задавался на глубине 15 км. На нижней панели рис. 3 представлен результат узкополосной фильтрации полученной реализации модельных данных (Р=0.004-0.005 Гц). Из рисунка видно, что, несмотря на то, что модельное поле возникло как результат действия большого количества коротких электромагнитных импульсов длительностью 1 мкс, полученная реализация модельных данных имеет весьма широкий спектр частот.

На рис. 4 построены зависимости фазовой скорости, градиента и амплитуды полученной реализации модельных данных от периода вариаций. Координатные оси представлены в логарифмическом масштабе. Подобные зависимости этих величин соответствуют экспериментальным данным [8].

3.0 •

2.0 •

1.0 ■

0.0

-1.0 ■

-2.0 •

-3.0 0.2

0.1

0.0

-0.1

-0.2

ч I I I I Г

1 г-1 г

1 I г

I I I

' I ' I ' I ' I ' I ' I 1 I ' I ' I ' I 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000

Время (с)

Z

Рис. 3. Реализация нефильтрованных (сверху) и фильтрованных (снизу) модельных

данных землетрясения

Рис. 4. Зависимости фазовой скорости, градиента и амплитуды реализации модельных

данных от периода вариаций

Заключение

В работе показано, что большое количество излучателей, генерирующих короткие электромагнитные импульсы длительностью ~1 мкс, могут создать на земной поверхности магнитное поле, имеющее весьма широкий спектр; в том числе возникают ультранизкочастотные магнитные вариации, наблюдаемые в большом числе экспериментов. По-видимому, теория образования микротрещин, возникающих в земной коре перед и в течение сейсмоактивного периода, вполне адекватно описывает процесс генерации УНЧ электромагнитных возмущений в земной коре.

Литература

1. Mogi K. Earthquake predictions. Academic Press Japan, 1985, 166 p.

2. Ismaguilov V.S., Kopytenko Yu.A., Hattori K., Voronov P.M., Molchanov O.A., Haya-kawa M. ULF Magnetic Emissions Connected with Under Sea Bottom Earthquakes // Natural Hazards and Earth Sys. Sci. 2001. V.1. Р. 1-9.

3. Scholz C.H., Molnar P., Johnson T. Detailed studies of frictional sliding of granite and implications for earthquake mechanism. // J. Geophys. Res. 1972. V.77. Р. 6392-6406, Cress G.O., Brady B.T., Rowell G.A. Sources of electromagnetic radiation from fracture of rock samples in laboratory. // Geophys. Res. Lett. 1987. V.14. Р. 331-334,

4. Molchanov O.A., Hayakawa M. Generation of ULF electromagnetic emissions by mi-crofracturing. // Geoph. Res. Lett. 1995. V. 22. P. 3091-3094.

5. Ковтун А.А. Использование естественного электромагнитного поля при изучении электропроводности Земли. Л.: Изд. ЛГУ, 1980. 195 с.

6. Семенов А.А. Теория электромагнитных волн. М.: Изд. МГУ, 1968. 316 с.

7. Ismaguilov V.S., Kopytenko Yu.A., Hattori K., Hayakawa M. Variations of phase velocity and gradient values of ULF geomagnetic disturbances connected with the Izu strong earthquakes. // Natural Hazards and Earth Sys. Sci. 2002. V.20. P. 1-5.