Opis metody liczenia magnitudy

Moduł Magnituda Mw liczy magnitudy stacyjne. Magnituda Mw liczona jest ze skalara momentu sejsmicznego M0 ,  który jest liczony ze spektrum przemieszczenia fali P lub S w ognisku metodą Andrew (Andrews 1986, Snoke 1987). Sygnał liczony jest osobno dla trzech składowych, po czym wartości omega dla poszczególnych składowych są sumowane i z sumy liczone jest M0 (Wiejacz i Wiszniowski 2006).

Metoda rozpoczyna się od obliczenia widma sygnału metodą multi-taper (Park,1987). Metoda ta wykorzystuje wielokrotne sumowanie widm cząstkowych, które są uzyskane dla różnych okien czasowych (Niewiadomski 1997). Skalowanie widma odbywa się w oparciu o równość Parsevala, zgodnie ze wzorem:


gdzie F(f) jest widmem obliczonym metodą multitaper, podczas gdy współczynnik skalujący Sc jest określony wzorem

gdzie

Za początek fali przyjmuje się początek okna czasowego (tP, tE). Długość okna czasowego jest dobierana tak, aby pasowała do kształtu fali i dlatego do pewnego stopnia jest to wynikiem arbitralnej oceny ludzkiego analityka. Domyślna długość okna czasowego może być ustawiona na każdej stacji dla każdego konkretnego regionu źródła zdarzenia. Jednak taka w pełni zautomatyzowana procedura może prowadzić do błędnych wyników w przypadku zmieniającego się poziomu hałasu sejsmicznego lub zbliżonych do siebie zdarzeń. Ważnym czynnikiem jest również wielkość wstrząsu; generalnie w przypadku większych zjawisk należy przyjąć dłuższe okno czasowe (jest to zasada metody długości trwania). W przypadku słabego zdarzenia, wielkość okna odpowiednia do większego wstrząsu może spowodować, że udział hałasu długookresowego może wpłynąć na podniesienie wielkości. W przeciwieństwie do tego, okno zbyt krótkie dla dużego zdarzenia powoduje, że jego długookresowe widma cząstkowe nie wchodzą w sumę i ostatecznie wielkość Mw może zostać błędnie obniżona.

Końcowe widmo przemieszczeń jest obliczane za pomocą wzoru uwzględniającego rozchodzenie się fal ciała w odległości R, tłumienie sprężyste zgodnie ze współczynnikiem jakości Q (Aki i Richards, 2002) oraz odpowiedź instrumentu GV(f):

Uwzględnianie odpowiedzi przyrządu jest ważne, ponieważ pozwala na porównywanie wyników i wykorzystanie danych zarejestrowanych przez różne czujniki i systemy akwizycji danych na różnych stacjach i/lub w różnym czasie po zmianie wyposażenia stacji.

Poprawka na tłumienie sprężyste w przypadku słabych zdarzeń polega w zasadzie na dużym wzmocnieniu wysokich częstotliwości w widmie, co wynika z postaci zależności tłumienia. W paśmie wysokich częstotliwości sygnał zdarzenia może być słabszy niż szum sejsmiczny, podczas gdy zarówno sygnał, jak i szum zostaną silnie wzmocnione przez korekcję. W takich przypadkach wartość Q jest podnoszona tak, aby kształt widma w danym paśmie częstotliwości pokrywał się z krzywą modelu Brune'a (1970). Wartości Ω0 i f0 wymagane przez model Brune'a oblicza się metodą Andrewsa (1986):

gdzie J i K są obliczane według wzorów Snoke'a (1987)

gdzie f1 i f2 to dolna i górna granica częstotliwości pasma, dla którego estymowany jest model Brune'a. Wyznaczenie granic pasma częstotliwości jest konieczne ze względu na hałas sejsmiczny. Wzory  Snoke'a są poprawne, jeśli

f1 << f0 << f2          (1)

ponieważ zakładają, że odpowiedź modelu Brune'a poza tym zakresem jest asymptotyczna. Wydaje się zatem, że wystarczające jest spełnienie

2f1 < f0 < ½f

co może być niemożliwe do spełnienia w przypadku słabego sygnału i dużych szumów. Dolna granica f1 dla większości zdarzeń wynosi 1 Hz. W szczególności dotyczy to słabych wstrząsów. Przy częstotliwościach poniżej 1 Hz szum o niskiej częstotliwości zaczyna odgrywać istotną rolę w widmach sygnału. Hałas ten staje się dominujący w przypadku słabych i odległych zdarzeń, dlatego należy odciąć niskie częstotliwości. W przypadku zdarzeń bliskich widma sygnału są zwykle większe niż szum i można zastosować dolną granicę pasma częstotliwości. Silne zdarzenia mają w swoim sygnale większą składową niskich częstotliwości, dlatego zgodnie z zależnością (1) trzeba zastosować dolną granicę pasma niskich częstotliwości mniejszą niż 1 Hz. Sygnał jest wtedy na tyle silny, że przeważa nad szumem w całym paśmie częstotliwości. Górna granica pasma częstotliwości f2 ma również na celu ograniczenie skutków szumu. Ograniczenie to ma znaczenie w przypadku słabych sygnałów, gdy korekcja tłumienia powoduje nadmierne wzmocnienie wysokich częstotliwości. W przypadku słabych zdarzeń konieczne staje się ustalenie górnej granicy pasma częstotliwości lub zwiększenie wartości Q. Ogólnie rzecz biorąc, efekty szumu uniemożliwiają poprawne określenie wielkości Mw poniżej wartości 1.7-2,0 na dowolnej stacji, ponieważ przy tej wielkości poziom widmowy jest następnie obliczany z widm szumu, a nie z widm sygnału zdarzenia.

Moment sejsmiczny jest obliczany z trzech składowych sygnału dwoma alternatywnymi metodami. W pierwszej metodzie wartości J, KΩ0 i f0 są obliczane dla każdego składowej osobno. Następnie obliczany jest moment sejsmiczny dla każdego kanału dla fali Pg lub Sg według wzoru:

gdzie i oznacza składową sygnału, a c0 jest odpowiednią prędkością fali (P lub S) u źródła. Prędkość fali S u źródła, c0 jest określana na podstawie głębokości i modelu prędkościowego 1D, natomiast prędkość fali P u źródła jest określona, albo z modelu, albo zależnością c0(P) = 1,73 c0(S). Przyjmuje się gęstość w źródła z modelu 1D lub domyślnie ρ0 = 2700 kg/m3. Średni współczynnik promieniowania Rc przyjmuje się dla fali P jako Rc(P)=0,52 i dla fali S Rc(S)=0,63 (Boore i Boatwright, 1984). Moment sejsmiczny jest obliczany jako pierwiastek kwadratowy z sumy kwadratów poszczególnych momentów sejsmicznych uzyskanych dla trzech składowych:

W drugiej metodzie wartości J i K są obliczane przy użyciu wszystkich trzech składowych sygnału:

Stosując tę metodę otrzymujemy tylko jeden zbiór wartości Ω0 i f0 , z których oblicza się M0. Metoda ta jest liczbowo lepsza niż pierwsza w sytuacji, gdy jakiś komponent ma stosunkowo bardzo słaby sygnał. W takim przypadku wartości J i K dla takiego składnika będą zarówno bardzo małe, a w pierwszej metodzie powstałby przypadek dzielenia dwóch bardzo małych liczb. Wynik takiego dzielenia byłby następnie sumowany na równych prawach z dwoma innymi niezerowymi składnikami. Jednak w praktyce Mw uzyskiwane dwoma wymienionymi metodami były zawsze prawie identyczne, co doprowadziło do stosowania głównie pierwszej metody.

Wielkość Mw jest obliczana z momentu sejsmicznego według wzoru Hanksa i Kanamoriego (1979). Wynik może być oparty na fali P lub S lub średniej z dwóch:

Domyślna długość okna przyjmowana dla liczenia spektrum wynoś 0.9 równicy czasów fali S i P dla fali P i 1.8 tej różnicy dla fali S. W przypadku braku piku P lub S wynosi róznica czasów wynosi (√3-1)TP lub (1=1/√3)TS, gdzie TP lub TS są różnicami między czasami P lub S a czasem w ognisku.

Instrukcja liczenia magnitudy Mw

Liczenie magnitudy Mw wywoływane jest z menu poleceniem Magnitude  Mw(SWIP4). Nazwa (SWIP4) informuje, że moduł ten jest programem zaadoptowanym z programu SWIP4. Po wywołaniu jawi się okno do liczenia magnitudy Mw (Rys. 17). U góry okna wypisywana są Origin (1), dla którego liczona będą magnitudy, oraz stacja sejsmiczna (2). Poniżej wyświetlany są sejsmogramy trzech składowych wybranej stacji (3), W zależności od wyboru opcji View seismogram na dole okna może być wyświetlany sejsmogram prędkościowy (opcja Vel.) lud przemieszczeniowy (opcja Disp.).Każdy kanał rysowany jest osobnym kolorem. Takim samym kolorem są rysowane wydma sygnałów. Czerwonymi  liniami zaznaczone są fale P i S a czarną linią rysowane jest w postaci prostokąta okno, z którego liczone jest widmo spektrum (4). Początek okna pokrywa się wstępnie z fazą P lub S, ale granice okien można zmienić ciągnąc je myszką z wciśniętym lewym przyciskiem myszy.

Rys. 17 Okno liczenia magnitudy Mw

Poniżej sejsmogramów rysowane są widma przemieszczeniowe poszczególnych składowych oraz krzywe teoretycznych widm modelu Bruna (5) estymowanych z sygnału. Z lewej strony rysowane są widma dla fali P, a z prawej strony dla fali S. Pogrubioną linią modelu zaznaczony jest przedział częstotliwości, dla którego estymowany jest model – wartości między dolnym a górnym limitem częstotliwości (6: Lower frequency [Hz] i Upper frequency [Hz]). Kolory oznaczają: niebieski – składowa Z, błękitny– składowa N, zielony – składowa E.

Poniżej z lewej strony znajduje się okno parametrów do liczenia Mw  (6) a z prawej strony wyświetlane są wartości M0, f0 i Mw dla fal P, S i całości oraz przyciski potwierdzające magnitudę (7). Ich naciśnięcie powoduje zapisanie danych Mw dla danej stacji. Można wybrać Mw z fali P, S lub średnią z obydwu fal..

Każda zmiana parametrów jak i okno do liczenia skutkuje przeliczeniem Mw i zmianą krzywych Brune'a. Wartości dolnej i górnej częstotliwości oraz współczynniki tłumienia Q są zapamiętywane w pliku konfiguracyjnym i przy powtórnym liczeniu przyjmowane jako domyślne, natomiast prędkości i gęstość są wyliczane z modelu 1D. Jeżeli w pliku konfiguracyjnym nie zdefiniowano jednej lub wielu stacji wcześniej i ustawiona jest w pliku konfiguracyjnym opcja Save, to pojawi się ostrzeżenie

a wartości dolnej i górnej częstotliwości oraz współczynniki tłumienia Q ustawią się na wartości domyślne. Koniecznie należy wtedy wpisać poprawne wartości szczególnie Q dla wszystkich brakujących stacji, gdyż komunikat ten już się nie powtórzy. Wartości dolnej i górnej częstotliwości trzeba praktycznie dobierać za każdym razem, tak aby krzywe Brune'a pasowały do widma sygnałów, natomiast Q raz ustalone nie powinno być zmieniane.

Na dole okna (8) znajdują się dodatkowe przyciski sterujące wyświetlaniem sejsmogramu i spektrogramów, oraz dodatkowe informacje o regionie i odległości od epicentrum i przycisk zamykający okno.