Все хотя бы раз слышали про землетрясения. Это опасное природное явление которое может привести к разрушению зданий, возникновению цунами и гибели людей. С помощью составления каталогов землетрясений и анализа сейсмической активности человечество пытается минимизировать риски от наступления сейсмических событий. Основной источник данных в сейсмологии - это непрерывные записи движения грунта, которые регистрируют с помощью сейсмических станций. Для того что бы составить представительный каталог землетрясений необходимо определить времена вступлений сейсмических волн в непрерывном потоке, рассчитать параметры очага и выполнить оценку магнитуды. Каждый этап рутинной обработки сейсмологических данных это предмет отдельной статьи, но мы с вами посмотрим на самое интересное - распознавание сейсмических волн методами машинного обучения. Определение времён вступлений сейсмических волн до сих пор выполняется (или корректируется) операторами-сейсмологами. Автоматизация этой сложной задачи позволит полностью исключить ручной труд при непрерывной обработке данных любой сейсмической сети. При этом, необходима такая модель, которая с одной стороны могла обеспечить точность отметок вступлений на уровне человека, а с другой была эффективна на этапе вывода (эксплуатация на CPU). Возможно ли это? Давайте посмотрим!

Как выглядят исходные данные?

Сейсмическая станция - это прибор который регистрирует движения грунта (скорость, ускорение, смещение). Данные представляют собой временной ряд время -> амплитуда. Типичная частота дискретизация в сейсмологии это 100Hz. Каждая сейсмическая станция обычно укомплектована тремя цифровыми каналами: два горизонтальных и один вертикальный.

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

Локальная сеть сейсмических станций, в зависимости от площади мониторинга, как правило, насчитывает 10 и более приборов. Все станции пишут непрерывные суточные архивы. Эти данные не идеальны. В непрерывном потоке могут быть пропуски из-за продолжительных обрывов связи. Датчики и регистраторы иногда выходят из строя и пишут совсем не то, что нужно, и так далее.

Итак, представим, что данные у нас есть. Но прежде чем делать что-то новое (да ещё и c ML) давайте посмотрим какие проблемы возникают у сейсмологов на этапе выделения вступлений сейсмических волн.

Проблема №1: рост числа сейсмических станций

Когда-то раньше сейсмические станции были дорогими и сложными в обслуживании. Сейчас приборы на базе raspberry может позволить себе не только любая геофизическая компания, но и самостоятельные исследователи или просто увлечённые этой темой люди. Любому инженеру понятно что чем больше приборов, тем больше данных, и, значит, точность расчётов будет лучше. Хорошо! А теперь давайте посмотрим, как выглядит 15 минут записи 28 станций.

Понятно, что какие то сильные землетрясения на такой развёртке будет видно, но слабые события тут уже разглядеть не получится. Нужно учитывать что оператор, который просматривает поток глазами, должен перекрывающими окнами по 15 минут просмотреть целые сутки. При 28 активных станциях для мониторинга Сахалина уже нужно делить сейсмометры на две группы (Север и Юг острова) и просматривать каждую группу отдельно. При этом, очевидно, что 28 станций это мало и через 3-5 лет их будет в два раза больше.

Проблема №2: пропуск слабых землетрясений

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

Мы видим, что масштаб амплитуды волн от землетрясения такой, что весь сейсмический шум перед вступлением "сжался" практически в одну линию. Нужно сказать, что это запись акселерометра, который предназначен для регистрации сильных движений. Если бы на его месте был бы велосиметр, то он, скорее всего, зашкалил и обрезал сигнал в области максимальных амплитуд.

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

Могу сказать, что только человек с большим опытом рутинной обработки может такое. Тут нужно применять техники фильтрации, иметь очень натренированный глаз и так далее. Человеческий фактор здесь решает всё - что-то оператор видит, а что-то пропускает.

Хорошо, а может этих слабых землетрясений мало? Давайте хотя бы примерно оценим где и сколько мы пропускаем. И тут на сцену выходит один из фундаментальных законов сейсмологии - закон Гутенберга-Рихтера. И нам этот закон говорит, что землетрясения распределены лог-нормально, и, значит, количество землетрясений растёт экспоненциально в область слабых магнитуд.

То есть, грубо говоря, кумулятивный график должен быть похож на прямую линию, но мы видим загиб ниже M=2.0, то есть здесь начинаются существенные пропуски событий. Конечно, в целом по Сахалину не каждый район достаточно плотно покрыт станциями чтобы уверенно регистрировать слабые события, но, очевидно, что из текущих данных мы можем выжать существенно больше вступлений за счёт поиска пропущенных оператором землетрясений.

Проблема №3: много импульсного шума

Фоновый шум мы уже разобрали, но дело в том, что это не весь шум который бывает. Если станция расположена вблизи населённого пункта, то мы видим на записях все проявления индустриальной активности: ходьба людей рядом со станцией, работа снегоуборочной техники, движение автотранспорта поблизости и так далее. Если у вас таких "шумных" станций много, то выглядит это примерно так:

На примере выше алгоритм STA/LTA (о нем чуть позже) ошибочно расставил отметки вступления. Эти отметки были поданы дальше в алгоритм локации, который по пороговому значению среднеквадратичной невязки посчитал что это реальное землетрясение. Но на самом деле нет. Это просто грузовики так одновременно проехали или что-то там ещё в близи станций случилось. Такой импульсный шум сильно усложняет обработку. Бывают всплески очень похожие на вступления сейсмической волны и оператор всегда тратит дополнительное время при появлении импульсного шума.

ОК! Проблемы мы как то обозначили, а теперь давайте посмотрим как выполняется автоматизация отметок фаз сейсмических волн не ML методами. Всё таки сейсмология не вчера появилась...

«Старый» STA/LTA

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

• Рассматривается только "энергетический" показатель превышения над уровнем шума, а не сама форма сигнала. Это значит, если вы прыгнули возле сейсмостанции или рядом проехала машина, то триггер сработает (см. Проблема №3).

• Сигнал должен иметь существенное превышение над фоновым значением. Это настраиваемый параметр, он задаётся отношением сигнал/шум и варьируется от 1,5 до 7,0. Что бы срабатывания не были слишком частыми в сейсмологической практике принято устанавливать значение 3,0 и выше. Таким образом у нас сразу выпадают слабые землетрясения сигнал от которых сопоставим с уровнем сейсмического шума (см. Проблема №2).

• Есть ещё такое явление как парные землетрясения, они происходят одно за другим с разницей в несколько секунд. На сейсмограмме это выглядит так, что сигнал от первых вступлений ещё не затух, как тут же приходят новые волны. В этом случае STA/LTA вторые вступления то же пропустит, так как для того что бы отметить второе вступление фоновый уровень должен придти в "норму", то есть опуститься до значений, которые были до первого землетрясения.

Кросс-корреляция и сопоставление шаблонов

Этот метод начал активно применяться с начала 2000-х когда вычислительные мощности позволили выполнить это. Алгоритм способен находить очень слабые землетрясения. Кросс-корреляция активно используется для восстановления афтершоковых последовательностей. Дело в том, что когда происходит сильное землетрясение то следом за главным событием обычно происходят тысячи более слабых землетрясений. Если в качестве мастер-события задать главный толчок, то тогда поиск других землетрясений из этой же очаговой зоны будет успешен. Однако, у этого метода есть главный недостаток: нет мастер-события - нет результата. Если землетрясение произойдет в другой очаговой зоне, то мы ничего не увидим. То есть применять этот алгоритм для решения наших проблем затруднительно, так как мы не знаем заранее где возникнет сейсмическое событие.

Deep learning

Итак, мы с вами добрались до самого интересного - это методы машинного обучения. Кроме очевидных преимуществ этого подхода, таких как распознавание именно формы сигнала и скорость вывода, стоит отметить открытые вопросы:

• Проблема обобщающей способности (Generalization). Если мы обучим нашу модель, например, на данных сейсмической сети Южной Калифорнии, будет ли эта модель успешно классифицировать вступления от сахалинских или итальянских землетрясений?

• Какую производительность будет показывать модель на других типах оборудования, например, на мобильных устройствах? Да, у нас с вами в кармане маленькая сейсмостанция, потому в современных телефонах есть акселерометр, который в основном используется для игр. И если модель будет работать для мобил, то мы можем использовать данные устройств, которые оказались вблизи сильного события, а это на несколько порядков больше стационарных станций.

Давайте заглянем в научные работы и посмотрим как там обстоят дела. Если исключить совсем пионерские работы, то основная масса научных статей на эту тему опубликована в последние пять лет. Нужно дополнительно очертить наше поле исследований: мы будем говорить о локальных землетрясениях, то есть о таких событиях, когда расстояние от очага до сейсмической станции не превышает 300 километров. Выделять мы будем P и S волны от локальных землетрясений. Сейчас можно уверенно говорить о существовании двух подходов к этой задаче.

Простой классификатор коротких записей

Самое простое и интуитивное - это нарезать волновой ряд на короткие отрезки по 4-6 секунд. Каждый отдельный отрезок рассматриваем как образец (sample) который принадлежат к одному из классов: P-волна, S-волна или шум (Noise). Так как сейсмические волны имеют определённую характерную форму, по которой их узнаёт человек, то мы можем классифицировать наши данные как картинки, использую всю мощь соответствующих архитектур (например CNN). На этапе вывода волновой поток от сейсмостанции "нарезается" перекрывающимися окнами (что бы сигнал полностью попал в какое-то окно) и там, где вероятность выше заданного порога ставится отметка вступления. Наиболее удачная модель в этой области - это GPD.

Выбор вступлений на записи всего землетрясения

Другой подход - это подавать в модель запись всего землетрясения (30-60 секунд) с применением более сложной системы триггеров. Наиболее значимый результат показывает модель EQT. Несмотря на то что в названии модели присутствует слово Transformer, эта архитектура больше относится к CNN-based за счет подавляющего числа соответствующих блоков. Данный метод показывает хорошие результаты при сканировании непрерывных архивов, но он не подходит для обработки данных в режиме реального времени, потому что нам необходимо подождать, пока вся запись землетрясения поступит в дата центр.

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

Данные, обучающая и тестовая выборки

Для того что бы правильно проверить обобщающую способность надо учить на одном наборе, а проверять на другом. Так и поступим. В качестве обучающей выборки возьмем самый большой на сегодняшний день датасэт волновых форм сейсмической сети Южной Калифорнии (размечен вручную). Архив содержит 4,5 миллиона сейсмограмм ровно разделённых по классам P/S-волн и Шум (P/S/N). Эти данные нормализованы и отфильтрованы, убраны отметки времени и станций, каждый сэмпл содержит три компоненты по 400 точек (4 секунды по 100Hz).

В качестве тестовых данных мы будем использовать записи Сахалинской и Дегстанской сейсмических сетей. То есть обучать мы будем на Калифорнии а проверять на регионах в другом конце Земли - всё по честному. Тестовых данных будет не так много, потому что сети сейсмических станций не такие плотные как в Калифорнии, но для проверки вполне достаточно. Приводим данные к унифицированному виду как в обучающей выборке и готово. Хотя, на практике это было не так быстро, и отняло у нас кучу времени.

Spectrogram-based models

Решая проблему обобщающей способности мы интуитивно понимаем, что форма сигнала от локальных землетрясений сохраняется похожей что для Калифорнии, что для Дагестана (или любого другого региона). Консервативные геофизики скажут, что сигнал искажается грунтовыми условиями под конкретной станцией. Но ведь если мы отправим сейсмолога из Махачкалы в Лос-Анджелес и дадим ему в руки волновой поток, то он с уверенностью отметит все вступления. Нам нужно каким то образом наиболее четко выделить форму сигнала независимо от производителя прибора и места его локации. И тут нам на помощь приходит STFT. У нас под рукой временной ряд, так почему же не сделать спектрограмму и распознавать вступления уже с неё. К тому же эта техника активно применяется в процессинге аудио.

Отлично! Теперь давайте сделаем модель CNN со спектрограммы. Если вы используете tensorflow то можете использовать готовое решение и положить вычисление STFT прямо в модель. Спектрограмма для низкочастотного "звука" сейсмических волн имеет свои характеристики, но в целом это не сильно отличается от примеров в документации.

Model: "model_3"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input_4 (InputLayer)         [(None, 400, 3)]          0         
_________________________________________________________________
stft_2 (STFT)                (None, 22, 33, 3)         0         
_________________________________________________________________
magnitude_2 (Magnitude)      (None, 22, 33, 3)         0         
_________________________________________________________________
magnitude_to_decibel_2 (Magn (None, 22, 33, 3)         0         
_________________________________________________________________
max_abs_scaler_2 (MaxABSScal (None, 22, 33, 3)         0         
_________________________________________________________________
conv2d (Conv2D)              (None, 22, 33, 32)        896       
_________________________________________________________________
batch_normalization_6 (Batch (None, 22, 33, 32)        128       
_________________________________________________________________
tf.nn.relu (TFOpLambda)      (None, 22, 33, 32)        0         
_________________________________________________________________
max_pooling2d (MaxPooling2D) (None, 11, 16, 32)        0         
_________________________________________________________________
conv2d_1 (Conv2D)            (None, 11, 16, 64)        18496     
_________________________________________________________________
batch_normalization_7 (Batch (None, 11, 16, 64)        256       
_________________________________________________________________
tf.nn.relu_1 (TFOpLambda)    (None, 11, 16, 64)        0         
_________________________________________________________________
max_pooling2d_1 (MaxPooling2 (None, 5, 8, 64)          0         
_________________________________________________________________
conv2d_2 (Conv2D)            (None, 5, 8, 128)         73856     
_________________________________________________________________
batch_normalization_8 (Batch (None, 5, 8, 128)         512       
_________________________________________________________________
tf.nn.relu_2 (TFOpLambda)    (None, 5, 8, 128)         0         
_________________________________________________________________
max_pooling2d_2 (MaxPooling2 (None, 2, 4, 128)         0         
_________________________________________________________________
flatten_1 (Flatten)          (None, 1024)              0         
_________________________________________________________________
dropout_14 (Dropout)         (None, 1024)              0         
_________________________________________________________________
dense_19 (Dense)             (None, 80)                82000     
_________________________________________________________________
dropout_15 (Dropout)         (None, 80)                0         
_________________________________________________________________
dense_20 (Dense)             (None, 3)                 243       
=================================================================
Total params: 176,387
Trainable params: 175,939
Non-trainable params: 448
_________________________________________________________________

Здесь, наверное, нет ничего супер инновационного. Производительность CNN мы примерно знаем (и чуть позже протестируем). Но у нас в активе есть ещё трансформеры. Все про них уже знают, это сегодняшний мейнстрим, во многом за счёт высокой обобщающей способности и легкой адаптации модели к различным приложениям. Однако, эта технология глубокого обучения имеет существенный недостаток – вычислительная сложность O(n^2) в механизме внимания. В 2019 году была представлена архитектура Linformer, способная выполнять операции с вниманием за линейное время. Этот метод обладает некоторыми недостатками, одним из которых является существенное снижение точности. Работа «Rethinking Attention with Performers» (статья, код-эталон, видео-разбор) это попытка подобрать оценку к исходной матрице внимания, которую можно вычислить за линейное время. Механизм называется FAVOR+, он полностью совместим с оригинальной архитектурой.

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

Model: "model_2"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input_3 (InputLayer)         [(None, 400, 3)]          0         
_________________________________________________________________
stft_1 (STFT)                (None, 22, 33, 3)         0         
_________________________________________________________________
magnitude_1 (Magnitude)      (None, 22, 33, 3)         0         
_________________________________________________________________
magnitude_to_decibel_1 (Magn (None, 22, 33, 3)         0         
_________________________________________________________________
max_abs_scaler_1 (MaxABSScal (None, 22, 33, 3)         0         
_________________________________________________________________
rearrange3d_1 (Rearrange3d)  (None, 11, 198)           0         
_________________________________________________________________
dense_11 (Dense)             (None, 11, 48)            9552      
_________________________________________________________________
cls_token_1 (ClsToken)       (None, 12, 48)            48        
_________________________________________________________________
pos_embeding2_1 (PosEmbeding (None, 12, 48)            576       
_________________________________________________________________
performer_block_2 (Performer (None, 12, 48)            18672     
_________________________________________________________________
performer_block_3 (Performer (None, 12, 48)            18672     
_________________________________________________________________
lambda_1 (Lambda)            (None, 48)                0         
_________________________________________________________________
layer_normalization_9 (Layer (None, 48)                96        
_________________________________________________________________
dropout_11 (Dropout)         (None, 48)                0         
_________________________________________________________________
dense_16 (Dense)             (None, 96)                4704      
_________________________________________________________________
dropout_12 (Dropout)         (None, 96)                0         
_________________________________________________________________
dense_17 (Dense)             (None, 48)                4656      
_________________________________________________________________
dropout_13 (Dropout)         (None, 48)                0         
_________________________________________________________________
dense_18 (Dense)             (None, 3)                 147       
=================================================================
Total params: 57,123
Trainable params: 57,123
Non-trainable params: 0
_________________________________________________________________

Тестирование производительности моделей

Теперь обучим наши модели. Для этого стандартно поделим обучающий датасэт на тренировочный и валидационный (80/20). Поставим раннюю остановку там, где сходимость по валидационным данным (validation loss) не улучшается в течении 5 эпох. Учить будем несколько раз что бы исключить случайно хороший результат. Поэтому после каждой итерации фиксируем результат для тестовых данных и всё заново.

%%time

interations = 5

callbacks = [
    keras.callbacks.EarlyStopping(monitor = 'val_loss', patience = 5, restore_best_weights = True, verbose = 1),
    keras.callbacks.ModelCheckpoint('w_performer_epoch_hpa{epoch}.h5', save_weights_only=True),
]

sparse_categorical_accuracy = []
loss = []
val_loss = []
val_sparse_categorical_accuracy = []

cal_top_1 = []
sakh_top_1 = []
dag_top_1 = []

for i in range(interations):
  print("Load dataset with random state {}".format(i)) 
  X_train, X_test = h5_tts('/content/scsn_ps_2000_2017_shuf.hdf5',\
                             batch_size = 480,\
                             test_size = 0.2, random_state = i, shuffle = False)

  print("compile the model")
  model_with_spec_cal = seismo_performer_with_spec(
                                        maxlen=400,
                                        nfft=64,
                                        hop_length=16,
                                        patch_size_1=22,
                                        patch_size_2=3,
                                        num_channels=3,
                                        num_patches=11,
                                        d_model=48,
                                        num_heads=2,
                                        ff_dim_factor=2,
                                        layers_depth=2,
                                        num_classes=3,
                                        drop_out_rate=0.1)



  LR = 0.0001
  model_with_spec_cal.compile(
      optimizer=keras.optimizers.Adam(learning_rate=LR),
      loss=keras.losses.SparseCategoricalCrossentropy(),
      metrics=[keras.metrics.SparseCategoricalAccuracy()],)

  print("Fit model on training data step {}".format(i+1))
  history = model_with_spec_cal.fit(
    X_train,
    #shuffle=False,
    epochs=400,
    validation_data=(X_test),
    callbacks = callbacks
    )
  
  #save convergance curve 
  sparse_categorical_accuracy.append(history.history['sparse_categorical_accuracy'])
  loss.append(history.history['loss'])
  val_sparse_categorical_accuracy.append(history.history['val_sparse_categorical_accuracy'])
  val_loss.append(history.history['val_loss'])

  print("Evaluate transformer based model on SAKHALIN test data")
  sakh_top_1.append(model_with_spec_cal.evaluate(X_test_sakh, y_test_sakh, batch_size=128)[1])

  print("Evaluate transformer based model on DAGESTAN test data")
  dag_top_1.append(model_with_spec_cal.evaluate(X_test_dag, y_test_dag, batch_size=128)[1])
  
  print("Evaluate transformer based model on CALIFORNIA test data")
  cal_top_1.append(model_with_spec_cal.evaluate(X_test)[1])

  model_with_spec_cal.save_weights('/content/drive/MyDrive/MODELS/w_model_performer_with_spec_seeds.hd5', save_format='h5')

Так как у нас три класса, то проверяем точность соответствующим образом. На этот счёт есть классная картинка от наших коллег из Калтэка.

Что бы протестировать производительность на этапе эксплуатации запустим наши модели на непрерывных данных. Поделим суточный архив одной сейсмической станции на перекрывающиеся окошки со сдвигом в 40 миллисекунд, получим 214К сэмплов, прогоним их на CPU и замерим время вывода. Для того что бы исключить накладные расходы на обращение к файловой системе (и много к чему ещё) измерять будем только время вывода самих моделей. А результаты сравним с топовой моделью GPD.

В таблице представлена общая точность для тестируемых моделей. Seismo-Performer и Spec-CNN это наши новые модели. GDP - это оригинальная модель, а GPD-fixed это переобученная и GPD на tf 2.5.0. При обучении на Калифорнийских данных (Panel A) мы видим лучший результат по точности как для тестовых так и для валидационных данных у архитектур на базе спектрограмм. Интересно, что если обучить все модели на специфичном Сахалинском наборе (Panel B), который содержит на несколько порядков меньше примеров, то наши новые модели существенно лучше распознают данные, которые они раньше не видели. В целом, можно сделать вывод, что и Seismo-Performer и Spec-CNN обладают лучшей обобщающей способностью, им нужно меньше данных для сходимости, так что наш скромный трюк со спектрами сработал.

По производительность на этапе вывода лидирует Seismo-Performer: в 2,5 раза быстрее оригинально GPD и примерно на 35% быстрее своего соседа Spec-CNN. Здесь у нас чистая архитектура трансформера с оптимизацией FAVОR+ и минимальное количество параметров - всего 57k. По точности в абсолютных цифрах показатели немного ниже Spec-CNN, но +35% к скорости это существенно, с учётом того, что при эксплуатации это будет не одна станция а, например, 20.

Если мы более подробно посмотрим на метрики производительности классификации, то увидим, что наши модели на основе спектров существенно лучше по показателю recall, а это значит, что они меньше "пропускают" полезный сигнал. Это очень важно на этапе эксплуатации, потому что ложные срабатывания всё же есть. Ложные срабатывания на реальных непрерывынх данных - это проблема всех существующих моделей. Дело в том при подготовке набора к обучению невозможно найти абсолютно все граничные случаи, когда сигнал является шумом, но он похож на реальное вступление. В качестве временного решения мы можем повысить порог значения SoftMax для класса при котором он будет считаться точным срабатыванием. Если у вас хороший recall, то вы можете сделать порог повыше и получить меньше пропусков и почти полное отсутсвие ложных срабатываний.

Это, можно сказать, конечный результат. Красные линии это отметки оператора. Звёздочки, это то что нашли модели (практически в одном месте для всех примеров). В качестве времени вступление выбирается "окошко" в котором вероятность на класс максимальна относительно соседних окон и выше заданного порога. Точность отметки в районе погрешности отметки оператора.

Заключение

Мы с вами посмотрели специфичную задачу из области наук о Земле. На выходе мы получили ML-модели которые отмечают вступления сейсмических волн на уровне оператора-сейсмолога. За счёт применения спектрограмм у нас получилось сделать сигнал более четким и улучшить обобщающую способность. А современные наработки из NLP позволили сделать очень быструю модель, которая может работать на CPU и готова к внедрению в сейсмологических дата центрах.

Оригинальная работа, репозиторий на GitHub c данными и плейбуками.

Спасибо, что дочитали до конца.