Tomografia komputerowa

1. Tomografia impedancyjna

Tomografia Impedancyjna TI (ang. Electrical Impedance Tomography - EIT) jest techniką obrazowania, w której wykorzystuje się właściwości elektryczne materiałów, w tym również tkanek biologicznych. W metodzie tej badany obiekt pobudzany jest ze źródła prądowego lub napięciowego a następnie obserwuje się powstały na jego brzegu rozkład napięć. Zebrane informacje przetwarzane są za pomocą algorytmu, który konstruuje obraz badanego obiektu.
Tomografia impedancyjna jest metodą stosunkowo nową, szybko rozwijającą się. W ostatnich latach znajduje zastosowanie w różnych dziedzinach, poczynając od przemysłu rafineryjnego, poprzez przemysł elektroniczny aż do biomedycyny.
Impedancyjna tomografia komputerowa znalazła zastosowania w:
  • medycynie,
  • defektoskopii wiroprądowej,
  • wykrywaniu i kontroli laminarnych i turbulętnych przepływów cieczy w kanałach zamkniętych (głównie przemysł rafineryjny),
  • badaniu stanu pni drzew (właściwości elektryczne zdrowego i przemarzniętego pnia drzewa są różne),
  • archeologii i geofizyce (skażenie wód gruntowych),
  • przemyśle hutniczym.
Zaletami impedancyjnej tomografii komputerowej są:
  • nieinwazyjne i nieniszczące badanie danego obiektu,
  • możliwość obrazowania małych zmian konduktywności,
  • niskie koszty urządzenia pomiarowego i jego eksploatacji,
  • możliwość przenoszenia aparatury.

W TI rozwiązanie zagadnienia prostego (ang. forward problem) polega na wyznaczeniu rozkładu potencjału wewnątrz obiektu, przy zadanych warunkach brzegowych i danych dotyczących rozpatrywanego obiektu. Z matematycznego punktu widzenia tomografia impedancyjna tomografia należy także do zagadnień odwrotnych pola elektromagnetycznego. Zagadnieniem odwrotnym pola elektromagnetycznego (ang. inverse problem) nazywamy proces identyfikacji, optymalizacji, bądź syntezy, w którym zmierza się do wyznaczenia parametrów opisujących dane pole, na podstawie posiadania niektórych informacji charakterystycznych dla tego pola. Zagadnienia odwrotne są trudne do analizy. Z reguły nie mają jednoznacznych rozwiązań i są źle uwarunkowane. Przyczyną tego jest często zbyt mała lub zbyt duża liczba informacji, które czasem są ze sobą sprzeczne bądź liniowo zależne.
Zadanie algorytmu rekonstrukcji (rys. 1) polega na wyznaczeniu takiego wektora σ (konduktywność), dla którego funkcja celu osiąga przynajmniej minimum lokalne. Wartość funkcji celu przy ustalonych danych pomiarowych zależy od macierzy stanu. Na podstawie precyzyjnie rozwiązanego zagadnienia prostego wyznacza się wartości potencjałów we wszystkich węzłach siatki przy założeniu, że znany jest rozkład konduktywności w obszarze. Kolejnym zadaniem jest analiza wrażliwościowa. Określa ona pochodną funkcji celu względem zmiennych decyzyjnych (względem wektora σ rozkładu konduktywności elektrycznej). Analiza wrażliwościowa w połączeniu z dowolną metodą optymalizacyjną daje odpowiedź na pytanie jakich należy dokonać poprawek zmiennych decyzyjnych, aby wartość funkcji celu uległa zmniejszeniu.


Rys. 1. Schemat algorytmu konstrukcji obrazu w TI,


W TI gęstość prądu elektrycznego, przepływającego przez dany ośrodek, rozkłada się w obszarze nierównomiernie i wybiera drogę o mniejszej rezystancji. Oznacza to, że nie płynie w jednej płaszczyźnie, ani nie płynie po liniach prostych. Przy rekonstrukcji obrazu analizowany jest cały obszar, bez możliwości ograniczenia strefy badań do fragmentu pewnej całości.

2. Algorytmy numeryczne

W impedancyjnej tomografii komputerowej można spotkać dwie grupy rozwiązania problemu rekonstrukcji obrazu. Do pierwszej z nich należą procedury jednoiteracyjne. Opierają się one na założeniu, że w przypadku niewielkich różnic między rozkładem rzeczywistym a rozkładem jednorodnym współczynników materiałowych ośrodka, transformację równania można aproksymować zależnością liniową. Wówczas transformacja odwrotna wymaga tylko jednego kroku obliczeniowego. Zaletą tego podejścia jest możliwość uzyskania obrazu w krótkim czasie. Do wad natomiast należy zaliczyć możliwość obrazowania tylko zmian konduktywności w stosunku do rozkładu jednorodnego. W drugiej grupie algorytmów rekonstrukcji obrazu znajdują się procedury iteracyjne. W każdym kroku procesu iteracyjnego równanie nieliniowe sprowadzane jest do problemu liniowego. Algorytmy te, w porównaniu do technik jednoiteracyjnych, pozwalają na uzyskanie lepszej rozdzielczości obrazu, jednak czas wykonywania obliczeń numerycznych jest znacznie dłuższy.
Metody rekonstrukcji obrazu stosowane w tomografii impedancyjnej można podzielić następująco:

A. Metody deterministyczne
  • metoda wstecznej projekcji,
  • metoda perturbacyjna,
  • metoda warstwowa,
  • algorytm Kaczmarza,
  • algorytm Wexlera,
  • gradientowe metody optymalizacyjne:
    • metoda Newtona-Rapsona,
    • metoda największego spadku,
    • metoda gradientów sprzężonych,
    • funkcje dzwonowe.
  • metody zbiorów poziomicowych,
B. Metody stochastyczne
  • algorytmy genetyczne (GA),
  • symulowane wyżarzanie (SA),
  • metoda Monte-Carlo (MC).
C. Sztuczne sieci neuronowe.
D. Hybrydowe metody optymalizacji.


Metody deterministyczne charakteryzują się tym, że poszukują minimum lokalnego funkcji celu, ciągłej w danym przedziale zmienności. Jeżeli mamy podane wartości początkowe zmiennych projektowych, obliczenia optymalizacyjne wprowadzają zmiany tych wartości, dążąc do znalezienia minimum lokalnego, a następnie minimum globalnego.
Metody stochastyczne korzystają z wyznaczonych funkcji celu tylko w testowych punktach badanego obszaru, bez nakładania dodatkowych wymagań na ciągłość i różniczkowalność funkcji. Dodatkowo poszukiwane jest minimum globalne funkcji celu w określonym zbiorze dopuszczalnych wartości zmiennych projektowych. Algorytmy numeryczne dla metod stochastycznych są względnie proste, ale wymagają bardzo dużej liczby obliczeń, a przez to wydłuża się czas konstrukcji obrazu.
Zastosowanie sztucznych sieci neuronowych ma na celu uzyskanie oczekiwanego rozwiązania w bardzo krótkim czasie, praktycznie „on-line”. Metoda ta polega na wyborze odpowiedniej sieci neuronowej i zgromadzeniu danych do jej trenowania. SSN zachowują się jak układy aproksymujące i mają pewną zdolność do uogólniania informacji, którymi są uczone. W procesie uczenia sieci przygotowywane są próbki do trenowania, obejmujące całą badaną powierzchnię. Zdolność zlokalizowania obiektu w miejscu, które było użyte w procesie trenowania sieci, jest miarą zdolności sieci do zapamiętywania danych uczących. Miarą zdolności generalizacji sieci jest możliwość rozpoznawania położenia obiektu, którego sieć wcześniej nie widziała. Zwraca się tutaj szczególną uwagę na to, aby sieć neuronowa nie dopasowała się za bardzo do próbek uczących, ponieważ mogłaby wtedy stracić właściwości generalizacji rozwiązań, co zwiększałoby błędy rekonstrukcji obrazu. Z drugiej strony nie można stosować zbyt ubogiej struktury sieci, aby aproksymacja nie była zbyt pobieżna, a przez to niedokładna (zjawisko niedouczenia sieci).
Metoda zbiorów poziomicowych (MZP) jest metodą umożliwiającą topologiczną zmianę kształtów figur oraz pozwalającą optymalizować wymiary i formę obiektów. Została ona przedstawiona przez Oshera i Sethiana. Podstawową zaletą metody jest możliwość wykonywania obliczeń numerycznych związanych z krzywymi lub płaszczyznami w układzie kartezjańskim bez konieczności ich parametryzowania. W rozważaniach metody zbiorów poziomicowych uwzględniany jest pewien obszar posiadający ruchomy brzeg Γ. Brzeg ten porusza się z wyznaczoną w kolejnych krokach prędkością , która może zależeć od pozycji, czasu, kształtu brzegu oraz warunków zewnętrznych. Idea zbiorów poziomicowych polega na zdefiniowaniu funkcji ϕ(x,t), która reprezentuje poruszający się brzeg. Metodę można z dużym powodzeniem zastosować do segmentacji obrazów otrzymanych za pomocą innych algorytmów i metod.
Hybrydowe metody optymalizacji łączą w sobie cechy metod deterministycznych i stochastycznych. Zazwyczaj w pierwszej fazie obliczeń algorytm stochastyczny lokalizuje położenie globalnego minimum funkcji celu. Gdy w pewnym momencie obliczeń efektywność tego algorytmu zaczyna spadać, wprowadza się algorytm deterministyczny, który w niewielkiej liczbie iteracji znajduje rozwiązanie. W metodach hybrydowych jako algorytmy stochastyczne stosuje się często algorytmy genetyczne oraz metodę symulowanego wyżarzania.
Do rozwiązanie zagadnienia prostego można użyć m.in.:
  • metodę elementów skończonych,
  • metodę granicy podobszarów,
  • metodę elementów brzegowych.
Metoda elementów skończonych MES (ang. Finite Element Method - FEM) jest obecnie najbardziej efektywną metodą uzyskiwania przybliżonego rozwiązania numerycznego, głównie dla zagadnień płaskich i przestrzennych o skomplikowanej geometrii, w których środowisko może wykazywać cechy nieliniowości. Rozwiązując takie zagadnienie analitycznie przyjmuje się wiele uproszczeń. W metodzie tej obszar dzielimy na skończoną liczbę elementów. Każdy element ma węzły, w których są szukane wielkości polowe. Węzły te są rozmieszczone najczęściej na bokach i narożnikach elementów w ten sposób, że dany węzeł, a z nim i jego wielkości polowe są wspólne dla dwóch lub większej liczby sąsiednich elementów. Zachowanie nieznanych wielkości pola nad każdym elementem jest aproksymowane funkcjami ciągłymi. Funkcje aproksymacji są zdefiniowane w odniesieniu do jednego elementu.
Metoda elementów brzegowych (MEB) (ang. Boundary Element Method - BEM) polega na wykorzystaniu rozwiązania fundamentalnego danego równania różniczkowego. Zaletą metody jest to, że wymaga ona aproksymacji zjawiska tylko na brzegu obszaru i tylko do tego ogranicza się poprowadzenie siatki węzłów i elementów. Prowadzi to do redukcji wymiaru zagadnienia o jeden. Ponadto uzyskanie dokładnych wartości poszukiwanego rozwiązania wewnątrz obszaru nie wymaga lokalnego zagęszczenia siatki na brzegu. Wadą metody jest znajomość i konieczność istnienia rozwiązania fundamentalnego równania różniczkowego. Kolejną niedogodnością jest generowanie przez MEB układu liniowych równań algebraicznych mających pełną i niesymetryczną macierz współczynników. Metoda elementów brzegowych wydaje się być obiecującą do rozwiązania zagadnienia prostego i odwrotnego tomografii impedancyjnej ze względu na mniejszą zajętość pamięci niż MES. Poza tym MEB używa mniejszej ilości elementów i węzłów dla osiągnięcia dokładności tego samego rzędu. Konieczność uwzględniania takich elementów modelu matematycznego, wymusza bardzo istotne zagęszczanie sieci elementów skończonych, co nawet przy obecnym rozwoju techniki komputerowej prowadzi do gigantycznych rozmiarów układów równań liniowych, które w procesie konstrukcji obrazu musza być rozwiązywane po wielokroć. Należy tu jeszcze wspomnieć o kłopotach z dyskretyzacją takich obszarów i uzyskaniem wysokiej jakości sieci elementów skończonych. Te problemy, choć nastąpił w tej dziedzinie ogromny postęp w minionym dziesięcioleciu, ciągle jeszcze są otwarte. Za największą zaletę w metodzie elementów brzegowych uznawane jest ograniczenie rozkładu sieci obiektu do dyskretyzacji jedynie na jego powierzchni i dlatego obserwuje się coraz większe zainteresowanie tą metodą w tomografii.
Metoda granicy podobszarów MGP (ang. Immersed Interface Method - IIM) jest metodą opartą na siatce kartezjańskiej. Główną jej ideą jest zastosowanie warunku skoku na brzegu rozdzielającym sąsiadujące obszary. Li i LeVeque określili lokalne właściwości rozwiązania i wyrazili je w warunku skoku w tej metodzie. Obliczyli i zmodyfikowali współczynniki różnic skończonych umieszczonych blisko brzegu rozpatrywanego obiektu. W metodzie granicy podobszarów zostało zaimplementowane równanie Poissona do całego badanego obszaru niezależnie od kształtu nieregularnego obiektu.

3. Metody pomiarowe

Opracowano dwie metody pozyskiwania informacji potrzebnych do budowy obrazu. Pierwsza z nich przedstawiona przez Barbera i Browna z Uniwersytetu w Sheffield polega na równomiernym rozmieszczeniu wokół badanego ciała 16 elektrod. Za pomocą pojedynczego generatora doprowadzone jest zasilanie do dwóch sąsiednich elektrod, natomiast na pozostałych elektrodach mierzone jest napięcie. Następnie napięcie podawane jest na kolejne elektrody i cykl pomiarowy zostaje powtórzony. Zaletą tej metody jest prosta konstrukcja i łatwy sposób akwizycji danych, wadą natomiast jest to, że nie można odwzorować bezwzględnych wartości parametrów elektrycznych, a jedynie zmiany tych parametrów. Metoda ta znalazła zastosowanie w gastroskopii do badania prędkości opróżniania treści żołądkowych, przepływu krwi w tułowiu, w głowie lub w ramieniu. Zespół inżynierów z Rensselaner Polytechnic Institute opracował dwa układy, 32 i 64 elektrodowy. Każdy z nich ma swój własny programowany generator. Mimo większego stopnia skomplikowania w stosunku do pierwszego rozwiązania, uzyskano lepszy stosunek sygnału do szumu oraz zmniejszenie wrażliwości na niedokładność położenia elektrod pomiarowych wokół badanego obiektu. Dzięki temu możliwe stały się badania tułowia oraz obserwacje zmian nowotworowych w piersiach. W obszarze 2D do symulacji komputerowej wykorzystuje się model zaproponowany przez Yorkeya. Jest to dwuwymiarowy obszar prostokątny, którego konduktywność tła wynosi 1 S/m, z umieszczonym wewnątrz obiektem o konduktywności równej 3 S/m. Wewnątrz obszaru wygenerowano regularną sieć elementów skończonych. W wyniku tego postępowania otrzymano obszar złożony z prostokątnych (w szczególnym przypadku kwadratowych) elementów o stosunkowo niewielkich rozmiarach. W jednakowych odstępach dookoła obszaru rozmieszczono n = 16 elektrod. Ich liczba uwarunkowana jest względami praktycznymi i sposobem akwizycji danych. Do dwóch (rozmieszczonych biegunowo lub sąsiednich) elektrod dołączane jest źródło pobudzające. Pozostałe elektrody służą do pomiarów napięcia na brzegu obszaru. Pomiary wykonywane są dla wszystkich możliwych sposobów podłączenia źródła zasilania, w celu zwiększenia liczby informacji o obiekcie oraz poprawy stosunku sygnału do szumu. Po wykonaniu pierwszej serii pomiarów następuje przełączenie układu pobudzającego na elektrody sąsiadujące. Proces pomiarowy powtarzany jest sekwencyjnie dla wszystkich możliwych układów połączeń źródła pobudzającego. W ten sposób dokonywane jest wielokrotne „prześwietlanie” badanego obiektu. Ze względu na symetrię układu dla n = 16 elektrod można uzyskać n / 2 = 8 niezależnych rozkładów prądowych. Niech każda konfiguracja źródła pobudzającego nosi nazwę kąta projekcji (ang. projection angle), w takim przypadku całkowita liczba tych kątów przy powyższych założeniach wynosi 8. Danymi wejściowymi dla algorytmu konstrukcji obrazu są pomiary napięć dokonane między sąsiednimi elektrodami. Pomiary wykonane na elektrodach z dołączonym źródłem pobudzającym są pomijane ze względu na nieznaczny spadek napięcia występujący między tymi elektrodami a badanym obszarem. Dla układu n = 16 elektrod oraz dowolnego kąta projekcji można otrzymać n − 3 = 13 niezależnych pomiarów. Stąd pełna liczba możliwych do uzyskania niezależnych pomiarów napięć pomiędzy sąsiednimi elektrodami napięciowymi przy n / 2 = 8 kątach wynosi: (n − 3)(n / 2) = 13 • 8 = 104.
Wartości potencjałów elektrod zależą od rozpływu prądu wewnątrz obszaru, a tym samym również od rozkładu konduktywności. Algorytm komputerowej rekonstrukcji obrazu poszukuje iteracyjnie takiego rozkładu konduktywności, dla którego obliczone wartości napięć międzyelektrodowych są możliwie najbliższe odpowiednim wartościom pomiarowym. Wartości początkowe konduktywności w punkcie wyjściowym procesu iteracyjnego są dobierane doświadczalnie, aby otrzymać możliwie najlepsze rozwiązanie. Zwykle przyjmuje się rozkład jednorodny. Systemy pomiarowe mają umożliwiać wykonywanie pomiarów z wykorzystaniem komputera i podłączonych do niego przyrządów oraz przetwarzanie danych pomiarowych i prezentowanie wyników w określonej postaci. System pomiarowy można przedstawić w postaci schematu blokowego, w którym operuje się pojęciami bloków funkcjonalnych realizujących określone zadania.

4. Badania i symulacje

Badanie stanu pni drzew

Do pomiarów wykorzystana została brzoza o średnicy 78 cm na wysokości 130 cm od ziemi (wysokość całkowita brzozy około 800cm). Na tej wysokości zostało umieszczonych równomiernie 16 elektrod w postaci pozłacanych gwoździ o średnicy 1,5mm i długości 15mm (10mm wbite w drzewo). Układ zasilany był napięciem sinusoidalnym 3,900V o częstotliwości 1,000kHz z generatora HP 33120A, pomiar precyzyjnym multimetrem firmy Fluke 8845A. Elektrody względem stron świata zostały rozmieszczone w taki sposób, że 7 była skierowana na kierunku północ, a 15 – na południe



Zdjęcia wykonywanych pomiarów podczas badań tomograficznych brzozy


Rozwiązując zadanie odwrotne należy wielokrotnie rozwiązać zadanie proste tak, aby w procesie iteracyjnym otrzymać rozkład napięć maksymalnie zbliżony do uzyskanego z pomiarów. Innymi słowy w procesie iteracyjnym szukane jest minimum funkcji celu określonej jako suma kwadratów różnic napięć obliczanych w kolejnych krokach oraz napięć pomierzonych.


Rozkład potencjałów dla różnych kątów projekcji



Wyniki rekonstrukcji obrazu z wykorzystaniem metody zbiorów poziomicowych


Diagnostyka medyczna

W badaniach medycznych zauważono, że między tkankami zdrowymi, a zmienionymi chorobowo występują różnice właściwości elektrycznych. W tabeli podano konduktywności i przenikalności elektryczne wybranych tkanek dla dwóch częstotliwości. Z porównania wynika, że różne tkanki mają inne konduktywności dla wybranych częstotliwości.

Tomografię impedancyjną można wykorzystać w następujących dziedzinach diagnostyki medycznej:
  • w gastroskopii,
  • w badaniach wentylacji płuc.
  • diagnostyka kobiet ciężarnych,
  • badanie aktywności kory mózgowej,
  • kardiografia impedancyjna,
  • diagnostyka oftalmologiczna,
  • badanie rozkładu temperatury wewnątrz ciała.


Tabela: Właściwości elektryczne wybranych tkanek


 


Symulacja rekonstrukcji obrazu w tomografii płuc za pomocą metody zbiorów poziomicowych


Zawilgocenie murów i wały przeciwpowodziowe

Skuteczność rekonstrukcji obrazu z wykorzystaniem metody zbiorów poziomicowych zostanie sprawdzona w metodzie badań nieniszczących stosowanych do oceny stanu zawilgocenia murów w budynkach. Omawiana metoda opiera się na teorii pola elektrycznego i pokazuje, jak szeroki zakres mogą mieć zastosowania elektrotechniki w różnych dziedzinach nauk. Podciągana kapilarnie z gruntu wilgoć jest problemem w budynkach, szczególnie starszych, z niewystarczającą izolacją poziomą i pionową ścian fundamentowych. Wilgoć jest podciągana z miejsc bezpośredniego styku miedzy murem a glebą. Budowlane materiały porowate (np. cegła czy beton), zarówno naturalne, jak i produkowane, mają pory, podobnie jak gąbka, i wilgoć może być podciągana wbrew grawitacji nawet na wysokość 3 metrów. Dodatkowo niszczące działanie wilgoci wzmagają zawarte w wodzie gruntowej szkodliwe dla murów zanieczyszczenia: roztwory soli, zasad czy kwasy organiczne. Spada wytrzymałość muru, zawarta w nim woda zamarzając zimą i zwiększając swoją objętość może prowadzić do pękania murów, podobnie jak krystalizująca się sól. Wilgoć stwarza niebezpieczeństwo nie tylko dla murów, ale również dla zdrowia ludzi. Po pierwsze sprzyja ona powstawaniu schorzeń reumatycznych, po drugie – tworzące się na powierzchni wilgotnej ściany grzyby – mogą wywoływać alergie oraz wiele innych chorób. Istnieje wiele różnych systemów osuszania. Bardzo ważne jest ciągłe kontrolowanie wilgotności muru w trakcie procesu osuszania. Można je przeprowadzić stosując metody elektryczne, ponieważ konduktywność muru zmienia się wraz z jego wilgotnością i zasoleniem, na przykład beton zmienia swą rezystywność od 0,03 do 0,06 Ωm dla 100 % nasycenia wodą aż po 0,1 do 60 Ωm dla 40 % nasycenia. Tomografia impedancyjna może być zastosowana do określenia rozkładu konduktywności wewnątrz badanego obiektu.


System pomiarowy – model rzeczywisty muru

Zadanie odwrotne pozwala na określenie położenia granicy pomiędzy obszarami na podstawie danych pomiarowych z elektrod. Przyjęto następujące wartości konduktywności do obliczeń: 10,6 mS/m — mokrego gazobetonu, 69 µS/m – suchego betonu. Na rysunku przedstawiono, na podstawie danych symulacyjnych, szukane położenie granicy między obszarami dla dwóch różnych sposobów przeprowadzania pomiarów.


Symulacja rekonstrukcji obrazu