How to correct the recording of a daamaged accelerometer in earthquake analysis

Hello everyone,

I have the recording of two accelerometers but it turned out that one of them did not work properly and the recording is messed up (See photo, The two recordings should be the same). However I think if processed well I can recover the lost acclerogram. Does anyone have an idea how I can proceed. I have attached the recordings.

Thank you in advance

[Sam Chak]

Hi @charbel

This is an attempt to plot the results. As I do not have expertise in seismology, I cannot definitively determine which output is preferable and which is less so. However, a typical accelerogram would be expected to resemble the first (blue) trace.

load('acceleration.mat');

t   = acceleration(:,1);    % time
dat1= acceleration(:,2);    % blue data
dat2= acceleration(:,3);    % red  data

subplot(211)
plot(t, dat1, 'color', "#0072BD"), grid on
xlim([t(1), t(end)]), ylim([-1, 2]), xlabel('Time'), ylabel('Acc / [cm/s^{2}]')
title('Accelerogram 1')

subplot(212)
plot(t, dat2, 'color', "#D95319"), grid on
xlim([t(1), t(end)]), ylim([-1, 2]), xlabel('Time'), ylabel('Acc / [cm/s^{2}]')
title('Accelerogram 2')

[charbel]

Thank you for your time to answer

Answer 1

[Star Strider]

If you have the Signal Ptocessing Toolbox, an alternative approach would be to use the highpass (or bandpass to also eliminiate the high-frequency noise spike at about 16.5 Hz) function to filter out the low-frequency components from the ‘broken’ accerometer.

That would be something like this —

load('acceleration.mat')
whos('file','acceleration')

  Name                  Size             Bytes  Class     Attributes

  acceleration      24000x3             576000  double              

acceleration

acceleration = 24000x3

0 0.0861 0.0806 0.0052 0.0868 0.0834 0.0104 0.0687 0.0651 0.0156 0.0393 0.0353 0.0208 0.0038 -0.0016 0.0260 -0.0243 -0.0321 0.0312 -0.0427 -0.0498 0.0365 -0.0429 -0.0493 0.0417 -0.0240 -0.0272 0.0469 0.0044 0.0015

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
plot(acceleration(:,1), acceleration(:,[2 3]))
grid
xlabel('Time')
ylabel('Amplitude')
title('Original Time Domain')
legend('Acceleromentr_1','Accelerometer_2', 'Location','best')

[FTs1,Fv] = FFT1(acceleration(:,[2 3]), acceleration(:,1));

figure
plot(Fv, abs(FTs1)*2)
grid
% set(gca, 'XScale','log')
xlabel('Frequency')
ylabel('Magnitude')
title('Fourier Transform')
xlim([0 20])

figure
plot(Fv, abs(FTs1)*2)
grid
% set(gca, 'XScale','log')
xlabel('Frequency')
ylabel('Magnitude')
title('Fourier Transform (Zoomed)')
xlim([0 1])
xline(0.45, '--k', 'Highpass Filter:  F_{co} = 0.5 Hz')

Fs = 1/mean(diff(acceleration(:,1)))                        % Sampling Frequency

Fs = 192

Fstd = std(diff(acceleration(:,1)))                         % Sampling Frequency Variation

Fstd = 5.0181e-15

acc_hp_filt = highpass(acceleration(:,[2 3]), 0.5, Fs, 'ImpulseResponse','iir');

figure
plot(acceleration(:,1), acc_hp_filt)
grid
xlabel('Time')
ylabel('Amplitude')
title('Highpass Filtered Time Domain')

figure
plot(acc_hp_filt(:,1), acc_hp_filt(:,2))
grid
xlabel('Accelerometer_1 (Highpass Filtered)')
ylabel('Accelerometer_2 (Highpass Filtered)')
axis('equal')

acc_bp_filt = bandpass(acceleration(:,[2 3]), [0.5 7.5], Fs, 'ImpulseResponse','iir');

figure
plot(acceleration(:,1), acc_bp_filt)
grid
xlabel('Time')
ylabel('Amplitude')
title('Bandpass Filtered Time Domain')

figure
plot(acc_bp_filt(:,1), acc_bp_filt(:,2))
grid
xlabel('Accelerometer_1 (Bandpass Filtered)')
ylabel('Accelerometer_2 (Bandpass Filtered)')
axis('equal')

function [FTs1,Fv] = FFT1(s,t)
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv,:);
end

It is necessary to use the Fourier transform in order to design the filtter passband frequencies correctly. This is overly detailed to show the relative effects of the two filter approaches. I filtered both signals with the same filters.

.

[charbel]

Thank you very much for this great response

[Star Strider]

As always, my pleasure!

Answer 2

[charan]

Hello,

A possible solution can be to remove the mean over a small interval from the existing data to change the data to be centered around zero or around the mean of the correct data. It can be done similar to the code shown below:

load("acceleration.mat");
plot(acceleration(:,1),acceleration(:,2));
hold on
plot(acceleration(:,1),acceleration(:,3));
hold off

n=acceleration(:,3);
for i=1:10:length(acceleration(:,3))-9
    n(i:i+9)=acceleration(i:i+9,3)-mean(acceleration(i:i+9,3))+mean(acceleration(:,2));
end
plot(acceleration(:,1),acceleration(:,2));
hold on
plot(acceleration(:,1),n,'g');
hold off
legend(["Accelogram 1","Corrected Accelogram 2"]);

The number of datapoints to be taken at a time for "mean" can be changed based on the requirement.

[Image Analyst]

To avoid changing the (probably) "good" data before the sensor went bad, you can start the correction just where it starts to go bad, like around x=50 or so, instead of at the first index.

[charbel]

Thank you for your great response