СПОСОБ КРАТКОСРОЧНОГО ПРОГНОЗИРОВАНИЯ ЗЕМЛЕТРЯСЕНИЙ
METHOD FOR SHORT-TERM EARTHQUAKE PREDICTION

Изобретение относится к области сейсмологии и может найти применение для прогнозирования землетрясений в национальных системах наблюдения и обработки данных геофизических измерений.
Для прогнозирования землетрясений используется диагностика предвестниковых аномальных возмущений состояния геофизических полей. Перед землетрясением наблюдается увеличение выхода различных газов: водорода, гелия, метана, радона. Факт выхода газов из зоны сейсмотектонической аномалии (СТА) в атмосферу перед землетрясением подтверждается резким [Патенты RU №2204852, 2003 г, №2275659, 2006 г, №2302020, 2007 г] увеличением фонового уровня радиации, ионизацией молекул воздуха, вплоть до свечения атмосферы [Тертышников А.В. Сейсмоозонные эффекты и проблема прогнозирования землетрясения, СПБ, 1999. 196 с.].
Известен «Способ краткосрочного прогнозирования землетрясений» -Патент RU №2423729, 2011 г. В способе-аналоге размещают измеритель на космическом носителе с регулируемым интервалом длительности и скважности дискретных отсчетов вдоль трассы полета носителя, измеряют свечение атмосферных газов в ультрафиолетовой полосе в диапазоне 100…300 нм, численным интегрированием функций получаемых регистрограмм рассчитывают их длину L, координаты гипоцентра очага отождествляют с максимумом регистрограмм, рассчитывают постоянную времени сейсмического процесса T из соотношения:
<mrow><mi>T</mi><mo>=</mo><mfrac><mrow><mi>Δ</mi><mi>t</mi></mrow><mrow><mi>ln</mi><mfrac><mrow><mrow><mo>(</mo><mrow><msub><mi>L</mi><mn>0</mn></msub><mo>−</mo><msub><mi>L</mi><mn>1</mn></msub></mrow><mo>)</mo></mrow></mrow><mrow><mrow><mo>(</mo><mrow><msub><mi>L</mi><mn>0</mn></msub><mo>−</mo><msub><mi>L</mi><mn>2</mn></msub></mrow><mo>)</mo></mrow></mrow></mfrac></mrow></mfrac></mrow>
,
y
прогнозируют время удара t[сут]=0,77M-4,4, где L0 - предельная длина дуги функции регистрограммы, равная:
<mrow><msub><mi>L</mi><mn>0</mn></msub><mo>=</mo><mfrac><mrow><msubsup><mi>L</mi><mn>2</mn><mn>2</mn></msubsup><mo>−</mo><msub><mi>L</mi><mn>1</mn></msub><msub><mi>L</mi><mn>3</mn></msub></mrow><mrow><mn>2</mn><msub><mi>L</mi><mn>2</mn></msub><mo>−</mo><msub><mi>L</mi><mn>1</mn></msub><mo>−</mo><msub><mi>L</mi><mn>3</mn></msub></mrow></mfrac></mrow>
,
где L1, L2, L3 - длины дуг функций регистрограмм в моменты измерений t1, t2, t3; Δt=(t2-t1)=(t3-t2) - межвитковый интервал времени космического носителя, на котором проводят измерения.
Недостатками аналога следует считать:
- неопределенность координат СТА, что усложняет планирование работы бортового измерителя над заданным районом;
- при скорости полета космического носителя ~10 км/с и размерах СТА ~300 км, максимальное время измерений не превышает 30 секунд, при возможности повторных измерений через 14 витков, что обуславливает значительный риск пропуска СТА.
Ионизация молекул воздуха над СТА приводит к образованию нескомпенсированного электрического заряда в виде «купола» над СТА. [Тертышников А.В Основы мониторинга чрезвычайных ситуаций. Учебное пособие - Химки: АГЗ МЧС РФ, 2010. 266 с]. Электростатическое поле влияет на ориентацию дипольных молекул водяного пара в воздухе. Установлено «Явление возникновения поляризационных аномалий электромагнитного поля над очагом землетрясения», см. Научное открытие №336, 2007 г. Изменение поляризации электромагнитного поля является одним из признаков-предвестников готовящегося землетрясения, однако зарегистрировать этот эффект сложно.
Ближайшим аналогом к заявляемому техническому решению является «Поляризационный датчик предвестника землетрясений», Патент RU №2343507, G.01.V, 9/00, 2009 г.
Поляризационный датчик предвестника землетрясений, установленный на космическом носителе, содержит канал приема поляризованного светового потока, отраженного от подстилающей поверхности, в которой последовательно включены объектив, деполяризатор, усилитель, пиковый детектор, пороговое устройство, аналогово-цифровой преобразователь, блок буферной памяти, а также отдельный генератор пилообразного напряжения, подключенный к металлизированным обкладкам деполяризатора.
К недостаткам ближайшего аналога следует отнести:
- малое время наблюдений, ~30 сек., при большой скважности измерений из-за смещения орбиты космического аппарата от зоны СТА;
- локальность измерений по пространству и времени, что снижает статистическую устойчивость результата;
- неопределенность времени землетрясения относительно выявленных предвестников.
Задача, решаемая изобретением, состоит в регистрации аномальных вариаций спектральной мощности регистрируемого солнечного потока на ряде длин волн, прошедшего всю толщу атмосферы.
тек
этал
Технический результат достигается тем, что способ краткосрочного прогнозирования землетрясений, включает создание в сейсмоопасном регионе системы измерений оптической плотности атмосферы из групп фотометров, разнесенных с постоянным шагом по координатам x, y на расстояние не более размеров зоны подготавливаемого землетрясения, построение гистограммы коэффициента пропускания атмосферы в дискретных интервалах длин волн фотометров: 340, 380, 440, 500, 675, 800 и 1020 нм, определение разницы (Δ) между средневзвешенной длиной волны текущего (λ) и эталонного (по Планку) солнечного спектра (λ), отождествление изменений Δ(t) с началом сейсмического процесса и расчет постоянной времени T функции сигнала Δ(t), формирование регистрограмм измерений для каждой из групп раздельно по координатам Δ(x, t) и Δ(y, t), расчет гипотетического центра очага как точки пересечения векторов, направляющие косинусов которых вычисляют через их проекции на осях крестообразных групп:
<mrow><mi>cos</mi><mi>α</mi><mo>=</mo><mfrac><mrow><mi>Δ</mi><mrow><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mrow><msqrt><mrow><mi>Δ</mi><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>+</mo><mi>Δ</mi><msup><mrow><mrow><mo>(</mo><mrow><mi>y</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msqrt></mrow></mfrac></mrow>
<mrow><mi>cos</mi><mi>β</mi><mo>=</mo><mfrac><mrow><mi>Δ</mi><mrow><mo>(</mo><mrow><mi>y</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mrow><msqrt><mrow><mi>Δ</mi><msup><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>+</mo><mi>Δ</mi><msup><mrow><mrow><mo>(</mo><mrow><mi>y</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msqrt></mrow></mfrac></mrow>
y
y
прогнозирование времени сейсмического удара t≈4,7T и магнитуды удара из соотношения lg t≈0,77M-4,4.
Изобретение поясняется чертежами, где
фиг.1 - одна из реализаций эманации радона в атмосферу накануне землетрясения;
фиг.2 - динамика электростатической напряженности над «куполом» подготавливаемого землетрясения;
фиг.3 - изменение оптической плотности атмосферы, измеряемое фотометрами системы «AERONET»;
фиг.4 - нормированные значения эталонного (по Планку) солнечного спектра и гистограмма его текущих значений в спектральных полосах фотометра;
фиг.5 - динамика изменения средневзвешенной длины волны солнечного спектра в сейсмическом процессе;
фиг.6 - периодограмма изменений оптической плотности атмосферы по серии состоявшихся землетрясений;
фиг.7 - функциональная схема устройства.
2
э
э
Техническая сущность изобретения состоит в следующем. Накануне сейсмического удара наблюдается активная эманация радона из земной коры в атмосферу. Одна из реализаций этого процесса иллюстрируется фиг.1. В результате радиоактивного распада радона и ионизации молекул аэрозолей в атмосфере возникает электростатическое поле, иллюстрируемое фиг.2. Молекула воды (HO) обладает постоянным некомпенсированным дипольным моментом. Электрический дипольный момент единицы объема воздуха с молекулами воды: P=а·N·E, где а - коэффициент электрической поляризуемости молекул; N - число частиц в единице объема, Е[в/м] - величина электрической напряженности возникающего поля. Между поляризуемостью молекул и диэлектрической проницаемостью (ε) вещества существует зависимость:
<mrow><mrow><mo>(</mo><mrow><mi>ε</mi><mo>−</mo><mn>1</mn></mrow><mo>)</mo></mrow><mo>⋅</mo><msub><mi>ε</mi><mn>0</mn></msub><mo>≅</mo><mi>N</mi><mo stretchy="false">(</mo><msub><mi>a</mi><mi>э</mi></msub><mo>+</mo><msub><mi>a</mi><mi>n</mi></msub><mo>+</mo><mfrac><mrow><msubsup><mi>M</mi><mrow><mo>∂</mo><mi>p</mi></mrow><mn>2</mn></msubsup></mrow><mrow><mn>3</mn><mi>K</mi><msup><mi>T</mi><mi>O</mi></msup></mrow></mfrac></mrow>
n
др
1
2
отр
2
1
2
1
O
<mrow><mi>n</mi><mo>=</mo><msqrt><mrow><mi>ε</mi><mi>μ</mi></mrow></msqrt></mrow>
<mrow><mi>n</mi><mo>=</mo><msqrt><mi>ε</mi></msqrt></mrow>
где a - коэффициент ионной поляризуемости, М - дипольный момент молекул, KT - тепловая энергия молекул. В оптическом диапазоне, на границе раздела двух сред, с коэффициентом преломления n и n, происходит отражение светового потока: K=(n-n)/(n-n). Коэффициент отражения является функцией коэффициента преломления, зависящего от первичных параметров среды:
, для воздуха µ=1, поэтому коэффициент преломления практически равен
. Согласно классическим представлениям [см., например, Преломление света, Физический энциклопедический словарь под редакцией A.M.Прохорова, изд. Сов. Энциклопедия, М, 1983 г., с.168] электроны и атомы вещества под действием световой волны совершают вынужденные колебания. Наличие собственной частоты атомов и молекул приводит к зависимостям коэффициента преломления n от концентрации вторичных излучателей в воздухе:
<mrow><mi>n</mi><mo>≈</mo><mn>1</mn><mo>+</mo><mi>Σ</mi><mfrac><mrow><msub><mi>N</mi><mi>i</mi></msub><mo>⋅</mo><mi>e</mi></mrow><mrow><mi>π</mi><mo>⋅</mo><msub><mi>m</mi><mi>i</mi></msub><msup><mi>c</mi><mn>2</mn></msup></mrow></mfrac><mo>∗</mo><mfrac><mn>1</mn><mrow><mn>1</mn><mo>/</mo><msubsup><mi>λ</mi><mrow><mi>o</mi><mi>i</mi></mrow><mn>2</mn></msubsup><mo>−</mo><mn>1</mn><mo>/</mo><msup><mi>λ</mi><mn>2</mn></msup></mrow></mfrac></mrow>
i
i
oi
где N - концентрация в воздухе вибраторов i-го сорта; ε - заряд электрона; m - масса вибратора i-го сорта; λ - собственная длина волны излучения вибратора i-го сорта; λ - текущая длина волны падающего светового потока. Итак, в области электростатической напряженности «купола» подготавливаемого землетрясения изменяются коэффициенты преломления и отражения падающего светового потока. В результате изменяется оптическая плотность атмосферы. По определению [см., например, Преломление света, Физический энциклопедический словарь под редакцией A.M. Прохорова, изд. Сов. Энциклопедия, М, 1983 г., с.944] под оптической плотностью атмосферы понимается десятичный логарифм отношения падающего потока излучения к потоку прошедшему через всю толщу атмосферы:
0
Д=lg П/П.
Контроль оптической плотности атмосферы, обусловленный изменениями радиационного баланса сейсмотектонической аномалии приземного электрического поля является одним из направлений мониторинга землетрясений. Для этого создана международная глобальная сеть наблюдений за характеристиками прозрачности атмосферы «AERONET», включающая порядка 500 станций [см. Интернет http://aeronet.gsfc.nasa.gov/new web/system _descriptions.html]. В работе сети используют фотометры с измерениями оптической плотности в спектральных участках с длиной волны: 340, 380, 440, 500, 675, 870, 1020 нм. Информация о наиболее значимых землетрясениях в различных регионах Земного шара и значениях оптической плотности аэрозоля на ближайших станциях существуют на сайте Геологического общества США [.] По результатам анализа данных сайта, за несколько суток до землетрясения наблюдается аномальное увеличение оптической плотности атмосферы. Одна из реализаций динамики изменения оптической плотности накануне сейсмического удара иллюстрируется графиком фиг.3. Однако, при экранировании солнечного потока облачностью, отмечается большое число пропусков измерений, т.е. параметр оптической плотности существенно зависит от климатических условий. В заявленном способе в качестве признака-предвестника землетрясения выбрана динамика изменения разности средневзвешенной длины волны текущего солнечного спектра и эталонного (по Планку) солнечного спектра:
тек
этал
Δ(t)=λ-λ
Для чего, предварительно, оптическую плотность атмосферы в каждом из дискретных каналов пересчитывают в коэффициент пропускания, из условия обратной пропорциональности. На рисунке фиг.4 иллюстрируются графики (1) эталонного (по Планку) нормированного солнечного спектра и гистограмма текущего (2) зарегистрированного солнечного спектра (одна из реализаций). Средневзвешенное значение длины волны делит площадь под графиками (1, 2) фиг.4 пополам и определяется из соотношения:
<mrow><mstyle displaystyle="true"><mrow><munderover><mo>∫</mo><mrow><mi>λ</mi><mo>=</mo><mn>340</mn></mrow><mrow><msub><mi>λ</mi><mrow><mi>с</mi><mi>р</mi></mrow></msub></mrow></munderover><mrow><mi>П</mi><mrow><mo>(</mo><mi>λ</mi><mo>)</mo></mrow><mi>d</mi><mi>λ</mi><mo>=</mo><mstyle displaystyle="true"><mrow><munderover><mo>∫</mo><mrow><msub><mi>λ</mi><mrow><mi>с</mi><mi>р</mi></mrow></msub></mrow><mrow><mn>1020</mn></mrow></munderover><mrow><mi>П</mi><mrow><mo>(</mo><mi>λ</mi><mo>)</mo></mrow><mi>d</mi><mi>λ</mi></mrow></mrow></mstyle></mrow></mrow></mstyle></mrow>
этал
тек
Из графиков фиг.4: λ=560 нм, λ=640 нм, Δ=80 нм.
y
y
Как отмечалось выше, коэффициент преломления атмосферы зависит от степени ионизации, а последняя - от объема эманации радона и размеров зоны подготавливаемого землетрясения. Дисперсия коэффициента преломления приводит к возрастанию средневзвешенной длины волны текущего спектра солнечного потока, что адекватно отражает сейсмический процесс. Динамика изменения параметра (Δ) содержит скрытую информацию о характеристиках предстоящего сейсмического удара: времени удара (t), магнитуде удара (M) и гипотетическом центре очага. Известно уравнение Гутенберга-Рихтера, связывающее время существования сигнала предвестника с магнитудой сейсмического удара: lg t=0,77M-4,4 [см. аналог]. Из математики известно [см., например, Н.С. Пискунов, Дифференциальное и интегральное исчисления для ВТУЗов, учебник том.1, изд. Пятое, Наука, М, 1964 г., стр.457-458], что сама величина и скорость ее изменения связаны дифференциальным уравнением первого порядка, общим решением которого является экспонента. Экспонента обладает тем свойством, что по трем ее дискретным отсчетам может быть восстановлена вся функция. Решение дифференциального уравнения для функции [Δ(t)] иллюстрируется графиком фиг.5. Постоянную времени сейсмического процесса определяют из соотношения:
<mrow><mi>T</mi><mo>=</mo><mfrac><mrow><mi>Δ</mi><mi>t</mi></mrow><mrow><mi>ln</mi><mfrac><mrow><mrow><mo>(</mo><mrow><msub><mi>Δ</mi><mn>0</mn></msub><mo>−</mo><msub><mi>Δ</mi><mn>1</mn></msub></mrow><mo>)</mo></mrow></mrow><mrow><mrow><mo>(</mo><mrow><msub><mi>Δ</mi><mn>0</mn></msub><mo>−</mo><msub><mi>Δ</mi><mn>2</mn></msub></mrow><mo>)</mo></mrow></mrow></mfrac></mrow></mfrac></mrow>
0
Установившееся значение Δ вычисляют как:
<mrow><msub><mi>Δ</mi><mn>0</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>Δ</mi><mn>1</mn></msub><msub><mi>Δ</mi><mn>3</mn></msub><mo>−</mo><msubsup><mi>Δ</mi><mn>2</mn><mn>2</mn></msubsup></mrow><mrow><mn>2</mn><msub><mi>Δ</mi><mn>2</mn></msub><mo>−</mo><msub><mi>Δ</mi><mn>1</mn></msub><mo>−</mo><msub><mi>Δ</mi><mn>3</mn></msub></mrow></mfrac></mrow>
1
2
3
1
2
3
где Δt - интервал времени между отсчетами функции сигнала Δ(t), Δ, Δ, Δ - значения функции сигнала в моменты отсчетов t, t, t. Экспонента достигает установившегося значения с вероятностью 0,999 при значении t=4,7T, за которое принимают время ожидаемого сейсмического удара.
Проводился статистический анализ времени существования сигнала предвестника. Данные о 30 состоявшихся землетрясениях представлялись в виде временных рядов отсчетов оптической плотности атмосферы за 30 суток до события и 30 суток после события. Для анализа использовалось Фурье-преобразование, которое позволило выявить регулярные составляющие временных рядов. Результат спектрального Фурье-анализа иллюстрируется графиком фиг.6. Из графиков фиг.6 следует, что на уровне значимости 0,05, время существования сигнала выявленного предвестника занимает интервал 1…2 суток.
Гипотетический центр очага землетрясения определяют путем пеленгации сигнала с двух, разнесенных по пространству, точек. Для чего формируют две крестообразные группы измерителей из фотометров системы «AERONET», измеряющие проекции сигнала на крестообразные оси координат. Известно, что положение радиус-вектора в пространстве определяется косинус направляющими [см. Г. Корн, «Справочник по математике для научных работников и инженеров», перев с англ, М, Наука, 1971 г., раздел. «Аналитическая геометрия», стр.73-74]. В прямоугольной системе координат косинус направляющая вектора равна отношению его проекции на данную ось к длине вектора. Длина вектора находится как корень квадратный из суммы квадратов его проекций. Точку пересечения радиус-векторов крестообразных групп отождествляют с гипоцентром очага землетрясения.
Пример реализации способа.
Заявленный способ может быть реализован по схеме фиг.7. Функциональная схема устройства фиг.7 содержит две крестообразные линейные группы (1, 2) измерителей, разнесенные в пространстве на пеленгационной базе (3). Группы образуют из единичных элементов-фотометров (4), осуществляющих измерение оптической плотности атмосферы, включенных в глобальную систему «AERONET» (5). Каждый из единичных элементов (4) подключен к канальному коммутатору (6), осуществляющему их циклический опрос, за интервал времени, задаваемый программируемой схемой выборки измерений (7). Одновременно, схема (7) осуществляет синхронизацию работы буферного запоминающего устройства (8) и устройства ввода данных (9) в компьютер (10), в составе элементов: процессора (11), винчестера (12), оперативного запоминающего устройства (13), дисплея (14), принтера (15), клавиатуры (16). Результаты обработки массива данных выводят на сервер (17) сети Интернет.
1
1
2
2
<mrow><mfrac><mrow><mi>Δ</mi><mrow><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mrow><mi>Δ</mi><mi>x</mi></mrow></mfrac></mrow>
<mrow><mfrac><mrow><mi>Δ</mi><mrow><mo>(</mo><mrow><mi>y</mi><mo>,</mo><mi>t</mi></mrow><mo>)</mo></mrow></mrow><mrow><mi>Δ</mi><mi>y</mi></mrow></mfrac></mrow>
Устройство работает следующим образом. Программы опроса линейных крестообразных групп (1, 2) формируют на ПЭВМ (10) и записывают на винчестер (12). Текущую программу опроса пересылают в программируемую схему выборки измерений (7), которая в режиме дежурного слежения циклически воспроизводится с установленной скважностью. Канальный коммутатор (6), в соответствии с текущей программой, осуществляет циклический опрос единичных элементов-фотометров (4), путем их подключения к буферному запоминающему устройству (8) на установленный дискретный интервал Δ(t). Каждый отсчет записывается в буферное ЗУ со своим адресом. После опроса всех датчиков одной группы формируется регистрограмма измерений по одной из координат Δ(x, t), которая через устройство ввода (9) пересылается в ОЗУ (13) для последующей обработки. Аналогично получают регистрограммы измерений других линейных групп: Δ(y, t), Δ(x, t), Δ(у, t). Скорость изменения функции сигнала Δ(x, y, t) в пространстве и времени содержит информацию о гипоцентре очага и моменте удара. Направление на центр купола очага определяют путем расчета градиента функции регистрограмм. Поскольку единичные измерители равномерно разнесены по координатам х, y с постоянным шагом, то производная по направлению представляется отношением конечных разностей
,
Гипоцентр определяют как точку пересечения векторов (градиентов) двух крестообразных групп. Поскольку крестообразные группы ориентированы по сторонам света, то направление градиента совпадает с азимутом. При известных координатах крестообразных групп гипоцентр очага определяют по известным азимутам на карте. Приведенные выше соотношения справедливы для любого момента времени наблюдений. Поэтому гипоцентр рассчитывают на начальном этапе наблюдений, что обеспечивает резерв времени для оповещения населения. По динамике изменения функции сигнала во времени прогнозируют характеристики ожидаемого сейсмического удара. Для чего, по дискретным отсчетам, рассчитывают постоянную времени сейсмического процесса T. Проводилась апостериорная обработка зарегистрированных данных фотометров по землетрясению в Турции 23.10.2011 г. с координатами 38°,43′, 19,20′′ с.ш. и 43°,30′, 46,78″ в.д. Результаты обработки и расчета функции сигнала представлены следующим рядом:
1
2
3
Δ(t), нм Δ=32 Δ=60 Δ=71
1
2
3
t, час t=2 t=4 t=6
0
Установившееся значение функции сигнала Δ≈100 нм.
Постоянная времени Т=4 час.
y
Ожидаемое время удара: t=4,7T=18,7 час.
Расчетная магнитуда удара М=7,3 балла.
Расчетное значение совпало с истинным с точностью до второго знака.
Устройство реализуют на существующей технической базе. В качестве единичного измерителя используют фотометры системы «AERONET» [System Description-Aerosol Robotic Network (AERONET) Homepage] Эффективность способа определяется такими показаниями, как адекватность, достоверность, устойчивость. Чтобы не пропустить ожидаемого события, устройство должно работать круглосуточно в дежурном режиме. Устройство, реализующее способ, имеет возможность адаптации к измеряемому процессу путем изменения программ, закладываемых в программируемую схему выборки измерений.