Третьоктавный анализ для железнодорожного транспорта

Третьоктавный анализ — это метод частотного анализа сигналов, при котором спектр делится на полосы, каждая из которых составляет 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

cumtrapz (generic function with 1 method)

Применим пользовательскую функцию к вектору 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))

3472

Зададим параметры для анализа спектра, которые содержались в техническом задании. В частности: размер "окна" для входа алгоритма БПФ, количество нулей для интерполяции графика спектральной плотности мощности (СПМ), перекрытие в 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)

1.5004321521175445

Значение амплитуды мы получаем из значения мощностей спектра с учётом параметра 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 для проведения подобного анализа.