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

2 次查看（过去 30 天）
charbel
charbel 2024-6-11
评论： Star Strider 2024-6-11
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
  2 个评论
Sam Chak
Sam Chak 2024-6-11
Ran in:
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')

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2024-6-11
Ran in:
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.
.

更多回答（1 个）

sai charan sampara
sai charan sampara 2024-6-11
Ran in:
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.
  2 个评论
Image Analyst
Image Analyst 2024-6-11
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Earthquake Engineering 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by