Цифровая обработка в ZETLAB

Землетрясение — сильное колебание земной коры, вызываемое причинами вулканического, тектонического характера [1] или искусственными процессами: взрывы, заполнение водохранилищ, обрушением подземных полостей горных выработок [2].

В России зоны повышенной сейсмической опасности (от 6 баллов и выше, с периодом повторяемости 500 лет) занимают около 40% от общей площади, в том числе 9% территории относится к 8-9 балльным зонам, карта сейсмической активности, в некотором приближении, приведена на рисунке 1. Стоит отметить, что в сейсмически активных зонах проживает более 20 млн.человек. Ежегодно на территории России происходят десятки и даже сотни землетрясений, большинство из них протекают незаметно, оставаясь таковыми как для рядовой общественности, так и для исследователей, изучающих протекающие в земле тектонические процессы. Однако, в новейшей истории России, три были сильнейшими: Шикотанское (Курилы), октябрь 1994 года, Нефтегорское, май 1995 года, Кроноцкое, декабрь 1997 года [3].

Рис.1. Карта сейсмической активности на территории РФ

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

Причина затруднений практической реализации мониторинга заключается в том, что малые сейсмические воздействия лишь немного превышают уровень микросейсм, что так же осложняется наличием помех, электрических наводок. На основе экспериментальных исследований в США разработаны нормы для размещения сейсмостанций относительно окружающих источников естественных техногенных шумов [GSE\WG2/2], эти нормы приведены в таблице [4]:

Нормы удаленности сейсмостанции от источников шума

Источник шумов	Расстояние от источника, км
1. Основные линии электропередач	10 — 15
2. Основные высокоскоростные дороги	5 — 10
3. Второстепенные дороги	1 — 3
4. Трубопроводы, насосные станции и т.д.	3 — 5
5. Дороги к лесозаготовкам и шахтам	1 — 3
6. Карьеры и большие шахты	10 — 20
7. Большие реки	3 — 5
8. Большие озера, океаны.	50 — 100
9. Бурение глубоких скважин	5 — 10
10. Города и небольшие деревни	10 — 15

Все приведенные в таблице данные справедливы и применяются до сих пор. Стоит отметить, что в этот перечень потенциально «вредных» для качества сигнала объектов абсолютно отнесли источники техногенного шума: 1-6 , 9, 10. Но так же и природные источники шумов — 7, 8.

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

Система предназначена для определения самого факта произошедшего землетрясения, в том числе и минимальной магнитуды, и самое главное — поиск начала события с точностью до нескольких отсчетов АЦП, для последующей передачи этой информации в алгоритм пеленгации эпицентра землетрясения, описание которого в данной статье рассматриваться не будет. Отметим, что система требует точного обнаружения времени начала события, в том числе и в условиях сезонных изменений естественного фона шумов, поскольку шумы варьируют не только в пространстве, но во времени т.к. любая станция подвержена суточным и сезонным вариациям [4], что гарантирует повышение точности в задаче поиска эпицентра землетрясения. Это свойство наиболее важно при анализе данных несколькими группами сейсмостанций одновременно, когда две или более сейсмостанций работают совместно в целях обнаружения эпицентра.

Помехи и электрические наводки, являясь шумами техногенного характера, зачастую, детерминированы, либо выходят за границы частотного диапазона полезного сигнала, это говорит о том, что от них, сравнительно легко, можно избавиться различными методами фильтрации: полосовая фильтрация, а так же методы адаптивной фильтрации. Некоторые алгоритмы адаптивной фильтрации очень неплохо описаны в статье А.Б. Сергиенко «Алгоритмы адаптивной фильтрации: особенности реализации в MATLAB» [5], однако, практическая реализация таких систем довольно затруднительна.

Для решения поставленной задачи будем оперировать реально записанным материалом, а именно, будем анализировать запись землетрясения, произошедшего в районе оз. Байкал 20 августа 2010 года (рисунок 2):

Рис.2. Зона расположения сейсмостанций (выделена цветом)

Сигнал был зарегистрирован тремя сейсмостанциями, установленными в районе оз. Байкал. Приблизительная схема регистрации приведена на рисунке 3:

На схеме регистрации (рисунок 3) исходные данные обрабатываются и архивируются на сейсмостанции, передавая на сервер только измеренные параметры: наличие события и время начала события. Для работы системы выбран трехкомпонентный (X,Y,Z) датчик, синхронно регистрирующий отрабатывающий сигнал по всем измерительным осям, что при взаимном анализе показаний по каждой из осей датчика, дает возможность определить вектор распространения сейсмической волны, возникшей в результате землетрясения. Оцифрованные данные пересылаются на верхний уровень по технологии OPC (стандартный протокол, основу которого составляет DCOM, предназначен для обмена данных по промышленной сети Ethernet, используется большинством производителей, как аппаратуры, так и программ для обеспечения стандартного обмена данными между элементами системы).

Рис.3. Схема регистрации сейсмического сигнала и передачи на вычислитель

Итак, мы имеем первичный преобразователь (трехкомпонентный датчик), сейсмостанцию (в приближении — это плата АЦП + промышленный компьютер, обеспечивающий математическую обработку и передачу данных на верхний уровень по технологии OPC). Задача состоит в том, что б обработать сигнал в реальном времени или записать исходный сигнал, проанализировать его и однозначно определить наличие события (землетрясение), если событие обнаружено, то необходимо указать время его начала с минимальной погрешностью но так же необходимо указать и длительность события.

Составим структурную схему (рисунок 5, сверху) потоков данных на нижнем уровне — уровень сейсмостанции включает в себя датчик, спутниковую систему синхронизации данных (основа — GPS или ГЛОНАСС), модуль АЦП. Информация в реальном масштабе времени поступает с датчика на модуль АЦП и далее передается на промышленный компьютер, где и обрабатывается проектируемым алгоритмом.

Рис.5. Схема потоков данных. Источник информации — сейсмодатчик. Антенна ГЛОНАСС в составе спутниковой системы необходима для синхронизации данных между разными сейсмостанциями.

Для первоначальной оценки параметров сигнала воспользуемся методом визуального анализом временной реализации данного сигнала, т.е. сейсмограммой (рисунок 6) сто секундной записи. На рисунке 6 представлены записи с трех компонент сверху вниз по порядку — ось X, ось Y, Ось Z. Видно, что приблизительно на 40 секунде виден всплеск сигнала, (для справки, соответствует московскому времени 17:31 20 августа 2010 года). В качестве рабочей записи примем этот временной диапазон, и все математические операции по идентификации сейсмического сигнала будем проводить именно с этим сигналом. Отметим, что в окнах сейсмограмм по оси абсцисс используется относительное время, а не абсолютное.

Рис.6. Сейсмограммы. Компоненты одного из сейсмодатчиков: X-верхний, Y-средний, Z-нижний графики. Показания до полосовой фильтрации. Данные представлены в относительном времени

Проанализируем шумовую составляющую записи, то есть ту часть, где заведомо не было полезного сигнала. На рисунке 7 представлен спектральный анализ сигнала — график спектральной мощности шумов в децибелах (относительно опорного значения для расчета дБ, принятого в случае единиц измерения ускорения в м/с², что составляет 0.0003 м/с²) в зависимости от частоты.

Параметры расчета:

• метод расчета — дискретное преобразование Фурье;
• от 0.01 Гц до 40 Гц;
• частотное разрешение 0.01 Гц;
• изображен максимальный график за 100 секунд усредненных по 5 секунд различных реализаций сигнала;

На рисунке 8 представлен более подробный график (за счет выбора иного частотного диапазона — в 10 раз меньше) со следующими параметрами:

• метод расчета — дискретное преобразование Фурье;
• от 0.01 Гц до 4 Гц;
• частотное разрешение 0.01 Гц;
• изображен максимальный график за 100 секунд усредненных по 5 секунд различных реализаций сигнала;

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

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

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

• ФВЧ 10 Гц;
• ФНЧ 1 Гц;

Фильтр реализуем с применением программы ZETFormula, входщей в состав ПО ZETLab, текст формулы, представлен на рисунке 16. Результат работы полосового фильтра приведен на рисунке 9 — спектральная мощность сигнала.

Рис.7. Выделен высокочастотный участок (относительно полезной частоты в нашей задаче), в данной задаче является шумами

Рис.8. Выделен низкочастотный участок (относительно полезной частоты в нашей задаче), в данной задаче является шумами, это микросейсмы (шумы океанов, крупных озер и т.д. [4])

Рис.9. Спектр шумовой составляющей сигнала по оси X (до события)

Сравнительным анализом временной реализации данных по осям X,Y и Z до и после проведения цифровой фильтрации в полосе 1-10 Гц (на рисунках 6 и 10) виден существенный выигрыш — теперь факт наличия события отчетливо виден на сейсмограммах (рисунок 10). Однако, момент начала события можно определить лишь с большим приближением (рисунок 11):

• ± 1-3 секунд для показаний по оси X;
• ± 1-5 секунд для показаний по оси Y;
• по оси Z факт наличия события невозможно определить (либо с погрешностью более 10 секунд);

Рис.10. Сейсмограммы. Компоненты одного из сейсмодатчиков: X-верхний, Y-средний, Z-нижний графики. Показания после полосовой фильтрации. Данные представлены в относительном времени

Рис.11. Компоненты одного из сейсмодатчиков: X-верхний, Y-средний, Z-нижний графики. Увеличение масштаба по оси времени. Данные представлены в относительном времени. Зона неопределенности выделена голубым цветом.

После полосовой фильтрации соотношение сигнал\шум увеличилось, но все же начало события с точностью до нескольких отсчетов АЦП определить невозможно. Далее, предварительно отфильтрованный сигнал, обработаем методом анализа отношения амплитуд в коротком и длинном временных окнах — метод STA\LTA (STA/LTA — Short Time Average to Long Time Average). Впервые данный метод был описан в теоретической работе Фрейбергера. Эти методики требуют небольшого объема вычислений, что является существенным аргументом использования их в системах реального времени [6]. Результат работы алгоритма из статьи [6] приведен на рисунке 12. Результат собственной реализации алгоритма приведен на следующем рисунке 13.

Приведем формулу расчета критерия STA\LTA:

где t₁<2, MX — математическое ожидание амплитуды сигнала за разное время наблюдения, например STA на участке 1 секунда, LTA на участке 10 секунд.

где t = ih (i=1,2,…) — дискретное время, 1/h — частота дискретизации, NS и NL- число отсчетов в коротком и длинном временных окнах.

Временная сложность детектора STA/LTA равна O(N), где N — число отчетов анализируемой записи. Столь малая временная сложность объясняет широкое применение детектора в геофизических системах, функционирующих в реальном времени [6].

Рис.12. Наглядная демонстрация алгоритма. Данные представлены в относительном времени

Рис.13. Результат работы программы STA\LTA. Данные представлены в относительном времени. Формула расчета выполнена на программном обеспечении ZETLab

Как видно из рисунка 13, четкость в определении времени начала достигнута. Теперь необходимо добиться автоматического детектирования начала события. В работе Фрейбергера сказано, что порог детектирования задается как константа.

Очевидны существенные недостатки такого алгоритма:

• возможны ложные срабатывания при низком значении пороговой линии для повышения точности работы;
• влияние межсезонных и суточных вариаций [4], тем самым может изменяться точность алгоритма во времени;
• при заранее высоком уровне порога снижение точности метода;

Как решение, устраняющее все вышеперечисленные недостатки, предлагаю использование адаптивного алгоритма обнаружения, когда порог подстраивается автоматически, в зависимости от статистических характеристик сигнала за время наблюдения NKL, где время предыстории сигнала NKL>NL в К раз, то есть NKL=NLxK.

Рис.14. Переходный процесс во время накопления статистики NKL. Зеленый график — порог, красный график — критерий STA\LTA. Данные представлены в относительном времени.

Рис.15. Пересечение графика адаптивного порога срабатывания с графиком критерия STA\LTA. Данные представлены в относительном времени.

Рис.16. Формула, реализующая алгоритм адаптации порога срабатывания. Формула расчета выполнена на программном обеспечении ZETLab

В качестве заключения отметим, что предложенный метод адаптации детектора STA/LTA позволяет автоматически определять участок сейсмограммы, что обеспечивает высокую точность в определении момента вступления волн. Математический аппарат ZETLab позволяет настраивать систему в кратчайшие сроки, позволяя изменять различные параметры алгоритма во время системы, без реконфигурации работы аппаратной части и остальных программ комплекса программ. Стоит отметить, что адаптация алгоритма STA\LTA с высокой точностью результата может применяться в условиях плавно нарастающего фронта сигнала и в условиях сезонных вариаций фонового шума. Так же, данный метод может быть применен и в прикладных областях, где необходимо определение момента вступления волн, например, при работе с датчиками акустической эмиссии, в методах неразрушающего контроля.

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

Автор статьи: Красовский А.А.
Статья опубликована в журнале «Цифровая обработка сигналов» № 3/2010.

Список использованной литературы:

1. Д.Н.Ушаков. Толковый словарь русского языка. Москва : Гос. ин-т «Сов. энцикл.»; ОГИЗ; Гос. изд-во иностр. и нац. слов., 1935-1940.
2. Интернет-энциклопедия. Землетрясение. ru.wikipedia.org. [В Интернете] [Цитировано: 18 08 2010 r.] //ru.wikipedia.org/wiki/Землетрясение.
3. МЧС. www.spas-extreme.ru. Интернет-портал МЧС. [В Интернете] [Цитировано: 18 08 2010 r.] //www.spas-extreme.ru/schoolsafety/prirodchs/zemletrysenie/russia.php.
4. О.К.Кедров Сейсмические методы контроля ядерных испытаний. Москва, Саранск : Институт физики Земли им. О.Ю. Шмидта; Рос. Акад. Наук, 2005.
5. статья «Алгоритмы адаптивной фильтрации: особенности реализации в MATLAB». А.Б.Сергиенко б.м. : Математика в приложениях, № 1 2003 r.
6. Автоматическое определение длительности сейсмического события в режиме реального времени. С.В.Баранов Москва : Сборник статей аспирантов, соискателей, докторантов и научных работников, № 3 2004 r.