Третьоктавный анализ для железнодорожного транспорта
Третьоктавный анализ — это метод частотного анализа сигналов, при котором спектр делится на полосы, каждая из которых составляет 1/3 октавы. Применение третьоктавного анализа для поездов позволяет исследовать шум и вибрации, оценивать их воздействие на окружающую среду, пассажиров и инфраструктуру. Этот метод особенно важен для:
-
Нормирования шума (соответствие ГОСТ, ISO, EN).
-
Диагностики состояния железнодорожного полотна и подвижного состава.
-
Оптимизации конструкций (например, снижение шума от колесных пар или рельсов).
Часто для анализа применяются акселерометры, установленные на колесных парах, раме вагона, рельсах, грунте (для оценки распространения вибраций в зданиях). Для железнодорожного транспорта ключевые третьоктавные полосы обычно лежат в диапазоне 1–1000 Гц (вибрация) и 20–5000 Гц (шум). Примеры важных частот:
-
Низкие (1–80 Гц): Вибрации от ударов колес о стыки рельсов.
-
Средние (80–500 Гц): Шум качения, гул двигателей.
-
Высокие (500–5000 Гц): Скрип тормозов, аэродинамический шум.
Построение третьоктавного спектра позволяет проанализировать пики в определенных полосах, которые могут указывать на:
-
Износ колес (рост уровня в полосах 63–250 Гц).
-
Дефекты рельсов (резонансы на 31.5–125 Гц).
-
Проблемы с подвеской (вибрации 2–20 Гц).
Рассмотрим задачу анализа реальных данных с акселерометра на ЖД путях и построения третьоктавного спектра.
Импорт и предобработка данных
Подключение необходимых библиотек:
using DSP, DelimitedFiles, Statistics, FFTW, Polynomials
Данные с акселерометра хранятся в файле data.txt. Считаем их в виде матрицы в переменную datamat:
datamat = readdlm("data.txt", ',');
Выделим из неё отдельные вектора для времени и "сырых" значений виброускорения (в м/c^2), отобразим на графике во временной области:
timevec = datamat[:,1];
data1 = datamat[:,2];
Plots.plot(timevec,data1,
title = "Импортированные данные виброускорений",
xlabel = "Время, с",
ylabel = "Виброускорение, м/с^2",
legend = false)
Мы видим, что регистрация значений виброускорений началась не сразу в начальный момент времени, поэтому логично отрезать часть сигнала близкого к нулю уровня. Для этого мы оставим только индексы вектора, соответствующие значениям сигнала, превышающим некий порог (в нашем случае - 1):
indices = findall(x -> x > 1, data1);
accel = data1[indices[1]:indices[end]];
t = timevec[1:length(accel)];
Plots.plot(t, accel,
title = "Вырезанные значения виброускорений",
xlabel = "Время, с",
ylabel = "Виброускорение, м/с^2",
legend = false)
Переход к виброскоростям
Для визуализации виброскоростей нам необходимо проинтегрировать вектор значений акселераций accel. Напишем пользовательскую функцию куммулятивного трапециевидного численного интегрирования cumtrapz, аналог той, что используется в среде MATLAB:
function cumtrapz(X::T, Y::T) where {T <: AbstractVector}
@assert length(X) == length(Y)
out = similar(X)
out[1] = 0
for i in 2:length(X)
out[i] = out[i-1] + 0.5*(X[i] - X[i-1])*(Y[i] + Y[i-1])
end
out
end
Применим пользовательскую функцию к вектору accel, а затем выделим постоянную составляющую (или тренд) в сигнале виброскоростей, чтобы при необходимости её удалить из сигнала поэлементным вычитанием. Для этого приблизим отсчёты сигнала виброскоростей полиномом 12-го порядка и оценим его значения в те же моменты времени, что и у сигнала скоростей:
velocity_trend = cumtrapz(t,accel);
poly = fit(t, velocity_trend, 12);
trend = poly.(t);
Plots.plot(t, velocity_trend, label = "Виброскорость",
title = "Виброскорость и её тренд",
xlabel = "Время, с")
Plots.plot!(t, trend, linewidth = 3, label = "Постоянная составляющая",
ylabel = "Виброскорость, м/с")
Спектральный анализ методом Уэлча
Проведём спектральный анализ виброускорений стандартными пакетами библиотеки DSP.jl. Для начала определим значение частоты дискретизации fs исходного сигнала в Гц, проанализировав шаг вектора времени dt:
dt = mean(diff(timevec));
fs = Int(floor(1/dt))
Зададим параметры для анализа спектра, которые содержались в техническом задании. В частности: размер "окна" для входа алгоритма БПФ, количество нулей для интерполяции графика спектральной плотности мощности (СПМ), перекрытие в 50% для метода Уэлча, выходное количество отсчётов БПФ (размер "окна" + кол-во нулей).
window_length = fs;
nzeros = 100000;
overlap = 0.5;
nFFT = window_length + nzeros;
noverlap = Int(floor(overlap*window_length));
stepsize = window_length - noverlap;
Стоит отметить, что количество точек входного "окна" равно частоте дискретизации. Это позволяет провести оценку спектра с разрешением в 1 Гц.
Также по заданию анализ проводится только в нижней половине спектра (первой зоны Найквиста) положительных частот. Для этого сигнал пропускается через фильтр нижних частот (ФНЧ) Баттерворта 12-го порядка, а его нормированная частота среза составляет 0.5. Результат фильтрации записывается в переменную filtered:
Wp = 0.5;
myfilter = digitalfilter(Lowpass(Wp), Butterworth(12));
filtered = filtfilt(myfilter,accel);
Теперь воспользуемся функцией welch_pgram для оценки спектральной плотности мощности методом Уэлча. Параметры перекрытия, размера входного "окна" и выходного количества точек БПФ соответствуют установленным ранее. Отобразим графики СПМ в дБ:
pxx = welch_pgram(filtered, n = stepsize, noverlap = div(stepsize, 2);
onesided = true,
nfft = nFFT,
fs = fs,
window = hanning);
pxx_nofilt = welch_pgram(accel, n = stepsize, noverlap = div(stepsize, 2);
onesided = true,
nfft = nFFT,
fs = fs,
window = hanning);
Plots.plot(freq(pxx), pow2db.(power(pxx)), label = "После фильтра",
title = "СПМ сигнала виброускорений",
xlabel = "Частота, Гц")
Plots.plot!(freq(pxx_nofilt), pow2db.(power(pxx_nofilt)), label = "До фильтра",
ylabel = "Спектральная плотность мощности, дБ/Гц")
Также отобразим амплитудный спектр в единицах измерения сигнала. Для этого нам необходимо определить параметр эквивалентной шумовой полосы для оконной функции Хэмминга, которую мы применяли в методе Уэлча:
mywin = hanning(window_length);
ENBW = fs * (sum(abs.(mywin).^2)) / (abs.(sum(mywin)).^2)
Значение амплитуды мы получаем из значения мощностей спектра с учётом параметра ENBW:
amplitude = sqrt.(2*power(pxx)*ENBW);
Plots.plot(freq(pxx), amplitude, legend = false,
title = "Амплитудный спектр сигнала виброускорений",
ylabel = "Амплитуда, м/с^2",
xlabel = "Частота, Гц",
xlim = (0,1000))
Третьоктавный спектр
Для построения третьоктавного спектра необходимо сперва определить границы частотных диапазонов и их количество count:
nstart = -20;
ncount = nstart;
count = 0;
fmax = 0;
while fmax*(2^(1/6)) <= fs/2
fmax = (1000).*((2^(1/3)).^ncount);
ncount = ncount + 1;
count = count + 1;
end
Подсчитаем среднюю мощность в каждом из диапазонов (вектора ff) в цикле и запишем результаты в вектор octave_power:
ff = (1000).*((2^(1/3)).^(nstart:(ncount-1)));
octave_power = zeros(count);
psdx = power(pxx);
for i = 1:count
f_low = ff[i]/(2^(1/6));
f_high = ff[i]*(2^(1/6));
idx = .&(freq(pxx) .> f_low, freq(pxx) .<= f_high);
idx = vec(idx);
octave_power[i] = mean(psdx[idx])
end
Перейдём к значениям средней амплитуды в полосах и отобразим третьоктавный спектр функцией bar:
octave_amplitude = sqrt.(2*octave_power*ENBW);
freq_ticks = round.(ff, digits = 0);
BW = exp.(0.2:0.24:0.24*count);
bar(freq_ticks, octave_amplitude, bar_width = BW, xaxis=:log,
xticks = (ff, freq_ticks), color = :green,
title = "Третьоктавный спектр виброускорений",
xlabel = "Частота, Гц",
ylabel = "Амплитуда, м/с^2",
legend = false)
Заключение
Третьоктавнй анализ - распространённый способ диагностики полотна или колёсных пар подвижного состава, позволяющий обнаружить, к примеру, неравномерный износ рельсов или проблемы с контактом "колесо–рельс". Анализ спектра в подобном диапазоне прописан в ГОСТ 31319-2009 и ГОСТ 31296.1-2019.
Мы рассмотрели возможности доступных инструментов среды Engee для проведения подобного анализа.