Index of /rozprawy2/10619

Pełen tekst

(1)AKADEMIA GÓRNICZO-HUTNICZA IM. STANISŁAWA STASZICA W KRAKOWIE. WYDZIAŁ GEOLOGII, GEOFIZYKI I OCHRONY ŚRODOWISKA KATEDRA GEOFIZYKI. ROZPRAWA DOKTORSKA ZWIĘKSZANIE MOŻLIWOŚCI WYZNACZANIA PARAMETRÓW ZBIORNIKOWYCH SKAŁ Z WYKORZYSTANIEM RENTGENOWSKIEJ MIKROTOMOGRAFII KOMPUTEROWEJ. mgr inż. Marek Dohnalik. PROMOTOR: prof. dr hab. inż. Jadwiga Jarzyna. KRAKÓW 2013.

(2) Pragnę złożyć szczególne podziękowania Pani prof. dr hab. inż. Jadwidze Jarzynie za wiele cennych rad i wskazówek podczas realizacji tej pracy. Dziękuję także za ogromną cierpliwość, wyrozumiałość oraz szczególną życzliwość.. Podziękowania kieruję również do Pani prof. dr hab. inż. Marii Ciechanowskiej za nieustanne mobilizowanie mnie do realizacji niniejszej pracy.. Bardzo dziękuję Pani mgr inż. Jadwidze Zalewskiej za umożliwienie mi uczestnictwa w wielu zagranicznych konferencjach, bez których nie zdobyłbym niezbędnej wiedzy i doświadczenia do realizacji tej pracy.. Dziękuję wszystkim obecnym i byłym pracownikom Zakładu Geofizyki Wiertniczej INiG za niesienie pomocy w zrozumieniu budowy skał oraz zjawisk w nich zachodzących. Dziękuję także mgr Konradowi Ziemianinowi z Zakładu Geologii i Geochemii INiG za wykonanie badań petrologicznych.. Arkadiuszowi Buniakowi dziękuję za pomysł na przedmiot badań niniejszej pracy.. Żonie dziękuję za wsparcie, wyrozumiałość i wiarę we mnie szczególnie podczas intensywnego okresu realizacji tej pracy..

(3) SPIS TREŚCI WSTĘP ....................................................................................................................................... 4 ROZDZIAŁ I – Opis rentgenowskiej mikrotomografii komputerowej ..................................... 6 I.1. Podstawy mikrotomografii rentgenowskiej ....................................................................... 6 I.1.1 Podstawowe własności promieniowania rentgenowskiego ........................................ 6 I.1.2. Oddziaływanie promieniowania rentgenowskiego z materią .................................... 8 I.1.3. Budowa rentgenowskiego mikrotomografu komputerowego .................................. 11 I.2. Otrzymywanie danych za pomocą mikrotomografu ........................................................ 14 I.2.1. Przebieg pomiaru ..................................................................................................... 15 I.3. Rekonstrukcja obrazów.................................................................................................... 18 I.3.1. Sztuczne elementy (artefakty) rekonstrukcji ........................................................... 21 ROZDZIAŁ II – Analiza przestrzennych obrazów struktury porów ....................................... 26 II.1. Przetwarzanie obrazu...................................................................................................... 26 II.2. Analiza liczbowa struktury porowej ............................................................................... 29 II.3. Wybór reprezentatywnego fragmentu próbki ................................................................. 42 II.3.1 Przykład zastosowania............................................................................................. 46 ROZDZIAŁ III – Opis obszaru badań ...................................................................................... 51 III.1. Obszar badań ................................................................................................................. 51 III.2. Charakterystyka obszaru badań ..................................................................................... 57 III.3. Zespoły facjalne górnego czerwonego spągowca ......................................................... 72 III.4. Charakterystyka petrograficzna piaskowców eolicznych ............................................. 74 III.4.1. Wnioski ................................................................................................................. 79 III.5. Wyniki wstępnych badań laboratoryjnych .................................................................... 80 ROZDZIAŁ IV- Analiza uzyskanych wyników ...................................................................... 93 IV.1. Analiza niepewności wyników pomiarów .................................................................... 93 IV.2. Wyniki badań laboratoryjnych ...................................................................................... 98 IV.3. Korelacja uzyskanych danych ..................................................................................... 125 IV.4. Zestawienie zależności w celu oceny możliwości określania predykcji parametrów zbiornikowych. ................................................................................................................ 133 PODSUMOWANIE ............................................................................................................... 148 WNIOSKI ............................................................................................................................... 150 SPIS TABEL .......................................................................................................................... 152 SPIS RYSUNKÓW ................................................................................................................ 154 LITERATURA ....................................................................................................................... 159. 3.

(4) WSTĘP Na terenie Polski większość złóż węglowodorów łatwo dostępnych oraz występujących w skałach charakteryzujących się dobrymi właściwościami zbiornikowymi została już odkryta. Od kilku lat prowadzone są poszukiwania gazu i ropy w formacjach, z których eksploatacja odbywa się w sposób niekonwencjonalny. Są to skały o słabych właściwościach zbiornikowych, gdzie węglowodory są uwięzione w porach o małych średnicach. Wydobywanie złóż w nowych warunkach wymusza także rozwój laboratoryjnych metod badawczych, w celu lepszego rozpoznania właściwości skał. Niniejsza praca obejmuje ocenę możliwości zastosowania nowej metody badawczej, jaką jest rentgenowska mikrotomografia komputerowa, do opisu struktury porowej skał. Na podstawie dostępnych opublikowanych opracowań, dotychczas przeprowadzonych własnych badań i doświadczenia autora rozprawy, przed realizacją pracy postawiono następujące tezy: 1.. Poszukiwania węglowodorów i wody w trudnych warunkach geologicznych wymuszają doskonalenie metod opisu przestrzeni porowej skał.. 2.. Współczynniki uzyskane na podstawie badania skał metodą mikrotomografii komputerowej umożliwiają sparametryzowanie przestrzeni porowej.. 3.. Opis. przestrzeni. porowej,. uzyskany. z. wykorzystaniem. mikrotomografii. komputerowej, jest bogatszy w porównaniu z wynikami z porozymetrii. 4.. Wyniki uzyskane metodą mikrotomografii komputerowej i NMR pozwalają na połączenie informacji o wielkości i krętości kanałów porowych z nasyceniem wodą swobodną i kapilarną.. Rentgenowska mikrotomografia komputerowa (micro-CT) jest stosunkowo nowym narzędziem służącym badaniu wewnętrznej struktury materiałów. Metoda ta opiera się na takich samych podstawach matematycznych jak tomografia komputerowa, szeroko stosowana w medycynie. Główną różnicą pomiędzy nimi, jest zastosowanie w mikrotomografii źródła rentgenowskiego, które emituje promieniowanie z bardzo małego obszaru anody, dzięki czemu możliwe jest zrekonstruowanie obrazu o wysokiej rozdzielczości. Ta różnica sprawia, że mikrotomografię można stosować do badania struktury porowej skał.. 4.

(5) Głównym celem niniejszej pracy było sprawdzenie możliwości wykorzystania nowej, w skali kraju metody, do badania struktury porowej skał. Do realizacji badań wykorzystano dostępny w Zakładzie Geofizyki Wiertniczej INiG w Krakowie, mikrotomograf Benchtop CT160Xi. Do binaryzacji i wizualizacji obrazów zastosowano oprogramowanie ImageJ oraz AVIZO. Zdecydowaną większość parametrów uzyskanych z analizy obrazu obliczono wykorzystując oprogramowanie MAVI. Zastosowane podczas realizacji pracy urządzenie pozwoliło na uzyskanie obrazów wewnętrznej struktury skał z rozdzielczością 5,8µm. Należy mieć jednak na uwadze, że pory występujące w skałach mają wielkość średnicy od kilku nanometrów do kilku milimetrów. Informacja ta sugeruje, że pewna, a w wielu przypadkach znaczna, część struktury porowej nie mogła zostać uwzględniona w obrazowaniu mikrotomograficznym oraz w późniejszej analizie. W celu weryfikacji danych otrzymanych z przetwarzania obrazu, uzyskane wartości odniesiono do innych parametrów, otrzymanych z wykorzystaniem sprawdzonych metod, jak porozymetria rtęciowa i jądrowy rezonans magnetyczny (NMR). Badania wykonane podczas realizacji pracy zostały przeprowadzone na próbkach piaskowca czerwonego spągowca facji rdzenia wydmy. Próbki pochodziły z trzech rejonów: Czarna Wieś – Parzęczewo, Siekierki – Miłosław i Środa Wielkopolska – Kromolice. Z badań petrograficznych wynikało, że próbki różniły się pomiędzy sobą przede wszystkim stopniem diagenezy. Następnie za pomocą badań laboratoryjnych stwierdzono występowanie znacznych różnic we właściwościach zbiornikowych, pomiędzy próbkami pochodzącymi z poszczególnych rejonów. Kolejnym celem pracy było sprawdzenie, czy parametry uzyskane z analizy obrazu także wykażą zróżnicowanie wartości pomiędzy próbkami pobranymi z różnych rejonów. Wszystkie otrzymane wyniki zostały przedstawione w formie tabelarycznej. Dla uzyskanych parametrów, wyznaczono macierze współczynników korelacji. Sprawdzono, w jakim stopniu wartości otrzymane z analizy obrazu korelują z parametrami uzyskanymi ze sprawdzonych metod badawczych. Następnie wykonano ocenę możliwości predykcji parametrów zbiornikowych na podstawie danych otrzymanych z analizy obrazu. Postawione we wstępie tezy zostały udowodnione podczas realizacji badań opisanych w niniejszej pracy.. 5.

(6) ROZDZIAŁ I – Opis rentgenowskiej mikrotomografii komputerowej I.1. Podstawy mikrotomografii rentgenowskiej I.1.1 Podstawowe własności promieniowania rentgenowskiego Wytwarzanie promieniowania X Promieniowanie rentgenowskie zostało odkryte w 1895 roku przez niemieckiego uczonego Wilhelma Röentgena. Promieniowanie to zostało odkryte przypadkowo, podczas badań związanych z promieniami katodowymi. Roentgen zauważył, że nieznane dotąd promienie pozwalają na prześwietlanie niektórych przedmiotów. Ponieważ nie był w stanie wyjaśnić, w jaki sposób powstaje to promieniowanie, nadał mu nazwę promieniowania X, które często nazywane jest też promieniowaniem rentgenowskim [Czerwiński, 1998]. Odkrycie to zostało bardzo entuzjastycznie przyjęte przez świat medyczny, gdyż promienie te między innymi umożliwiły prześwietlanie ludzkiego ciała.. Rys. I.1. Widmo promieniowania elektromagnetycznego z zaznaczonym podziałem fal ze względu na długość [Cnudde, 2008]. 6.

(7) Promienie X są zlokalizowane na wysokoczęstotliwościowym końcu widma promieniowania elektromagnetycznego. Ich długość fali mieści się w zakresie od 10 pm do 10 nm. Zakres promieniowania rentgenowskiego znajduje się pomiędzy ultrafioletem i promieniowaniem gamma (rys. I.1). W praktyce promieniowanie przenikliwe uzyskuje się w lampie rentgenowskiej (rys. I.2), gdzie źródłem promieniowania jest anoda, w którą uderza i na której zostaje zatrzymana wiązka elektronów, przyspieszana za pomocą różnicy potencjałów rzędu kilku do kilkuset tysięcy woltów.. Rys. I.2. Schemat lampy rentgenowskiej: A – anoda, K – katoda, V – różnica potencjałów, X – foton promieni rentgenowskich. Fotony emitowane jako promieniowanie hamowania, mają ciągłe widmo energii od zera do wartości maksymalnej, równej energii kinetycznej elektronu przyspieszonego w polu między katodą a anodą. Zdecydowana większość promieniowania niskoenergetycznego absorbowana jest w materiale anody oraz w wewnętrznym i zewnętrznych filtrach lampy. W związku z tym widmo emitowanej wiązki promieniowania hamowania ma postać krzywej posiadającej maksimum dla pewnej określonej energii (rys. I.3) [Skrzyński, 2004].. Rys. I.3. Przykładowe widmo promieniowania emitowanego z lampy rentgenowskiej tomografu (lampa OPTI 155 zamontowana w tomografie Siemens Somatom DR) [Skrzyński, 2004]. 7.

(8) Wybicie elektronów z powłoki wewnętrznej atomu, przez elektron z wiązki bombardującej próbkę, powoduje emisję kwantu charakterystycznego promieniowania rentgenowskiego. Przebieg tego zjawiska przedstawiono schematycznie na rysunku I.4. Powstała w rezultacie jonizacji powłoki (np. K) atomu luka zostaje zapełniona w wyniku przejścia na opróżnione miejsce elektronu z powłoki bardziej oddalonej od jądra (np. LIII), a różnica energii między tymi dwoma poziomami energetycznymi atomu zostaje wypromieniowana w postaci kwantu promieniowania X (promieniowanie charakterystyczne) [Szummer i in., 1994]. Generowanie promieniowania rentgenowskiego jest procesem bardzo niewydajnym energetycznie. Podczas tego procesu jedynie niewielka część energii elektronów zamieniana jest na promieniowanie X, większość energii zostaje zużyta w postaci wydzielonego ciepła [Stock, 2009].. Rys. I.4. Schemat ilustrujący emisję kwantu charakterystycznego promieniowania X [Szummer i in., 1994]. I.1.2. Oddziaływanie promieniowania rentgenowskiego z materią Absorpcja promieniowania [Gostkowska i Zajdel, 1977] Podczas penetracji badanego obiektu promieniowanie ulega osłabieniu. Spowodowane jest to całkowitą absorpcją (efekt fotoelektryczny) niektórych fotonów oraz odbiciem i rozproszeniem części pozostałych fotonów wiązki. Fotony mogą ulegać rozproszeniu koherentnemu (efekt Rayleigh’a ) oraz niekoherentnemu (efekt Comptona).. 8.

(9) Efekt Rayleigh’a [Cnudde, 2005] Efekt ten zachodzi, gdy foton promieniowania X o bardzo niskiej energii oddziałuje na mocno związany elektron. Jeżeli elektron zostanie wprowadzony w stan wibracji o częstotliwości równej częstotliwości promieniowania, to wyemituje promieniowanie o tej samej długości fali w przypadkowym kierunku (rys. I.5), podczas powrotu do stanu spoczynku. Ponieważ proces ten jest sprężysty, praktycznie nie następuje tu utrata energii. Dla fotonów o energii mniejszej niż 50 keV, sprężyste (koherentne) rozpraszanie jest dominującym procesem.. Rys. I.5. Schemat rozpraszania sprężystego (koherentnego) [Cnudde, 2005]. Efekt fotoelektryczny [Gostkowska i Zajdel, 1977; Ketcham i Carlson, 2001] Zjawisko to polega na przekazaniu całkowitej energii fotonu elektronowi swobodnemu lub związanemu na orbicie atomu. Zachodzi najczęściej, gdy energia fotonu jest niewiele większa od energii wiązania elektronu.. Rys. I.6. Ilustracja efektu fotoelektrycznego [wikipedia: efekt fotoelektryczny]. 9.

(10) W rezultacie tego zjawiska znika foton, a elektron uzyskujący energię kinetyczną równą pierwotnej energii fotonu zostaje wyemitowany (rys. I.6). Zjawisko to jest przeważające przy badaniu próbek skał, gdyż dominuje w przedziale energii od 50-100 keV. Efekt Comptona [Gostkowska i Zajdel, 1977] Nazywamy tak nieelastyczne rozpraszanie fotonów na elektronach. W jego wyniku foton, który zderzy się z elektronem, wybija go i zmienia kierunek ruchu (rys. I.7). Energia fotonu maleje, natomiast elektron uzyskuje energię kinetyczną, równą różnicy energii fotonu przed i po zderzeniu. Zmniejszeniu energii fotonu odpowiada zmniejszenie częstości fali elektromagnetycznej. Strata energii fotonu zależy od kąta rozproszenia (odchylenia fotonu od pierwotnego kierunku) i jest największa przy rozproszeniu wstecznym, pod kątem 180o. Efekt ten dominuje dla promieniowania o energii z zakresu 5-10 MeV. W mikrotomografii rentgenowskiej fotony odbite, nie przekazują użytecznych informacji, natomiast. są. zazwyczaj. źródłem. szumu,. który. negatywnie. wpływa. na. jakość. rekonstruowanych obrazów.. Rys. I.7. Schemat rozpraszania Comptonowskiego (niekoherentnego) [Cnudde, 2005]. Dla tak opisanych zjawisk wiązka promieniowania jest osłabiana wg prawa Lamberta-Beer’a (I.1) [Stock, 2009]. Zmniejszenie natężenia promieniowania (-ΔI) w wyniku przejścia wiązki przez elementarną grubość absorbenta jest proporcjonalne do natężenia I oraz do grubości warstwy Δx:.  I    I  x. 10. (I.1).

(11) Dla absorbenta o grubości x, po scałkowaniu otrzymuje się wykładniczą postać prawa pochłaniania:. I  I 0 e  x. (I.2). gdzie: I – natężenie wiązki po przejściu przez ośrodek, Io – oznacza natężenie wiązki padającej, µ – całkowity współczynnik absorpcji (pochłaniania), zwany również współczynnikiem tłumienia liniowego, x – grubość absorbenta.. Detekcja promieniowania W. przypadku. systemów. mikrotomograficznych,. promieniowanie. przenikliwe. jest. rejestrowane przez detektor, którego funkcją jest konwersja promieniowania rentgenowskiego na sygnał elektryczny. Zazwyczaj detektor tworzy układ w postaci warstwy scyntylacyjnej i elementu fotopowielającego.. Scyntylator. absorbuje. promieniowanie. X. i. przetwarza. je. na. promieniowanie światła widzialnego, które zamieniane jest na sygnał elektryczny, proporcjonalny do natężenia fotonów. Uzyskany sygnał elektryczny jest następnie zapisywany w postaci obrazu, którego piksele tworzą macierz współczynników osłabienia promieniowania. Efektywność pracy scyntylatora zależy między innymi od [Stock, 2009]: - długości fali emitowanego promieniowania, - rodzaju materiału absorbującego (dla wysokich energii promieniowania stosowane są warstwy scyntylacyjne, wykonane z materiału o wysokiej liczbie atomowej), - sprawności emisji światła.. I.1.3. Budowa rentgenowskiego mikrotomografu komputerowego Badania wykonane. podczas. realizacji. pracy zostały przeprowadzone z. użyciem. mikrotomografu Benchtop CT160Xi wyprodukowanego przez firmę X-tek (Nikon). Poniżej zostanie omówiona budowa tego urządzenia oraz podstawowe parametry głównych jego elementów. 11.

(12) Najważniejszymi elementami budowy mikrotomografu rentgenowskiego są (rys. I.8): . źródło promieniowania rentgenowskiego,. . manipulator,. . detektor promieniowania rentgenowskiego.. Rys. I.8. Zdjęcie komory pomiarowej mikrotomografu Benchtop CT160Xi. Źródło promieniowania: Źródłem promieniowania w aparacie jest działo rentgenowskie posiadające wolframową katodę i anodę, o napięciu przyspieszającym z zakresu 40-160 kV, natężeniu 0-500μA i maksymalnej mocy 60W. Zarówno napięcie jak i natężenie prądu źródła są płynnie regulowane w swoich zakresach, z dokładnością do jedności. Przy optymalnych parametrach pracy możliwe jest uzyskanie wiązki o rozdzielczości ok. 4µm (rys. I.9).. Okno wzorca. Zbliżenie fragmentu z liniami 4µm. Rys. I.9. Radiogram wzorca do testowania rozdzielczości wiązki. 12.

(13) Uchwyt na próbki Manipulator znajdujący się pomiędzy lampą rentgenowską a detektorem, umożliwia przesuw próbki w 3 prostopadłych osiach X, Y, Z z dokładnością do 0,01mm. Wykorzystanie pełnego zakresu przesuwu na linii źródło – detektor, umożliwia osiągnięcie 100. krotnego powiększenia geometrycznego. Manipulator sterowany jest za pomocą trzech joystick’ów, które oprócz przesuwu uchwytu na próbkę w trzech prostopadłych osiach, umożliwiają także jego pełny obrót z dokładnością do 1/1000 stopnia. Manipulator, na którym znajduje się uchwyt na próbki, jest w stanie przesuwać próbki o masie do 5kg.. Detektor Aparat Benchtop CT160Xi wyposażony jest w płaski detektor PaxScan 2520V, firmy Varian. Matryca aktywna detektora składa się z 1900 (poziomo) x 1516 (pionowo) pikseli, co przy liniowym rozmiarze piksela 127µm daje powierzchnię 24,2cm  19,3cm. Dzięki zastosowaniu wysokowydajnych podzespołów elektronicznych urządzenia oraz gigabitowej karty sieciowej, detektor jest w stanie rejestrować nawet 30 klatek (projekcji) na sekundę, które na bieżąco mogą być zapisywane na twardym dysku komputera. 14. bitowy przetwornik analogowo-cyfrowy (A/D) zapewnia uzyskanie obrazów o dużym kontraście - posiadających 8192 odcieni w skali szarości.. 13.

(14) I.2. Otrzymywanie danych za pomocą mikrotomografu. W dalszej części rozdziału bardzo często będą pojawiać się pojęcia dotyczące opisu otrzymywanych obrazów. Zatem, w celu łatwiejszego zrozumienia występujących pojęć na wstępie objaśniam najczęściej używane zwroty. Skala szarości – przedział odcieni od koloru czarnego (niskie wartości w skali szarości) do koloru białego (wysokie wartości w skali szarości). Projekcja – (rys. I.10a) w najprostszym ujęciu projekcja jest cyfrowym zdjęciem rentgenowskim przedstawionym w skali odcieni szarości. Skala odcieni szarości obrazu jest proporcjonalna do wartości tłumienia promieniowania rentgenowskiego. W przypadku urządzenia Benchtop CT160Xi, kolory ciemne, aż do czarnego odpowiadają wysokim wartościom absorpcji, natomiast kolory jasne, aż do białego odpowiadają niskim wartościom lub praktycznie zerowej absorpcji. Przekrój zrekonstruowanego obrazu – (rys. I.10b) w przypadku obrazów zrekonstruowanych stosowana jest następująca skala odcieni szarości: kolory jasne, aż do białego odpowiadają wysokim wartościom współczynnika promieniowania X (kolorom ciemnym na projekcji), a kolory ciemne odpowiadają niskim wartościom współczynnika promieniowania X (kolorom jasnym na projekcji). Woksel – podstawowy element (najczęściej sześcian) w grafice trójwymiarowej, jest odpowiednikiem piksela w grafice dwuwymiarowej.. a). b). Rys. I.10. Projekcja próbki i jej zrekonstruowany przekrój. 14.

(15) Aby uzyskać informacje dotyczące badanego obiektu, należy wykonać trzy etapy pomiaru i przetwarzania (rys. I.11): . zarejestrować dane,. . zrekonstruować obraz,. . przeanalizować (przetworzyć) uzyskany obraz.. Rys. I.11. Schemat podstawowych etapów uzyskiwania danych. I.2.1. Przebieg pomiaru Zasada działania mikrotomografu jest analogiczna do działania medycznych tomografów, służących do badania ludzi. Główną różnicą pomiędzy tymi urządzeniami jest konfiguracja źródło – manipulator – detektor oraz rozdzielczość wiązki. W przypadku większości systemów mikrotomograficznych ( poza systemami ‘in - vivo’ do badań żywych organizmów) podczas wykonywania pomiaru, próbka wykonuje obrót o 360°, zatrzymując się co pewien krok kątowy, natomiast źródło promieniowania i detektor pozostają bez ruchu. W systemach medycznych źródło wraz z detektorami obracają się wokół pacjenta.. 15.

(16) Drugą różnicą jest rozdzielczość; systemy medyczne osiągają rozdzielczość woksela ok. 0,50,50,5mm3, natomiast mikrotomograficzne do 111µm3, a nanotomograficzne nawet do 505050nm3. Aby poprawnie przeprowadzić pomiar, konieczne jest wykonanie kilku podstawowych czynności. Należy umieścić próbkę na manipulatorze oraz ustawić odpowiednie powiększenie, a następnie ustawić w lampie rentgenowskiej napięcie przyspieszające i sprawdzić, czy próbka jest w pełni penetrowana przez promienie rentgenowskie. W przypadku tomografu Benchtop CT160Xi sygnał docierający do detektora, po przebyciu najbardziej absorbującej drogi promieniowania, musi być o ok. 2000 punktów (3% skali całego spektrum sygnału), silniejszy od wartości obrazu dark current (sygnału pochodzącego od termicznej generacji elektronów detektora, rejestrowanego przy wyłączonej lampie rentgenowskiej). Gwarantuje to odpowiednią różnicę pomiędzy sygnałem i szumem. W razie niespełnienia tego warunku należy zwiększyć napięcie przyspieszające promieniowania (wszystkie skanery umożliwiają także płynną regulację natężenia fotonów promieniowania X oraz zmianę czasu ekspozycji dla każdej projekcji, w celu uzyskania jak najlepszego kontrastu obrazu). Kolejno należy upewnić się, że podczas wykonywania obrotu o 360° żaden fragment obiektu nie jest rzutowany poza powierzchnię detektora. Niektóre algorytmy programów rekonstruujących dane umożliwiają wykonanie pomiaru bez pełnego obrotu próbki lub w przypadku, gdy rzut obiektu wychodzi poza obszar detektora (pomiar z offset’em) [UGCT] (rys. I.12).. Rys. I.12. Przykład skanowania bez i z offset’em [UGCT]. Podczas wykonywania pomiaru aparat zarejestruje zadaną liczbę obrazów uzyskanych co zadany krok kątowy (360°/liczba projekcji). 16.

(17) Uzyskiwane podczas pomiaru obrazy są znormalizowane, poprzez zastosowanie korekcji niejednorodności czułości detektora (flat filed) oraz zminimalizowanie wpływu termicznej generacji elektronów w matrycy detektora (dark current) [Bielecki i in., 2009]. W przypadku systemu Benchtop CT160Xi, po uruchomieniu kreatora pomiaru, system rejestruje sygnał docierający do detektora przy wyłączonym promieniowaniu rentgenowskim (dark current – rys. I.13b), a następnie sygnał docierającego promieniowania do detektora, przy próbce wyciągniętej z pola widzenia (flat field – rys. I.13c). Wszystkie kolejno uzyskiwane obrazy są normalizowane, poprzez odjęcie od uzyskiwanej projekcji wartości obrazu dark current i podzielenie przez wartości obrazu flat field [UGCT]. Na końcu tej operacji obraz zostaje przeskalowany w celu uzyskania najlepszego kontrastu danych. Obraz znormalizowany (rys. I.13d) wyraźnie różni się od obrazu wejściowego (rys. I.13a).. a). b). c). d). Rys. I.13. Przykład normalizacji projekcji mikrotomograficznej. W celu uzyskania przestrzennego obrazu wewnętrznej struktury należy przejść do kolejnego etapu, jakim jest rekonstrukcja. 17.

(18) I.3. Rekonstrukcja obrazów Podczas procesu projekcji wstecznej zostaje odtworzona przestrzeń struktury badanej próbki na podstawie uzyskanych dwuwymiarowych projekcji. Bardzo dobrze ten proces ilustrują rysunki I.14 i I.15 z pracy Stock’a [2009], które schematycznie przedstawiają proces zapisu projekcji (rys. I.14) oraz proces projekcji wstecznej (rys. I.15). Rysunek I.14 przedstawia układ kół o jednakowym współczynniku pochłaniania promieniowania rentgenowskiego. Histogramy oznaczone literami i, ii oraz iii przedstawiają wartości absorpcji zarejestrowane na projekcjach dla trzech pozycji kątowych. Wyższym słupkom, odpowiada większa absorpcja promieniowania X. W przypadku rzeczywistym wartości pochłaniania na krawędziach prześwietlanych obiektów byłyby znacznie niższe niż na środku, natomiast w celu uproszczenia opisu wybrano przykład monotonicznego rozkładu absorpcji.. Rys. I.14. Zapis projekcji układu kół,. Rys. I.15. Proces projekcji wstecznej. dla 3 pozycji kątowych [Stock, 2009]. [Stock, 2009]. Na podstawie uzyskanych profili można zademonstrować, w jaki sposób zostaje wykonana projekcja wsteczna. Rysunek I.15a przedstawia projekcję wsteczną wzdłuż kierunku, w którym projekcja została uzyskana (rys. I.14i). Na rysunku I.15b dodano projekcję, z kolejnej pozycji kątowej (rys. I.14ii), a na obrazie I.15c przedstawiono przecinające się linie z 3 pozycji. Wszystkie punkty przecięcia się linii z 3 projekcji określają możliwe miejsca występowania kół. Na rysunku I.15d ciemne koła przedstawiają odwzorowany układ z rysunku I.14. Położenie kół wypełnionych czerwonym kolorem określa dodatkowe pozycje, w których, na podstawie obliczeń projekcji wstecznej tylko dla trzech pozycji, mogą występować badane obiekty. Błędne pozycje (czerwone koła) można wyeliminować dopiero po naniesieniu informacji z projekcji dla dodatkowych pozycji kątowych.. 18.

(19) Wadą stosowania samej projekcji wstecznej, jest nanoszenie na obraz wartości pochłaniania w miejscu, gdzie podczas prowadzenia pomiaru nie znajdował się żaden obiekt. Dopiero zastosowanie filtrowanej projekcji wstecznej powoduje minimalizację tego efektu. Filtracja uzyskanych projekcji ma na celu wygenerowanie ujemnych wartości pochłaniania, które mają zniwelować wartości pochłaniania powstałe na obrazie uzyskanym po samej projekcji wstecznej [Stock, 2009]. Poniżej, na rzeczywistym przykładzie, przedstawiono wpływ liczby projekcji na jakość zrekonstruowanego obrazu. Badanym obiektem był układ zilustrowany na rysunku I.16 składający się z ołówka, polimeru i próbki skały. Celowo wybrano obiekty o różnych kształtach i różnych wartościach pochłaniania promieniowania rentgenowskiego, aby móc je łatwiej rozróżnić na zrekonstruowanych przekrojach. Projekcja badanego układu (rys. I.16b) pokazuje, że największą sumaryczną wartość pochłaniania promieniowania rentgenowskiego (największe zaczernienie) posiada próbka skały, kolejno grafit ołówka i próbka polimeru.. a). b). Rys. I.16. Przykładowy badany model oraz jego projekcja. Rysunek I.17 przedstawia wpływ liczby projekcji, z których została zrekonstruowana struktura obiektu, na jakość uzyskanego obrazu. Poszczególne obrazy uzyskano z 3, 5, 25, 125, 625 i 3125 projekcji. Przedstawione przekroje poprzeczne wykonano na tej samej wysokości układu, oznaczonej na rysunku I.16b czerwoną linią. Każdy z obrazów osobno znormalizowano [Wojnar i in., 2002], w celu uzyskania jak najlepszego kontrastu i wyraźniejszego przedstawienia różnicy wpływu liczby projekcji, na jakość rekonstrukcji.. 19.

(20) 20 Ilosć projekcji - 625. Ilosć projekcji - 5. Rys.I.17. Przykład wpływu liczby projekcji na jakość rekonstruowanego obrazu. Ilosć projekcji - 125. Ilosć projekcji - 3. Ilosć projekcji - 3125. Ilosć projekcji - 25.

(21) Dla rekonstrukcji wykonanej z 3 i 5 projekcji praktycznie nie można określić kształtu ani liczby skanowanych obiektów, dopiero na obrazie uzyskanym po rekonstrukcji z 25 projekcji wyraźnie zauważalne stają się kontury badanych obiektów. Analizując rys. I.17 można zauważyć, że wraz ze wzrostem liczby projekcji poprawia się jakość obrazu oraz liczba dostrzegalnych szczegółów do tego stopnia, że na przekroju otrzymanym z rekonstrukcji 3125 projekcji można dostrzec wewnętrzną mikrostrukturę drewna, z którego wykonany jest ołówek. Linie, wyraźnie widoczne na pierwszych czterech rysunkach, zanikają na dalszych, gdyż zwiększa się stosunek sygnału do szumu. Ma także miejsce efekt, o którym wcześniej wspomniano: gdy linie powstałe podczas wykonywania projekcji wstecznej przecinają się, następuje ich wzmocnienie, jeżeli nie przecinają się, to osiągane w tym miejscu wartości są niskie i podobne do wartości tła.. I.3.1. Sztuczne elementy (artefakty) rekonstrukcji Metoda mikrotomografii, jak każda metoda pomiarowa, posiada swoje ograniczenia oraz generuje sztuczne elementy (artefakty) w uzyskiwanych obrazach. Najczęściej spotykanymi artefaktami są: . koncentryczne pierścienie (ring artifact),. . twardnienie wiązki (beam hardening),.  smugi (streak artifact). Koncentryczne pierścienie (ang. ring artifact) Ten rodzaj artefaktu, pojawia się na zrekonstruowanych obrazach jako pełne lub niepełne okręgi (rys. I.18), których środkiem jest środek obrotu manipulatora. Powodem powstawania tego rodzaju artefaktów są niedokładności detektora, które zazwyczaj mogą być powodowane wypalonymi pikselami detektora lub błędnym odczytem intensywności promieniowania [Sijbers i Postnov, 2004]. Błąd ten można minimalizować poprzez ruchy skanowanej próbki w różnych kierunkach pomiędzy kolejnymi pozycjami projekcji. Jednak taka korekcja wydłuża czas pomiaru i nie eliminuje w pełni artefaktu.. 21.

(22) Rys. I.18. Przykład artefaktu pierścienia obrazu próbki skały. Twardnienie wiązki (ang. beam hardening) Lampa rentgenowska w mikrotomografie rentgenowskim emituje polichromatyczną wiązkę promieniowania, przez co może powstawać artefakt twardnienia wiązki. Efekt ten powstaje, gdyż wraz ze wzrostem drogi penetracji przez promieniowanie X, średnia energia wiązki rośnie. Dzieje się, tak ponieważ fotony o niższej energii są znacznie łatwiej absorbowane na powierzchni próbki [Van de Casteele, 2004]. Artefakt ten jest bardzo dobrze widoczny, przy tworzeniu liniowego przekroju wartości pochłaniania przez próbkę. Na rys. I.19 czerwonymi strzałkami zaznaczono zawyżone wartości pochłaniania promieniowania na krawędzi próbki.. Rys. I.19. Wykres wartości pochłaniania w przekroju liniowym próbki. (Oś odciętych przedstawia długość przekroju, a oś pionowa, wartości pochłaniania). Artefakt ten można stosunkowo łatwo minimalizować, poprzez stosowanie jednorodnych filtrów metalowych, których zadaniem jest absorbowanie niskoenergetycznego spektrum wiązki. Rysunek I.20 przedstawia zmianę widma energii przy wzroście napięcia 22.

(23) przyspieszającego (linia przerywana B), ciągła linia F na tym rysunku przedstawia efekt zastosowania filtru metalowego eliminującego promieniowanie niskoenergetyczne.. Rys. I.20. Wykresy widma energii fotonów dla różnych napięć przyspieszających (I – natężenie fotonów, E – energia fotonów, kVp – zmiana potencjału) [Stock, 2009]. Na rysunku I.21 przedstawiono efekt zastosowania filtru o grubości 0,4mm wykonanego z miedzi, kompensującego efekt twardnienia wiązki. Wartości pochłaniania uzyskiwane na krawędzi próbki są zbliżone do wartości z wewnątrz analizowanego przekroju.. Rys. I.21. Kompensacja artefaktu twardnienia wiązki przy zastosowaniu filtra miedzianego (Oś odciętych przedstawia długość przekroju, a oś pionowa, wartości pochłaniania). 23.

(24) Smugi (ang. streak artifact) Artefakt ten występuje, gdy w skanowanym obiekcie znajdują się fragmenty o bardzo dużym kontraście w tłumieniu promieniowania. Spowodowany jest wytworzeniem wtórnego, rozproszonego promieniowania podczas penetracji gęstego obszaru. Promieniowanie to na rekonstruowanym obrazie tworzy smugi o obniżonych wartościach absorpcji [Cnudde, 2005]. Artefakt ten jest bardzo popularny w medycznych zastosowaniach tomografii komputerowej. Występuje zazwyczaj, gdy w ciele pacjenta znajdują się bardzo gęste obiekty, takie jak amalgamat wypełniający zęby [Różyło – Kalinowska i in., 2006], czy części metalu służące lepszemu zrastaniu się kości. Podczas badania próbek skał, artefakt ten pojawia się w otoczeniu minerałów o wysokiej gęstości, np. pirytu. Pochłanianie obserwowane wokół ziaren pirytu jest porównywalne z pochłanianiem porów (rys. I.22). W przypadku próbek posiadających niski współczynnik porowatości, występowanie tego artefaktu może uniemożliwić analizę struktury porowej próbki, gdyż występowanie na przekroju sztucznie generowanych niskich wartości absorpcji może znacznie zawyżyć obliczony współczynnik porowatości.. Rys. I.22. Przykład występowania artefaktów pasmowych w przekroju próbki piaskowca. Efekt objętości częściowej (ang. partial volume effect) Zaczernienie każdego woksela jest wprost proporcjonalne do współczynnika tłumienia fragmentu próbki, który dany woksel reprezentuje (a przez to – do jego gęstości). Należy jednak pamiętać o efekcie objętości częściowej (ang. partial volume effect). Przejawia się on 24.

(25) tym, że wartość zaczernienia woksela, który pokrywa fizycznie materiał niejednorodny (np. granicę pomiędzy porem a skałą), będzie średnią ważoną z wartości zaczernień, które dają poszczególne fazy będące w obrębie woksela (rys. I.23) [Zalewska i in., 2010]. Na poniższym rysunku (rys. I.23a) po lewej stronie przedstawiony jest schemat odzwierciedlający rzeczywisty kształt obiektu o dwóch fazach gęstościowych. Natomiast po prawej stronie, obraz tego obiektu złożony z dziewięciu pikseli. Piksele na rysunku I.23b posiadają więcej niż dwie wartości, będące wynikiem uśrednienia wartości faz będących w obrębie jednego piksela.. a). b). Rys. I.23. Wpływ efektu objętości częściowej na obraz badanego obiektu [Kobayashi Y. i in. 2010]. 25.

(26) ROZDZIAŁ II – Analiza przestrzennych obrazów struktury porów II.1. Przetwarzanie obrazu Ostatnim krokiem na drodze do otrzymania wyników liczbowych, charakteryzujących badaną próbkę, jest przetwarzanie i analiza obrazu. Ze względu na duże obciążenie komputera analizą trójwymiarowych obrazów, zazwyczaj istnieje konieczność ograniczenia się do analizy mniejszego fragmentu niż obraz całej badanej próbki. Pierwszym etapem przetwarzania obrazu jest filtracja – etap przygotowujący do segmentacji. Filtracja stosowana jest w celu poprawy jakości obrazu, zwiększenia stosunku sygnału do szumu oraz uwydatnienia różnicy wartości na krawędziach pomiędzy fazami o różnych współczynnikach tłumienia. Następnie, obraz należy poddać segmentacji, czyli podzielić go na obszary jednorodne pod względem wartości tłumienia promieniowania X. Najprostszym przypadkiem segmentacji jest binaryzacja. Proces ten polega na sprowadzeniu obrazu o szerokim zakresie wartości w skali szarości, do obrazu dwuwartościowego – binarnego (zazwyczaj czarno-białego). W przypadku próbek skał binaryzacja polega na wydzieleniu z przestrzeni próbki wokseli, reprezentujących szkielet skalny oraz strukturę porową (rys. II.1). Rysunek II.1a przedstawia w odcieniach szarości przekrój przez próbkę piaskowca (kolor szary, reprezentuje szkielet skalny, natomiast kolor czarny odpowiada strukturze porowej), natomiast rysunek II.1b jest binarnym obrazem II.1a, gdzie kolor biały odpowiada przestrzeni porowej, a kolor czarny szkieletowi skalnemu.. a). b). Rys. II.1. Przykład binaryzacji struktury piaskowca. 26.

(27) Binaryzacja Najbardziej popularnym sposobem przeprowadzenia binaryzacji jest progowanie (ang. thresholding) z wykorzystaniem histogramu rozkładu stopni szarości. Proces ten polega na wyznaczeniu. granicznej. wartości. zaczernienia. (progu). charakterystycznego. dla. analizowanych faz. Po wyznaczeniu wartości progowych wartość zaczernienia (GV) każdego woksela obrazu jest porównywana z wartością progową (Th), po czym każdemu wokselowi przypisywana jest nowa wartość zaczenienia (0 lub 1):. gdzie: GV (x,y,z) – wartość poziomu szarości danego punktu w obrazie źródłowym, GV’ (x,y,z) – wartość odpowiedniego punktu w obrazie binarnym, Th – próg binaryzacji. Binaryzacja. obrazów. zawartych. w. oprogramowania ImageJ oraz AVIZO.. pracy. została. wykonana. z. wykorzystaniem. Zastosowano metodę binaryzacji wzdłuż granic. [Kaczmarczyk i in., 2010]. W opinii autora ustalanie wartości progowej (rozgraniczającej na obrazie strukturę porową od szkieletu) metodą wzdłuż granic (ang. thresholding along boundaries) jest dobrą metodą podziału obrazu mikrotomograficznego skały na pory i szkielet skalny. Metoda jest generalnie prosta do zastosowania w przypadku większości próbek i daje wiarygodne wyniki [Dohnalik i in., 2009]. Pierwszym etapem metody jest obliczenie histogramu wartości zaczernienia obrazu (rys. II.2). Histogram przestrzennych obrazów skał jest wykresem bimodalnym, tzn. składającym się z dwóch maksimów. Pierwsze - o niższych wartościach – odpowiada powietrzu, a drugie, o większych wartościach, odpowiada matrycy skalnej (rys. II.2a). W tej metodzie zakłada się, że wartość rozdzielająca powietrze od szkieletu skały leży w okolicy minimum (Th min) pomiędzy dwoma szczytami. Jednak ze względu na efekt objętości częściowej (partial volume effect), wartość rozgraniczająca pory od szkieletu jest większa od wartości Thmin [Youssef i in., 2007]. Metodyka poszukiwania wartości progowej polega na zaznaczeniu wokseli o wartości zaczernienia Thmin5. Efektem tej czynności jest przybliżone wyznaczenie granicy pomiędzy porami a szkieletem skalnym.. 27.

(28) Pośród zaznaczonych wokseli oraz wokseli sąsiadujących z nimi (przy sąsiedztwie 4) należy znaleźć odpowiednią ilość par ‘woksel zaznaczony’ – ‘woksel sąsiadujący’, w których woksel sąsiadujący ma wartość zaczernienia o około 20 wyższą od woksela zaznaczonego. Zakłada się, że przy takiej różnicy woksel sąsiadujący odpowiada szkieletowi skalnemu. Wartości punktów sąsiadujących należy mierzyć w miejscach jednorodnych i reprezentatywnych dla próbki. Należy unikać obszarów o podwyższonych wartościach tłumienia mogących sztucznie zawyżyć wartość rozgraniczającą i wpłynąć na podwyższenie objętości struktury porowej (przy zawyżonej wartości, do grupy ‘pory’ zostaną przypisane woksele należące do matrycy skalnej). W następnym kroku należy obliczyć średnie zaczernienie wybranych wokseli zaznaczonych oraz wokseli sąsiadujących. Średnia pomiędzy tymi dwoma wartościami da poszukiwaną wartość progową, przyjętą do binaryzacji.. a). b). c). Rys. II.2. Binaryzacja obrazu metodą progowania wzdłuż granic (a) histogram całego obrazu z zaznaczoną wartością Thmin, (b) zaznaczone woksele graniczne pomiędzy skałą a porami, (c) obraz binarny przy wyliczonej wartości progowej (kolorem czarnym zaznaczono szkielet skalny, kolorem białym – pory). 28.

(29) II.2. Analiza liczbowa struktury porowej Dysponując obrazem binarnym można przystąpić do liczbowego parametryzowania struktury badanego obiektu. Pierwszym, podstawowym etapem obliczeń jest wyliczenie współczynnika porowatości próbki, uzyskanego z równania (II.1): (II.1). gdzie: Kp – współczynnik porowatości analizowanej próbki, Vp – objętość wokseli reprezentujących pory, Vsz – objętość wokseli reprezentujących szkielet.. Wydzielanie obiektów niepołączonych W celu uzyskania dodatkowych informacji o obrazie struktury porowej, można wydzielić wszystkie odizolowane obiekty poprzez analizę sąsiedztwa punktów. Liczbowa procedura tej analizy dla obrazu 2D, dla sąsiedztwa „4” (analizującego połączenia w prostopadłych kierunkach) przedstawiona jest poniższą zależnością (II.2): N4(p)=[(x+1,y),(x-1,y),(x,y+1), (x,y-1)],. (II.2). gdzie : N4(p) – zbiór punktów sąsiadujących z punktem p, p – jest pikselem, którego otoczenie jest badane, x, y – współrzędne. Piksel q należy do sąsiedztwa piksela p, jeżeli należy do zbioru punktów N4(p) [Fisher i in., 2003]. Oprogramowanie, zaczynając od wybranego początkowego punktu, należącego do warstwy porów, analizuje jego sąsiedztwo, a następnie analizuje sąsiedztwo kolejnych punktów spełniających warunek (II.2). Na rysunku II.3 przyjęto wartość 1, dla pikseli należących do porów, a wartość 0 dla pikseli należących do matrycy skalnej. Rysunek II.3 ilustruje sytuację, gdy w analizowanym obrazie 29.

(30) występują dwa obiekty niepołączone między sobą. Pierwszy z nich składa się z 8 punktów (lewa strona rysunku), a drugi z 6 punktów (prawa strona rysunku). Rys. II.3. Schematyczny opis określania sąsiedztwa dla obrazu 2D [Fisher i in., 2003]. Graficzna ilustracja analizy sąsiedztwa punktów przedstawiona jest na rysunku II.4. Rysunek II.4a przedstawia obraz binarny, na którym pory oznaczone są kolorem białym, natomiast rysunek II.4b przedstawia wydzielonych 11 obiektów oznaczonych różnymi kolorami.. a). b). Rys. II.4. Ilustracja graficzna wydzielania obiektów niepołączonych. Analogiczna sytuacja zachodzi dla wokseli obrazu przestrzennego. W niniejszej pracy w celu wydzielenia. odseparowanych. obiektów. analizowane. jest. sąsiedztwo. 26. punktów. sąsiadujących, do woksela początkowego oznaczonego kolorem niebieskim (rys. II.5). Rysunek II.6b przedstawia wydzielone obiekty w przestrzeni struktury porowej oznaczonej kolorem zielonym na rysunku II.6a. 30.

(31) Rys. II.5. Schemat sąsiedztwa 26 punktów [MAVI 1.3.1 documentation, Neighbourhood]. a). b). Rys. II.6. Obraz struktury porowej (a) oraz wydzielonych obiektów niepołączonych (b). Klasyfikacja objętościowa Każdy z wydzielonych obiektów otrzymuje swój numer porządkowy oraz jest oznaczony kolorem, lub odcieniem szarości. Podczas dalszej analizy możliwe jest uzyskanie informacji o objętości każdego z obiektów (tab. II.1). Oprogramowanie MAVI umożliwia także filtrację obrazu 3D ze względu na objętość poszczególnych obiektów i pozwala wyświetlić podtypy porów z wybranego zakresu ich objętości. Funkcja ta będzie stosowana w pracy w celu klasyfikowania obiektów ze względu na ich objętość oraz wizualizowania ich dystrybucji przestrzennej.. 31.

(32) Tabela II.1. Zestawienie objętości poszczególnych wydzielonych obiektów. Nr obiektu. Objętość [µm3]. 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20. 32203872 481248 432 6912 648 648 216 15120 216 216 432 648 648 432 1080 5832 49032 648 648 432 1296. Krętość kanalików porowych Krętość jest cechą geometryczną ścieżki połączonych punktów (pikseli, wokseli) zależną od ilości zakrętów na ścieżce (rys. II.7). Krętość jest liczona jako stosunek długości ścieżki pomiędzy dwoma punktami do odległości pomiędzy tymi punktami. Parametr ten określony jest wzorem (II.3).  ij . lj li. gdzie, lj - długość ścieżki pomiędzy dwoma punktami, li – odległość pomiędzy dwoma punktami.. 32. (II.3).

(33) lj li. Rys. II.7. Schemat krętego obiektu. W programie MAVI funkcja Geometric Tortuosity pozwala na zbadanie krętości kanałów porowych obecnych w próbce. Aby program MAVI w danym kierunku obliczył krętość, próbka musi posiadać następujące cechy:  obecność co najmniej jednego woksela reprezentującego warstwę porów na przeciwległych ścianach próbki w określonym kierunku (x, y, z),  obecność ciągłej ścieżki pomiędzy tymi punktami, utworzonej przez woksele sąsiadujące ze sobą. Jeżeli spełnione są podane wyżej warunki, dla każdej z j ścieżek program oblicza krętość τi dla kierunku i-tego według wzoru. (II.3). [MAVI 1.3.1 documentation, Geometric. Tortuosity]. Wynikiem tych obliczeń jest podanie maksymalnej, minimalnej i średniej wartości krętości w badanym kierunku. Krętość, zgodnie z równaniem (II.3) może przyjmować wartości większe lub równe 1. Wartość minimalna jest przypisywana ścieżkom prostym - takim, w których wszystkie woksele połączone między sobą tworzą ścieżkę linii prostej równoległej do osi, w której mierzona jest krętość. Im wyższa wartość krętości, tym bardziej przebieg badanej ścieżki odbiega od prostej. W analizie obrazu mikrotomograficznego skał, analiza krętości geometrycznej pozwala wyciągać wnioski na temat własności transportowych badanej próbki.. 33.

(34) Badanie niejednorodności rozkładu struktury porowej W celu określenia jednorodności przestrzennego rozkładu struktury porowej skorzystano z modułu SUBFIELD FEATURES oprogramowania MAVI. Moduł ten pozwala zmierzyć tzw. porowatość lokalną. Metoda pomiaru polega na podzieleniu obrazu wejściowego na mniejsze podpróbki (sześciany) i wyliczeniu dla każdej z nich współczynnika porowatości. Analizę porowatości lokalnej wejściowych obrazów o rozmiarach 900 × 900 × 400 wokseli wykonano dla podpróbek o wymiarze 1003 wokseli, czyli dla każdej próbki uzyskano 324 wartości porowatości lokalnej. Jako kryterium oceny jednorodności struktury porowej przyjęto wartość względnego odchylenia standardowego (σr), wyrażonego wzorem (II.4):. (II.4) gdzie, σ100 – wartość odchylenia standardowego współczynnika porowatości lokalnej dla podpóbki o rozmiarze 100 x 100 x 100 woskeli, KpCT – średnia wartość porowatości dla obrazu 900 x 900 x 400 wokseli. W celu określenia wartości granicznej, poniżej której porowatość uznawana jest za równomiernie rozłożoną, przeprowadzono badania testowe na dwóch próbkach białego piaskowca szydłowieckiego. Jest on uważany za materiał, którego struktura porowa jest jednorodnie wykształcona. Stwierdzenie o jednorodnym rozkładzie porów poparte jest różnorodnymi badaniami parametrów petrofizycznych na próbkach. W przypadku tych dwóch próbek wartość względnego odchylenia standardowego wyniosła 0,16. W pracy przyjęto, że dla próbek piaskowca czerwonego spągowca, których wartość względnego odchylenia standardowego nie przekracza 0,3, struktura porowa będzie określona jako jednorodnie wykształcona. Wyliczenie wymiaru pudełkowego struktury porowej W celu wyliczenia wymiaru pudełkowego zastosowano moduł FractalCount (3D) w oprogramowaniu ImageJ [ImageJ], który został opracowany przez Pera Christiana Hendena i Jensa Bachego-Wiiga [Henden, Bache – Wiig]. Analizie poddane zostały obrazy o wymiarze 900 × 900 × 400 wokseli. W pojawiającym się oknie dialogowym modułu, użytkownik musi podać wartość progową (threshold) dla warstwy obrazu, dla której obliczenia będą. 34.

(35) prowadzone oraz początkowy rozmiar pudełka (start box size), minimalny rozmiar pudełka (min box size) i współczynnik redukcji (box division factor). Analiza wykonywana przez algorytm jest analogiczna do tradycyjnego sposobu pomiaru wymiaru pudełkowego na obrazach dwuwymiarowych [Peitgen i in., 1995]. Na przestrzenny obraz wejściowy nakładana jest siatka o wielkości zadanego przez użytkownika pudełka (sześcianu). Następnie badana jest liczba sześcianów zawierających w sobie woksele należące do warstwy porów. W kolejnych etapach rozmiar pudełka siatki jest zmniejszany o zadany przez użytkownika stopień redukcji. Algorytm kończy obliczenia po osiągnięciu rozmiaru siatki określonego przez użytkownika jako minimalny rozmiar pudełka. Wynikiem obliczeń jest prosta log(liczba zliczeń) od - log(rozmiar pudełka), której wartość współczynnika kierunkowego (slope) jest wymiarem pudełkowym (rys. II.8).. Rys. II.8. Wynik obliczeń algorytmu FractalCount 3D w postaci funkcji log(liczba zliczeń) od log(rozmiar pudełka), wartość współczynnika kierunkowego prostej jest wymiarem pudełkowym analizowanej struktury. Średnia długość cięciwy poprowadzonej przez warstwę porów Za kolejny interesujący parametr, możliwy do wyliczenia z obrazów mikrotomograficznych, uznano średnią długość cięciw (ang. mean chord length). Parametr ten określa średnią długość cięciw poprowadzonych, przez warstwę porów [Osher, Schladitz 2009, Młynarczuk, 2004]. W poniższej tabeli (II.2) przedstawiono zależność pomiędzy rozmiarem wielkości podsystemów porów w analizowanym obrazie, a średnią wartością długości cięciw dla badanej próbki.. 35.

(36) Tabela II.2. Zależność długości średniej cięciwy od objętości analizowanych porów. A. B. Objętość analizowanych porów [woksel] < 10. Długość średniej Objętość cięciwy analizowanych porów [woksel] [woksel] 1,8 < 10. Długość średniej cięciwy [woksel] 1,8. 10 – 100. 2,1. <100. 2,0. 100 – 1 000. 3,2. < 1 000. 2,26. 1 000 – 10 000. 4,6. < 10 000. 2,34. 10 000 <. 5,3. 10 000 <. 4,99. W części A tabeli przedstawiono zmiany średniej długości cięciw w zakresie poszczególnych klas objętości, natomiast w części B zakresy zmian dla kolejnych klas wraz z mniejszymi. Zauważalny jest wyraźny wzrost długości średniej cięciwy podczas zwiększania się objętości analizowanych podsystemów porów. Można zatem spodziewać się, że próbki o dobrze rozbudowanej strukturze porowej będą posiadały wysoką wartość tego parametru.. Charakterystyka Eulera Hans Jorg – Vogel już w artykule z 2002 roku [Vogel, 2002] zwraca uwagę na wpływ charakterystyki Eulera na właściwości analizowanej struktury. Odnosi się do próbki gleby i rozważa występujące w niej bardzo nieregularne kształty (szczeliny, rozprzestrzenianie się korzeni czy kanały wytworzone przez robaki). Autor zaznacza, iż bardzo ważnym parametrem w ocenie struktury porowej jest charakterystyka Eulera (ang. Euler number) [Toriwaki, Yonekura, 2002]. Podkreśla, że sposób, w jaki pory są połączone ma większy wpływ na właściwości struktury, niż ich ilość, czy rozmiar. Charakterystyka Eulera może być określona następującym wzorem [Vogel, 2002] (II.5): χ=N−C+H gdzie: N – liczba izolowanych porów, C – ilość połączeń, H – liczba ziaren szkieletu całkowicie otoczonych pustką.. 36. (II.5).

(37) Wartość współczynnika H dla obrazów skał jest praktycznie pomijalna, gdyż w naturze nie są spotykane cząstki szkieletu skalnego całkowicie otoczone powietrzem. Dla powyższych założeń, wartość charakterystyki Eulera może być traktowana jako miara połączenia struktury porowej. Gdy otrzymywane są wartości dodatnie (N > C), oznacza to, że struktura porowa jest słabo połączona, w przypadku przeciwnym (C > N) struktura porowa posiada wiele połączeń pomiędzy porami. Tabela II.3 przedstawia zmiany wartości charakterystyki Eulera dla poszczególnych zakresów objętości systemów porów występujących w obrazie. Porównaniu został poddany obraz o rozmiarze 200 x 200 x 200 wokseli, o wyliczonym współczynniku porowatości 22%. Pierwsza kolumna przedstawia zakres objętości systemów porów pozostających w obrębie analizowanego obrazu, druga – ilość obiektów poddanych analizie, a w ostatniej kolumnie podano wartości charakterystyki Eulera. Tabela II.3. Zależność wartości charakterystyki Eulera od objętości analizowanych porów. Objętość analizowanych porów [woskel]. Ilość analizowanych obiektów. Wartość charakterystyki Eulera. do 10. 3562. 3425. do 100. 6439. 6015. do 1000. 6467. 5969. do 10000. 6468. 5962. wszystkie. 6469. -5903. powyżej 10. 3306. -8944. pow. 100. 32. -11916. pow. 1000. 2. -11872. pow. 10000. 1. -11865. Powyższe porównanie wykazało wyraźną zależność pomiędzy ilością izolowanych obiektów w przestrzeni porowej, a wartością charakterystyki Eulera. Widoczny jest wzrost wartości charakterystyki Eulera dla obiektów o małej objętości 10 i 100 wokseli, następnie wartości zaczynają maleć, gdyż w analizowanym obszarze pojawiają się połączone systemy porów. Trend ten utrzymuje się, gdy z całego obrazu stopniowo usuwane są izolowane obiekty -. 37.

(38) wartość charakterystyki Eulera maleje od -5903 do -11916. Następnie zaczyna nieznacznie rosnąć, gdyż z obrazu usuwane są dobrze połączone podsystemy porów, osiągając wartość -11865 dla jednego bardzo dobrze rozbudowanego systemu porów V klasy objętości. W tabelach wynikowych (Rozdział IV, tabele IV.2 - IV.4) przedstawiono znormalizowane wartości charakterystyki Eulera otrzymane poprzez podzielenie wyliczonej wartości tego parametru przez liczbę wszystkich obiektów wydzielonych podczas analizy obrazu.. Granulometry – informacja o średnicy struktury porowej W tym punkcie omówiony zostanie sposób, w jaki za pomocą analizy obrazu można uzyskać informację o rozkładzie średnic struktury porowej. Przedstawione dane otrzymano z wykorzystaniem funkcji Granulometry programu MAVI [MAVI 1.4.0 Documentation, Granulometry]. Algorytm ten dla każdego punktu obrazu należącego do warstwy porów wyznacza maksymalną średnicę kuli, która całkowicie zawiera się w analizowanej warstwie obrazu oraz obejmuje cały woksel. W ten sposób użytkownik uzyskuje informację o lokalnej grubości czy średnicy pustki. Po zakończeniu omawianego procesu obliczeniowego program w postaci tabelarycznej podaje informację, ile kul o danej średnicy można wpisać w analizowany obraz (histogram, rys. II.9). Środki opisanych kul leżą na linii wyznaczającej oś środkową analizowanej struktury. Pliki wynikowe pomijają średnice o długości poniżej 4 wokseli, ponieważ zazwyczaj odnoszą się one do powierzchni chropowatej struktury. Na rysunku II.9 przedstawiono przykładowy wynik obliczeń z wykorzystaniem algorytmu Granulometry. W tabeli zestawiono dane otrzymane dla dwóch próbek (9513 i 7483) o zupełnie odmiennych właściwościach zbiornikowych. W kolumnie drugiej przedstawiony jest udział średnic o poszczególnych długościach. Natomiast w kolumnie trzeciej, znajdują się krzywe kumulacyjne objętości kul o poszczególnej średnicy, wraz z mniejszymi. Wartości wyrażone na osiach odciętych, podane są w µm. Kolejną zaletą tego modułu jest możliwość zobrazowania fragmentu próbki w którym występują pory o największej średnicy. Rysunek II.10. przedstawia pory o największych średnicach (zakres 110 - 144µm) występujące w próbce 9513.. 38.

(39) 0,05. 6E+09. Sumaryczna objętość kul wraz z mniejszymi [µm3 ]. 0,04 0,03 0,02. 0,01. 144. 139. 135. 130. 126. 122. 118. 113. 109. 99. 104. 94. 87. 81. 76. 69. 60. 51. 40. 0 24. 9513. Krzywa kumulacyjna objętości kolejnych kul. Histogram udziału kul o danej średnicy Udział kul o danej średnicy [%]. Nr próbki. 5E+09 4E+09 3E+09. 2E+09 1E+09 0 200. 150. 0,1. 0,08 0,06 0,04. 0,02 66. 62. 60. 56. 54. 51. 48. 43. 40. 36. 29. 0. Sumaryczna objętość kul wraz z mniejszymi [µm3 ]. 0,12. 24. Udział kul o danej średnicy [%]. 50. 0. Średnica kuli [µm]. Średnica kuli [µm]. 7483. 100. 2,5E+07 2,0E+07. 1,5E+07 1,0E+07 5,0E+06 0,0E+00 80. 60. 40. 20. 0. Średnica kuli [µm]. Średnica kuli [µm]. Rys. II.9. Histogram występowania kul o poszczególnej średnicy oraz ich krzywe kumulacyjne. Rys. II.10. Wizualizacja porów o największych średnicach ( z zakresu 110 - 144µm), występujących w próbce 9513. 39.

(40) Analiza szkieletu (wyliczenie liczby koordynacyjnej struktury porowej) Szkielet jest obrazem odwzorowującym charakterystykę geometryczną i topologiczną struktury [Yang, 2005; MAVI 1.3.1 Documentation, Skeleton]. Transformacja obrazu w szkielet, zwana także transformacją do osi środkowej (ang. medial axis transformation), polega na izotropowym wycinaniu („erozji”) białych pikseli na całej powierzchni elementów obrazu [Dong, 2008]. Punkty, w których następuje zetknięcie się płaszczyzn wycinania, zostają dołączone do obrazu szkieletu próbki. Wynik przykładowej transformacji zaprezentowano na rysunku II.11.. (a). (b). Rys. II.11. Obraz binarny (a) oraz efekt jego transformacji do osi środkowej (b). W szkielecie obrazu można wyróżnić następujące elementy:  zakończenia - piksele umieszczone na końcu pojedynczej gałęzi,  węzły - piksele stanowiące skrzyżowania dwóch gałęzi,  gałęzie - każdy ciąg pikseli, zakończony skrzyżowaniem lub węzłem,  pętle - układ gałęzi tworzący figurę zamkniętą. Określenie tych elementów następuje na podstawie analizy sąsiedztwa badanego punktu. Po zastosowaniu funkcji analizy stworzonego szkieletu (Skeleton Analyzer), generowany jest obraz, na którym zaznaczone są poszczególne elementy topologiczne szkieletu przy pomocy kodu kolorów przedstawionych w tabeli II.4 [MAVI 1.3.1 Documentation, Skeleton Analyzer]. Wynik przykładowej analizy szkieletu przedstawiono na rysunku II.12.. 40.

(41) Tabela II.4. Kod kolorów używany przez program MAVI przy wyświetleniu wyniku analizy szkieletu [MAVI 1.3.1 Documentation, Skeleton Analyzer]. klasa punktu 1 2 3 4 5 6 7 8. opis. kolor. punkt izolowany (I) punkt wewnątrz łuku (C) końcowy punkt łuku (CE) punkt skrzyżowania łuków (CC) punkt skrzyżowania płaszczyzny i łuku (SC) punkt wewnątrz płaszczyzny (S) punkt na krawędzi płaszczyzny (SE) punkt skrzyżowania płaszczyzn (SS). ciemnozielony granatowy oliwkowy fioletowy morski czerwony zielony niebieski. wartość zaczernienia 1 2 3 4 5 6 7 8. Rys. II.12. Analiza szkieletu zaprezentowanego na rysunku II.11. (b). Funkcja Skeleton Analyzer pozwala na uzyskanie informacji o budowie sieci porowej w postaci podania liczby występujących:  węzłów – punktów leżących na skrzyżowaniu łuków,  połączeń pomiędzy porami (gardzieli) – liczba pikseli stanowiących punkty wewnętrzne łuków,  zakończeń systemu porów – liczba pikseli będących zakończeniem łuku lub płaszczyzny. W przypadku analizy obrazu mikrotomograficznego skał, funkcja transformacji obrazu do szkieletu pozwala uzyskać informacje o przebiegu osi stanowiących środki przestrzeni porowych (zarówno porów jak i kanalików porowych). Na szkielecie obrazu można także łatwo zidentyfikować połączenia porów oraz wydzielić pory izolowane. Ze względu na bardzo duże obciążenie stacji roboczych przez algorytm tworzący drzewo topologiczne struktury porowej, dla tych obliczeń zdecydowano się zmniejszyć rozmiar analizowanego obrazu. 41.

(42) Tabela II.5 pokazuje zależność czasu pracy algorytmu CenterLineTree programu AVIZO (opis w rozdziale II.3.1) tworzącego drzewo topologiczne struktury porowej od rozmiaru obrazu. Tabela II.5. Wpływ rozmiaru obrazu na czas jego analizy. Rozmiar obrazu. Czas analizy. 900x900x400. Oprogramowanie zawiesiło się.. 600x600x400. 1,5 h. 300x300x400. 45 min.. 200x200x200. 25 min.. 100x100x100. 11 min.. II.3. Wybór reprezentatywnego fragmentu próbki Podstawowym problemem, tego etapu obliczeń, było wyznaczenie odpowiedniego, reprezentatywnego fragmentu (RF) obrazu do wyznaczenia przebiegu osi środkowych. W przypadku mikrotomografii konieczne jest, aby analizowany obraz był na tyle mały, by wystarczyło zasobów komputera do jego przetworzenia, a zarazem na tyle duży by reprezentował sobą cechy całego obiektu [Brun i in., 2010]. Do dalszych obliczeń zdecydowano się wybrać wycinek o rozmiarze 200 x 200 x 200 (2003) wokseli (1,12 x 1,12 x 1,12mm), w dalszej części opisano metodę wyboru takiego fragmentu. Dobór reprezentatywnego fragmentu jest problemem od dawna pojawiającym się w szeregu metod badawczych, nie tylko w tomografii. Jednym z pierwszych autorów rozważających wybór odpowiedniego rozmiaru dla RF w celu badania przepływu płynów w materiałach porowatych był Bear [1988]. Badania przeprowadzone przez Piller i in. [2009] polegały na określeniu wpływu rozmiaru obrazu na uzyskane wyniki symulacji przepuszczalności skał. Badaniu poddano wycinek próbki skały piaskowca o średnicy 5mm. Następnie dla kilku rozmiarów podpróbek – 3003 (1,35mm)3, 5003. (2,25mm)3. i. 6003 (2,7mm)3 obliczono wartości. współczynnika. przepuszczalności. Uzyskane wyniki wyniosły odpowiednio: 409, 250 i 260mD. Otrzymane dane pokazały zdecydowaną różnicę pomiędzy wartościami dla próbek o rozmiarze 3003 i 5003, natomiast wartości dla obrazów 5003 i 6003 były już do siebie zbliżone. Jak podkreślili autorzy znaczne różnice w otrzymanych wartościach mogły wynikać, albo ze zbyt małej 42.

(43) wielkości reprezentatywnego fragmentu lub dużej nierównomierności rozkładu porowatości badanej próbki. Interesującą metodykę służącą wyznaczaniu rozmiaru reprezentatywnego fragmentu proponują Okabe i Oseto [2006]. W opisywanym przypadku, w celu określenia wpływu rozmiaru RF oraz niejednorodności porowatości na otrzymane wyniki badań, zastosowali numeryczną symulację przepływów metodą LBM (Lattice-Boltzman Method). Sam opis metody, która nie była przedmiotem pracy, można znaleźć w następującej literaturze: [Chen, Doolen, 1998 i Rothman, Zaleski, 1997]. W opisywanym przypadku wejściowy obraz posiadał rozmiar 2403 wokseli, co przy rozdzielczości pomiaru 5,67um dało rozmiar rzeczywisty próbki 1,36 x 1,36 x 1,36mm3. Na poniższym rysunku przedstawiono zakres uzyskiwanych wartości współczynnika porowatości (II.13a) oraz przepuszczalności (II.13b) w zależności od zastosowanego rozmiaru RF.. a). b). Rys. II.13. Wartości współczynnika porowatości (a) i przepuszczalności (b) uzyskane dla podpróbek o różnych rozmiarach [Okabe, Oseto, 2006]. W doświadczeniu analizowane były sześciany o boku 60 (340,2um), 120 (680,4um) i 180 (1002,6um) wokseli. Powyższa ilustracja otrzymanych zależności wyraźnie pokazuje, że wraz ze wzrostem rozmiaru obiektu maleje rozrzut otrzymywanych wyników. Powyższe rozważania dają wyraźną i spodziewaną odpowiedź, że im większy rozmiar reprezentatywnego fragmentu, tym dokładniejsze wartości otrzymywanych wyników. Zastanawiający jest jedynie fakt, szczególnie w przypadku niedawnych publikacji [Piller i in., 2009, Al-Raoush, Papadopoulos, 2010], dlaczego autorzy nie rozważają wpływu geometrii 43.

(44) struktury porowej, przy wyborze rozmiaru RF. Nałożenie pewnych warunków, jak chociażby wybranie wycinka o zbliżonej porowatości do porowatości całej próbki, powinno pozwolić uzyskać bardziej zbliżone wyniki, przy zachowaniu zdecydowanie mniejszego rozmiaru RF. Ponieważ podczas realizacji tej pracy nie istniała możliwość przeprowadzenia tak skomplikowanych obliczeń (symulacji) w celu wybrania RF, w opisywanych badaniach zdecydowano się skorzystać jedynie z parametrów obrazu, charakteryzujących geometrię struktury porowej. Jako najważniejszy parametr zdecydowano się traktować współczynnik porowatości (Kp), następnie, zgodnie ze spostrzeżeniami Hansa Vogla [Vogel, 2002], wartość charakterystyki Eulera (Eu) oraz długość średniej cięciwy (Ch). Każdemu z parametrów przypisano odpowiednio wagi 0,5, 0,3 i 0,2. Wyznaczenie odpowiedniego fragmentu próbki przebiegało w następujący sposób: dla całego obrazu struktury porowej (900 x 900 x 400 wokseli) wyznaczono sumę według poniższego wzoru: S = 0,5·Kppróbki + 0,3·Eupróbki + 0,2·Chpróbki. (II.6). analogiczne obliczenia wykonano dla każdej z analizowanej podpróbki o rozmiarze (200 x 200 x 200). Za reprezentatywny fragment uznano ten, dla którego różnica pomiędzy wartością wyliczoną dla całego obrazu oraz podpróbki była najmniejsza. Poniżej (tabela II.6) przedstawiono przykład takich obliczeń dla próbki numer 12915. Z przeprowadzonych porównań, wynika, że najbardziej reprezentatywnym fragmentem dla próbki 12915 jest sześcian o rozmiarze 200 x 200 x 200 wokseli, którego początek w układzie przestrzennym leży w punkcie o współrzędnych (400,000,000).. 44.

(45) Tabela II.6. Przedstawienie obliczonych wartości sum wg wzoru (II.6) i różnic w celu wyznaczenia reprezentatywnego fragmentu dla próbki 12915. badany obszar cała próbka (000,000,000) (000,000,200) (000,200,000) (000,200,200) (000,400,000) (000,400,200) (000,600,000) (000,600,200) (200,000,000) (200,000,200) (200,200,000) (200,200,200) (200,400,000) (200,400,200) (200,600,000) (200,600,200) (400,000,000) (400,000,200) (400,200,000) (400,200,200) (400,400,000) (400,400,200) (400,600,000) (400,600,200) (600,000,000) (600,000,200) (600,200,000) (600,200,200) (600,400,000) (600,400,200) (600,600,000) (600,600,200). wartość sumy 4,5196 4,2704 4,1754 4,4586 4,1065 4,5815 4,1641 4,4181 4,3350 4,3029 4,2764 4,8296 4,6689 4,7123 4,7306 4,5503 4,4535 4,4972 4,3423 4,8437 4,6148 4,8352 4,7676 4,6315 4,4307 4,5635 4,4900 4,5767 4,5733 4,6017 4,4279 4,7270 4,5716. 45. różnica sum 0,2492 0,3442 0,0610 0,4132 0,0619 0,3555 0,1015 0,1846 0,2167 0,2432 0,3099 0,1493 0,1926 0,2110 0,0307 0,0661 0,0224 0,1773 0,3241 0,0952 0,3156 0,2480 0,1119 0,0889 0,0439 0,0297 0,0571 0,0537 0,0821 0,0917 0,2074 0,0520.

(46) II.3.1 Przykład zastosowania W niniejszej pracy do uzyskania szkieletu struktury porowej zastosowano moduł CenterLineTree programu AVIZO [Avizo Help, CenterLineTree], który oparty jest na algorytmie TEASAR (tree-structure extraction algorithm for accurate and robust skeletons) [Sato i in, 2000]. Moduł ten generuje drzewo topologiczne analizowanej struktury. Od typowej transformacji obrazu na szkielet, z użyciem osi środkowej, różni się tym, iż nie generuje pętli. Bezpośrednim wynikiem pracy modułu CenterLineTree jest uzyskanie obrazu typu SpatialGraph (rys. II.14), składającego się z linii, zakończeń oraz węzłów. Moduł ten zachowuje informacje o odległości każdego punktu otrzymanego drzewa topologicznego od krawędzi obiektu w obrazie źródłowym. Informacja ta może zostać wykorzystana do określenia lokalnej grubości struktury porowej, którą można przedstawić za pomocą skali kolorów (rys. II.15). W tym przypadku zmiana koloru od fioletowego do pomarańczowego odpowiada względnemu wzrostowi średnicy kanalików porowych.. Rys. II.14. Drzewo topologiczne struktury porów. 46.

(47) min. max. Rys. II.15. Drzewo topologiczne struktury porów przedstawiające informację o grubości kanalików porowych. Po wczytaniu do programu MAVI wygenerowanego drzewa punktów można poddać go analizie w module Skeleton Analyzer. W wyniku obliczeń uzyskamy informacje (rys. II.16), które punkty należą do gałęzi (kolor niebieski), węzłów (kolor fioletowy) lub zakończeń porów (kolor oliwkowy). Rysunek II.16 przedstawia przekrój przez próbkę z zaznaczonymi poszczególnymi punktami, których kolorystyka również zdefiniowana jest w tabeli II.4.. 47.

(48) Rys. II.16. Przekrój obrazu po analizie szkieletu. Po wczytaniu obrazu, przeanalizowanego za pomocą modułu Skeleton Analyzer, do programu AVIZO (rys.II.17) można przedstawić rozkład odcinków z zachowaniem wcześniej zdefiniowanych kolorów.. Rys. II.17. Przestrzenna dystrybucja zakończeń, węzłów i kanalików. 48.

(49) Po usunięciu ze struktury przedstawionej na rysunku II.17 wszystkich węzłów, pozostaną w niej same niepołączone odcinki, odwzorowujące długość kanałów porowych. Rysunek II.18 przedstawia histogram długości kanałów porowych występujących w analizowanej próbce.. Rys. II.18. Histogram długości kanalików porowych. Dla analizowanej próbki otrzymano następujące dane: . liczba kanalików: 131,. . liczba zakończeń: 100,. . liczba węzłów: 66.. Sumaryczna długość kanalików wyniosła 2113 jednostek (j), co daje średnią długość kanalików na poziomie 16,1j, przy boku woksela 5,8 μm rzeczywista długość wynosi 93,4μm. Dla tak uzyskanych danych możliwe jest także obliczenie średniej liczby koordynacyjnej (mówiącej ile średnio kanałów łączy się w węźle) według teoretycznego wzoru:. LK . 2  Lk  LZ LW. gdzie, LK - liczba kanalików, LZ - liczba zakończeń, LW - liczba węzłów. 49. (II.7).

(50) Dla powyższych danych : LK . 2 131  100  2,45 66. co daje średnią liczbę koordynacyjną 2,45. Ridgway i Tarbuck [1967] w swojej pracy podają empiryczny wzór opisujący relację pomiędzy liczbą koordynacyjną, a porowatością dla układu szklanych kulek: 1 – Kp = 1,072 – 0,11937z + 0,0043z2. (II.8). gdzie, Kp – porowatość badanej próbki, z – średnia liczba koordynacyjna. Natomiast Dullien [1992], prowadząc badania dla próbki piaskowca Berea określił wartość liczby koordynacyjnej na poziomie 2,8. Wartość ta okazała się być bardzo zbliżoną do wartości wyliczonej dla tej próbki ze wzoru (II.8), która wyniosła 2,9. Z powyższych badań wynika, iż zależność (II.8) może też mieć zastosowanie dla badania struktury porowej niektórych skał.. 50.

(51) ROZDZIAŁ III – Opis obszaru badań III.1. Obszar badań Polski basen permski (czerwony spągowiec i cechsztyn) stanowił wschodnią część południowego basenu permskiego obejmującego swym zasięgiem zachodnią i centralną Europę [Karnkowski, 1999]. Na rysunku III.1.A. przedstawiono przykład profilu poprowadzonego przez utwory czerwonego spągowca od obszarów centralnej Polski, przez Niemcy, Holandię i na terytorium Wielkiej Brytanii kończąc [Doornenbal, Stevenson, 2010]. Profil przedstawia miąższości utworów czerwonego spągowca występujące na terytorium poszczególnych krajów, a kolorami oznaczono jego poszczególne formacje. W czerwonej ramce po prawej stronie rysunku zaznaczono część profilu poprowadzonego przez terytorium Polski, którego powiększenie pokazane jest na rysunku III.1.B. Basen czerwonego spągowca był kontynentalnym basenem aluwialnym o asymetrycznym rozkładzie facji osadowych. Na kształt basenu zasadniczy wpływ miał ekstensyjny reżim tektoniczny oraz budowa geologiczna podłoża podpermskiego [Pokorski, 1997]. Basen ma cechy półrowu tektonicznego z największą subsydencją wzdłuż północno-wschodniej krawędzi.. 51.

(52) 52 Rys. III.1. A. Profil czerwonego spągowca [Doornenbal, Stevenson, 2010].

(53) 53 Rys. III.1.B. Fragment profilu czerwonego spągowca poprowadzonego przez terytorium Polski [Doornenbal, Stevenson, 2010].