Способ предсказания землетрясений

 

Использование: геофизика, в национальных системах сейсмического контроля для прогнозирования землетрясений. Сущность изобретения: синтезированное изображение подстилающей поверхности получают попиксельным расчетом отношения квантованных уровней цифровых изображений. Осуществляют методами пространственного дифференцирования и фрактального анализа выделение контура очага на синтезированном изображении. Для участка изображения внутри контура проводят линеаментный анализ. По серии последовательных во времени синтезированных изображений очага отслеживают изменение средневзвешенной суммы азимутов и рассчитывают период этой функции. По расчетному периоду прогнозируют магнитуду и время ожидаемого удара наблюдаемого очага. Технический результат: повышение чувствительности. 5 ил.

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

Известны два класса признаков предстоящего землетрясения: статические и динамические. Статические признаки появляются за несколько лет до землетрясения, характеризуются большими диапазонами изменения наблюдаемых параметров, но не позволяют точно предсказать время удара и его магнитуду. Наиболее полный перечень статических признаков землетрясения составлен Т. Рикитаке Предсказание землетрясений, перевод с англ., М.: Мир, 1979 г., стр.314-333, табл.15.13. Данные о предвестниках.

Другой класс - динамические признаки-предвестники. Они появляются за несколько суток (часов) до удара, но в силу своей "инерционности", скрытности не могут быть зарегистрированы существующими техническими средствами. Среди динамических признаков-предвестников известны такие, как: - раскачка очага землетрясения, сопровождаемая распространением от него сверхнизких литосферных волн с периодами 3,5... 4,5 час (см., например, Давыдов В. Ф., Землетрясения. Телеметрия предвестников, М.: изд-во МГУЛ, 2001 г., стр.19); - появление над эпицетральной областью в атмосфере сильного вертикального электростатического поля величиной до нескольких кВ/м; - изменение над проекцией очага на ионосферу электронной концентрации (до 20%) в слоях ионосферы (см., например, Краткосрочный прогноз катастрофических землетрясений с помощью радиофизических наземно-космических методов, Доклады конференции, РАН, ОИФЗ им. О.Ю. Шмидта, М., 1998 г., стр.27-31).

Очаг землетрясения аккумулирует огромную энергию тектонических напряжений. В потенциальном поле механических напряжений очага наблюдаются аномалии других физических полей, захватывающих литосферу, ионосферу и магнитосферу Земли. В частности, при сжатии земной коры наблюдается анизотропия пород, изменение проводимости грунтов, их излучательной способности и, как следствие, изменение мощности сигнала собственного, восходящего излучения подстилающей поверхности над очагом (см., например, Давыдов В.Ф., Землетрясения. Телеметрия предвестников, М.: изд-во МГУЛ, 2001 г., стр.55).

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

Известен "Способ контроля лесопожарной опасности", патент РФ 2147253, кл. А 62 С 3/02, 2000 г. - аналог. В способе-аналоге осуществляют регистрацию собственного излучения подстилающей поверхности, преобразуют зарегистрированную функцию электрического сигнала в цифровые матрицы отсчетов зависимости амплитуды от координат, калибруют тракт зондирования по измерениям эталонных участков, выделяют методами пространственного дифференцирования контуры на двумерных изображениях лесных массивов, а величину влажности лесных горючих материалов внутри контуров рассчитывают по регрессивной зависимости влажности от коэффициента вариации сигнала (отношения мощности процесса и мощности переменной составляющей).

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

Ближайшим аналогом по технической сущности к заявляемому решению является "Способ обнаружения аномалий подстилающей поверхности", патент РФ 2160912, кл. G 01 V 8/00, 9/00, 2000 г. В способе ближайшего аналога получают изображение подстилающей поверхности в виде цифровой матрицы функции яркости А(х, у) от пространственных координат, разбивают изображение на мозаику сравнительно однородных по яркости участков на основе априорных данных, вычисляют функции фрактальной размерности каждого участка, составляют матрицу эталонов из коэффициентов фрактальной размерности характерных участков, а аномалию фиксируют по выходу разницы между текущим и эталонным значением фрактальной размерности за пороговый уровень для анализируемого участка.

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

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

Решение поставленной задачи обеспечивается тем, что в способе предсказания землетрясений, при котором осуществляют регистрацию собственного излучения подстилающей поверхности в виде зависимости амплитуды сигнала А(х, у) от пространственных координат, методами пространственного дифференцирования выделяют контуры на изображениях, проводят фрактальный анализ фрагментов изображений внутри выделенных контуров, дополнительно собственное излучение принимают антенной с линейной поляризацией по двум разнесенным по диапазону частот каналам приема, формируют синтезированную матрицу изображения из попиксельных отношений амплитуд сигналов этих каналов, контур очага выделяют на синтезированном изображении, вычисляют азимуты линеаментов поля градиентов изображения внутри контура очага и их среднюю сумму, рассчитывают по серии получаемых изображений период Т_о (час) изменения средней суммы азимутов, прогнозируют время удара, отсчитываемое от начала периодических изменений средней суммы азимутов и его магнитуду по регрессионным зависимостям t_у≅4,7Т_о, М≅110/Т_о ².

Изобретение поясняется чертежами, где: фиг. 1. Вариации направлений осей сжатия очага накануне землетрясения по данным наземных измерений; фиг.2. Визуализированное изображение синтезированной матрицы;
фиг.3. Линеаментное поле градиентов синтезированной матрицы;
фиг.4. Функция изменения результирующего азимута линеаментов;
фиг.5. Функциональная схема устройства, реализующего способ.

Введенные операции, образующие способ как совокупность существенных признаков, обеспечивают достижение таких качественных свойств, как:
- статистическая устойчивость метода за счет использования при обработке дополнительных независимых признаков синтезированного изображения;
- достоверность предсказания за счет прямого расчета параметров предстоящего удара: магнитуда и времени.

Это позволяет судить о соответствии заявляемого технического решения критерию "изобретательский уровень".

Техническая сущность изобретения заключается в следующем. Как отмечалось выше, в области очага землетрясения наблюдаются аномалии различных физических полей: гравитационного, магнитного, электростатического, механических напряжений земной коры. При механических напряжениях среда становится анизотропной (см., например, Фотоупругость, Физический энциклопедический словарь под редакцией А.М. Прохорова, М.: Сов. энциклопедия, 1983 г., стр.827). Величина анизотропии пропорциональна механическим напряжениям. Накануне землетрясения, вследствие раскачки очага, изменяется картина поля механических напряжений и, как результат, величина анизотропии пород. На фиг.1 (I-VI) представлена иллюстрация вариации направлений осей сжатия для роя землетрясений Гармского района, СССР, 1966 г., по данным наземных измерений (см., например, Т. Рикитаке Предсказание землетрясений, 1976 г., стр.309). Вследствие анизотропии изменяется поляризация собственного восходящего излучения поверхности над очагом. Преимущественно линейная поляризация приурочена к участкам поверхности, расположенным по направлениям оcей сжатия очага. Для других участков по периферии осей сжатия поляризация излучения менее выражена.

Известен также эффект Фарадея, заключающийся во вращении плоскости поляризации линейно поляризованных электромагнитных волн при их распространении в среде, находящейся в постоянном магнитном поле. В частности, эффект Фарадея наблюдается при распространении электромагнитных волн через плазменные образования, например ионосферу, находящуюся под воздействием постоянного магнитного поля Земли (см., например, Фарадея эффект, Советский энциклопедический словарь под редакцией А.М. Прохорова, 4-е изд., M.: Сов. энциклопедия, 1989 г., стр.1414). Кроме того, имеет место и вращательная дисперсия, состоящая в зависимости угла поворота плоскости поляризации электромагнитных волн от волнового числа к = 2π/λ, т.е. Δϕ = ϕ(2π/λ) (см., например, А.Л. Драбкин, В.Л. Зузенко Антенно-фидерные устройства, M.: Сов. радио, 1961 г., стр. 740). Таким образом, поляризация собственного восходящего излучения поверхности над очагом связана как с вариациями поля механических напряжений, так и с диапазоном частот регистрируемого излучения.

Следует ожидать, что закрутка вектора поляризации собственного излучения подстилающей поверхности будет существенно отличаться в двух каналах приема, если эти каналы разнести по диапазону частот. Например, для радиометра спутника NОАА, США, регистрирующего собственное излучение подстилающей поверхности в диапазоне от 2 до 10 мкм, по пяти независимым частотным каналам наибольшая разница поляризации принимаемого сигнала будет в первом (2 мкм) и пятом (10 мкм) каналах. Кроме того, вследствие анизотропии сжатой породы над очагом поляризация восходящего излучения изменяется от участка к участку или, применительно к получаемым изображениям зондируемой поверхности, от пиксела к пикселу. При двухканальном приеме на общую антенну с линейной поляризацией разница в поляризации сигналов преобразуется в разницу амплитуд пикселей идентичных участков изображений. Для подчеркивания контраста синтезированного изображения вычисляют попиксельные отношения амплитуд сигналов полученных изображений. Синтезированное изображение очага, полученное постфактум (т. е. после состоявшегося землетрясения по его известным координатам из архива запрашивались изображения подстилающей поверхности за несколько суток до удара), иллюстрируется фиг.2. Синтезированное изображение очага, опосредованно отражающее вариации поляризаций сигнала, содержит информацию об ориентации поля его механических напряжений.

Следующей задачей является извлечение скрытой информации из синтезированного изображения. Поскольку амплитуды пикселов синтезированного изображения пропорциональны действующим механическим напряжениям, следовательно, производные по координатам функции яркости изображения отражают ориентацию этих напряжений (ориентацию осей сжатия) (см., например, Производная по направлению в книге Н.С. Пискунов Дифференциальное и интегральное исчисления, учебник для ВТУЗов, 5-е изд., т.1, М.: Наука, 1964 г., стр.264-266). Производные по направлению функции изображения А(х, у) задают поле градиентов. По определению векторное поле градиентов скалярной функции синтезированного изображения А(х, у) для любой точки вычисляют как

На фиг. 3 представлено линеаментное поле градиентов синтезированной матрицы изображения внутри выделенного контура очага при длине штриха линеамента 10 пикселей. Программа обработки для получения линеаментного поля градиентов синтезированного изображения представлена ниже, в примере реализации. Длина штриха каждого линеамента векторного поля устанавливается во входных параметрах программы. При известных координатах верхнего и нижнего концов линеамента (х_i, у_i) и (х_j, у_j) длина линеамента определяется по формуле
.

Азимут линеамента

Средневзвешенная сумма азимутов линеаментов векторного поля определена как функция

Возникающий накануне удара переходный колебательный процесс (раскачка очага) изменяет картину механических напряжений по пространству очага. Последнее отражается в периодических изменениях функции А(х, у) синтезированного изображения и средневзвешенной суммы азимутов линеаментов. Получая серию последовательных во времени снимков контролируемого региона, отслеживают динамику изменения функции совпадающую по периоду с колебаниями очага. Рассчитывают период изменения функции Для этого формируют последовательный ряд значений функции В соответствии с теоремой отсчетов Котельникова-Шеннона непрерывная функция однозначно определяется своими дискретными значениями, которые отсчитываются через интервалы времени Δt=1/2 F_max, где F_max - максимальная частота спектра функции, при этом среднеквадратическая ошибка восстановления функции по дискретным отсчетам зависит от интервала дискретизации (см., например, Марпл С.П. Цифровой спектральный анализ и его приложения, перевод с англ., М.: Мир, 1990 г., стр.181). При известном (см. , например, Давыдов В.Ф. Землетрясения. Телеметрия предвестников, M. : изд-во МГУЛ, 2001 г., стр.30, рис.17) периоде раскачки очага (3,2. .. 4,5 час) и допустимой среднеквадратической ошибке дискретизации 0,1 максимальный интервал Δt обновления изображений контролируемого очага составит ≈ 750 сек. Чтобы обеспечить запас времени для оповещения населения о предстоящем ударе, период функции рассчитывают на раннем этапе наблюдений. Аппроксимируя в первом приближении периодическую функцию синусоидой вида α₀•sinωt, вычисляют конечные разности первого порядка:

Составляют ряд отношений:

который представляет собой последовательность производных от функции (α₀sinωt)ʹ = α₀•ω•cosωt в моменты времени t₁, t₂, . . . , t_n. Поскольку ω = 2π/T, то, имея последовательный ряд дискретных отсчетов функции и ее производных во времени, программным методом вычисляют период функции Т_о. Восстановление гармонических функций, заданных последовательностью отсчетов на некотором участке, представляется стандартной математической операцией, входящей в комплект специализированного программного обеспечения (см., например, Дьяконов А. П., Справочник по МАТН САД, М.: Ск "Пресс", 1998 г., Использование функции предсказания Predikt, с.242, п.11.15). Восстановленная по дискретным отсчетам периодическая функция изменения результирующего азимута линеаментов иллюстрируется фиг.4. Частота колебаний механической системы задается выражением где с - жесткость упругого элемента (породы), m - колебательная масса. Параметры с, m являются характеристиками очага землетрясения, а их соотношение с/m определяет магнитуду М. Магнитуда удара обратно пропорциональна квадрату периода М=М(1/Т²). Регрессивная зависимость магнитуды от периода (см. патент РФ 2170447, 2001 г.) М= 110/Т². Раскачка очага как переходный процесс подчиняется общим законам генерирования и установления колебаний (см., например, Зернов Н.В., Карпов В.Г. Теория радиотехнических цепей, учебник, М.-Л.: Энергия, 1963 г., стр.407) и задается формулой t_y = π/dekr•T₀ где dekr - логарифмический декримент затухания, равный натуральному логарифму отношения амплитуд колебаний двух смежных периодов. Переходный процесс считают законченным, когда амплитуда колебаний составляет порядка 1% от максимальной. Тогда t_у=4,7Т_о.

Пример реализации способа.

Заявляемый способ может быть реализован на базе устройства по схеме фиг. 5. Функциональная схема устройства фиг. 5 содержит систему орбитальных спутников наблюдения 1 типа NOAA (США) с установленными на них сканирующими радиометрами 2 типа AVNRR. Радиометры осуществляют прием восходящего излучения подстилающей поверхности в полосе сканирования 3 по пяти независимым каналам одновременно.

Зарегистрированный радиометром сигнал с сопутствующей информацией (время приема, эфемериды) передается по открытой автономной радиолинии 4 и принимается в режиме открытого доступа наземными пунктами приема 5, расположенными в г. Красноярск, г. Москва (ИКИ), где записывается на видеомагнитофон 6 типа "Арктур". Из принятых изображений подстилающей поверхности в пункте приема 5 на основе сопутствующей информации формируется база данных изображений регионов, которая в виде файлов выводится на сервер 7 сети "Интернет" (адрес сервера www http: \\ SmiS. IkI. RSSI.RU). По запросу потребителей информация в виде кадров перекачивается через "Интернет" и автономный сервер 8 в индивидуальный архив накопления, хранения и выдачи информации 9 Центра мониторинга МЧС 10. Архив 9 может быть организован на базе стриммеров типа FT-120. Федеральной службой мониторинга 10 МЧС ведется база всех сейсмоопасных районов страны, которая позволяет выборочно отслеживать наиболее опасные регионы и оперативно формировать последовательность отснятых радиометрами AVNRR кадров изображений этих регионов. Обработку накопленных массивов данных изображений сейсмоопасных регионов осуществляют на персональных ПЭВМ 11 типа IBM в стандартном наборе элементов: процессора-вычислителя 12, винчестера 13, оперативного ЗУ 14, принтера 15, дисплея 16, клавиатуры 17. Предварительно специализированный комплекс программ обработки изображений типа МАТН САД 7.0 PLUS записывают на винчестер 13. Попиксельное вычисление отношения амплитуд сигналов в двух разнесенных по диапазону волн изображений осуществляют стандартной процедурой программного расчета (см., например, Векторизация элементов матрицы, МАТН САД, 7.0 PLUS, издание 3-е стереотипное, М.: Информ. -издат. дом Филинъ, 1998 г., стр.211). Синтезированное изображение подстилающей поверхности, полученное попиксельным расчетом отношения квантованных уровней цифровых изображений I и V каналов радиометра AVNRR, иллюстрируется фиг.2. Далее осуществляют, по операциям ближайшего аналога, методами пространственного дифференцирования и фрактального анализа выделение контура очага на синтезированном изображении. Для участка изображения внутри контура проводят линеаментный анализ, т.е. программным методом рассчитывают направления линеаментов поля градиентов. Расчет реализуется программой, представленной в конце описания.

По серии последовательных во времени синтезированных изображений очага отслеживают изменение средневзвешенной суммы азимутов и рассчитывают период этой функции. По расчетному периоду прогнозируют магнитуду и время ожидаемого удара наблюдаемого очага. Известные способы предсказания землетрясений (патенты РФ 2170446, 2170447, 2001 г.) имеют существенную зону нечувствительности, ограниченную среднеквадратической ошибкой измерения амплитуды раскачки очага в несколько см. Заявленный способ отличается высокой чувствительностью, поскольку динамическая раскачка очага даже на несколько мм сразу изменяет статику механических напряжений по пространству очага и связанную с этим средневзвешенную сумму азимутов линеаментов.

Формула изобретения

Способ предсказания землетрясений, при котором осуществляют регистрацию собственного излучения подстилающей поверхности в виде зависимости амплитуды сигнала А(х, у) от пространственных координат, методами пространственного дифференцирования выделяют контуры на изображениях, проводят фрактальный анализ фрагментов изображений внутри выделенных контуров, отличающийся тем, что собственное излучение принимают антенной с линейной поляризацией по двум разнесенным по диапазону частот каналам приема, формируют синтезированную матрицу изображения из попиксельных отношений амплитуд сигналов этих каналов, контур очага выделяют на синтезированном изображении, вычисляют азимуты линеаментов поля градиентов изображения внутри контура очага и их среднюю сумму, рассчитывают по серии получаемых изображений период Т_о (ч) изменения средней суммы азимутов, прогнозируют время удара (t_у), отсчитываемое от начала периодических изменений средней суммы азимутов и его магнитуду (М) по регрессионным зависимостям
t_у≅4,7 Т_о,
М≅110/Т_о ².

РИСУНКИ

Рисунок 1, Рисунок 2, Рисунок 3, Рисунок 4, Рисунок 5, Рисунок 6, Рисунок 7, Рисунок 8