Analiza zmian w czasie dwuwymiarowego rozkładu prawdopodobieństwa położenia epicentrum

w czasie dwuwymiarowego rozkładu prawdopodobieństwa położenia epicentrum wstrząsów

Charakterystyki probabilistyczne badanej próby statystycznej najczęściej określa się za pomocą me-tod parametrycznych (SILVERMAN, 1986). Zakłada się, że próba statystyczna ma znany rozkład okre-ślony funkcją analityczną, dla której znajduje się parametry najlepiej dostosowujące model do próby.

Metody parametryczne sprawdzają się w sytua-cjach, w których mamy dobrze znany model roz-kładu zmiennej. Bardziej skomplikowane rozkłady, np. rozkłady wielomodalne, są trudne do analizy;

w tym przypadku metody parametryczne nie dają pożądanych rezultatów (KAY, 1993). W takich sy-tuacjach można zastosować jądrowe estymatory nieparametryczne (równanie 4.3). Aby można było zastosować funkcję jako jądro, musi ona spełniać następujący warunek:

K x dx( ) =

-¥

ò

Zakłada się także, że funkcja jądrowa jest syme-tryczna względem zera i ma w tym punkcie maksi-mum.

Jakość estymacji oceniana jest na podstawie błędu średniokwadratowego (MSE), zdefiniowane-go jako wartość oczekiwana kwadratu błędu esty-macji w punkcie x:

MSE_x=Eéf x -f x gdzie f^{^} jest estymatorem funkcji rzeczywistej gę-stości prawdopodobieństwa f^{^}. Aby uzyskać całko-wity wskaźnik jakości estymacji, należy scałkować MSE_x po całej przestrzeni zmiennej losowej. War-tość takiej całki określana jest jako scałkowany błąd średniokwadratowy:

Korzystając z wyrażenia (4.5), można wykazać, że najbardziej efektywne jest jądro Epanecznikowa (SILVERMAN, 1986). Inne jądra są również stosowa-ne, mimo że mają mniejszą efektywność. Spadek efektywności związany z zastosowanym jądrem w stosunku do jądra Epanecznikowa zwykle nie jest duży, np. w przypadku najprostszego jądra jed-nostajnego efektywność jest tylko o 7% mniejsza, a dla najczęściej stosowanego jądra normalnego —

o 5%. Porównywalne wyniki (w stosunku do jądra Epanecznikowa) uzyskuje się zatem przez zwięk-szenie liczebności próby, w przypadku wymienio-nych jąder stosownie o 7% i 5%. Jak widać, zwięk-szenie liczby prób nie jest bardzo duże, zatem wybierając jądro, należy się raczej kierować właściwościami funkcji jądrowej, właściwościami uzyskiwanego estymatora i np. oczekiwaną wydaj-nością obliczeń niż efektywwydaj-nością. Najczęściej sto-sowane jest jądro normalne z uwagi na właściwo-ści funkcji e^x.

Rozkład epicentrów wstrząsów jest dwuwymia-rowy i taka też musi być funkcja gęstości prawdo-podobieństwa. Próba, która jest katalogiem sej-smicznym epicentrów wstrząsów, opisuje rozkład wstrząsów [X_i Y_i], gdzie X i Y są współrzędnymi epicentrów (MIREK, MIREK, 2007). Dwuwymiaro-we, radialne jądro normalne w przestrzeni dwuwy-miarowej przyjmie postać:

a estymator jądrowy gęstości prawdopodobieństwa w takiej przestrzeni jest zdefiniowany przez:

f x y

Jądro K jest nazywane jądrem radialnym, gdyż jego wartość jest stała dla stałej odległości od punktu próby losowej, dla którego wyznaczamy wartość funkcji f^{^}. Poziomice jądra (4.6) tworzą na płaszczyźnie okręgi. Można zwiększyć efektyw-ność jądra przez dostosowanie jego kształtu do kształtu rozkładu funkcji prawdopodobieństwa, któ-rej kształt może być skorelowany np. z kształtem pola eksploatacji lub lokalną tektoniką. Do modyfi-kacji kształtu funkcji jądra (4.6) stosuje się trans-formację liniową:

gdzie R jest macierzą kowariancji zmiennych X i Y.

Dwuwymiarowe jądro normalne z zastosowaniem transformacji liniowej przyjmuje postać:

K x y

Macierz R na diagonalnej ma wartości wariancji zmiennych X i Y. W przypadku zmiennych

nieza-4.5. Badanie kierunków migracji ognisk wstrząsów „tektonicznych” w celu wyjaśnienia uwarunkowań... 67

9*

leżnych macierz R przyjmie poza diagonalną war-tości zerowe. Jądro (4.6) jest szczególnym przy-padkiem jądra z transformacją liniową, jeśli za ma-cierz R podstawimy mama-cierz jednostkową.

Bardzo ważny jest parametr wygładzania h. De-cyduje on o kształcie funkcji jądrowej. Duża jego wartość spowoduje rozciągnięcie i spłaszczenie funkcji jądrowej, natomiast mała wartość wysmukla tę funkcję. Wartość parametru h, a co za tym idzie

— kształt funkcji jądrowej wpływa na jakość esty-matora jądrowego. Parametr ten jest dobierany tak, aby zminimalizować scałkowany błąd średniokwa-dratowy MISE (4.5) (GROENBORN, WALLNER, 1992).

Problem polega jednak na tym, że nie znamy rze-czywistej funkcji gęstości prawdopodobieństwa f.

Można zatem założyć a priori kształt estymowanej funkcji i na tej podstawie minimalizować wyraże-nie (4.5). Jeśli założymy, że rozkład w analizowa-nym obszarze jest normalny, to można wykazać, że optymalna wartość parametru h wynosi:

h n X Y gdzie n jest liczebnością próby. Należy pamiętać, że wartość ta jest optymalna tylko dla przypadku, w którym zmienne losowe X i Y mają rozkład nor-malny. Taka sytuacja może mieć miejsce np. w wy-dzielonych strefach czasoprzestrzennych, jeśli taki warunek był brany pod uwagę w procesie wydziela-nia stref. W przypadku gdy analizowany obszar jest większy, a funkcja gęstości prawdopodobieństwa ma wiele lokalnych maksimów, wartość (4.10) moż-na traktować jako pierwsze przybliżenie, którego wielkość może być szybko wyznaczona. Współ-czynnik h można dokładniej określić za pomocą metody uwiarygodnienia krzyżowego (cross-valida-tion), która także jest oparta na minimalizacji scałkowanego błędu średniokwadratowego MISE.

Tak wyznaczony parametr h może być stosowa-ny w całej przestrzeni jako wartość stała. Badastosowa-ny obszar cechuje się jednak nierównomiernym roz-kładem, co oznacza, że są miejsca o większym skupieniu danych i miejsca, w których danych tych jest mniej. Nasuwa się wniosek, że w miejscach o większym skupieniu można zmniejszyć wartość parametru wygładzania, aby zwiększyć rozdziel-czość metody i znaleźć ewentualne lokalne zmiany w funkcji gęstości prawdopodobieństwa. Zwiększe-nie parametru wygładzania w obszarach o mZwiększe-niej- mniej-szej liczbie danych spowoduje większe uogólnienie funkcji w tym obszarze, a także podniesie jakość estymacji. Można zatem modyfikować wartość pa-rametru wygładzania, biorąc pod uwagę liczbę da-nych w badanym obszarze. Ponieważ wartość funk-cji gęstości prawdopodobieństwa zależny od liczby

danych w badanym rejonie, estymatora jądrowego f^{^} wyznaczonego dla stałego parametru h można użyć jako pierwszego etapu obliczeń, a następnie na jego podstawie modyfikować wartość parametru h i ponownie wyznaczyć wartość f^{^}. Wartość współczynnika modyfikującego s można określić jako:

gdzie c jest dodatnią stałą, a ~s — średnią geome-tryczną wartości estymowanych dla stałego h funk-cji f^{^} w punktach określonych przez dane wejścio-we:

W ostatecznej formie estymator funkcji gęstości prawdopodobieństwa określa wyrażenie:

Stała c we wzorze (4.11) decyduje o intensyw-ności modyfikacji i zawiera się w przedziale [0; 1].

Może być wyznaczana metodą uwiarygodnienia krzyżowego. Często w zastosowaniach aplikacyj-nych wartość tego parametru ustala się na 0,5.

Opierając się na metodzie estymacji nieparame-trycznej funkcji gęstości prawdopodobieństwa, ob-liczono rozkłady na terenie GZW. Na rys.

4.19—4.36 zaprezentowano zmiany wartości funk-cji gęstości prawdopodobieństwa obserwowane w latach 1987—2004. Do analizy brano dane w oknach czasowych o długości jednego roku i energii³ 1e5 J. Analiza wykazała wyraźne zmia-ny funkcji gęstości prawdopodobieństwa wystąpie-nia wstrząsu w czasie, a tym samym czasową zmia-nę aktywności sejsmicznej w tym rejonie.

Literatura

GROENBORN P., WELLNER J.A., 1992: Information Bounds and Nonparametric Maximum Likehood Estimation. Berlin, Birkhäuser Verlag.

IDZIAKA., LASOCKIS., 1997: Badania struktury seryjnej indu-kowanych wstrząsów sejsmicznych z obszaru Górnośląskie-go Zagłębia WęgloweGórnośląskie-go. In: „Results from recent study in seismology and engineering geophysics”. Proc. Regional Conference with International Participation, Ostrava, 8—9

68 4. Sejsmiczność obszaru Górnośląskiego Zagłębia Węglowego

4.5. Badanie kierunków migracji ognisk wstrząsów „tektonicznych” w celu wyjaśnienia uwarunkowań... 69

Rys. 4.19. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1987 r.

Rys. 4.20. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1988 r.

70 4. Sejsmiczność obszaru Górnośląskiego Zagłębia Węglowego

Rys. 4.21. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1989 r.

Rys. 4.22. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1990 r.

4.5. Badanie kierunków migracji ognisk wstrząsów „tektonicznych” w celu wyjaśnienia uwarunkowań... 71

Rys. 4.23. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1991 r.

Rys. 4.24. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1992 r.

72 4. Sejsmiczność obszaru Górnośląskiego Zagłębia Węglowego

Rys. 4.25. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1993 r.

Rys. 4.26. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1994 r.

4.5. Badanie kierunków migracji ognisk wstrząsów „tektonicznych” w celu wyjaśnienia uwarunkowań... 73

Rys. 4.27. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1995 r.

Rys. 4.28. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1996 r.

10 — Geneza...

74 4. Sejsmiczność obszaru Górnośląskiego Zagłębia Węglowego

Rys. 4.29. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1997 r.

Rys. 4.30. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1998 r.

4.5. Badanie kierunków migracji ognisk wstrząsów „tektonicznych” w celu wyjaśnienia uwarunkowań... 75

Rys. 4.31. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 1999 r.

Rys. 4.32. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 2000 r.

10*

76 4. Sejsmiczność obszaru Górnośląskiego Zagłębia Węglowego

Rys. 4.33. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 2001 r.

Rys. 4.34. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 2002 r.

April 1997. Ed. Z. KALAB. Academy of Sciences of Czech Republic, Institute of Geonics, s. 151—158.

KAY S.M., 1993: Fundamentals of Statistical Signal Pro-cessing: Estimation Theory. New Jersey, PTR Prentice Hall.

KUSTOWSKIB., 2000: Time changes in dominant trends of epi-centre migration for local seismicity from a mine. Acta Geophys. Polon., 48, s. 455—464.

LASOCKIS., IDZIAKA., 1998: Dominant directions of epicenter distribution of regional mining-induced seismicity series in Upper Silesian Coal Basin in Poland. Pure Appl. Geo-phys., 153, s. 21—40.

LASOCKIS., KUSTOWSKI B., 1998: Hazard sejsmiczny i domi-nujące kierunki migracji ognisk podczas eksploatacji ścia-nowej: Przykłady z KWK „Wujek” (Seismic hazard and dominant migration directions of seismic sources in mining stopes: Examples from Wujek coal mine). Pr. GIG. Seria:

Konferencje: „Tąpania ’98: Bezpieczne prowadzenie robót górniczych”, s. 81—89.

LASOCKI S., KUSTOWSKI B., 2002: Analysis of deflection for identification of epicenter migration directions of swarm activity from Hida Mountains in Japan. Acta Geophys. Po-lon., 50, s. 23—34.

4.5. Badanie kierunków migracji ognisk wstrząsów „tektonicznych” w celu wyjaśnienia uwarunkowań... 77

Rys. 4.35. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 2003 r.

Rys. 4.36. Rozkład gęstości prawdopodobieństwa wystąpienia wstrząsu w 2004 r.

LASOCKIS., WĘGLARCZYKS., GIBOWICZS.J., 1997: A new met-hod to estimate directional character of mining-induced se-ismicity: application to the data from Wujek coal mine, Poland. Rockbursts and Seismicity in Mines. In: Proc. 4th Int. Symp. on Rockbursts and Seismicity in Mines, Kraków, Poland, 11—14 Aug. 1997. Eds. S.J. GIBOWICZ, S. LASOCKI. Rotterdam, Balkema, s. 207—211.

MIREKJ., MIREKK., 2007: Nieparametryczna estymacja funk-cji gęstości prawdopodobieństwa położenia źródła sej-smicznego. Warsztaty Górnicze 2007 z cyklu „Zagrożenia naturalne w górnictwie”, Ślesin 2007, s. 331—337.

NADARAYAE.A., 1989: Nonparametric Estimation of Probabi-lity Densities and Regresion Curves. Dordrecht, Kluwer.

SILVERMAN B.W., 1986: Density Estimation for Statistic and Data Analysis. London, Chapman and Hall.

Janusz Mirek, Beata Orlecka-Sikora, Stanisław Lasocki

4.6. Estymacja górnego ograniczenia