___________________________________________________________________________
Nowe możliwości obrazowania emisji sejsmicznej
indukowanej podziemnymi pracami górniczymi
dla kopalń LGOM
Andrzej Leśniak1), Anna Pięta1), Kamil Szostek1), Eugeniusz Koziarz2) 1)
AGH Akademia Górniczo-Hutnicza, lesniak@agh.edu.pl 2)
KGHM Polska Miedź S.A. Oddział Zakłady Górnicze „Rudna”, eugeniusz.koziarz@kghm.com
Streszczenie
Prezentowane rezultaty dotyczą nowoczesnych sposobów obrazowania źródeł emisji sejs-micznej indukowanej podziemnymi robotami górniczymi. Badania dotyczą rozwiązań, które mogą być osiągnięte dla poprawy dokładności lokalizacji oraz automatycznej oceny błędu lokalizacji dla standardowych konfiguracji sieci pomiarowej, typowej dla kopalń miedzi LGOM oraz sieci czujników jednoskładowych, wzbogaconej opcjonalnie o czujniki trójskładowe. W artykule przedyskutowano sposób wizualizacji struktury prędkościowej ośrodka geologicz-nego, wpływającej na efektywność klasycznej procedury lokalizacji dla różnych typów sieci pomiarowej oraz dla różnych rozmieszczeń czujników, w oparciu o zaimplementowane w kopalni „Rudna” programy „Litos-Wstrząsy” i zintegrowany z nim „Lokalizator 3C”. W dys-kusji uwzględniono model jednorodny i izotropowy oraz model niejednorodny (warstwowy) o budowie odzwierciedlającej budowę ośrodka geologicznego w rejonie oddziałów wydobyw-czych ZG „Rudna”. Z uwagi na typową dla polskich kopalń rud miedzi płaską konfigurację sieci pomiarowej oraz fakt, że dla klasycznych metod lokalizacji źródeł emisji występuje znaczny błąd określania głębokości źródeł wstrząsów sejsmicznych, zaprezentowano sposób wizualizacji dostosowany do tego typu warunków. W dalszej części przedstawiono nowocze-sną metodę „czasu odwróconego”, potencjalnie pozwalającą na znaczne poprawienie do-kładności lokalizacji, nawet w skrajnie trudnych przypadkach sieci płaskiej, złożonej z czujni-ków jednoskładowych i nakładających się w czasie źródeł emisji sejsmicznej. Metoda ta po-zwala również na oszacowanie wielkości błędów lokalizacji źródeł emisji wraz z przybliżoną oceną wieloznaczności sejsmicznego zagadnienia odwrotnego. Dyskutowany sposób obra-zowania zjawiska emisji sejsmicznej został dopasowany do charakteru metody i poparty przy-kładami dla danych syntetycznych.
Słowa kluczowe: źródła sejsmiczne, sejsmiczność indukowana, wizualizacja
Modern methods of visualization of seismic emission induced
by underground mining works for LGOM mines
Abstract
The modern methods of visualization of seismic emission process induced by underground mining works are presented in the article. The solutions used in mining practice for quality enhancement of seismic source location and automatic evaluation of errors of seismic source location are discussed. In the analysis not only the standard configuration of seismic sensors network used in copper mines of LGOM is analyzed but also the standard network supported by additional 3C geophones. The presented results were obtained using two
inte-grated computer programs Litos-Wstrząsy and Lokalizator 3C, implemented in “Rudna” mine. The detailed discussion of velocity model of geological medium and its influence on effectiveness of standard source location method for different configurations and types of seismic networks is presented. The two standard velocity models are discussed – the homo-geneous and heterohomo-geneous (stratified) ones. Their geometrical and physical parameters of the velocity structure are typical for one of mining panel from “Rudna” copper mine. The most important consequence of flat seismic network used in that mine is relatively large error of z-component of localized seismic sources when the standard location methods are used. In the next part of the paper the modern “time reversal” method is presented. It allows improvement of seismic source location quality even for flat networks and superimposed seismic events. The methods allow also the source location error estimation together with the ambiguity of the source location inverse problem. The discussed results are illustrated with the synthetic examples.
Key words: seismic sources, induced seismicity, visualization
Wstęp
Współczesny rozwój technik informatycznych przyczynia się do zwiększenia możli-wości wizualizacji złożonych procesów, mających miejsce w kopalniach podziem-nych. Nie chodzi wyłącznie o wizualizację procesów udostępniania i eksploatacji złóż, lecz również o różne zjawiska towarzyszące pracom podziemnym. Jednym z istotnych i skomplikowanych zjawisk jest emisja sejsmiczna. Wizualizacja tego procesu dynamicznego jest z jednej strony sporym wyzwaniem dla projektantów systemów komputerowych, zaś z drugiej umożliwia szczegółowe ewidencjonowanie i analizowanie tego zjawiska w czasie bieżącym. Z uwagi na kluczowe znaczenie emisji w ocenie zagrożenia ludzi i sprzętu oraz możliwości jej wykorzystania do ana-lizy powiązanych z nią zjawisk (zmiany naprężeń, powstawanie stref spękań itd.) wszelkie metody, przyczyniające się do pełnego poznania i zrozumienia tych zja-wisk, są wysoce pożądane. Wizualizacja procesów wielokrotnie przyczynia się do intuicyjnego i głębokiego wejrzenia w naturę tych procesów i zwiększa efektywność prowadzonych prac górniczych.
Artykuł obejmuje wybrane i ważne z punktu widzenia zastosowań w praktyce górniczej zagadnienia wizualizacji emisji sejsmicznej. Przedstawione przykłady i zamieszczone rezultaty są dziełem autorów artykułu, a część z nich została zaim-plementowana w systemach komputerowych Litos-Wstrząsy i Lokalizator 3C, pracu-jących w kopalni „Rudna”.
1. Sposoby obrazowania zjawisk sejsmicznych
Podstawowym sposobem obrazowania zjawisk sejsmicznych, zarówno indukowa-nych, jak i tych naturalindukowa-nych, są zapisy sejsmogramów. Pozwalają one na uzyskanie istotnych informacji na temat samego zjawiska. Na rys. 1 przedstawiono przykład sejsmogramu, który prezentowany jest przez program Litos-Wstrząsy po zarejestro-waniu zjawiska sejsmicznego.
Rys. 1. Przykładowe sejsmogramy prezentowane w programie Litos-Wstrząsy Należy zauważyć, że w przypadku sieci czujników jednoskładowych, zapis sejsmo-gramu dany jest jedynie w postaci pojedynczego szeregu czasowego. Natomiast w przypadku sieci czujników trójskładowych każde zjawisko jest opisane trzema szeregami, odpowiednio dla osi X, Y oraz Z. Dzięki temu możliwa jest ich wizualiza-cja w postaci tzw. hodografów, tak jak na rys. 2. Przedstawiają one drganie czujnika w wybranych dwóch płaszczyznach XY, XZ lub YZ. Taka wizualizacja pozwala ana-lizować charakter zjawiska i badać jego polaryzację.
Rys. 2. Przykład reprezentacji sejsmogramu w postaci trzech hodografów w różnych płaszczyznach
Wraz z rozwojem grafiki komputerowej możliwe stało się prezentowanie zjawisk sejsmicznych w bardziej przystępny sposób. Za pomocą uniwersalnych i darmowych bibliotek graficznych, takich jak OpenGL [7], program Litos-Wstrząsy, działający w kopalni O/ZG „Rudna”, został wzbogacony o wyświetlanie hodografów trójwymia-rowych, które prezentują rzeczywisty ruch czujnika pomiarowego. Dodatek do pro-gramu, nazwany „Lokalizator 3C”, pozwala między innymi śledzić ruch czujnika w sposób interaktywny, a przez to pomaga w interpretacji zarejestrowanego zjawi-ska. Możliwe jest na przykład wyświetlanie tylko części zarejestrowanych sejsmo-gramów, a więc wyświetlanie fragmentu, który jest najbardziej charakterystyczny i najbardziej interesujący dla użytkownika. Istotnym atutem programu są jego niskie wymagania sprzętowe oraz możliwość działania niezależnie od wspomnianego pro-gramu Litos-Wstrząsy (dla celów wizualizacyjnych).
2. Model geologiczny w wizualizacji
W obrazowaniu emisji sejsmicznej nieodłącznym elementem jest prezentacja mode-lu geologicznego. W zależności od potrzeb i preferencji użytkownika oraz sposobu wykorzystywania modelu (np. w przypadku techniki śledzenia promienia), może być on skonstruowany w różny sposób. W programie Lokalizator 3C wybrano metodę, polegającą na konstrukcji granic warstw o różnej prędkości, ponieważ algorytmy używane do lokalizacji źródła zjawiska sejsmicznego opierają się właśnie na tej cesze ośrodka. Rys. 3 przedstawia stworzony i wykorzystywany model geologiczny kopalni O/ZG „Rudna”. Model składa się z granic, które reprezentują charaktery-styczne zmiany w prędkości propagacji fali P. Dla warstwy najniższej, poniżej grani-cy czerwonej, prędkość fali ustalono na 2250 m/s, następnej na 4750 m/s i kolejnej na 4700 m/s. Obszar nad granicą żółtą nie jest wykorzystywany w obliczeniach. Należy zaznaczyć, że rysunek przedstawia model z opcjonalnym przewyższeniem, co poprawia jego czytelność. W rzeczywistości, obszar kopalni o powierzchni prawie 100 km2 ma niewielkie różnice wysokości. Model taki użytkownik może oglądać z dowolnej strony, wykorzystując mysz do obrotu wokół wybranego punktu obser-wacji oraz do przesuwania wirtualnej przestrzeni.
Rys. 3. Model geologiczny dla obszaru LGOM z dziesięciokrotnym przewyższeniem skali pionowej. Przedstawione granice odpowiadają granicom warstw o różnej prędkości
propagacji fali P. Kolory dobrano arbitralnie
Istotnym dodatkiem do wyświetlanego modelu jest umieszczenie w oknie wizua-lizacyjnym siatki kilometrowej wraz ze współrzędnymi. Dodatkowo wizualizację można wzbogacić o mapę kopalni dowolnego typu, wyświetlaną np. pod modelem geologicznym. Może być ona widoczna w tle modelu dzięki odpowiedniemu zasto-sowaniu przezroczystości. Poza mapą, możliwe jest również prezentowanie konturów obszarów kopalni (rys. 4). Pozwala to w jeszcze wygodniejszy sposób na obserwo-wanie zarówno modelu geologicznego w korelacji z planem kopalni, jak i lokalizację
zjawisk sejsmicznych. Ponadto, niejednokrotnie konieczne jest zwiększenie czytel-ności czy konieczność zbliżenia do wybranego obszaru. Dzięki interaktywczytel-ności wi-zualizacji, nie stanowi to problemu i użytkownik może w intuicyjny sposób wybrać interesujący go obszar do szczegółowej analizy. W aplikacji zaimplementowano również skróty klawiszowe, które pozwalają na przedstawienie modelu w jednym ze standardowych widoków, takich jak rzut ortogonalny lub perspektywiczny, z dołu lub z dowolnej strony modelu. Kolejnym usprawnieniem jest możliwość wycięcia i wy-świetlenie tylko tego fragmentu modelu, który jest najbardziej istotny dla prezentacji zjawiska sejsmicznego i tras sejsmicznych (rys. 5). Funkcjonalności programu ob-sługiwane są za pomocą licznych skrótów klawiszowych, których listę można w do-wolnej chwili wyświetlić na ekranie.
Rys. 4. Przykład prezentowania konturów obszarów kopalni oraz mapy kopalni w oknie wizualizacji
Rys. 5. Przykład detali, które zwiększają czytelność interaktywnej wizualizacji: model został ograniczony do wąskiego pasa. Dodatkowo widoczne są szczegóły związane z hipocentrum
3. Sposób prezentacji sieci pomiarowej i zjawisk sejsmicznych
W wizualizacji, poza wspomnianymi elementami modelu geologicznego oraz kopalni (w tym jej mapy), możliwe i pożądane jest wyświetlanie lokalizacji czujników sej-smicznych, zarówno tych jedno- jak i trójskładowych. W programie wizualizacyjnym Lokalizator 3C czujniki przedstawiane są za pomocą stożków o wierzchołkach skie-rowanych w dół. Kolory tych stożków mogą odzwierciedlać wybraną cechę czujnika. Przykładowo, w programie kolor błękitny reprezentuje czujniki jednoskładowe, na-tomiast różowy – trójskładowe. W pobliżu czujników wyświetlane są również ich nazwy (rys. 6).
W związku z głównym zadaniem programu Lokalizator 3C, którym jest naniesie-nie poprawki na wartość głębokości zlokalizowanego źródła zjawiska sejsmicznego (korekta składowej „z” hipocentrum obliczonego programem Litos-Wstrząsy), możli-we jest również wyświetlenie tras sejsmicznych od źródła zjawisk do wybranych czujników trójskładowych. Pojawiają się one jako czerwone krzywe łamane, ponie-waż na granicy warstw o różnych prędkościach promienie sejsmiczne zmieniają kierunek. Dodatkowo, nowo określone hipocentrum jest przedstawiane w postaci białej kuli. Błąd lokalizacji również jest rysowany, w tym wypadku za pomocą prze-zroczystego walca (rys. 6), którego odpowiednie wymiary prezentują błąd lokalizacji, również liczony w programie.
Rys. 6. Wizualizacja wyświetla również trasy sejsmiczne, czujniki oraz sejsmogram dla wybranego czujnika trójskładowego
4. Obrazowanie z użyciem metody czasu odwróconego
Modelowanie za pomocą metody czasu odwróconego (ang. time reversal imaging) jest nową, ale już często wykorzystywaną metodą do obrazowania źródeł wstrząsów sejsmicznych. Jej podstawowe założenia oraz praktyczne wykorzystanie w zagad-nieniach związanych z medycyną czy naukami o Ziemi zostało przedstawione w Finka [1]. Metoda wykorzystuje symetrię równania falowego względem czasu. Praktyczna realizacja opiera się na odwróceniu sygnału sejsmicznego zarejestrowa-nego przez odbiorniki względem czasu i prowadzenia modelowania, tak jakby sy-gnał ten był sysy-gnałem emitowanym ze źródła, zlokalizowanego w miejscu położenia odbiorników sejsmicznych (rys. 7). Podczas realizowanego w dalszym etapie mode-lowania numerycznego zarejestrowany sygnał propaguje przez ośrodek geologiczny, intensywnie interferując w obszarach odpowiadających ogniskom fali sejsmicznej.
Rys. 7. Schematyczne przedstawienie założeń metody time reversal imaging
W związku z faktem, że całość zarejestrowanego sygnału sejsmicznego jest uży-ta w modelowaniu odwrotnym, metoda TRI jest niezależna od błędu związanego z określaniem czasu przyjścia fali sejsmicznej, zarejestrowanym na sejsmogramie [2]. Co więcej, poprzez modelowanie wstecznej propagacji fali sejsmicznej mogą być zlokalizowane również zarejestrowane słabe i wzajemnie nakładające się zda-rzenia sejsmiczne.
Ograniczenia stosowania metody związane z niedokładnością określenia modelu prędkościowego ośrodka oraz z rozmieszczeniem czujników, w których został zare-jestrowany sygnał sejsmiczny, zostały opisane w [6] oraz [5].
Poniżej przedstawiono praktyczny przykład, wykonany dla danych syntetycz-nych, ilustrujący możliwości zastosowania metody TRI do lokalizacji źródeł wstrzą-sów sejsmicznych. Przykład pokazuje zastosowanie metody TRI dla przypadku, gdy zarejestrowany przez czujniki sygnał jest złożeniem fal sejsmicznych, wygenerowa-nych przez kilka różwygenerowa-nych ognisk sejsmiczwygenerowa-nych.
Lokalizacja położenia źródła sygnału sejsmicznego została przeprowadzona dla modelu trójwarstwowego, ze zmiennym nachyleniem warstw. Modelowania propa-gacji fali ciśnieniowej realizowano dla ośrodka dwuwymiarowego. Obliczenia prowa-dzono z zastosowaniem metody różnic skończonych z wykorzystaniem schematu obliczeniowego opartego o siatki przesunięte w czasie i przestrzeni [4].
Założono istnienie pięciu ognisk fali sejsmicznej, oraz sieci sejsmologicznej zło-żonej z siedmiu czujników sejsmicznych. W pierwszym etapie modelowań przepro-wadzono modelowania numeryczne symulujące rejestracje zjawisk sejsmicznych przez siedem założonych stacji. Założono punktowy, eksplozywny mechanizm źró-dła. Symulacje zdarzenia sejsmicznego przeprowadzono, przykładając we wcześniej ustalonych punktach siatki obliczeniowej sygnał Rickera. Symulowano zdarzenia sejsmiczne generujące fale, których widmo przesunięte jest w stronę wyższych czę-stotliwości, pozwalające na uzyskanie wysokiej rozdzielczości poziomej i pionowej pomiarów. W opisywanym przykładzie częstotliwość dominująca sygnału Rickera zawarta była w przedziale od 130-170 Hz. W celu uniknięcia błędów numerycznych, sygnał źródłowy, po wyemitowaniu 99% energii, był wyciszany odpowiednio dobraną procedurą numeryczną. Modelowanie przeprowadzono dla ośrodka geologicznego o rozmiarach 1700 na 2500 metrów. Ilość punktów budujących siatkę obliczeniową dobrano w ten sposób, aby spełniała warunki stabilności rozwiązania numeryczne-go. Propagację fali akustycznej dla wybranych chwil czasowych przedstawiono na rys. 8. Na rys. 8a widoczne są zmiany ciśnienia, obserwowane dla pięciu zdarzeń sejsmicznych, wygenerowanych w tym samym czasie. Rys. 8b przedstawia kolejny etap modelowań, w którym dochodzi do interferencji fal sejsmicznych, ich odbicia i załamania. Rys. 8c przedstawia wymodelowany rozkład ciśnienia po odbiciu fali od powierzchni terenu.
Rys. 8. Przedstawienie pierwszego etapu modelowania – propagacji fali sejsmicznej z pięciu źródeł, w których zjawisko zaistniało jednocześnie
W drugiej części obliczeń przeprowadzono powtórne modelowanie propagacji fali sejsmicznej, w którym odwrócone zapisy sejsmiczne imitowały sygnał źródłowy. Na rys. 9 przedstawiono wizualizacje fazy modelowania, przeprowadzonej w celu okre-ślenia obszarów, gdzie następuje zogniskowanie emitowanej energii sejsmicznej. Przedstawione na rys. 9 zmiany ciśnienia modelowane były dla określonego w pierwszej fazie modelowania dwuwymiarowego ośrodka geologicznego.
Rys. 9. Wizualizacja modelowania wstecznego, zarejestrowanego przez odbiorniki sygnału sejsmicznego w głąb modelu na podstawie wartości ciśnienia. Na rysunku b) widać wyraźnie
zaznaczoną lokalizację obszarów źródeł z pierwszego etapu modelowania. Obszary te za-znaczają się tylko dla chwili czasowej, w której doszło do wyemitowania energii ze źródła
– rys. c)
Analiza rozkładu wartości ciśnienia, pozwala na wydzielenie, w pewnej chwili czasowej obszarów, w których obserwowana jest podwyższona wartość ciśnienia. Obszary te zaznaczają się na rys. 9b jako pięć plam barwnych o podwyższonej in-tensywności. Wydzielone obszary odpowiadają położeniom źródeł sejsmicznych, założonych w pierwszym etapie modelowań. W trakcie kolejnych kroków oblicze-niowych obszary o podwyższonej wartości ciśnienia ulegają rozproszeniu (rys. 9c). Z punktu widzenia określenia lokalizacji źródeł wstrząsu sejsmicznego istotny jest wybór parametru zapewniającego optymalny sposób wizualizacji obszarów, odpo-wiadających ogniskom fali sejsmicznej [1].
Wizualizację drugiej fazy modelowania może być przeprowadzona w oparciu o zaproponowany w tym celu współczynnik PAPR (ang. peak-to-average power
ratio), obliczany na podstawie maksymalnej wartości ciśnienia uzyskanego podczas
obliczeń oraz wartości średniej ciśnienia. Obrazowania wykorzystujące współczyn-nik PAPR pozwalają na zachowanie obszarów o wzmożonych koncentracjach ci-śnienia, odpowiadające lokalizacjom źródeł fali sejsmicznej z pierwszego etapu modelowań (rys. 10). Jeśli porównamy te wyniki z wizualizacją tradycyjną (rys. 9), łatwo zauważyć znaczny wzrost rozdzielczości obrazu, co umożliwia bardziej wiary-godną lokalizację źródeł emisji i ocenę błędu ich lokalizacji.
Rys. 10 Wizualizacja ogniskowania źródeł fali sejsmicznej na podstawie wartości parametru PAPR. Mapy wartości parametru PAPR obliczono dla chwil czasowych, przedstawionych
Podsumowanie – inne możliwości wizualizacji zjawisk sejsmicznych
W artykule zaprezentowano wybrane możliwości wizualizacji procesu emisji sej-smicznej, począwszy od graficznej prezentacji niejednorodnego ośrodka geologicz-nego, poprzez zobrazowanie sejsmicznej sieci pomiarowej i zarejestrowanych zja-wisk, aż po wizualizację procedury lokalizacji i graficznego zobrazowania błędów lokalizacji źródeł emisji. Zaprezentowane wyniki zostały wygenerowane w oparciu o autorski program Lokalizator 3C, sprzęgnięty ze używanym w kopalni „Rudna” programem Litos-Wstrząsy. Jako uzupełnienie możliwości tych programów przed-stawiono rezultaty obrazowania procesu propagacji fal sejsmicznych z użyciem no-wej metody czasu odwróconego. Syntetyczne wyniki zaprezentowane w artykule wskazują na duży potencjał metody, szczególnie w wypadku interpretacji słabych sygnałów sejsmicznych o charakterze interferencyjnym.
Wizualizacje zjawisk sejsmicznych mogą zostać wzbogacone o inne, ważne elementy. W większości przypadków, przedstawiane elementy wymagają jednak, aby odpowiednie dane zostały przeliczone za pomocą innych algorytmów. Niemniej w miarę rozwoju technologii coraz więcej obliczeń może być wykonywanych w cza-sie rzeczywistym, np. z wykorzystaniem procesorów kart graficznych. Pozwala to na interaktywną prezentację wyników nawet bardzo złożonych obliczeń, dynamicznie zmieniając poszczególne parametry.
Przykładem ciekawego rozwinięcia programu wizualizacyjnego jest zaimplemen-towanie wyświetlania mechaniki wstrząsu oraz wizualizacji dynamicznego spadku naprężeń. Jest to zagadnienie niezwykle istotne dla poprawnej interpretacji zjawi-ska, a atrakcyjne i przejrzyste przedstawienie go na wizualizacji podnosi znacząco jej wartość.
Badania zostały sfinansowane przez KGHM Polska Miedź S.A. w ramach zlecenia KGHM-KR-U-0191-2014 oraz częściowo z badań statutowych KGIS, GGiOŚ, AGH w ramach zada-nia 11.11.140.613.
Bibliografia
[1] Fink M., 1999, Time-reversed acoustics, Scientific American; 555:67-73.
[2] Gajewski D., Tessmer E., 2005, Reverse modelling for seismic event characterization: Geophysical Journal International, 163, 276-284.
[3] Kremers S., Fichtner A., Brietzke G., Igel H., Larmat C., Huang L., Käser M., 2011, Ex-ploring the potentials and limitations of the time-reversal imaging of finite seismic sources: Solid Earth.
[4] Levander A., 1988: Fourth-order finite-difference P-SV seismograms, Geophysics, 53-11, 1425-1436.
[5] Saenger E.H., 2011, Time reverse characterization of sources in heterogeneous media, NDT & E International 44, 751-759.
[6] Steiner, B., Saenger E.H., Schmalholz S.M., 2011, Time reverse imaging with limited s-wave velocity model information, Geophysics 76 (5), MA33-MA40.
[7] Wright, R.S. Jr., Lipchak B., Haemel N., 2007, OpenGL Superbible: comprehensive tuto-rial and reference, Fourth Edition.
