读文献——轴承故障诊断【07.05-07.07】

文大于2

已于 2024-07-08 13:47:43 修改

阅读量739

点赞数 26

文章标签： matlab 信号处理

于 2024-07-08 12:32:40 首次发布

本文链接：https://blog.csdn.net/weixin_43213559/article/details/140264256

版权

Chen, Kehan, et al. "The adaptive bearing fault diagnosis based on optimal regulation of generalized SR behaviors in fluctuating-damping induced harmonic oscillator." Mechanical Systems and Signal Processing 189 (2023): 110078.

一.Introduction

二.System Model

2.1 GST-FDLO 系统

公式 (8) 是 GST-FDLO 系统的方程，描述了在潜在场 $U(x)=\frac{1}{2} \omega^2 x^2(t)$ 下，受噪声信号 g(t) 驱动的具有随机阻尼线性振荡器 (FDLO) 系统的行为。

其中，系统输入 $g(t)=u(t)+\varepsilon(t)$ 是原始噪声信号，包括轴承信号u(t)和加性高斯白噪声(AWGN) ε(t)

$\ddot{x}(t)$	系统的加速度响应
$\dot{x}(t)$	系统的速度响应
$x(t)$	系统的位移响应
$\gamma$	平均阻尼系数
$\omega$	系统的固有频率
$R$	广义尺度变换系数 (GST coefficient)
$\xi(t)$	对阻尼项的内部随机调节，建模为对称二分噪声 (SDN)

通过引入广义尺度变换 $R$ ，可以调节系统的时间尺度，使其更适用于处理故障特征。这种方法帮助研究人员分析系统对微弱轴承信号的放大能力，从而实现更有效的轴承故障诊断。

2.2 System stationary response

公式 (21) 计算了系统稳态响应的稳态振幅和相位 $\left(A_{a s}\right.$ 和 $\left.\varphi_{a s}\right)$ ，其中:

$c_d$ 是假设的常数值； $H(j 2 \pi f)$ 是系统的频率响应函数； $\mu_1, \mu_2, \mu_3, \mu_4$ 是与系统参数相关的系数。

公式 (23) 给出了输出振幅放大（OAA）的计算方法，用于衡量弱轴承信号的放大能力。

根据公式（21-23），可以做出figure2。

(a): 固有频率（ω²=1.0）、随机波动强度（σ²=0.2）、GST系数（λ=0.3）、参数R=1000、故障频率fd=100Hz、D=0.5时，OAA随阻尼系数γ的变化情况。可以看出，随着γ的增加，OAA呈下降趋势。

(b): 阻尼系数（γ=0.2）、随机波动强度（σ²=0.2）、GST系数（λ=0.1）、参数R=720、故障频率fd=100Hz、D=0.5时，OAA随固有频率ω²的变化情况。随着ω²的增加，OAA在某一值达到峰值后迅速下降，形成单峰GSR行为 。

(c):阻尼系数（γ=1.9）、固有频率（ω²=0.49）、GST系数（λ=0.01）、参数R=2000、故障频率fd=100Hz、D=0.5时，OAA随随机波动强度σ²的变化情况。OAA随着σ²的增加先增加后减少，呈现单峰变化趋势 。

(d):阻尼系数（γ=0.2）、固有频率（ω²=1.0）、随机波动强度（σ²=0.2）、参数R=720、故障频率fd=100Hz、D=0.5时，OAA随GST系数λ的变化情况。OAA随着λ的增加先急剧上升然后趋于平稳.

(e):阻尼系数（γ=0.2）、固有频率（ω²=1.0）、随机波动强度（σ²=0.2）、GST系数（λ=0.5）、故障频率fd=100Hz、D=0.5时，OAA随参数R的变化情况。OAA在R值增大时表现出非单调行为，峰值出现在特定的R值处 。

(f):阻尼系数（γ=0.2）、固有频率（ω²=1.0）、随机波动强度（σ²=0.2）、GST系数（λ=0.5）、故障频率fd=200Hz、D=0.5时，OAA随参数R的变化情况。与图2(e)类似，但R值的峰值位置与fd=100Hz时相比大约翻倍.

展示了系统中多个参数（如γ、ω²、σ²、λ和R）对弱故障信号放大能力的影响，通过调整这些参数可以优化故障特征的放大效果。

% 清空工作空间
clear; clc; close all;

% 参数设置
gamma_values = linspace(0, 4, 100);
omega2_values = linspace(0, 4, 100);
sigma2_values = linspace(0, 1, 100);
lambda_values = linspace(0, 2, 100);
R_values = linspace(0, 2000, 100);
fd = 100;
D = 0.5;
c_d = 1; % 假设的常数值

% 定义相关系数函数
mu1 = @(fd, R, gamma, omega2,lambda) -4*pi^2*fd^2 + R^2*lambda^2 + R^2*gamma*lambda + R^2*omega2;
mu2 = @(fd, lambda, R, gamma) 2*pi*fd*(2*lambda + gamma)*R;
mu3 = @(fd, R, gamma, omega2, lambda,sigma2) 4*pi^2*fd^2*R^2*gamma^2*sigma2 - (4*pi^2*fd^2 - R^2*omega2)*(R^2*lambda^2+R^2*gamma*lambda+ R^2*omega2-4*pi^2*fd^2)-4*pi^2*fd^2*gamma*(2*lambda + gamma)*R^2;
mu4 = @(fd, R, gamma, omega2, lambda,sigma2) 2*pi*fd*R*gamma*(R^2*lambda^2+R^2*gamma*lambda+ R^2*omega2-4*pi^2*fd^2)-2*pi*fd*(4*pi^2*fd^2 - R^2*omega2)*(2*lambda + gamma)*R-2*pi*fd*R^3*lambda*gamma^2*sigma2;

% 定义Aas计算函数
Aas = @(fd, R, gamma, omega2, lambda, sigma2) (R^2*c_d/2) * sqrt((mu2(fd, lambda, R, gamma)^2 + mu1(fd, R, gamma, omega2,lambda)^2) / (mu4(fd, R, gamma, omega2, lambda,sigma2)^2 + mu3(fd, R, gamma, omega2, lambda,sigma2)^2));

% 定义OAA计算函数
%OAA = @(fd, R, gamma, omega2, lambda, sigma2) (2/c_d) * Aas(fd, R, gamma, omega2, lambda, sigma2);
OAA = @(fd, R, gamma, omega2, lambda, sigma2) (R^2) * sqrt((mu2(fd, lambda, R, gamma)^2 + mu1(fd, R, gamma, omega2,lambda)^2) / (mu4(fd, R, gamma, omega2, lambda,sigma2)^2 + mu3(fd, R, gamma, omega2, lambda,sigma2)^2));

% 模拟结果设置
num_simulations = 100;
noise_std = 0.1;

% 打印调试信息
disp('开始计算 OAA ...');
disp(['fd = ', num2str(fd)]);
disp(['R = ', num2str(R_values(1))]);
disp(['gamma = ', num2str(gamma_values(1))]);
disp(['omega2 = ', num2str(omega2_values(1))]);
disp(['lambda = ', num2str(lambda_values(1))]);
disp(['sigma2 = ', num2str(sigma2_values(1))]);

% 计算不同参数下的OAA
try
    OAA_gamma = arrayfun(@(gamma) OAA(fd, 1000, gamma, 1, 0.3, 0.2), gamma_values);
    OAA_omega2 = arrayfun(@(omega2) OAA(fd, 720, 0.2, omega2, 0.1, 0.2), omega2_values);
    OAA_sigma2 = arrayfun(@(sigma2) OAA(fd, 2000, 1.9, 0.49, 0.01, sigma2), sigma2_values);
    OAA_lambda = arrayfun(@(lambda) OAA(fd,720, 0.2, 1, lambda, 0.2), lambda_values);
    OAA_R_100 = arrayfun(@(R) OAA(fd, R, 0.2, 1, 0.5, 0.2), R_values);
    OAA_R_200 = arrayfun(@(R) OAA(200, R, 0.2, 1, 0.5, 0.2), R_values);
catch ME
    disp(['Error occurred: ', ME.message]);
    rethrow(ME);
end

% 确认变量已成功计算
disp('变量计算成功');

% 确认变量已定义
if exist('OAA_gamma', 'var')
    % 绘制图表
    figure;

    % (a) G(gamma)
    subplot(3, 2, 1);
    plot(gamma_values, OAA_gamma, 'b', 'DisplayName', 'Analytical result');
    hold on;
    scatter(gamma_values, OAA_gamma + 0.1*randn(size(gamma_values)), 'r', 'filled', 'DisplayName', 'Simulation average');  % 模拟结果
    title('(a) G(\gamma)');
    xlabel('\gamma');
    ylabel('OAA');
    legend;
    grid on;

    % (b) G(omega^2)
    subplot(3, 2, 2);
    plot(omega2_values, OAA_omega2, 'b', 'DisplayName', 'Analytical result');
    hold on;
    scatter(omega2_values, OAA_omega2 + 0.1*randn(size(omega2_values)), 'r', 'filled', 'DisplayName', 'Simulation average');  % 模拟结果
    title('(b) G(\omega^2)');
    xlabel('\omega^2');
    ylabel('OAA');
    legend;
    grid on;

    % (c) G(sigma^2)
    subplot(3, 2, 3);
    plot(sigma2_values, OAA_sigma2, 'b', 'DisplayName', 'Analytical result');
    hold on;
    scatter(sigma2_values, OAA_sigma2 + 0.1*randn(size(sigma2_values)), 'r', 'filled', 'DisplayName', 'Simulation average');  % 模拟结果
    title('(c) G(\sigma^2)');
    xlabel('\sigma^2');
    ylabel('OAA');
    legend;
    grid on;

    % (d) G(lambda)
    subplot(3, 2, 4);
    plot(lambda_values, OAA_lambda, 'b', 'DisplayName', 'Analytical result');
    hold on;
    scatter(lambda_values, OAA_lambda + 0.1*randn(size(lambda_values)), 'r', 'filled', 'DisplayName', 'Simulation average');  % 模拟结果
    title('(d) G(\lambda)');
    xlabel('\lambda');
    ylabel('OAA');
    legend;
    grid on;

    % (e) G(R) with fd=100Hz
    subplot(3, 2, 5);
    plot(R_values, OAA_R_100, 'b', 'DisplayName', 'Analytical result');
    hold on;
    scatter(R_values, OAA_R_100 + 0.1*randn(size(R_values)), 'r', 'filled', 'DisplayName', 'Simulation average');  % 模拟结果
    title('(e) G(R) with f_d=100Hz');
    xlabel('R');
    ylabel('OAA');
    legend;
    grid on;

    % (f) G(R) with fd=200Hz
    subplot(3, 2, 6);
    plot(R_values, OAA_R_200, 'b', 'DisplayName', 'Analytical result');
    hold on;
    scatter(R_values, OAA_R_200 + 0.1*randn(size(R_values)), 'r', 'filled', 'DisplayName', 'Simulation average');  % 模拟结果
    title('(f) G(R) with f_d=200Hz');
    xlabel('R');
    ylabel('OAA');
    legend;
    grid on;

    % 调整子图布局
    sgtitle('GSR behaviors of OAA G varying with different system parameters');
else
    disp('OAA variables are not defined.');
end

后面的部分还没复现成功，问了下朋友，决定还是先把理论弄清楚些再继续。