Engee 文档
Notebook

铁路运输的第三个八度分析

第三倍频程分析是一种对信号进行频率分析的方法，其中频谱被划分为频带，每个频带是1/3倍频程。 对列车使用第三个八度音程分析可以研究噪声和振动，并评估它们对环境、乘客和基础设施的影响。 这种方法对于:

*噪声归一化（符合GOST，ISO，EN）。

*铁路轨道和机车车辆状况的诊断。

*优化结构(例如减少轮组或导轨的噪音)。

安装在轮组、车架、轨道和地面上的加速度计通常用于分析（评估建筑物振动的传播）。 对于铁路运输，关键的第三个八度频带通常位于1-1000Hz（振动）和20-5000Hz（噪声）的范围内。 重要频率的例子:

*低（1-80Hz）：车轮撞击轨道接头产生的振动。

*中等（80-500赫兹）：滚动噪音，发动机嗡嗡声。

*高（500-5000赫兹）：吱吱作响的刹车，空气动力学噪音。

A三分之一倍频程谱的构建允许您分析某些波段的峰值，这些峰值可能表明:

*车轮磨损（在63-250赫兹波段增加）。

*钢轨缺陷（31.5–125Hz时的谐振）。

*悬挂问题（2-20Hz振动）。

让我们考虑分析来自铁路轨道上加速度计的真实数据并构建三分之一倍频程频谱的问题。

数据导入和预处理

连接必要的库:

In [ ]: 
using DSP, DelimitedFiles, Statistics, FFTW, Polynomials

加速度计数据存储在文件中 data.txt. 我们把它们算作一个变量中的矩阵 datamat:

In [ ]: 
datamat = readdlm("data.txt", ',');

让我们从中提取时间和"原始"振动加速度值的单独矢量（以m/s^2为单位），并将它们显示在时域的图形上。:

In [ ]: 
timevec = datamat[:,1];
data1 = datamat[:,2];
Plots.plot(timevec,data1,
            title = "Импортированные данные виброускорений",
            xlabel = "Время, с",
            ylabel = "Виброускорение, м/с^2",
            legend = false)
Out[0]:

我们看到振动加速度值的记录在初始时刻没有立即开始，因此切断接近于零的部分信号电平是合乎逻辑的。 要做到这一点，我们将只留下对应于超过某个阈值（在我们的例子中，1）的信号值的矢量索引。:

In [ ]: 
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)
Out[0]:

过渡到振动速度

为了可视化振动速度，我们需要对加速度值的矢量进行积分。 accel. 让我们编写一个自定义累积梯形数值积分函数。 cumtrapz，环境中使用的类似物 MATLAB:

In [ ]: 
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
Out[0]:
cumtrapz (generic function with 1 method)

让我们将自定义函数应用于矢量 accel，然后我们将在振动速度信号中选择一个恒定的分量（或趋势），以便在必要时通过元素相减将其从信号中删除。 为此，我们用12阶多项式近似振动速度信号，并在与速度信号相同的时间点估计其值。:

In [ ]: 
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 = "Виброскорость, м/с")
Out[0]:

Welch方法的光谱分析

我们将使用标准库包对振动加速度进行光谱分析。 DSP.jl. 首先，让我们确定采样频率的值 fs 以Hz为单位的原始信号，分析时间矢量的步 dt:

In [ ]: 
dt = mean(diff(timevec));
fs = Int(floor(1/dt))
Out[0]:
3472

让我们设置职权范围中包含的频谱分析参数。 特别是：用于FFT算法的输入的"窗口"的大小，用于内插功率谱密度图（SPM)的零的数目，用于Welch方法的50%的重叠，FFT样本的输出数目（"窗口"的大小+零的数目）。

In [ ]: 
window_length = fs;
nzeros = 100000;
overlap = 0.5;
nFFT = window_length + nzeros;
noverlap = Int(floor(overlap*window_length));
stepsize = window_length - noverlap;

值得注意的是，输入"窗口"的点数等于采样频率。 这使得能够以1Hz的分辨率估计频谱。

此外，根据分配，分析仅在正频率的频谱的下半部分（第一奈奎斯特区）进行。 为此，信号通过12阶巴特沃斯低通滤波器，其归一化截止频率为0.5。 过滤结果记录在变量中 filtered:

In [ ]: 
Wp = 0.5;
myfilter = digitalfilter(Lowpass(Wp), Butterworth(12));
filtered = filtfilt(myfilter,accel);

现在让我们使用该函数 welch_pgram 用Welch方法来估计频谱功率密度。 重叠参数，输入"窗口"的大小和FFT点的输出数量对应于先前设置的那些。 让我们在dB中显示SPM图表:

In [ ]: 
 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 = "Спектральная плотность мощности, дБ/Гц")
Out[0]:

我们还将以信号的测量单位显示振幅谱。 为此，我们需要确定汉明窗函数的等效噪声带的参数，我们在Welch方法中使用了该参数。:

In [ ]: 
mywin = hanning(window_length);
ENBW = fs * (sum(abs.(mywin).^2)) / (abs.(sum(mywin)).^2)
Out[0]:
1.5004321521175445

我们从频谱的功率值中获得振幅值，同时考虑到参数 ENBW:

In [ ]: 
amplitude = sqrt.(2*power(pxx)*ENBW);
Plots.plot(freq(pxx), amplitude, legend = false, 
            title = "Амплитудный спектр сигнала виброускорений",
            ylabel = "Амплитуда, м/с^2",
            xlabel = "Частота, Гц",
            xlim = (0,1000))
Out[0]:

第三个八度频谱

要构建三分之一倍频程频谱，首先需要确定频率范围的边界及其数量。 count:

In [ ]: 
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:

In [ ]: 
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:

In [ ]: 
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)
Out[0]:

结论

第三世纪分析是诊断机车车辆腹板或轮组的常用方法，它允许检测例如钢轨的不均匀磨损或轮轨接触问题。 Gost31319-2009和GOST31296.1-2019中规定了类似范围的频谱分析。

我们已经审查了可用的工具进行这种分析的可能性。