Geodynamika 2011/2012
Ćwiczenie dodaktowe
Ostatnia aktualizacja 11 kwietnia 2012
1
.
Zadanie
Nale ˙zy przeprowadzi´c
analiz˛e widmow ˛a
(zwan ˛
a te ˙z analiz ˛a harmoniczn ˛a, analiz ˛a
Fouriera) dla danych grawimetrycznych z Józefosławia – grawimetr spr ˛e ˙zynowy
LC&R ET26. Prosz ˛e poda´c cz ˛estotliwo´sci wyró ˙znionych fal i ich amplitudy. Pro-
sz ˛e równie ˙z zamie´sci´c wykres z odpowiednio opisanymi osiami i warto´sciami.
Uzupełnienie 1.
Analiza widmowa
pozwala znale´z´c składowe harmoniczne (szereg
sinusoid, ich
amplitudy
i
cz ˛estotliwo´sci
) w analizowanym sygnale. Jest to nadzwy-
czaj u˙zyteczna technika, nie tylko w pływach lecz w całej geodezji jak i w innych
dziedzinach nauki i ˙zycia. Po szczegóły odsyłam do literatury i internetu. Krótko
mówi ˛ac transformacja i odwrotna transformacja Fouriera pozwalaj ˛a na zmian˛e sygnału
z dziedziny czasu na dziedzin˛e cz˛estotliwo´sci i odwrotnie.
1
.1.
Dane
Dane do zadania w postaci pliku tekstowego znajduj ˛
a si ˛e na stronie
www.geo.
republika.pl
. Fragment danych przedstawiony jest poni ˙zej – dane s ˛
a podane
w formacie yyyy mm dd
hh mm ss g
, co oznacza odpowiednio rok, miesi ˛
ac,
dzie ´n, godzin ˛e, minut ˛e i sekund ˛e oraz warto´s´c przyspieszenia siły ci ˛e ˙zko´sci
w nm/s
2
.
Wykaz 1. Dane do zadania (przykładowy fragment)
2007
1
2
16
0
0
-1059.343
2007
1
2
17
0
0
-1374.666
2007
1
2
18
0
0
-1783.991
2007
1
2
19
0
0
-2227.802
Uzupełnienie 2. Zał ˛aczone dane s ˛a wolne od „dziur”. W przypadku szybkiej trans-
formacji Fouriera nierównomierno´s´c próbkowania sygnału jest du˙zym problemem. S ˛a
sposoby, ˙zeby z tym sobie radzi´c ale to jest poza tematem tego ´cwiczenia. W ka˙zdym
razie nie mo˙zna stosowa´c tych algorytmów dla nierównomiernie rozło˙zonych danych
w dziedzinie cz˛estotliwo´sci.
Zestawy
Z danych nale ˙zy wybra´c swój podzbiór okre´slony datami (poni ˙zej)
w zale ˙zno´sci od numeru na li´scie (numery równie ˙z podane s ˛
a na stronie
internetowej.
Wykaz 2. Zestawy
Numer
0:
od 2007
1
2
16
0
do 2007
3
6
3
0
Numer
1:
od 2007
1 12
22
0
do 2007
3 16
9
0
Numer
2:
od 2007
1 23
4
0
do 2007
3 26
15
0
Numer
3:
od 2007
2
2
10
0
do 2007
4
5
21
0
Numer
4:
od 2007
2 12
16
0
do 2007
4 16
3
0
Numer
5:
od 2007
2 22
22
0
do 2007
4 26
9
0
Numer
6:
od 2007
3
5
4
0
do 2007
5
6
15
0
Numer
7:
od 2007
3 15
10
0
do 2007
5 16
21
0
Numer
8:
od 2007
3 25
16
0
do 2007
5 27
3
0
Numer
9:
od 2007
4
4
22
0
do 2007
6
6
9
0
Numer
10:
od 2007
4 15
4
0
do 2007
6 16
15
0
Numer
11:
od 2007
4 25
10
0
do 2007
6 26
21
0
Numer
12:
od 2007
5
5
16
0
do 2007
7
7
3
0
Numer
13:
od 2007
5 15
22
0
do 2007
7 17
9
0
Numer
14:
od 2007
5 26
4
0
do 2007
7 27
15
0
Numer
15:
od 2007
6
5
10
0
do 2007
8
6
21
0
Numer
16:
od 2007
6 15
16
0
do 2007
8 17
3
0
Numer
17:
od 2007
6 25
22
0
do 2007
8 27
9
0
Numer
18:
od 2007
7
6
4
0
do 2007
9
6
15
0
Numer
19:
od 2007
7 16
10
0
do 2007
9 18
17
0
Numer
84:
od 2009
5 16
15
0
do 2009
7 18
2
0
Numer
85:
od 2009
5 26
21
0
do 2009
7 28
8
0
Numer
86:
od 2009
6
6
3
0
do 2009
8
7
14
0
Numer
87:
od 2009
6 16
9
0
do 2009
8 17
20
0
Numer
88:
od 2009
6 26
15
0
do 2009
8 28
2
0
Numer
89:
od 2009
7
6
21
0
do 2009
9
7
8
0
Numer
90:
od 2009
7 17
3
0
do 2009
9 17
14
0
Numer
91:
od 2009
7 27
9
0
do 2009
9 27
20
0
Numer
92:
od 2009
8
6
15
0
do 2009 10
8
2
0
Numer
93:
od 2009
8 16
21
0
do 2009 10 18
8
0
Numer
94:
od 2009
8 27
3
0
do 2009 10 28
14
0
Numer
95:
od 2009
9
6
9
0
do 2009 11
7
20
0
Numer
96:
od 2009
9 16
15
0
do 2009 11 18
2
0
Numer
97:
od 2009
9 26
21
0
do 2009 11 28
8
0
Numer
98:
od 2009 10
7
3
0
do 2009 12
8
14
0
Numer
99:
od 2009 10 17
9
0
do 2009 12 18
20
0
Numer
100:
od 2009 10 27
15
0
do 2009 12 29
2
0
2
.
Rozwiązanie
Do oblicze ´n mo ˙zna wykorzysta´c istniej ˛
ace oprogramowanie. Spo´sród wielu
(bardzo wielu) mo ˙zliwo´sci mo ˙zna,
wykorzysta´c Matlaba lub jego darmowy (!) odpowiednik
1
Octave’a
Zach ˛ecam do u ˙zywania tych dedykowanego do oblicze ´n matematycznych
programów. Na zach ˛et ˛e poni ˙zej podaje przykładowy skrypt rozwi ˛
azuj ˛
acy
zadanie w Octave’ie (linie komentarza zaczynaj ˛
a si ˛e od znaku %).
% Wczytanie danych z pliku do macierzy A
A =
load
("./zestawy/dane0.dat");
% Ile pozycji jest w danych
L=
length
(A);
% Wykorzystanie wbudowanej funkcji fft.
% Szybka transformata fouriera wymaga,
% zeby liczba danych byla potega dwojki
% dlatego wykorzystujemy funkcje
% ’nextpow2’ (patrz wyjasnienia w manualu.
FFT=
fft
( A(:,7) , 2^
nextpow2
(L) ) / L ;
% Czestotliwosc i amplituda
freq = 1/2 *
linspace
(0,1,2^
nextpow2
(L) / 2+1 )’;
ampl = 2 *
abs
( FFT( 1:2^
nextpow2
(L) / 2+1 ) );
% Mielismy dane godzinne a wyniki chcemy w cpd
freq = freq * 24;
% I zapis pliku wynikowego
PLIK=
fopen
(’wyniki.dat’,’w’);
for
i=1:
length
(freq)
fprintf
(PLIK,"%10.3f %10.3f\n",freq(i),ampl(i))
end
wykorzysta´c Excela – jest mo ˙zliwo´s´c przeprowadzenia analizy Fouriera
w tym pakiecie – ale tutaj nic nie podpowiem gdy ˙z sam tego nie robiłem.
Z tego co wiem open office nie oferuje jeszcze takich oblicze ´n.
wykorzysta´c program Tsoft
http://seismologie.oma.be/TSOFT/tsoft.
html
Doskonałe i łatwe w u ˙zyciu (okienka, menu) oprogramowanie do analizy
1
prawie odpowiednik
3
szeregów czasowych. Ma te ˙z moduł do oblicze ´n teoretycznych pływów
ziemskich. Obszerna dokumentacja i pliki pomocy pozwol ˛
a Pa ´nstwu na
stosunkowo łatwe rozwi ˛
azanie problemu. Poni ˙zej pokazany jest „zrzut
ekranu” z programu Tsoft z rozwi ˛
azaniem dla zestawu nr 0.
znale´z´c sobie co´s „pod siebie” – ogromny wybór.
3
.
Przykład
Rozwi ˛
azanie dla zestawu nr 0.
-4000
-3500
-3000
-2500
-2000
-1500
-1000
-500
2007
-01-09
2007
-01-24
2007
-02-08
2007
-02-23
∆
g
[
n
m
/
s
2
]
Dane pomiarowe
Rozwi ˛
azanie wykorzystuj ˛
ace kod Octave’a. Wyniki w pliku,
Wykaz 3. Wyniki działania skryptu (fragment)
0.926
375.653
0.938
252.204
0.949
84.242
0.961
50.482
4
0.973
63.132
0.984
76.388
0.996
455.769
1.008
459.650
1.020
72.500
1.031
61.801
1.043
45.080
1.055
53.428
1.066
55.850
0
500
1000
1500
2000
2500
3000
3500
0
2
4
6
8
10
12
∆
g
[
n
m
/
s
2
]
Cz ˛estotliwo´s´c
[
cpd
]
Widmo
Powy ˙zszy rysunek jest nieczytelny ze wzgl ˛edu na wielkie warto´sci (wła-
´sciwo´s´c metody) dla bardzo małych cz ˛estotliwo´sci. Ograniczmy si ˛e do
najciekawszego obszaru (w kontek´scie pływów). Oczywi´scie pływy długo-
okresowe s ˛
a nadzwyczaj interesuj ˛
ace jednak w naszych danych niemo ˙zliwe
do wyró ˙znienia ze wzgl ˛edu na dryft grawimetru.
0
50
100
150
200
250
300
350
400
450
500
1
1
.5
2
2
.5
∆
g
[
n
m
/
s
2
]
Cz ˛estotliwo´s´c
[
cpd
]
Widmo
5
I wystarczyłoby teraz z pliku lub wykresu (z mała dokładno´sci ˛
a) odczy-
ta´c cz ˛estotliwo´sci i amplitudy głównych pików w widmie (
fal pływowych
).
Mo ˙zna si ˛e równie ˙z pokusi´c o ich nazwanie.
Uzupełnienie 3. Rozdzielczo´s´c widma silnie zale˙zy od liczby danych pomiarowych.
Im dłu˙zszy ci ˛ag obserwacyjny, tym wi˛ecej szczegółów (tym samym ró˙znych fal) mo˙zemy
w widmie wyró˙zni´c. Przykład widma z pełnego roku obserwacji znajd ˛a pa ´nstwo
w „podlinkowanej” pracy
[pdf]
.
Uzupełnienie 4. Na koniec zwracam uwag˛e, ˙ze cz˛esto przy analizie harmonicznej
wykorzystuje si˛e
moc widmow ˛
a
. W przypadku naszych danych jednostk ˛a mocy
widmowej były by
(
nm
·
s
−
2
)
2
·
Hz
−
1
. W przypadku
widma amplitudowego
mamy
wprost amplitudy fal pływowych – nm
·
s
−
2
.
Powodzenia!
6
