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

Полный текст

Введение

Решение задачи поиска предсказания землетрясений и цунами на основании изучения изменений, наблюдаемых физических полей Земли, сегодня активно развивающаяся область исследований в геофизике[1–3]. В частности, изучение аномалий акустических полей в верхних слоях литосферы одно из направлений в этой области, весьма перспективно [4–6]. Исследовательской базой этого направления служит анализ динамических характеристик измеряемых параметров акустических шумов, образующихся под действием сил, порождаемых внутренними процессами Земли, которые проецируются на ее поверхность в виде продольных и поперечных волн. Возникающая при этом деформация в почвенном слое вызывает образование высокочастотного акустического сигнала, который получил название геоакустической эмиссии (ГАЭ) [7]. Сигналы ГАЭ опосредовано отражают процессы, порождающие сейсмические явления и поэтому в определенной степени могут служить базой для их описания. Преобразованные датчиками в электрические сигналы, ГАЭ представляет собой импульсный поток в диапазоне частот от единиц герц до сотен килогерц. Поскольку природа известных сейсмических явлений описывается нелинейными динамическими процессами [8], то генерируемые при этом поверхностные сейсмоакустические сигналы носят нестационарный стохастический характер. Этот факт накладывает известные ограничения на стремление получить устойчивые признаки возникающих сейсмических явлений с использованием статистических наблюдений, которые описаны моделями нелинейной динамики [9–11]. Кроме того, на результаты анализа наблюдений существенное влияние оказывают климатические и погодные условия, накладывающие на сигнал дополнительный акустический шум в точке наблюдения.

В настоящее время в Институте космофизических исследований и распространения радиоволн дальневосточного отделения Российской академии наук (ИКИР ДВО РАН) накоплена база измерений ГАЭ. На основании проводимых теоретических исследований в период с 2001 года по настоящее время сформулирована и проверяется гипотеза о существовании связи обнаруженной изменчивости характеристик сигналов ГАЭ с фиксируемыми изменениями состояния приповерхностного слоя литосферы [12–15]. Поиск подтверждающих фактов, указывающих на правомерность этой гипотезы, является основной целью проводимых исследований. Предполагается, что в силу влияния сейсмических волн на деформационные процессы в неоднородной среде литосферного слоя, в нем происходит растрескивание среды, порождающее акустические импульсы. Это явление отражается на интенсивности генерации потока импульсов акустической эмиссии определенным образом, связанным с характером сейсмических волн. В результате сейсмическое воздействие меняет состояние слоя, что отражается в динамике характеристик сигналов ГАЭ (интенсивность потока, амплитудно-фазовые показатели импульсов, скважность, частота следования, частота модуляции импульсов). Наблюдая изменение характеристик, можно попытаться выявить определенные закономерности изменений этих показателей, и связать их с определяющим фактором, в интересуемом случае - с параметрами сейсмических волн. Изменения характеристик сигнала ГАЭ с 2001 года фиксируются в Лаборатории акустических исследований ИКИР ДВО РАН в ходе наблюдений, как изменение статистических распределений амплитуды, межимпульсных интервалов и концентрации импульсов во времени. В общем случае задача исследования формулируется, как стремление качественно описать наблюдаемую изменчивость свойств среды под действием деформаций, вызванных воздействием внешних факторов. Изменения свойств среды фиксируются опосредовано через обнаружение аномального поведения динамики характеристик сигналов ГАЭ.

Пусть за некоторое время $\Delta t$ в результате измерений получаем вектор значений определенного параметра сигнала ГАЭ (например, серии межимпульсных интервалов, амплитуд, значений частотной или амплитудной модуляции, или их производных) $X=(X_1,X_2,\dots ,X_u)$. Обозначим место вектора, как точку в u-мерном фазовом пространстве. Получим n последовательных измерений и разместим их в том же фазовом пространстве. В результате получим фазовую траекторию, описывающую динамику выбранного параметра за время $T=n\Delta t$. Очевидно, что траектория однозначно описывается матрицей вида

 

$U=\left\{\begin{array}{cccc}{X}_{\mathrm{1,}u}& X_{\mathrm{2,}u}& \cdots & X_{n,u}\\ ⋮& ⋮& ⋮& ⋮\\ X_{\mathrm{1,2}}& X_{\mathrm{2,2}}& \cdots & X_{n\mathrm{,2}}\\ X_{\mathrm{1,1}}& X_{\mathrm{2,1}}& \cdots & X_{n\mathrm{,1}}\end{array}\right\}.$ (1)

Будем называть (1) отображением динамики выбранного параметра сигнала ГАЭ с дискретностью $\Delta t$ в масштабе измерения T, или матрицей-отображений динамики параметра, или просто — отображением. В ходе наблюдений в фазовом пространстве складывается непрерывная картина динамики характеристик ГАЭ, которую можно представить как последовательное соединение фрагментов траектории вида (1). Согласно законам нелинейной динамики, состояния открытых систем отображаются в фазовом пространстве некоторыми петлями сложных траекторий [16–18]. При достаточно длительном наблюдении за изменениями состоянием системы, следует ожидать появление самоподобных траекторий, имеющими сходное описание. Этот факт создает прецедент для поиска закономерности смены состояния системы и, как следствие, возможность обнаружения закономерностей в сигнале на уровне описания динамики его характеристик через отображения. Степень сходства описания можно определить через сравнение отображений вида (1). Однако в данном случае при выборе метода сравнения отображений следует учитывать специфику их построения, связанную с порядком их составления из векторов параметров.

$U=\left\{\begin{array}{c}{a}_1,a_2,a_3,\cdots a_n\end{array}\right\},$ (2)

где $a_i=\left\{X_{i,j}\right\}$, $j=\mathrm{1,}u$. То есть, существует определяющая порядок построения отображения временная ось, и поэтому корректное сравнение должно осуществляться синхронным попарным сопоставлением векторов сравниваемых отображений вдоль временной оси. Примем следующее определение степени сходства отображений: два равно размерных отображения A и B подобны до степени h, равного сумме коэффициентов корреляции Пирсона попарно выбираемых векторов сравниваемых матриц $a_i$ и $b_i$ вида (2):

 $h(A,B)=\sum _{i=1}^{n}r(a_i,b_i)=\sum _{i=1}^n\frac{\mathrm{cov}(a_i,b_i)}{{\sigma }_{a_i}{\sigma }_{b_i}}$ (3)

где ${\sigma }_a,{\sigma }_b$ — дисперсии векторов a и b, соответственно. Из свойства корреляции вытекает, что максимальное значение степени сходства отображения шириной в n векторов равно $h_{\max }=\left|n\times h\right|$. Если $h(A,B)=-h_{\max }$, то каждая из матриц является негативным отображением по отношению к другой. Если $h(A,B)=0$, это означает, что отображения полностью не совпадают.

Разбивая эпоху наблюдения сигнала ГАЭ на последовательные фрагменты равной длительности Т, сформируем множество отображений $\left(U_m\right)$, отражающих динамику выбранного параметра сигнала ГАЭ за эпоху анализа длительностью $T\times M$. Сравнивая между собой попарно отображения полученного множества фрагментов по критерию (3) получим симметричную матрицу значений $r_k(k=\mathrm{1,}M)$ степеней совпадений множества отображений $\Omega$ (4):

$\Omega =\left\{\begin{array}{ccccc}1& r_{\mathrm{2,1}}& \cdots & r_{n-\mathrm{1,1}}& r_{k\mathrm{,1}}\\ r_{\mathrm{1,2}}& 1& \cdots & r_{n-\mathrm{1,2}}& r_{k\mathrm{,2}}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ r_{\mathrm{1,}k-2}& r_{\mathrm{2,}k-2}& \cdots & r_{k-\mathrm{1,}k-2}& r_{k,k-2}\\ r_{\mathrm{1,}k-1}& r_{\mathrm{2,}k-1}& \cdots & 1& r_{k,k-1}\\ r_{\mathrm{1,}k}& r_{\mathrm{2,}k}& \cdots & r_{k-\mathrm{1,}k}& 1\end{array}\right\}.$ (4)

Зададимся значением порога сходства отображений ${\lambda }_0$, и применим правило отбора сравниваемых значений 5 построчно:

$v=\left\{\begin{array}{cc}\mathrm{1,}& \lambda >{\lambda }_0\\ 0& \end{array}\right..$ (5)

В результате получим бинарную матрицу, где любая строка, содержащая хотя бы одну <1> вне главной диагонали, будет обозначать кластер, объединяющий все отображения с номерами элементов этой строки, в которых стоят <1>. Представленная последовательность вычислений по существу представляет алгоритм кластеризации, в котором критерием близости выступает корреляционная связь отображений вида (3). Формирование матрицы (4) позволяет выделить кластеры отображений с характерными особенностями динамики характеристик ГАЭ для их дальнейшего сопоставления с зафиксированными ранее сейсмическими событиями.

Таким образом, в настоящем исследовании задача обнаружения смены состояния среды фиксируется на основании распознавания отображений динамики характеристик наблюдаемого сигнала ГАЭ вида (4). Предшествующий распознаванию отображений этап их разбиения на классы носит трудоемкий методический и вычислительный характер. Среди эффективных алгоритмов кластеризации в современной обработке отображений выделяются инструменты, разработанные в рамках искусственного интеллекта. Достижения в этой области на основе структур нейронных сетей показали высокую эффективность их внедрения в практику цифровой обработки сигналов и распознавания образов. Что и явилось мотивацией к использованию нейронных сетей в настоящем исследовании на этапах препроцессинга сигналов, формирования тензоров отображений и их кластеризации.

Методика эксперимента

Особенности приема сигнала ГАЭ вблизи поверхности почвенного слоя земли, описанные в ряде более ранних работ [19, 20], состоят в том, что расчет вычисления статистических распределений амплитудно-частотных характеристик импульсов заменен на распределение характеристик, полученных по результатам сравнения структуры этих импульсов [21]. В проводимых нами исследованиях на этапе предобработки импульсного сигнала ГАЭ амплитудно-фазовая структура каждого составляющего сигнал импульса претерпевает инвариантное преобразование, позволяющее описать его структуру в виде бинарной матрицы. Это преобразование позволяет провести сравнительный анализ структур последовательности импульсов в сигнале ГАЭ и их кластеризацию по степени сходства матриц. В результате на фрагменте сигнала заданной длительности удается описать и определить численно разнообразие структур импульсов, а далее сформировать множество, составляющее классы - алфавиты. Каждый элемент алфавита связан со статистикой его встречаемости на фрагменте сигнала и представляет объект, включающий амплитудно-фазовую характеристику импульса. Далее выполняется ранжирование элементов алфавита по значению размерности символов и степени сходства их матриц [21]. Весь анализируемый сигнал разбивается на равные фрагменты, для которых вычисляются отдельные алфавиты и выстраивается их распределение по времени. Каждое полученное распределение элементов алфавита выкладывается на двумерную сетку, где по одной оси вносятся элементы алфавита по их размерности, а по другой выстраиваются алфавиты, полученные на последовательно нарезанных фрагментах сигнала. В итоге формируется трехмерное отображение изменчивости состава алфавитов по времени, где третья ось представлена значениями частоты встречаемости элементов определенного алфавита для каждого фрагмента сигнала. Таким образом сформированная матрица представляет отображение и позволяет наблюдать динамику характеристик сигнала ГАЭ. В данном случае в символы алфавита "упакованы" амплитудно-фазовые характеристики сигнала ГАЭ, которые последовательно отражены в алфавитах от одного временного фрагмента к другому.

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

Первый этап методики предусматривает подготовку данных и начинается с применения структурно-лингвистической обработки временных файлов, содержащих сигналы ГАЭ, которые были получены в ходе наблюдений на Камчатском геодинамическом полигоне в рамках выполнения государственного задания ИКИР ДВО РАН с 2016 по 2021 годы [22]. Каждый файл это отрезок сигнала, представленный 15 минутным его фрагментом. Длительность фрагмента принята в лаборатории акустических исследований ИКИР ДВО РАН за условный стандарт — минимально достаточный блок для обработки данных выбранными методами. Файл содержит от нескольких десятков до нескольких тысяч импульсов ГАЭ. В результате структурно-лингвистической обработки из каждого файла выделяется алфавит. После сортировки по значениям размерности элементов каждый алфавит представляется как вектор. Далее составляется ряд векторов, полученных из последовательно следующих друг за другом фрагментов сигнала. В памяти компьютера формируется блок размерностью 512 на 2880 точек, который визуально может быть представлен трехмерным графиком, составленным из 2880 векторов. Пример трех сформированных отображений представлен на рис. 1.

 

Рис. 1. Контроль результата работы функции предобработки отображений

[Figure 1. The сontrols result of the image preprocessing function]

 

На отображениях снизу вверх размещены элементы алфавитов по значениям их размерностей, слева направо приведена временная шкала, цветом отражены значения частоты появления элементов в алфавитах. Выбор ряда из 2880 векторов определен средне месячным периодом анализа наблюдений (4 файла * 24 часа * 30 суток) [24].

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

Общее число подготовленных таким образом данных за период наблюдений с 2017 по 2021 годы составил 60 графиков (12 месяцев * 5 лет). Практика исследований показала, что внешнее представление полученных графиков выглядит как отображение со сложным цветным узором. Визуальный анализ полученного материала практически затруднителен в виду большого разнообразия трудно подающихся описанию цветных паттернов в поле отображениях. В обход трудностей зрительного анализа было принято собрать тензор данных для обработки отображений нейросетью. Путем интерполяции подготовленные отображения преобразуются в растровое представление в формате Portable Network Graphics (PNG), масштабируются до размеров 64х64 пикселя. Формат новых отображений выбирается исходя из подгонки их к доступным параметрам входного слоя используемой архитектуры нейронной сети. Размер получившегося массива отображений недостаточен для качественной обработки данных нейронной сетью (60 отображений). По этой причине использовался метод генерации предыскаженных клонов отображений на базе существующих отображений. Варианты предыскажений отображений ограничены из-за вероятности потерять полезную информацию. которая привязана к масштабу и векторному построению исходных отображений. Использовались преобразования и трансформации, которые только смещали отображение по оси времени, а также производилась линейная корректировка яркости и контрастности, поскольку такие изменения допустимы. В этом состояла одна из особенностей обработки исходной базы отображений. Для каждого исходного подготовленного отображения была создана серия из 50 дубликатов отображений с предыскажениями, заложенными в генераторе аугментации. В результате выполненной предобработки итоговое число нового массива отображений составило 3000. Исходные 60 неискаженных отображений используются только на заключительном этапе распознавания отображений по классам.

Второй этап обработки данных направлен непосредственно на определение особенностей динамики характеристик сигналов ГАЭ путем выделение специфики паттернов полученных отображений, передаваемой изменениями палитры рисунка, и дальнейшем разбиении отображений на группы по степени сходства пространства признаков применением алгоритмов кластерного анализа.

Отображения относятся к образам так называемых «тепловых карт». Распознавание таких отображений имеет свои особенности, связанные с отсутствием на них определенных объектов. Место объектов занимают градиенты цветности или уровней серого, которые несут основную информационную нагрузку. Для выполнения кластеризации было принято решение провести воспользоваться алгоритмом Ллойда — (метод K-средних), использующий минимум-дисперсионный подход. Для возможности сравнения распределений по кластерам принудительное разбиение осуществлялось от 2 до 6 кластеров. Результат выполнения кластеризации представлен в табл. 1.

 

Таблица 1

Результаты работы алгоритма кластеризации [The results of the clustering algorithm]

Заданное число кластеров	Полученное число элементов в кластерах
2	929            2071
3	808        907        1285
4	800     863     523      814
5	516    863    271    788     562
6	861   336   561   463    538    268

 

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

 

Таблица 2

Результаты оценки эффективности кластерного анализа [The results of evaluating the effectiveness of cluster analysis]  

Число полученных кластеров	Значение Контурного коэффициента
2	0,722281
3	0,953454
4	0,859461
5	0,756904
6	0,745786

 

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

Следующий этап представляемой методики состоит в обучении нейронной сети на основании выделенных групп отображений по лучшему результату оценки качества кластеризации. Затем выполняется проверка результата обучения с использованием выборки 60 исходных неискаженных отображений на предмет качества классификации используемой нейронной сети. В качестве обучающей выборки данных был выбран результат кластеризации с 3 выделенными группами, файлам которых были присвоены соответствующие метки со значениями классов 0, 1, 2.

Для проведения эксперимента была разработана нейронная сеть с архитектурой типа ResNet. Для эффективной работы нейросети были проведены многочисленные эксперименты с изменением глобальных параметров, численностью сверточных слоев, числа остаточных связей, числа фильтров. Найденная оптимальная архитектура нейросети, на которой были получены лучшие результаты обучения и классификации представлена на рис. 2.

 

Рис. 2. Описание архитектуры нейросети

[Figure 2. Description of the neural network architecture]

 

Результаты

Результаты тестирования нейросети следующие: ошибка классификации на тестовой выборке составила 0,3426; точность классификации на тестовой выборке составила 0,9133. Для всех классов нейросеть показала уверенное распознавание отображений по классам обучения с малыми значениями ошибок классификации уже на 40 эпохе обучения. Нормированная матрица ошибок представлена на рис. 3.

 

Рис. 3. Матрица ошибок классификации

[Figure 3. Classification error matrix]

 

Нейросеть верно распознала 89% отображений в Классе 0, 98% отображений в Классе 1 и 87% отображений в Классе 2. Таким образом средняя точность распознавания составила 91 %. Графики обучения модели на рис. 4 и 5 показывают ровную тенденцию без всплесков переобучения.

 

Рис. 4. Графики, показывающие долю верных ответов в процессе обучения модели на 64 эпохах

[Figure 4. Graphs showing the proportion of correct answers in the learning process of the model at 64 epochs]

 

Рис. 5. Графики ошибки в процессе обучения модели на 64 эпохах

[Figure 5. Error graphs in the learning process of the model at 64 epochs]

 

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

 

Рис. 6. Примеры отображений, принадлежащих определенным классам

[Figure 6. Examples of images belonging to different classes]

 

На рис. 6 представлены ряды по 10 экземпляров отображений, представляющих определенный класс. Можно видеть достаточно заметные отличия паттернов отображений для каждого из классов. Однако сделать убедительные выводы по состоянию поверхностного слоя литосферы на основании только наблюдаемой динамики характеристик сигнала ГАЭ по паттернам отображений и времени проявления паттернов оказалось недостаточно. Здесь в явном виде помимо нелинейности самих сейсмических явлений, формирующих специфические паттерны на отображениях, значительно влияют погодные и сезонные природные влияния в виде осадков, температурных градиентов и перепадов атмосферного давления, которые были ранее зафиксированы в наших исследованиях (например, [24]).

Заключение

Проведенное исследование показало, что разработанная методика позволяет выйти на новый уровень оценки динамики характеристик сигналов ГАЭ, представленных в виде отображений. Анализ отображений с использованием алгоритмов кластеризации и классификации, построенных на базе нейронных сетей, дает возможность уверенно обнаружить общие и отличительные черты и закономерности, скрытые в импульсном потоке сигналов ГАЭ, несмотря на то, что исходные обрабатываемые отображения не имеют строго обозначенных контурных очертаний, а это в значительной степени затрудняет формирование пространства отличительных признаков, на основании которых формируются кластеры. К сожалению малое число исходных отображений не позволяет провести достаточно глубокий анализ, как с позиции применяемых методов предобработки данных, на основе использования нейросетевого подхода, так и с позиции применения различных методов кластеризации. Однако разработанная методика, и разработанные в ходе ее реализации компьютерные программы, позволяют проводить дальнейшие эксперименты в автоматическом режиме на всех этапах от формирования датасета и до вывода результатов распознавания и классификации. Разработанный классификатор на базе архитектуры ResNet успешно справляется с поставленной технической целью и обеспечивает перспективу развития архитектуры в сторону ее адаптации для решения задачи поиска признаков изменчивости состояния приповерхностного слоя литосферы.

Аббревиатуры

ГАЭ	Геоакустическая эмиссия

Об авторах

Юрий Игоревич Сенкевич

Институт космофизических исследований и распространения радиоволн ДВО РАН

Автор, ответственный за переписку.
Email: senkevich@ikir.ru
ORCID iD: 0000-0003-0875-6112

доктор технических наук, доцент, ведущий научный сотрудник, Лаборатории акустических исследований

Россия, 684034, Камчатский край, Елизовский район, с. Паратунка, ул. Мира, 7

Список литературы

Дополнительные файлы