nass-micro-station-measurem.../ground-motion/index.org

371 lines
13 KiB
Org Mode

#+TITLE: Ground Motion Measurements
:DRAWER:
#+STARTUP: overview
#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas
#+HTML_LINK_HOME: ../index.html
#+HTML_LINK_UP: ../index.html
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/readtheorg.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/zenburn.css"/>
#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.stickytableheaders.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/readtheorg.js"></script>
#+HTML_MATHJAX: align: center tagside: right font: TeX
#+PROPERTY: header-args:matlab :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :output-dir figs
#+PROPERTY: header-args:shell :eval no-export
:END:
* Experimental Setup
The goal here is to compare the ground motion at the location of the micro-station (ID31 beamline at ESRF) with other measurements of the ground motion.
This will also permit to confirm that the device use for the measurement (geophones, amplifiers, ADC) are working well and that the data processing is correct.
One L22 geophone is put on the ID31 floor.
The signal is then filtered with a first order low pass filter with a cut-off frequency of $1kHz$.
Then the signal is amplified by a Voltage Amplifier with the following settings:
- AC/DC option set to DC
- Amplification of 60dB (1000)
- Low pass filter at the output with a cut-off frequency of $1kHz$
On figure [[fig:ground_motion_measurements]] is an overview of multiple measurements made at famous location.
#+name: fig:ground_motion_measurements
#+caption: Power spectral density of ground vibration in the vertical direction at multiple locations ([[https://vibration.desy.de/][source]])
#+attr_html: :width 800px
[[file:./img/ground_motion_measurements.png]]
* Measurement Analysis
:PROPERTIES:
:header-args:matlab+: :tangle matlab/ground_meas_id31.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:ground_meas_id31>>
** ZIP file containing the data and matlab files :ignore:
#+begin_src bash :exports none :results none
if [ matlab/ground_meas_id31.m -nt data/ground_meas_id31.zip ]; then
cp matlab/ground_meas_id31.m ground_meas_id31.m;
zip data/ground_meas_id31 \
mat/data_028.mat \
mat/id09_floor_september2018.mat \
mat/ground_motion_dist.mat \
ground_meas_id31.m
rm ground_meas_id31.m;
fi
#+end_src
#+begin_note
All the files (data and Matlab scripts) are accessible [[file:data/ground_meas_id31.zip][here]].
#+end_note
** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
<<matlab-dir>>
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+end_src
** Load data
#+begin_src matlab
data = load('mat/data_028.mat', 'data'); data = data.data;
#+end_src
** Time domain plots of the measured voltage
#+begin_src matlab
figure;
hold on;
plot(data(:, 3), data(:, 1));
hold off;
xlabel('Time [s]'); ylabel('Voltage [V]');
xlim([0, 100]);
#+end_src
#+NAME: fig:ground_motion_id31_time
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_id31_time.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:ground_motion_id31_time
#+CAPTION: Measurement of the ground motion - Time domain
#+RESULTS: fig:ground_motion_id31_time
[[file:figs/ground_motion_id31_time.png]]
** Computation of the ASD of the measured voltage
#+begin_src matlab :results none
dt = data(2, 3) - data(1, 3);
Fs = 1/dt;
win = hanning(ceil(10*Fs));
#+end_src
#+begin_src matlab :results none
[px_dc, f] = pwelch(data(:, 1), win, [], [], Fs);
#+end_src
#+begin_src matlab :results none
figure;
hold on;
plot(f, sqrt(px_dc));
hold off;
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude Spectral Density $\left[\frac{V}{\sqrt{Hz}}\right]$')
xlim([0.1, 500]);
#+end_src
#+NAME: fig:ground_motion_id31_asd_volt
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_id31_asd_volt.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:ground_motion_id31_asd_volt
#+CAPTION: Amplitude Spectral Density of the measured Voltage
#+RESULTS: fig:ground_motion_id31_asd_volt
[[file:figs/ground_motion_id31_asd_volt.png]]
** Scaling to take into account the sensibility of the geophone and the voltage amplifier
The Geophone used are L22. Their sensibility is shown on figure [[fig:geophone_sensibility]].
#+begin_src matlab :results none
S0 = 88; % Sensitivity [V/(m/s)]
f0 = 2; % Cut-off frequency [Hz]
S = S0*(s/2/pi/f0)/(1+s/2/pi/f0);
#+end_src
#+begin_src matlab :results none
freqs = logspace(-1, 2, 1000);
figure;
hold on;
plot(f, abs(squeeze(freqresp(S, f, 'Hz'))));
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude $\left[\frac{V}{m/s}\right]$');
xlim([0.1, 100]);
#+end_src
#+NAME: fig:geophone_sensibility
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/geophone_sensibility.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:geophone_sensibility
#+CAPTION: Sensibility of the Geophone
#+RESULTS: fig:geophone_sensibility
[[file:figs/geophone_sensibility.png]]
We also take into account the gain of the electronics which is here set to be $60dB$.
#+begin_src matlab :results none
G0_db = 60; % [dB]
G0 = 10^(G0_db/20); % [abs]
#+end_src
We divide the PSD measured (in $\text{V^2}/\sqrt{Hz}$) by the square of the gain of the voltage amplifier to obtain the PSD of the voltage across the geophone.
We further divide the result by the square of the magnitude of sensibility of the Geophone to obtain the PSD of the velocity in $(m/s)^2/Hz$.
#+begin_src matlab :results none
psd_gv = px_dc./abs(squeeze(freqresp(G0*S, f, 'Hz'))).^2;
#+end_src
Finally, we obtain the PSD of the ground motion in $m^2/Hz$ by dividing by the square of the frequency in $rad/s$.
#+begin_src matlab
psd_gm = psd_gv./(2*pi*f).^2;
#+end_src
** Time domain plots of the ground motion
We can inverse the dynamics of the geophone to convert the measured voltage into the estimated ground motion.
#+begin_src matlab
est_vel = lsim(inv(G0*S)*(s/2/pi)/(1+s/2/pi), data(:, 1), data(:, 3)); % Estimated velocity above 1Hz
est_vel = est_vel - mean(est_vel(data(:,3)>10)); % The mean value of the velocity if removed
est_dis = lsim(1/(1+s/2/pi), est_vel, data(:, 3)); % The velocity is integrated above 1Hz
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
plot(data(:, 3), est_vel);
hold off;
xlabel('Time [s]'); ylabel('Velocity [m/s]');
xlim([10, 100]);
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/time_domain_velocity.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:time_domain_velocity
#+CAPTION: Time domain velocity ([[./figs/time_domain_velocity.png][png]], [[./figs/time_domain_velocity.pdf][pdf]])
[[file:figs/time_domain_velocity.png]]
#+begin_src matlab :exports none
figure;
hold on;
plot(data(:, 3), est_dis);
hold off;
xlabel('Time [s]'); ylabel('Displacement [m]');
xlim([10, 100]);
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/time_domain_displacement.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:time_domain_displacement
#+CAPTION: Time domain absolute displacement ([[./figs/time_domain_displacement.png][png]], [[./figs/time_domain_displacement.pdf][pdf]])
[[file:figs/time_domain_displacement.png]]
** Computation of the ASD of the velocity and displacement
The ASD of the measured velocity is shown on figure [[fig:ground_motion_id31_asd_velocity]].
#+begin_src matlab :results none
figure;
hold on;
plot(f, sqrt(psd_gv));
hold off;
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD of the measured Velocity $\left[\frac{m/s}{\sqrt{Hz}}\right]$')
xlim([0.1, 500]);
#+end_src
#+NAME: fig:ground_motion_id31_asd_velocity
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_id31_asd_velocity.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:ground_motion_id31_asd_velocity
#+CAPTION: Amplitude Spectral Density of the Velocity
#+RESULTS: fig:ground_motion_id31_asd_velocity
[[file:figs/ground_motion_id31_asd_velocity.png]]
We also plot the ASD in displacement (figure [[fig:ground_motion_id31_asd_displacement]]);
#+begin_src matlab :results none
figure;
hold on;
plot(f, sqrt(psd_gm));
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD of the displacement $\left[\frac{m}{\sqrt{Hz}}\right]$')
xlim([0.1, 500]);
#+end_src
#+NAME: fig:ground_motion_id31_asd_displacement
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_id31_asd_displacement.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:ground_motion_id31_asd_displacement
#+CAPTION: Amplitude Spectral Density of the Displacement
#+RESULTS: fig:ground_motion_id31_asd_displacement
[[file:figs/ground_motion_id31_asd_displacement.png]]
And we also plot the PSD of the displacement in $\frac{{\mu u}^2}{Hz}$ as it is a usual unit used (figure [[fig:ground_motion_id31_psd_displacement]]).
One can then compare this curve with the figure [[fig:ground_motion_measurements]].
#+begin_src matlab :results none
figure;
hold on;
plot(f, psd_gm.*1e12);
hold off;
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('PSD of the measured displacement $\left[\frac{{ \mu m }^2}{Hz}\right]$')
xlim([0.1, 500]); ylim([1e-13, 1e3]);
#+end_src
#+NAME: fig:ground_motion_id31_psd_displacement
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_id31_psd_displacement.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:ground_motion_id31_psd_displacement
#+CAPTION: Power Spectral Density of the measured displacement
#+RESULTS: fig:ground_motion_id31_psd_displacement
[[file:figs/ground_motion_id31_psd_displacement.png]]
** Save
We save the PSD of the ground motion for further analysis.
#+begin_src matlab
save('./mat/psd_gm.mat', 'f', 'psd_gm', 'psd_gv');
#+end_src
** Comparison with other measurements of ground motion
Now we will compare with other measurements.
*** Load the measurement data
First we load the measurement data.
Here we have one measurement of the floor motion made at the ESRF in 2018, and one measurement made at CERN.
#+begin_src matlab
id09 = load('./mat/id09_floor_september2018.mat');
cern = load('./mat/ground_motion_dist.mat');
#+end_src
*** Compute PSD of the measurements
We compute the Power Spectral Densities of the measurements.
#+begin_src matlab
Fs_id09 = 1/(id09.time3(2)-id09.time3(1));
win_id09 = hanning(ceil(10*Fs_id09));
[id09_pxx, id09_f] = pwelch(1e-6*id09.x_y_z(:, 3), win_id09, [], [], Fs_id09);
#+end_src
#+begin_src matlab
Fs_cern = 1/(cern.gm.time(2)-cern.gm.time(1));
win_cern = hanning(ceil(10*Fs_cern));
[cern_pxx, cern_f] = pwelch(cern.gm.signal, win_cern, [], [], Fs_cern);
#+end_src
*** Compare PSD of Cern, ID09 and ID31
And we compare all the measurements (figure [[fig:ground_motion_compare]]).
#+begin_src matlab :results silent
figure;
hold on;
plot(id09_f, id09_pxx, 'DisplayName', 'ID09');
plot(cern_f, cern_pxx, 'DisplayName', 'CERN');
plot(f, psd_gm, 'k', 'DisplayName', 'ID31');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('PSD [$m^2/Hz$]');
legend('Location', 'northeast');
xlim([0.1, 500]);
#+end_src
#+NAME: fig:ground_motion_compare
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_compare.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:ground_motion_compare
#+CAPTION: Comparison of the PSD of the ground motion measured at different location
#+RESULTS: fig:ground_motion_compare
[[file:figs/ground_motion_compare.png]]