Создание моделей акустической диагностики оборудования

Введение

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

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

Метод основан на смешивании сигналов. Мы разбиваем каждый эталонный сигнал на короткие сегменты, а затем создаём примеры с разной долей дефектной составляющей по формуле смешивания. Из каждого сегмента извлекаются диагностические признаки: статистики временной области, спектральные характеристики, энергии в частотных полосах и признаки огибающей. На этих признаках обучаются три модели случайного леса, каждая предсказывает износ своего элемента. В результате модель сохраняется в файл и готова к применению на новых сигналах, на основании которых будет выполняться диагностика ресурса каждого элемента подшипника.

Присоединение функций и библиотек

Подключим библиотеки, необходимые для работы с данными и обучения модели:

MAT – чтение акустических сигналов из файлов формата MATLAB.
WAV – запись сигналов в аудиофайлы и чтение WAV-файлов для диагностики.
Statistics – базовые статистические функции: среднее, медиана, стандартное отклонение, дисперсия.
StatsBase – расширенные статистики: асимметрия, эксцесс, квантили.
DSP – цифровая обработка сигналов, в частности преобразование Гильберта для вычисления огибающей.
FFTW – быстрое преобразование Фурье для спектрального анализа.
DecisionTree – реализация случайного леса для регрессии.
Random – фиксация генератора случайных чисел для воспроизводимости.
Serialization – сохранение обученной модели в файл JLS.

for pkg in (["MAT", "WAV", "Statistics", "DSP", "FFTW", "DecisionTree", "StatsBase", "Random", "Serialization"])
    try
        eval(Meta.parse("using $pkg"))
    catch
        Pkg.add(pkg)
        eval(Meta.parse("using $pkg"))
    end
end

Загрузим вспомогательные функции из внешних файлов:

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

include("media_player.jl")
include("вычислить_признаки.jl")
include("выборка.jl")

выборка (generic function with 1 method)

Исходные данные

Прочитаем файл Bearing.mat с четырьмя временными последовательностями акустических сигналов подшипника. Результатом станет словарь с ключами "normal", "roller", "inner", "outer", где каждому ключу соответствует матрица из 120000 отсчётов по 12 секунд. Эти данные послужат эталонными сигналами для обучения модели диагностики:

normal — подшипник без дефектов, эталонное состояние;
roller — подшипник с неисправностью тела качения;
inner — подшипник с неисправностью на внутреннем кольце;
outer — подшипник с неисправностью на наружном кольце.

подшипник = matread("Bearing.mat")

Dict{String, Any} with 4 entries:
  "normal" => [0.0886615 0.0586209 … -0.04068 0.0400542]
  "roller" => [-0.0027614 -0.096324 … 0.172831 -0.123776]
  "inner"  => [-0.0830044 -0.195734 … -0.243328 0.111755]
  "outer"  => [0.00852784 0.42355 … 0.41421 0.0966489]

Преобразуем двумерные матрицы из MAT-файла в одномерные векторы длиной 120000 отсчётов для удобства обработки. Функция vec() преобразует матрицу в вектор, что необходимо для поэлементных операций при смешивании сигналов. Каждый вектор соответствует акустической записи подшипника с определённым типом дефекта.

норма = vec(подшипник["normal"]);
ролик = vec(подшипник["roller"]);
внутр = vec(подшипник["inner"]);
внеш = vec(подшипник["outer"]);

Выведем первые и последние значения нормального сигнала для визуального контроля загруженных данных.

display(норма)

120000-element Vector{Float64}:
  0.08866153846153846
  0.05862092307692307
 -0.056951999999999996
 -0.05862092307692307
 -0.04985907692307692
 -0.015646153846153844
  0.092208
  0.06049846153846154
 -0.0002086153846153846
  0.026285538461538457
 -0.04109723076923077
 -0.032544
  0.02858030769230769
  ⋮
  0.09429415384615383
  0.08991323076923076
  0.021070153846153846
 -0.055491692307692304
 -0.016272
  0.05027630769230769
  0.09095630769230768
  0.08156861538461539
 -0.03588184615384615
 -0.12308307692307692
 -0.040679999999999994
  0.040054153846153843

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

Зададим частоту дискретизации сигналов. Нам известно, что записи длятся 12 секунд и содержат по 120 тысяч отсчётов, поэтому частота составляет 10 кГц. Это значение потребуется при вычислении спектральных признаков и частотных характеристик сигналов.

частота_дискретизации = 10000 # Гц

10000

Воспроизведение и визуализация сигналов

Сохраним четыре эталонных сигнала в аудиофайлы формата WAV для возможности прослушивания. Каждый сигнал нормализуется делением на свой максимум. Частота дискретизации 10000 Гц записывается в заголовок WAV-файла.

wavwrite(норма/maximum(abs.(норма)), "Норма.wav", Fs=частота_дискретизации)
wavwrite(ролик/maximum(abs.(ролик)), "Ролик.wav", Fs=частота_дискретизации)
wavwrite(внутр/maximum(abs.(внутр)), "Внутр.wav", Fs=частота_дискретизации)
wavwrite(внеш/maximum(abs.(внеш)), "Внеш.wav", Fs=частота_дискретизации)

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

Воспроизведём аудиосигнал подшипника в нормальном состоянии.

media_player("$(@__DIR__)/Норма.wav", mode="audio")

Воспроизведём аудиосигнал подшипника с полностью неисправным телом качения.

media_player("$(@__DIR__)/Ролик.wav", mode="audio")

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

гр_ролик = plot(1:2000, ролик[1:2000], xlabel="Отсчёты времени", ylabel="Амплитуда", 
title="Нормальное состояние и дефект ролика", label="Дефект тела качения", linecolor=:orangered, linewidth=0.5)
plot!(гр_ролик, 1:2000, норма[1:2000], label="Нормальное состояние", linecolor=:blue, linewidth=0.7)
display(гр_ролик)

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

media_player("$(@__DIR__)/Внутр.wav", mode="audio")

гр_внутр = plot(1:2000, внутр[1:2000], xlabel="Отсчёты времени", ylabel="Амплитуда", 
title="Нормальное состояние и дефект внутреннего кольца", label="Дефект внутреннего кольца", linecolor=:orangered, linewidth=0.5)
plot!(гр_внутр, 1:2000, норма[1:2000], label="Нормальное состояние", linecolor=:blue, linewidth=0.7)
display(гр_внутр)

media_player("$(@__DIR__)/Внеш.wav", mode="audio")

гр_внеш = plot(1:2000, внеш[1:2000], xlabel="Отсчёты времени", ylabel="Амплитуда", 
title="Нормальное состояние и дефект внешнего кольца", label="Дефект внешнего кольца", linecolor=:orangered, linewidth=0.5)
plot!(гр_внеш, 1:2000, норма[1:2000], label="Нормальное состояние", linecolor=:blue, linewidth=0.7)
display(гр_внеш)

Деление на сегменты

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

function разделить_на_сегменты(сигнал, длина_сегмента, перекрытие=0.5)
    шаг = Int(floor(длина_сегмента * (1 - перекрытие)))
    сегменты = []
    for i in 1:шаг:(length(сигнал) - длина_сегмента + 1)
        push!(сегменты, сигнал[i:i+длина_сегмента-1])
    end
    return сегменты
end

разделить_на_сегменты (generic function with 2 methods)

Зададим длину сегмента 4000 отсчётов (0.4 секунды) и перекрытие 75%. Такая длина достаточна для захвата нескольких периодов вращения и стабильного вычисления признаков. Перекрытие 75% позволяет получить 117 сегментов из сигнала длиной 120000 отсчётов вместо 30 без перекрытия.

длина_сегмента = 4000
перекрытие = 0.75

0.75

Разделим все четыре исходных сигнала на сегменты длиной 4000 отсчётов с перекрытием 75%. Каждая переменная содержит массив из 117 векторов, которые станут базовыми блоками для создания обучающей выборки. Это обеспечивает достаточное разнообразие данных для обучения модели.

сегменты_норма = разделить_на_сегменты(норма, длина_сегмента, перекрытие);
сегменты_ролик = разделить_на_сегменты(ролик, длина_сегмента, перекрытие);
сегменты_внутр = разделить_на_сегменты(внутр, длина_сегмента, перекрытие);
сегменты_внеш = разделить_на_сегменты(внеш, длина_сегмента, перекрытие);

Генерация обучающей выборки

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

Random.seed!(42)

TaskLocalRNG()

Сгенерируем полную обучающую выборку, вызвав функцию выборка(). На выходе получим матрицу признаков X размером 31850×53 и матрицу целевых значений y размером 31850×3. Каждая строка содержит 53 признака для одного сегмента, а соответствующая строка — процент износа трёх элементов подшипника.

X_матрица, y_матрица = выборка(сегменты_норма, сегменты_ролик, сегменты_внутр, сегменты_внеш);

Генерация обучающей выборки:
 ✅ Нормальные сигналы
 ✅ Одиночные дефекты
 ✅ Комбинированные дефекты

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

X_матрица

31850×53 Matrix{Float64}:
 0.0118985   0.0123083   0.0740049  …  0.0968883  0.0119781   3.52274
 0.0115287   0.0123083   0.0736846     0.0964796  0.0118108   3.50702
 0.0111265   0.0127255   0.0719866     0.0968276  0.0117078   3.56513
 0.0110415   0.0127255   0.0715529     0.0961194  0.0116846   3.5976
 0.0110879   0.010848    0.0715749     0.0949451  0.0112787   3.68794
 0.0115873   0.0106394   0.0715779  …  0.0973208  0.0119432   3.64734
 0.0115536   0.0106394   0.0723734     0.101552   0.0131384   3.53682
 0.0119747   0.0105351   0.0729215     0.101801   0.0130867   3.55802
 0.0123654   0.0114738   0.072629      0.101026   0.0129805   3.51131
 0.012056    0.0118911   0.0729242     0.0977965  0.0123715   3.55692
 0.0118416   0.0120997   0.0730542  …  0.0974264  0.0121984   3.58336
 0.0109197   0.0123083   0.0728902     0.0958466  0.0117823   3.63756
 0.00994772  0.0114738   0.0725086     0.100375   0.0127277   3.52092
 ⋮                                  ⋱  ⋮                      
 0.0211568   0.0176073   0.456967      0.0700357  0.00561789  4.14688
 0.0273427   0.0251264   0.630943      0.0478307  0.00235748  4.59753
 0.0147288   0.0135144   0.372571   …  0.0418331  0.00184276  4.70477
 0.0148811   0.00737998  0.247604      0.0730399  0.00620312  4.15096
 0.0120027   0.00836822  0.294217      0.0825136  0.00809944  3.77173
 0.0175897   0.0128769   0.369934      0.0660802  0.00521613  4.57255
 0.0135876   0.0131207   0.0825215     0.0877914  0.00952871  3.66911
 0.013041    0.00968947  0.167135   …  0.0708217  0.00573486  4.04914
 0.013665    0.0134183   0.0860728     0.0964847  0.0116604   3.64174
 0.0186058   0.0138228   0.670599      0.0488422  0.00257517  4.48511
 0.0186161   0.0152096   0.331693      0.0622088  0.00439671  4.29467
 0.0294753   0.0305285   0.590931      0.0517937  0.00289528  4.41478

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

y_матрица

31850×3 Matrix{Float64}:
  0.0        0.0       0.0
  0.0        0.0       0.0
  0.0        0.0       0.0
  0.0        0.0       0.0
  0.0        0.0       0.0
  0.0        0.0       0.0
  0.0        0.0       0.0
  0.0        0.0       0.0
  0.0        0.0       0.0
  0.0        0.0       0.0
  0.0        0.0       0.0
  0.0        0.0       0.0
  0.0        0.0       0.0
  ⋮                   
 83.0758    50.6422   81.4412
 35.4018    90.1001   51.4261
  6.37642   52.8442   27.5528
 75.4885    15.796     9.24868
 85.2361    18.1263   72.7981
 91.3287    37.5423   24.939
  0.984328   7.39619   0.0
 51.7633     8.79462  18.2312
  0.0        0.0      58.2913
 63.482     92.1392    0.0
 57.4886    42.6725   41.0776
 63.6515    82.9749    0.0

Нормализация

Выполним нормализацию: вычтем среднее и разделим на стандартное отклонение для каждого из 53 признаков. Это приведёт все признаки к единому масштабу с нулевым средним и единичной дисперсией. Нормализация необходима для корректной работы модели и последующего применения к новым данным.

средние = mean(X_матрица, dims=1)[:]
ско = std(X_матрица, dims=1)[:]
ско[ско .< 1e-10] .= 1.0
X_норм = (X_матрица .- средние') ./ ско'

31850×53 Matrix{Float64}:
 -1.16239   -0.602387   -1.37075    …   1.71603    1.979      -1.72866
 -1.24785   -0.602387   -1.37255        1.69161    1.92531    -1.77095
 -1.34081   -0.519822   -1.3821         1.7124     1.89225    -1.61459
 -1.36046   -0.519822   -1.38453        1.67008    1.88479    -1.52721
 -1.34974   -0.891368   -1.38441        1.5999     1.75453    -1.28409
 -1.2343    -0.932651   -1.38439    …   1.74188    1.96781    -1.39336
 -1.2421    -0.932651   -1.37992        1.99476    2.35134    -1.69077
 -1.14477   -0.953292   -1.37684        2.00962    2.33474    -1.63372
 -1.05447   -0.767519   -1.37848        1.9633     2.30069    -1.7594
 -1.12599   -0.684953   -1.37683        1.77031    2.10525    -1.63669
 -1.17553   -0.64367    -1.37609    …   1.74819    2.04968    -1.56553
 -1.38862   -0.602387   -1.37702        1.65378    1.91616    -1.41967
 -1.61325   -0.767519   -1.37916        1.92439    2.21955    -1.73356
  ⋮                                 ⋱   ⋮                     
  0.977417   0.446225    0.781769       0.111248  -0.0620721  -0.0490778
  2.4071     1.9342      1.75964       -1.21579   -1.10838     1.1636
 -0.508242  -0.363704    0.307404   …  -1.57422   -1.27356     1.4522
 -0.473032  -1.57765    -0.394997       0.290784   0.125736   -0.0381029
 -1.1383    -1.38209    -0.133          0.856963   0.734287   -1.05862
  0.152967  -0.489863    0.292581      -0.125143  -0.191001    1.09639
 -0.771998  -0.441628   -1.32288        1.17238    1.19296    -1.33477
 -0.898318  -1.12063    -0.847294   …   0.158218  -0.0245357  -0.312119
 -0.754119  -0.382733   -1.30292        1.69191    1.87702    -1.40841
  0.387814  -0.302685    1.98253       -1.15534   -1.03852     0.861098
  0.390198  -0.0282467   0.0776393     -0.356513  -0.453962    0.348613
  2.89999    3.00321     1.53474       -0.978947  -0.935788    0.671842

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

function нормализовать(пр)
    return (пр .- средние) ./ ско
end

нормализовать (generic function with 1 method)

Обучение моделей

Создадим три независимые модели случайного леса (Random Forest) для предсказания износа каждого элемента подшипника. Модели для колец используют 400 деревьев с глубиной 35, а для ролика 600 деревьев с глубиной 50, поскольку его дефект сложнее детектировать из-за меньшей амплитуды. Каждая модель будет обучаться на своём целевом столбце матрицы y.

модель_внутр = DecisionTree.RandomForestRegressor(
    n_trees=400, max_depth=35, min_samples_split=2, min_samples_leaf=1, rng=42)
модель_внеш = DecisionTree.RandomForestRegressor(
    n_trees=400, max_depth=35, min_samples_split=2, min_samples_leaf=1, rng=43)
модель_ролик = DecisionTree.RandomForestRegressor(
    n_trees=600, max_depth=50, min_samples_split=2, min_samples_leaf=1, rng=44)

RandomForestRegressor
n_trees:             600
n_subfeatures:       -1
partial_sampling:    0.7
max_depth:           50
min_samples_leaf:    1
min_samples_split:   2
min_purity_increase: 0.0
ensemble:            nothing

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

DecisionTree.fit!(модель_внутр, X_норм, y_матрица[:, 1])
DecisionTree.fit!(модель_внеш, X_норм, y_матрица[:, 2])
DecisionTree.fit!(модель_ролик, X_норм, y_матрица[:, 3])

RandomForestRegressor
n_trees:             600
n_subfeatures:       -1
partial_sampling:    0.7
max_depth:           50
min_samples_leaf:    1
min_samples_split:   2
min_purity_increase: 0.0
ensemble:            Ensemble of Decision Trees
Trees:      600
Avg Leaves: 14299.196666666667
Avg Depth:  36.45166666666667

Сохранение моделей

Вычислим эталонные среднее квадратичное значение (root mean square – RMS) и размах нормального сигнала для будущей нормализации амплитуды новых записей. При диагностике сигнал будет масштабироваться так, чтобы его RMS совпадал с эталонным значением 0.073. Это позволит корректно обрабатывать сигналы, записанные с разным усилением.

эталонный_rms = sqrt(mean((норма .- mean(норма)) .^ 2));
эталонный_размах = maximum(abs.(норма .- mean(норма)));

Сохраним обученные модели и все сопутствующие параметры в файл модель_подшипника.jls через стандартный механизм сериализации Julia. В файл записываются три модели Random Forest, векторы средних и СКО для нормализации признаков, параметры сегментации сигналов и эталонные значения амплитуды. Время выполнения данной кодовой ячейки может составить около 15 минут. При последующей диагностике новых сигналов этот файл будет загружаться одной командой deserialize, избавляя от необходимости повторного обучения и сериализации.

serialize("модель_подшипника.jls", (
    модель_внутр = модель_внутр,
    модель_внеш = модель_внеш,
    модель_ролик = модель_ролик,
    средние_признаков = средние,
    ско_признаков = ско,
    длина_сегмента = длина_сегмента,
    перекрытие = перекрытие,
    частота_дискретизации = частота_дискретизации,
    эталонный_rms = эталонный_rms,
    эталонный_размах = эталонный_размах))

Заключение

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

Ключевые инженерные решения, принятые в работе, стоит выделить. Разбиение сигналов на короткие перекрывающиеся сегменты дало 468 базовых блоков вместо четырёх целых сигналов, что критически важно для обучения. Формула смешивания норма + степень * (дефект - норма) позволяет создавать искусственные примеры с любой степенью износа, не имея реальных записей промежуточных состояний. Модель для тела качения использует 600 деревьев вместо 400, потому что его дефект слабее по амплитуде и труднее детектируется.

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

Модели акустической диагностики