Diagnóstico de Fallas en Rodamientos Mediante MATLAB: Extracción y Clasificación de Características

El diagnóstico de fallas en rodamientos es crucial para el mantenimiento predictivo y la fiabilidad de la maquinaria rotatoria. El proceso general se puede desglosar en cuatro etapas clave: adquisición de señales, preprocesamiento, extracción de características y reconocimiento de fallas. La extracción de características actúa como el nexo fundamental entre la señal bruta y la capacidad de diagnosticar fallas.

Preparación Inicial

1. Conjunto de Datos

Se empleará el Conjunto de Datos de Fallas en Rodamientos de Case Western Reserve University. Este conjunto de datos abarca:

  • 10 tipos de condiciones: Normal, Falla en Pista Interior, Falla en Pista Exterior, Falla en Cuerpo Rodante, entre otras.
  • Frecuencias de muestreo: 12 kHz (común) y 48 kHz (alta frecuencia).
  • 200 muestras por condición, cada una de longitud 1024.

2. Dependencias de Software

  • MATLAB Signal Processing Toolbox: Para el procesamiento de señales.
  • MATLAB Deep Learning Toolbox: Para aplicaciones de aprendizaje profundo.
  • Statistics and Machine Learning Toolbox: Para análisis estadístico y aprendizaje automático.

Pasos Centrales y Implementación en MATLAB

1. Preprocesamiento de Señales

Objetivo: Eliminar ruido y tendencias, realzando las características de la falla.

Métodos: Denoising mediante wavelets (Ondas) y eliminación de tendencias (Detrending).

Código MATLAB (Denoising y Detrending):


function senal_limpia = preprocesar_senal(senal_cruda)
    % Denoising con Wavelet (db4, nivel de descomposición 4)
    nombre_wavelet = 'db4';
    nivel = 4;
    [coeficientes, limites] = wavedec(senal_cruda, nivel, nombre_wavelet);
    
    % Umbralización (umbral blando)
    umbral = wdcbm2(coeficientes, limites, 'sqtwolog');
    tipo_umbral = 's'; % 's' para umbral blando
    mantener_aproximacion = 0; % No retener componente de aproximación
    senal_denoised = wdencmp('lvd', coeficientes, limites, nombre_wavelet, umbral, tipo_umbral, mantener_aproximacion);
    
    % Eliminación de tendencia (ajuste polinomial de orden 2)
    vector_tiempo = 1:length(senal_denoised);
    coef_polinomio = polyfit(vector_tiempo, senal_denoised, 2);
    tendencia = polyval(coef_polinomio, vector_tiempo);
    senal_limpia = senal_denoised - tendencia;
end

2. Extracción de Características

La extracción de características es fundamental. Se combinan características del dominio del tiempo, frecuencia y tiempo-frecuencia para capturar completamente la información de la falla.

(1) Características del Dominio del Tiempo

Estas características reflejan la distribución de amplitud y la correlación temporal de la señal, útiles para la detección temprana de fallas.

Características Comunes: Raíz cuadrática media (RMS), Pico (Peak), Pico a Pico (Peak-to-Peak), Curtosis, Asimetría (Skewness), Factor de Pulso (Pulse Factor).

Código MATLAB (Características Temporales):


function caracteristicas_tiempo = extraer_caracteristicas_tiempo(senal)
    % Raíz Cuadrática Media (RMS)
    caracteristicas_tiempo.rms = rms(senal);
    % Pico (Valor absoluto máximo)
    caracteristicas_tiempo.pico = max(abs(senal));
    % Pico a Pico
    caracteristicas_tiempo.pico_a_pico = max(senal) - min(senal);
    % Curtosis
    caracteristicas_tiempo.curtosis = kurtosis(senal);
    % Asimetría (Skewness)
    caracteristicas_tiempo.asimetria = skewness(senal);
    % Factor de Pulso
    caracteristicas_tiempo.factor_pulso = caracteristicas_tiempo.pico / caracteristicas_tiempo.rms;
end

(2) Características del Dominio de Frecuencia

Estas características muestran la distribución de frecuencias y son la base principal para el diagnóstico de fallas en rodamientos, ya que las frecuencias de falla se identifican en el análisis de frecuencia.

Características Comunes: Centroide Espectral (Spectral Centroid), Frecuencia Cuadrática Media (Mean Square Frequency), Varianza de Frecuencia (Frequency Variance), Curtosis Espectral (Spectral Kurtosis), Energía de Banda Espectral (Band Energy).

Código MATLAB (Características de Frecuencia):


function caracteristicas_frecuencia = extraer_caracteristicas_frecuencia(senal, fs)
    % Cálculo de la FFT
    N = length(senal);
    fft_datos = fft(senal);
    mag_fft = abs(fft_datos(1:floor(N/2)+1)); % Amplitud antes de normalizar
    mag_fft_norm = mag_fft / max(mag_fft); % Amplitud normalizada
    
    % Eje de Frecuencias
    eje_frecuencia = (0:floor(N/2)) * fs / N;
    
    % Cálculo de la Densidad Espectral de Potencia (PSD) usando pwelch
    [psd, frec_psd] = pwelch(senal, hanning(256), 128, 256, fs);
    
    % 1. Centroide Espectral
    caracteristicas_frecuencia.centroide = sum(eje_frecuencia .* mag_fft_norm) / sum(mag_fft_norm);
    
    % 2. Frecuencia Cuadrática Media
    caracteristicas_frecuencia.freq_cuadratica_media = sum(eje_frecuencia.^2 .* mag_fft_norm) / sum(mag_fft_norm);
    
    % 3. Varianza de Frecuencia
    caracteristicas_frecuencia.varianza_freq = sum((eje_frecuencia - caracteristicas_frecuencia.centroide).^2 .* mag_fft_norm) / sum(mag_fft_norm);
    
    % 4. Curtosis y Asimetría Espectral (calculados a partir de la distribución de probabilidad)
    pdf = mag_fft_norm / sum(mag_fft_norm); % Función de Densidad de Probabilidad
    denominador_asimetria = caracteristicas_frecuencia.varianza_freq^(3/2);
    if denominador_asimetria == 0, denominador_asimetria = eps; end % Evitar división por cero
    caracteristicas_frecuencia.asimetria_espectral = sum((eje_frecuencia - caracteristicas_frecuencia.centroide).^3 .* pdf) / denominador_asimetria;
    denominador_curtosis = caracteristicas_frecuencia.varianza_freq^2;
    if denominador_curtosis == 0, denominador_curtosis = eps; end % Evitar división por cero
    caracteristicas_frecuencia.curtosis_espectral = sum((eje_frecuencia - caracteristicas_frecuencia.centroide).^4 .* pdf) / denominador_curtosis;
    
    % 5. Energía de Banda Espectral: Dividido en 5 bandas (0-10%, 10-20%, etc. del ancho de banda Nyquist)
    bandas = [0, 0.1*fs/2; 0.1*fs/2, 0.2*fs/2; 0.2*fs/2, 0.3*fs/2; 0.3*fs/2, 0.4*fs/2; 0.4*fs/2, 0.5*fs/2];
    for i = 1:size(bandas, 1)
        mascara_banda = (eje_frecuencia >= bandas(i,1)) & (eje_frecuencia <= bandas(i,2));
        caracteristicas_frecuencia.(['energia_banda_', num2str(i)]) = sum(mag_fft_norm(mascara_banda));
    end
end

(3) Características del Dominio Tiempo-Frecuencia

Las señales de falla en rodamientos suelen ser no estacionarias (ej. fallas por impacto). Las características tiempo-frecuencia, como la energía de paquetes de ondas (Wavelet Packet Energy) o el espectro marginal de la Transformada de Hilbert-Huang (HHT), capturan información temporal y de frecuencia simultáneamente, siendo más adecuadas para señales no estacionarias.

Características Comunes: Energía de Paquetes de Ondas (Wavelet Packet Energy), Espectro Marginal de HHT.

Código MATLAB (Energía de Paquetes de Ondas):


function wt_features = extraer_wt_features(senal, fs)
    % Descomposición por Paquetes de Ondas (db4, nivel 4)
    nombre_wavelet = 'db4';
    nivel = 4;
    [wt, ~] = wpdec(senal, nivel, nombre_wavelet);
    
    % Cálculo de la energía de cada nodo
    num_nodos = 2^nivel;
    energia = zeros(1, num_nodos);
    for i = 0:num_nodos - 1
        coef_nodo = wpcoef(wt, i);
        energia(i+1) = sum(coef_nodo.^2);
    end
    
    % Normalización de la energía
    wt_features.energia_normalizada = energia / sum(energia);
    
    % Extracción de la Entropía de Energía (Energy Entropy)
    prob = wt_features.energia_normalizada / sum(wt_features.energia_normalizada + eps); % Añadir eps para evitar log(0)
    wt_features.entropia_energia = -sum(prob .* log2(prob));
end

(4) Cálculo de Frecuencias de Falla Característica

Las frecuencias de falla características son identificadores clave en el diagnóstico. Se calculan basándose en los parámetros del rodamiento (ej. frecuencia de falla en pista interior, exterior, cuerpo rodante).

Fórmulas (Ejemplo para Rodamientos de Bolas):

  • Frecuencia de Falla en Pista Interior (BPFI): \\(f\_{BPFI}=\\frac{n}{60}⋅f\_r⋅(1+\\frac{d}{D}cosϕ)\\)
  • Frecuencia de Falla en Pista Exterior (BPFO): \\(f\_{BPFO}=\\frac{n}{60}⋅f\_r⋅(1−\\frac{d}{D}cosϕ)\\)
  • Frecuencia de Falla en Cuerpo Rodante (BSF): \\(f\_{BSF}=\\frac{D}{2d}⋅f\_r⋅(1−(\\frac{d}{D}cosϕ)^2)\\)
  • Frecuencia de Falla en Jaula (FTF): \\(f\_{FTF}=\\frac{1}{2}⋅f\_r⋅(1−\\frac{d}{D}cosϕ)\\)

Parámetros:

  • \\(n\\): Número de cuerpos rodantes.
  • \\(f\_r\\): Frecuencia de rotación del eje (Hz, \\(f\_r=\\frac{rpm}{60}\\)).
  • \\(d\\): Diámetro del cuerpo rodante (mm).
  • \\(D\\): Diámetro del círculo de paso (mm).
  • \\(ϕ\\): Ángulo de contacto (°).

Código MATLAB (Cálculo de Frecuencias de Falla):


function frec_falla = calcular_frecuencia_falla(n, rpm, d, D, phi)
    % Frecuencia de Rotación del Eje (Hz)
    fr = rpm / 60;
    
    % Frecuencia de Falla en Pista Interior (BPFI)
    bpfi = fr * (1 + (d/D) * cosd(phi)) / 2;
    
    % Frecuencia de Falla en Pista Exterior (BPFO)
    bpfo = fr * (1 - (d/D) * cosd(phi)) / 2;
    
    % Frecuencia de Falla en Cuerpo Rodante (BSF)
    bsf = (D/(2*d)) * fr * (1 - (d/D * cosd(phi))^2);
    
    % Frecuencia de Falla en Jaula (FTF)
    ftf = fr * (1 - (d/D) * cosd(phi)) / 2;
    
    % Devolver frecuencias de falla
    frec_falla = struct('bpfi', bpfi, 'bpfo', bpfo, 'bsf', bsf, 'ftf', ftf);
end

3. Reconocimiento de Fallas

Objetivo: Utilizar las características extraídas para clasificar el tipo de falla del rodamiento (ej. normal, pista interior, pista exterior).

Métodos: Máquinas de Soporte Vectorial (SVM) (adecuadas para datos de alta dimensión y poca cantidad) o Redes Neuronales Convolucionales (CNN) (para características tiempo-frecuencia).

(1) Clasificación SVM

Código MATLAB (Entrenamiento SVM):


function modelo_svm = entrenar_svm(caracteristicas, etiquetas)
    % Normalización de Características (Z-score)
    caracteristicas_norm = zscore(caracteristicas);
    
    % División Entrenamiento/Prueba (70%/30%)
    rng(42); % Semilla aleatoria fija para reproducibilidad
    particion_cv = cvpartition(size(caracteristicas_norm, 1), 'HoldOut', 0.3);
    caracteristicas_entreno = caracteristicas_norm(training(particion_cv), :);
    etiquetas_entreno = etiquetas(training(particion_cv), :);
    caracteristicas_prueba = caracteristicas_norm(test(particion_cv), :);
    etiquetas_prueba = etiquetas(test(particion_cv), :);
    
    % Entrenamiento del Modelo SVM (Kernel RBF)
    modelo_svm = fitcsvm(caracteristicas_entreno, etiquetas_entreno, 'KernelFunction', 'rbf', 'BoxConstraint', 1, 'KernelScale', 'auto');
    
    % Evaluación del Modelo
    predicciones = predict(modelo_svm, caracteristicas_prueba);
    precision = sum(predicciones == etiquetas_prueba) / length(etiquetas_prueba);
    disp(['Precisión SVM: ', num2str(precision * 100), '%']);
end

(2) Clasificación CNN

Se presenta un ejemplo de CNN 1D para procesar características extraídas del dominio tiempo-frecuencia.

Código MATLAB (Entrenamiento CNN 1D):


function modelo_cnn = entrenar_cnn(caracteristicas, etiquetas)
    % Normalización de Características (Z-score)
    caracteristicas_norm = zscore(caracteristicas);
    
    % División Entrenamiento/Prueba (70%/30%)
    rng(42);
    particion_cv = cvpartition(size(caracteristicas_norm, 1), 'HoldOut', 0.3);
    caracteristicas_entreno = caracteristicas_norm(training(particion_cv), :);
    etiquetas_entreno = categorical(etiquetas(training(particion_cv), :));
    caracteristicas_prueba = caracteristicas_norm(test(particion_cv), :);
    etiquetas_prueba = categorical(etiquetas(test(particion_cv), :));
    
    % Definición de la Arquitectura 1DCNN
    num_clases = numel(categories(etiquetas_entreno));
    capas = [
        sequenceInputLayer(size(caracteristicas_entreno, 2)) % Capa de entrada (dimensión de características)
        convolution1dLayer(32, 5, 'Padding', 'same')      % Capa convolucional (32 filtros, kernel 5)
        batchNormalizationLayer()                         % Capa de normalización por lotes
        reluLayer()                                       % Capa de activación ReLU
        maxPooling1dLayer(2, 'Stride', 2)                 % Capa de max-pooling
        convolution1dLayer(64, 3, 'Padding', 'same')      % Capa convolucional (64 filtros, kernel 3)
        batchNormalizationLayer()
        reluLayer()
        maxPooling1dLayer(2, 'Stride', 2)
        fullyConnectedLayer(num_clases)                   % Capa totalmente conectada (número de clases)
        softmaxLayer()                                    % Capa Softmax
        classificationLayer()                             % Capa de clasificación
    ];
    
    % Opciones de Entrenamiento
    opciones = trainingOptions('adam', ...
        'MaxEpochs', 50, ...
        'MiniBatchSize', 32, ...
        'ValidationData', {caracteristicas_prueba, etiquetas_prueba}, ...
        'Plots', 'training-progress', ...
        'Verbose', 1);
    
    % Entrenamiento del Modelo CNN
    modelo_cnn = trainNetwork(caracteristicas_entreno, etiquetas_entreno, capas, opciones);
end

Ejemplo de Flujo Completo

1. Carga de Datos


% Cargar datos del dataset de Case Western Reserve University (Ejemplo: Falla Pista Interior)
% Supongamos que el archivo 'bearing_IR_fault.mat' contiene:
% vibracion (señal de vibración), rpm (velocidad), d (diámetro cuerpo rodante), 
% D (diámetro paso), phi (ángulo contacto), n (número cuerpos rodantes)
load('bearing_IR_fault.mat');

2. Preprocesamiento de Señal


senal_limpia = preprocesar_senal(vibracion);

3. Extracción de Características


% Características Temporales
caracteristicas_tiempo = extraer_caracteristicas_tiempo(senal_limpia);
% Características de Frecuencia (asumiendo fs=12000 Hz)
caracteristicas_frecuencia = extraer_caracteristicas_frecuencia(senal_limpia, 12000);
% Características Tiempo-Frecuencia (Energía de Paquetes de Ondas)
caracteristicas_wt = extraer_wt_features(senal_limpia, 12000);

% Combinar características (requiere manejo de estructuras/arrays)
% Ejemplo simplificado: concatenar valores numéricos
caracteristicas_combinadas = [struct2array(caracteristicas_tiempo), struct2array(caracteristicas_frecuencia), struct2array(caracteristicas_wt)];

% Calcular Frecuencias de Falla
frecuencias_falla = calcular_frecuencia_falla(n, rpm, d, D, phi);

4. Entrenamiento del Modelo SVM


% Asumimos que 'etiquetas' es un vector con los tipos de falla (ej. 1=Normal, 2=Pista Interior, etc.)
modelo_svm_entrenado = entrenar_svm(caracteristicas_combinadas, etiquetas);

5. Prueba del Modelo


% Simular una nueva señal de prueba (ej. Falla Pista Interior)
senal_prueba_cruda = ...; % Cargar o generar señal de prueba
fs_prueba = 12000;

% Preprocesar y extraer características de la señal de prueba
senal_prueba_limpia = preprocesar_senal(senal_prueba_cruda);
test_time = extraer_caracteristicas_tiempo(senal_prueba_limpia);
test_freq = extraer_caracteristicas_frecuencia(senal_prueba_limpia, fs_prueba);
test_wt = extraer_wt_features(senal_prueba_limpia, fs_prueba);
caracteristicas_prueba_combinadas = [struct2array(test_time), struct2array(test_freq), struct2array(test_wt)];

% Realizar predicción con el modelo SVM entrenado
prediccion = predict(modelo_svm_entrenado, caracteristicas_prueba_combinadas);
disp(['Tipo de falla predicho: ', num2str(prediccion)]);

Optimización del Rendimiento

  1. Selección de Características: Utilizar métodos como Laplacian Score o ReliefF para identificar las características más relevantes, reduciendo la dimensionalidad y mejorando la velocidad de entrenamiento.
  2. Aumento de Datos (Data Augmentation): Aplicar técnicas como adición de ruido, desplazamiento temporal o escalado de amplitud a los datos de entrenamiento para mejorar la robustez y generalización del modelo.
  3. Fusión de Modelos: Emplear aprendizaje por conjuntos (Ensemble Learning), como Random Forest o AdaBoost, combinando las prediccciones de múltiples modelos (ej. SVM y CNN) para aumentar la precisión.

Consideraciones Importantes

  • Calidad de Datos: Asegurar que las señales de vibración adquiridas estén libres de interferencias de ruido (utilizar filtros hardware o software de denoising).
  • Ajuste de Parámetros: Adaptar los parámetros del cálculo de frecuencias de falla según el tipo específico de rodamiento (ej. rodamiento rígido de bolas, rodamiento de rodillos cónicos).
  • Requerimientos de Tiempo Real: Para diagnósticos en tiempo real, considerar el uso de modelos ligeros (como SVM) o técnicas de cuantización de modelos (ej. TensorFlow Lite) para reducir la latencia computacional.

Aplicaciones Extensibles

  • Diagnóstico de Múltiples Fallas: Detectar simultáneamente varias fallas (ej. pista interior y cuerpo rodante), simplemente ampliando el conjunto de etiquetas.
  • Diagnóstico bajo Variación de Condiciones (Variable Operating Conditions): Aplicar aprendizaje por transferencia (Transfer Learning), como la adaptación de dominio, para migrar modelos entrenados en condiciones de laboratorio a entornos industriales con condiciones variables.
  • Alerta Temprana Inteligente: Integrar detección de anomalías (ej. One-Class SVM) para predecir fallas en sus etapas iniciales.

En resumen, el diagnóstico de fallas en rodamientos con MATLAB se centra en la extracción de características (dominio tiempo, frecuencia, tiempo-frecuencia) y la clasificación inteligente (SVM, CNN). Siguiendo este flujo, se puede lograr un diagnóstico de fallas de rodamientos con alta precisión y robustez, aplicable a la gestión de la salud de equipos industriales en sectores como la energía eólica, metalurgia y automoción.

Etiquetas: matlab Diagnóstico de Fallas Rodamientos procesamiento de señales Extracción de Características

Publicado el 10-27 12:06

Etiquetas populares

©2026 Friki Work