Level Set -menetelmä impedanssitomografiassa
Antti Voss Filosofian maisterin tutkielma Fysiikan koulutusohjelma Itä-Suomen yliopisto, Sovelletun fysiikan laitos
ITÄ-SUOMEN YLIOPISTO, Luonnontieteiden ja metsätieteiden tiedekunta Fysiikan koulutusohjelma, laskennallinen ja teknillinen fysiikka
Antti Voss
Filosofian maisterin tutkielma, 59 sivua
Tutkielman ohjaajat: Professori Marko Vauhkonen
Apulaisprofessori Ville Kolehmainen kesäkuu 2014
Tiivistelmä
Impedanssitomografiassa (Electrical Impedance Tomography, EIT) kuvannetta- vaan kohteeseen syötetään heikkoa vaihtovirtaa kohteen reunoille kiinnitettyjen elekt- rodien kautta. Virransyötön seurauksena syntyvät jännitteet mitataan käyttäen sa- moja elekrodeja ja kohteen sisäinen johtavuusjakauma voidaan rekonstruoida näiden reunalta tehtyjen jännitemittausten perusteella. Johtavuusjakauman rekonstruoimi- nen on matemaattisesti epälineaarinen ja huonostiasetettu käänteisongelma.
Tässä tutkielmassa impedanssitomografian rekonstruktio-ongelman ratkaisemi- seksi sovelletaan tasa-arvojoukkojen menetelmää eli niin kutsuttua level set - menetelmää. Rekonstruktio-ongelman ratkaisu perustuu yleiseen Tikhonov regula- risointimenetelmään. Tutkielman lähestymistavassa pääpiirre on, että kohteen johta- vuus oletetaan paloittain vakioksi sekä eri johtavuuksien alueet ja niiden väliset raja- pinnat esitetään level set -formuloinnin avulla. Numeeriset testit osoittavat, että le- vel set -lähestymistavalla kohteessa olevien geometrialtaan hyvinkin monimutkaisten epähomogeenisuuksien paikat ja muodot saadaan paikannettua hyvällä tarkkuudella.
Myös eri alueiden johtavuudet voidaan estimoida menetelmän avulla tarkasti. Mene- telmä on melko riippumaton asetettavasta alkuarvauksesta ja se sietää kohtuullisen hyvin häiriötä mittausdatassa.
Avainsanat: Tasa-arvojoukkojen (Level Set) menetelmä, Impedanssitomografia, Täy- dellinen elektrodimalli (CEM), Tikhonov regularisointi, Äärellisten elementtien me- netelmä (FEM), Gauss-Newton algoritmi
LYHENTEET
EIT Impedanssitomografia
LS Pienimmän neliösumman menetelmä CEM Täydellinen elekrodimalli
FEM Äärellisten elementtien menetelmä
SYMBOLIT
Ω Kuvattava kohde
∂Ω Kuvattavan kohteen reuna el Elektrodi l
D Sähkövuon tiheys B Magneettinen induktio E Sähkökenttä
H Magneettikenttä ρp Varaustiheys J Virrantiheys Js Lähdevirta Jo Ohminen virta
Sähköinen permittiivisyys µ Magneettinen permeabiliteetti σ Sähköjohtavuus
Il Elektrodille l syötettävä virta
zl Elektrodin l ja kohteen välinen kontakti-impedanssi Ul Elektrodilla l mitattu jännite
H1(Ω) Hilbert-avaruus
Qh Avaruuden H1(Ω) äärellisulotteinen aliavaruus uh FEM-approksimaatio kohteen potentiaalijakaumalle Uh FEM-approksimaatio elektrodijännitteille
φ Level set -funktio
Λ Level set -funktion määrittämän alue Γ Level set -funktion nollakäyrä
d(x) Etäisyysfunktio H(x) Askelfunktio δ(x) Deltafunktio
1,2 Askel- ja deltafunktion approksimaatioiden sileysparametrit
Sisältö
1 Johdanto 5
2 Impedanssitomografia (EIT) 8
2.1 Matemaattinen malli . . . 9
2.1.1 Täydellinen elektrodimalli . . . 9
2.2 Käänteisongelma . . . 15
2.2.1 Jacobin matriisin laskeminen . . . 16
3 Level Set- menetelmä 18 3.1 Level set -funktio ja Hamilton-Jacobi yhtälöt . . . 18
3.2 Merkki-etäisyysfunktio . . . 20
4 Materiaalit ja menetelmät 22 4.1 Level Set-lähestymistapa ja EIT . . . 22
4.1.1 Kaksi johtavuutta . . . 23
4.1.2 Kolme tai neljä johtavuutta . . . 25
4.1.3 Askelfunktion ja deltafunktion approksimaatiot . . . 26
4.1.4 Gauss-Newton algoritmi ja Matlab implementointi . . . 27
4.1.5 Hilat . . . 28
5 Tulokset 29 5.1 2D . . . 29
5.1.1 Yksi level set -funktio . . . 29
5.1.2 Kaksi level set -funktiota . . . 37
5.2 3D . . . 41
3
6 Yhteenveto & Pohdinta 44
Liitteet 47
Viitteet 51
Luku 1
Johdanto
Impedanssitomografia (Electrical Impedance Tomography, EIT) on kuvantamismo- daliteetti sähköä johtaville kohteille. Kuvantamismenetelmän idea on, että kohteen reunalle kiinnitetyillä elektrodeilla syötetään pienitaajuista vaihtovirtaa kuvannetta- vaan kohteeseen ja sen seurauksena syntyvät jännitteet elektrodeilla mitataan. Koh- teen sisäinen johtavuus- tai resistiivisyysjakauma voidaan rekonstruoida reunalta tehtyjen ei-invasiivisten jännitemittausten perusteella. Rekonstruoidun johtavuusja- kauman avulla saadaan merkittävää tietoa kohteen rakenteellisista ominaisuuksis- ta ja voidaan esimerkiksi tutkia kehon funktionaalisia prosesseja. Impedanssitomo- grafialla onkin lukuisia sovelluskohteita monilla tieteenaloilla kuten esimerkiksi lää- ketieteessä [17, 1, 54, 2, 33], teollisuudessa [34, 22, 36, 48, 18, 39, 40], geologiassa [27, 4, 19, 37, 12, 46] ja rakenteiden kunnon tutkimisessa [13, 5, 23].
Matemaattisessa mielessä johtavuuden (tai resistiivisyyden) rekonstruktio- ongelma on huonostiasetettu epälineaarinen käänteisongelma. Huonostiasetetulla on- gelmalla tarkoitetaan, että pienetkin virheet esimerkiksi mittausdatassa voivat ai- heuttaa suuria virheitä lopullisessa ratkaisussa/estimaatissa. Matemaattisen taustan EIT:n käänteisongelmalle kehitti Alberto Calderón [7] 1980-luvun alussa, jonka jäl- keen on kehitetty monenlaisia lähestymistapoja rekonstruktio-ongelman ratkaisemi- seksi. Yksi erittäin paljon käytetty tapa on käsitellä käänteisongelmaa epälineaari- sena LS-ongelmana lisättynä regularisointitermillä. Tässä lähestymistavassa etsitään johtavuudelle estimaattia, joka minimoi mitattujen jännitteiden ja suoran ongelman avulla laskettujen jännitteiden erotuksen normin neliön ja regularisointifunktionaalin summan.
Tässä tutkielmassa kuvannettavan kohteen sisällä olevien mahdollisten epäho-
5
mogeenisuuksien muodot, paikat ja johtavuudet estimoidaan soveltamalla tasa- arvojoukkojen eli niin sanottua level set -menetelmää. Level set -menetelmä on te- hokas ja mukautuva numeerinen menetelmä muotojen ja pintojen muutosten seu- raamiseen. Alkuperäisen level set -idean kehittivät Osher ja Sethian [32] 1980-luvun lopulla ja sittemmin siitä on tullut suosittu työkalu monilla tutkimuksen aloilla esi- merkiksi kuvankäsittelyssä, optimoinnissa ja virtauslaskennassa [14, 28]. Level set -menetelmien suuri etu numeerisessa laskennassa on, ettei kohteiden pintoja tai reu- noja tarvitse parametrisoida missään vaiheessa.
Level set -ideaa on käytetty onnistuneesti myös erilaisten käänteisongelmien rat- kaisemiseen [41, 6, 15]. EIT:n rekonstruktio-ongelma on myös ratkaistu soveltamal- la level set -menetelmää [11, 42, 43, 44, 24, 38]. Tämän tutkielman lähestymistapa EIT:n käänteisongelman ratkaisuun level set -menetelmällä on hyvin samanlainen kuin edellä mainituissa viitteissä. Lähestymistavan pääpiirre on esittää kohteen joh- tavuus ja kohteessa olevien mahdollisten eri johtavuuksien väliset rajapinnat level set -formuloinnin avulla. Johtavuus kohteessa oletetaan olevan paloittain vakio funk- tio. Tämän jälkeen kun johtavuudelle on olemassa esitys, niin impedanssitomografian käänteisongelma ratkaistaan epälineaarisena LS-ongelmana käyttäen yleistettyä Tik- honov regularisointia.
Tutkielman tavoitteet ja sisältö
Tämän tutkielman päämäärä on tutkia ja tarkastella level set -menetelmän toimivuut- ta EIT:n rekonstruktio-ongelman ratkaisussa sekä 2D- että 3D-geometriassa. Mene- telmää testataan monissa erilaisissa tilanteissa. Rekonstruktioita tehdään monille geo- metrialtaan erilaisille kohteille, jotka sisältävät kahta tai kolmea eri johtavuuden ar- voa. Menetelmää testataan tilanteissa, joissa johtavuuksien arvot tunnetaan tai niitä estimoidaan samanaikaisesti kohteen muodon estimoinnin yhteydessä. Myös menetel- män sietoa mittauskohinalle testataan. Numeerisissa testeissä tehtävät rekonstruktiot perustuvat simulointien avulla laskettuun dataan.
Tutkielma on jaettu kuuteen eri kappaleeseen. Kappaleessa 2 on esitetty EIT:n perusperiaate ja siihen liittyvä suora ongelma ja käänteisongelma. Kappaleessa 3 on käyty läpi level set -menetelmän perusperiaatteita. Kappaleessa 4 on johdettu level set -menetelmän matemaattinen formulointi EIT:n käänteisongelman ratkaisussa sekä
esitetty kaikki numeeriset testit yksityiskohtaisesti. Tulokset on esitetty kappaleessa 5 ja lopuksi viimeisessä kappaleessa on lyhyt yhteenveto tutkielmasta ja kommentit saaduista tuloksista.
Luku 2
Impedanssitomografia (EIT)
Impedanssitomografia on kuvantamismodaliteetti, jossa kuvattavan kohteen sisäinen johtavuusjakauma rekonstruoidaan kohteen reunoilta tehtyjen jännitemittausten pe- rusteella. Tyypillinen mittausasetelma on esitetty kuvassa (2.1). Joukko elektrodeja {el}on kiinnitetty kuvattavan kohteen Ω reunalle ∂Ω. Elektrodien kautta kohteeseen syötetään sähkövirtaa ja syntyvät jännitteet mitataan käyttäen samoja elektrodeja.
Yleensä syötettävät sähkövirrat ovat matalataajuista (10-100kHZ) ja pieniamplitu- dista (1-5mA) vaihtovirtaa.
Kuva 2.1: Esimerkki mittausasetelmasta.
EIT kuuluu diffuuseihin tomografiamodaliteetteihin, kun sähkövirta johdetaan 8
kohteeseen se leviää diffuusisti koko kohteen tilavuuteen. Tämän seurauksena tomo- grafiakuvan rekonstruointi on hankalampaa kuin esimerkiksi röntgentomografiassa, jossa röntgensäteet kulkevat suoraan ja yhden mittauksen informaatio tulee hyvin ra- jatulta alueelta. Sähkönjohtavuusjakauman estimointi reunalta tehtyjen mittausten perusteella on epälineaarinen, huonostiasetettu käänteisongelma. Tämä on tyypilli- nen piirre diffuuseille tomografioille. Käytännön määritelmä huonostiasetetulle kään- teisongelmalle on, että hyvinkin pienet virheet mittauksissa ja matemaattisissa mal- leissa voivat aiheuttaa suuria virheitä ratkaisussa. Toisin sanoen tämä tarkoittaa sitä, että klassiset (LS) ratkaisut ovat epästabiileja ja/tai ei-yksikäsitteisiä.
Sähkönjohtavuusjakauman rekonstruointiin on olemassa monenlaisia lähestymis- tapoja, mutta yleisin tapa on rakentaa matemaattinen malli (ns. elektrodimalli) suo- ralle ongelmalle. Elektrodimallin avulla elektrodijännitteet reunalla voidaan laskea, kun syötettävät virrat ja kohteen sähkönjohtavuus tunnetaan. Perinteinen lähestymis- tapa ratkaista EIT:n käänteisongelma elektrodimalliin perustuen on Tikhonov regula- risointi [49, 51]. Tikhonov regularisoinnissa klassiseen epälineaariseen LS-ongelmaan lisätään sakkofunktionaali, jolloin alkuperäinen huonostiasetettu ongelma korvataan hyvin käyttäytyvällä approksimaatiolla, jonka ratkaisu on (toivottavasti) lähellä al- kuperäisen ongelman ratkaisua.
2.1 Matemaattinen malli
EIT:n käänteisongelma on ratkaista kohteen sisäinen johtavuusjakauma kohteen reu- nalta tehtyjen jännitemittausten perusteella. Käänteisongelman ratkaisemiseksi tarvi- taan matemaattinen malli (t.s suora malli), jonka avulla lasketaan jännitteet reunalla, kun syötettävät sähkövirrat ja johtavuusjakauma tiedetään. Toistaiseksi tarkin malli on niin sanottu täydellinen elektrodimalli (Complete Electrode Model, CEM)[45, 10].
Seuraavaksi esitetty täydellisen elektrodimallin johtaminen ja numeerinen ratkaisu perustuu pääosin viitteisiin [50, 26, 52, 25, 45, 10, 9].
2.1.1 Täydellinen elektrodimalli
Sähkömagneettisten kenttien käyttäytyminen kohteen Ω sisällä voidaan mallintaa Maxwellin yhtälöillä
∇ ·D=ρc (2.1)
∇ ·B = 0 (2.2)
∇ ×E =−∂B
∂t (2.3)
∇ ×H =J+ ∂D
∂t , (2.4)
missä ∇× on roottorioperaattori, ∇· on divergenssioperaattori, ∂t∂ on osittaisderi- vaatta ajan suhteen, E on sähkökenttä, H on magneettikenttä, B on magneettinen induktio, D on sähkövuon tiheys (electric displacement), ρc on varaustiheys ja J on virrantiheys, J =Js+Jo= lähdevirta + ohminen virta. Jos oletetaan, että kohde Ω on lineaarinen ja isotrooppinen väliaine, voidaan kirjoittaa
D=E (2.5)
B =µH (2.6)
J =σE, (2.7)
missä on permittiivisyys,µ on permeabiliteetti jaσ =σ(x) ,x∈Ω, on johtavuus.
EIT:ssä tehdään yleensä ns. kvasi-staattinen approksimaatio, joka tarkoittaa aikariippuvuuden huomiotta jättämistä, vaikkakin käytetään vaihtovirtaa. Kvasi- staattinen approksimaatio on pätevä, jos vaihtovirran taajuus on tarpeeksi pieni. Ap- proksimaation johdosta aikaderivaatat yhtälöissä (2.3) ja (2.4) ovat nollia. Lisäksi lähdevirta Js kohteen sisällä voidaan olettaa nollaksi, koska Js = 0 EIT:ssä käytet- tävillä vaihtovirran taajuuksilla. Nyt koska ∇ ×E = 0 niin on olemassa sähköinen potentiaali u siten että
E =−∇u. (2.8)
Ottamalla divergenssi yhtälöstä (2.4) ja hyödyntämällä yhtälöitä (2.7) ja (2.8) saa- daan yhtälö, joka mallintaa sähköiset ominaisuudet kohteen sisällä
∇ ·σ∇u= 0, u∈Ω. (2.9)
Yhtälö (2.9) on malli sähköiselle potentiaalille vain kohteen sisällä. Jotta saataisiin tarkka matemaattinen malli koko ongelmalle, myös reunaehdot täytyy määrittää.
Reunaehdot
Reunaehdot sähköiselle potentiaalille täydellisessä elektrodimallissa ovat seuraavat
u+zlσ∂u
∂n =Ul, x∈el, l= 1,2, ..., L (2.10)
Z
el
σ∂u
∂ndS =Il, x∈el, l= 1,2, ..., L (2.11) σ∂u
∂n = 0, x∈∂Ω\ ∪Ll=1el, (2.12) missä Il on elektrodille el syötetty virta, L on elektrodien lukumäärä, zl on elektro- din el ja kohteen välinen kontakti-impedanssi jaUl on elektrodillael mitattu jännite.
Yhtälö (2.10) ottaa huomioon oikosulku-ilmiön elektrodeilla (ts. jännite Ul on vakio koko elektrodilla el) sekä kontakti-impedansseista zl johtuvan ylimääräisen jännite- häviön elektrodeilla. Kaksi alinta ehtoa (2.11) ja (2.12) tarkoittavat käytännössä, että virtaa syötetään vain elektrodeilta el. Yksityiskohtaisempi johtaminen reunaehdoille on esitetty esimerkiksi viitteissä [52, 50, 10].
Edellisten reunaehtojen lisäksi myös kaksi ehtoa syötettävälle virralle ja mitatta- ville jännitteille on oltava voimassa, jotta suoran ongelman ratkaisu olisi olemassa ja se olisi yksikäsitteinen [45, 35]. Nämä ehdot ovat
L
X
l=1
Il = 0 (2.13)
L
X
l=1
Ul = 0. (2.14)
Ylempi ehto tarkoittaa, että varauksen säilymislain pitää olla voimassa ja alempi yhtälö vaatii, että potentiaalimittauksille täytyy valita jokin referenssipiste. Koko- naisuudessaan täydellinen elekrodimalli koostuu yhtälöistä (2.9)-(2.14). Täydellisen elektrodimallin ratkaisun olemassaolo ja yksikäsitteisyys on todistettu viitteessä [45].
Suoran ongelman numeerinen ratkaisu
Täydellinen elektrodimalli (CEM) liittää yhteen elektrodijännitteet Ul, syötettävät virrat Il ja sähkönjohtavuusjakauman σ. EIT:n suora ongelma on ratkaista kohteen
sisäinen sähköinen potentiaali u ja elektrodipotentiaalit Ul, kun syötettävät virrat ja johtavuusjakauma tiedetään. Täydellisen elektrodimallin (suoran ongelman) ana- lyyttinen ratkaisu on olemassa vain hyvin yksinkertaisille geometrioille sen moni- mutkaisten reunaehtojen takia. Tämän takia EIT:n suoran ongelman ratkaisemiseksi käytetään numeerisia menetelmiä. Osittaisdifferentiaaliyhtälöiden ratkaisemiseksi on olemassa monia eri numeerisia menetelmiä, mutta tässä tutkielmassa EIT:n suoran ongelman ratkaisussa käytetään äärellisten elementtien menetelmää (Finite Element Method, FEM), koska se soveltuu hyvin osittaisdifferentiaaliyhtälöiden ratkaisuun monimutkaisissa geometrioissa ja hankalissa reunaehdoissa [3, 20, 29].
Ensimmäinen askel täydellisen elektrodimallin FEM-ratkaisussa on kirjoittaa rat- kaistava osittaisdifferentiaaliyhtälö ja reunaehdot niin sanotussa variationaalimuodos- sa (ts. heikossa muodossa). Tämän jälkeen variationaaliongelman ratkaisulle etsitään äärellisdimensionaalinen approksimaatio. Täydellisen elektrodimallin heikon muodon johtaminen yksityiskohtaisesti on esitetty viitteessä [45]. Täydellisen elektrodimallin variationaalimuodolla on heikko ratkaisu (u, V)
B((u, U),(v, V)) =
L
X
l=1
IlVl, (2.15)
missä v ∈H1(Ω) ja V ∈RL. Yhtälön (2.15) vasenta puolta
B((u, U),(v, V)) =
Z
Ω
σ∇u· ∇v dx+
L
X
l=1
Z
el
(u−Ul)(v−Vl)dS (2.16) kutsutaan bilineaarimuodoksiB : (H1(Ω)⊗RL)×(H1(Ω)⊗RL) →R.
Äärellistä ratkaisua variationaalimuodolle (2.15) etsittäessä kohde Ω täytyy ensin diskretisoida ja jakaa pieniin elementteihin. Tämän jälkeen tuntemattomille funktioil- le (ujaU) kirjoitetaan äärellisdimensionaaliset approksimaatiot. Potentiaalijakaumaa u kohteen sisällä approksimoidaan äärellisellä summalla
uh =
Nn
X
i=1
αiϕi (2.17)
ja samoin myös jännitteitäU elektrodeilla
Uh =
L−1
X
j=1
βjnj, (2.18)
missä funktiot ϕi muodostavat kannan äärellisulotteiselle aliavaruudelle Qh ⊂ H1(Ω) ja kantavektorit nj ∈ RL valitaan esimerkiksi n1 = (1,−1,0, ...,0)T, n2 = (1,0,−1, ...,0)T, ..., nL−1 = (1,0, ...,−1)T, jolloin ehto (2.14) täyttyy. Approksimaa- tioissa (2.17) ja (2.18)Nn on solmupisteiden määrä elementtikannassa, L on elektro- dien määrä ja kertoimetαi jaβjovat tuntemattomia, jotka täytyy määrittää. Nyt hyö- dyntämällä äärellisten elementtien teoriaa [3] ja sijoittamalla approksimaatiot (2.17) ja (2.18) variationaalimuotoon (2.15) saadaan matriisiyhtälö
Ab=f, (2.19)
missä vektori b ∈ RNn+L−1 on muotoa b=
α β
, (2.20)
missä α ∈ RNn ja β ∈ RL−1. Datavektori f on muotoa f =
0
PL
l=1Il(nj)l
=
0 Ie
, (2.21)
missä 0= (0, ...,0)T ∈ RNn jaIe= (I1−I2, I1−I3, ..., I1−IL)T ∈ RL−1. Matriisi A ∈ R(Nn+L−1)×(Nn+L−1) on harva blokkimatriisi ja se on muotoa
A=
B C
CT D
, (2.22)
missä osamatriisit ovat muotoa
B(i, j) =
Z
Ω
σ∇ϕi· ∇ϕjdx+
L
X
l=1
1 zl
Z
el
ϕiϕj dS ,
i, j = 1,2, ..., Nn (2.23)
C(i, j) = −
1 z1
Z
e1
ϕi dS− 1 zj+1
Z
ej+1
ϕidS
,
i= 1,2, ..., Nn, j = 1,2, ..., L−1 (2.24)
D(i, j) =
L
X
l=1
1 zl
Z
el
(ni)l(nj)ldS
=
|e1|
z1 , i6=j
|e1|
z1 +|ezj+1|
j+1 , i=j
, i, j = 1,2, ..., L−1. (2.25) Ratkaisemalla yhtälöstä (2.19) b =A−1f saadaan suoran ongelman approksimatiivi- nen ratkaisu. Ensimmäisten kertoimien αi, i = 1,2, ..., Nn avulla saadaan laskettua sisäinen potentiaalijakauma uh solmupisteissä ja jännitteet Uh elektrodeilla saadaan määritettyä lopuista βj, j = 1,2, ..., L−1 kertoimista. Potentiaalit elektrodeilla saa- daan laskettua kertoimien β ja yhtälön (2.18) avulla
Uh =
L−1
X
j=1
βjnj =
β1
−β1 0 0 0 ...
+
β2 0
−β2 0 0 ...
+· · ·+
βL−1 0 0 ... 0
−βL−1
=
PL−1
l=1 βl
−β1
−β2
−β3 ...
−βL−1
(2.26)
Tämä voidaan kirjoittaa lyhyesti matriisimuodossa
Uh =Cβ , (2.27)
missä β = (β1, β2, ..., βL−1)T ja C ∈ RL×(L−1) on matriisi, joka on muotoa
C=
1 1 1 · · · 1
−1 0 0 · · · 0 0 −1 0 · · · 0 0 0 −1 · · · 0 ... ... ... . .. ... 0 0 0 · · · −1
. (2.28)
Käytännössä varsinaiset EIT-mittaukset koostuvat jännitteistä, jotka on mitattu elektrodiparien välillä Ui = Ulh - Ukh, i = 1,2, ..., m, missä m mittausten lukumäärä.
Esimerkiksi jännitemittaukset voisivat muodostua potentiaalieroista U1 = U2h −U1h, U2 = U3h−U2h, jne. Suoran ongelman ratkaisuna saaduista elektrodijännitteistä Uh aitoja mittauksia vastaavat potentiaalierot voidaan määrittää seuraavasti
U =MUh =MCβ , (2.29) missä M ∈Rm×L on mittausmatriisi (measurement pattern), joka kuvaa miltä elekt- rodeilta jännitteet on mitattu.
2.2 Käänteisongelma
EIT:n käänteisongelma on määrittää johtavuusjakaumaσkohteen sisällä, kun syötet- tävät virrat ja niitä vastaavat mitatut jännitteet tiedetään. Fysikaalisen mallinnuksen (suora malli) ja FEM-diskretoinnin avulla voidaan kirjoittaa EIT-havaintomalli
V =U(σ;z, I) +e , (2.30)
missä e on havaintovirhe, V on mitatut jännitteet elektrodeilla ja U(σ;z, I) on dis- kreetti matemaattinen malli, joka yhdistää johtavuuden σ ja mittaukset V toisiinsa.
EIT-rekonstruktiomenetelmiä on olemassa monia erilaisia, mutta yksi yleinen mene- telmä on käsitellä EIT:n käänteisongelmaa epälineaarisena LS-ongelmana. Tyypilli- nen lähestymistapa on käyttää yleistettyä Tikhonov-regularisointia [49, 51], jolloin minimoitava LS-funktionaali on muotoa
F(σ) =kL1(V −U(σ))k2+αA(σ), (2.31) missä L1 on painotusmatriisi, α on regularisointiparametri jaA(σ) on regularisointi- funktionaali. Useasti regularisointifunktionaali on muotoa
A(σ) =kL2(σ−σ∗)k2, (2.32) missä L2 on regularisointimatriisi ja σ∗ on johtavuuden prioriestimaatti, joka voi ol- la esimerkiksi johtavuuden odotusarvo. Datan riippuvuus johtavuudesta σ on epä- lineaarinen, joten myös funktionaali (2.31) on epälineaarinen johtavuuden suhteen.
Funktionaalin minimoimiseksi ja ratkaisunσ löytämiseksi voidaan käyttää iteratiivi- sia menetelmiä. Tässä työssä epälineaarisen käänteisongelman ratkaisussa käytetään Gauss-Newton algoritmia, jolloin iteratiivinen algoritmi johtavuusestimaatille ˆσ saa- daan seuraavasti
ˆ
σk+1= ˆσk+hkJkTW1Jk+αW2−1JkTW1(V −U(σk))−αW2(σ−σ∗) , (2.33) missä Jk on jännitemittausten Jacobin matriisi sen hetkisen johtavuusestimaatin ˆσk suhteen, hk on askelpituus, W1 = LT1L1 ja W2 = LT2L2.
2.2.1 Jacobin matriisin laskeminen
Suoran ongelman ratkaisuna saatujen jännitemittausten U Jacobin matriisia johta- vuusjakauman suhteen tarvitaan usein iteratiivisissa menetelmissä ja se on muotoa
J = ∂U
∂σ1, . . . , ∂U
∂σNn
!
. (2.34)
Lasketut jännitteet elektrodeilla voidaan kirjoittaa U = MUh = MCβ = MCbe , missä b on sama kuin yhtälössä (2.20) ja matriisi Ce = ( 0 C ), 0 ∈ RL×Nn. Tällöin lasketut jännitteet U = MCβ = MCbe = fMb ja niiden osittaisderivaatta m:nnen johtavuusparametrin σm suhteen saadaan
∂U
∂σm = ∂fMb
∂σm =fM ∂b
∂σm . (2.35)
Termi ∂σ∂b
m voidaan kirjoittaa yhtälön (2.19) avulla muotoon
∂b
∂σm = ∂(A−1f)
∂σm . (2.36)
Yhtälön (2.36) oikea puoli voidaan muokata muotoon
∂(A−1f)
∂σm =−A−1 ∂A
∂σmA−1f =−A−1 ∂A
∂σmb . (2.37)
Nyt yhtälöiden (2.36) ja (2.37) avulla osittaisderivaatat (2.35) voidaan kirjoittaa
∂U
∂σm
=fM ∂b
∂σm
=−MAf −1 ∂A
∂σm
b
=−
(A−1)TfMT
T ∂A
∂σmb =−
A−1MfT
T ∂A
∂σmb
=−ˆbT ∂A
∂σmb , (2.38)
missä ˆb = A−1fMT. Matriisissa A ainoastaan osamatriisi B riippuu johtavuudesta σ, joten osittaisderivaatta ∂A/∂σm voidaan laskea seuraavasti
∂A
∂σm =
∂B
∂σm 0
0 0
. (2.39)
Johtavuudelle voidaan kirjoittaa äärellisulotteinen esitystapa muodossa σ(x) =
Nn
X
k=1
σkϕk(x). (2.40)
Tällöin matriisi ∂B/∂σm on muotoa
∂B(i, j)
∂σm = ∂
∂σm
Z
Ω
σ∇ϕi· ∇ϕj dx
=
Z
Ω
∂σ
∂σm∇ϕi· ∇ϕj dx
=
Z
Ω
ϕm∇ϕi· ∇ϕj dx . (2.41)
Luku 3
Level Set- menetelmä
Level set -menetelmät yleisesti ovat numeerisia menetelmiä tai tekniikoita, joiden avulla voidaan seurata muuttuvien tai etenevien rajapintojen ja niiden muotojen ke- hittymistä. Alkuperäisen level set -menetelmän tarkoitus oli analysoida ja seurata etenevien rajapintojen liikettä nopeuskentässä ja sen kehittivät Stanley Osher ja Ja- mes A. Sethian 1980-luvun lopulla. Viime aikoina level set -menetelmille on kehitetty monia eri sovelluskohteita esimerkiksi optimoinnissa ja kuvankäsittelyssä. [14, 28]
3.1 Level set -funktio ja Hamilton-Jacobi yhtälöt
Alkuperäinen idea Osherin ja Sethian [32] kehittämän level set -menetelmän taka- na on varsin yksinkertainen. Rajapinta Γ ympäröi aluetta Λ ja tavoitteena on ana- lysoida ja laskea tämän rajapinnan muutos tunnetun nopeuskentän v vaikutuksen alaisuudessa. Muutoksen aiheuttava kenttä voi riippua paikasta, ajasta, (raja)pinnan geometriasta (esim. normaalivektorin suunnasta) tai esimerkiksi ulkoisesta fysikaali- sesta virtauksesta. Pääidea level set -menetelmässä on määrittää sileä ja vähintään Lipschitz-jatkuva funktio φ(x, t), joka kuvaa rajapinnan Γ implisiittisesti joukoksi, missä funktio menee nollaan eli φ(x, t) = 0. Tämä on ydinidea kaikissa level set - menetelmissä sovelluksista riippumatta [30, 31, 32]. Esimerkki 2D-tilanteessa, missä level set -funktionφnollataso määrittää alueen Λxy-tasoon, on esitetty kuvassa (3.1).
18
Kuva 3.1: Esimerkki 2D-tapauksesta (a) Level set -funktio φ leikkaa nollatason φ(x, y) = 0 ja (b) Nollataso ja sen määrittämä rajapinta Γ eri kulmasta.
Level set -funktiolla voi olla hieman eriäviä määritelmiä riippuen käytettävästä so- velluksesta. Kuitenkin tyypillinen lähestymistapa on määritellä implisiittinen funktio, jolla on seuraavat ominaisuudet [30]
φ(x, t)<0, kunx∈Λ φ(x, t)>0, kunx /∈Λ
φ(x, t) = 0, kunx∈∂Λ = Γ(x, t)
(3.1)
Rajapinnan Γ(x, t) seuraamiseen ja kehittymiseen myöhemmillä ajanhetkillä t kun φ(x, t) muuttuu tarvitsee vain siis paikantaa kohdat missä φ(x, t) = 0. Tällä yksin- kertaisella idealla on erittäin suuri merkitys numeerisen laskennan kannalta, koska näin topologiset muutokset kuten alueiden Λ hajoaminen ja yhdistyminen ovat hy- vin määriteltyjä ja tapahtuvat itsestään menetelmän sisällä. Myöskään minkäänlais- ta rajapintojen parametrisointia ei tarvita eikä etukäteistietoa mahdollisten alueiden määrästä. [30, 31]
Alkuperäisen level set -menetelmän formuloinnissa rajapinnan Γ muutos saadaan aikaan, kun funktionφarvot muuttuvat nopeuskentänv vaikutuksesta. Tätä funktion φ muutosta voidaan kuvata matemaattisesti osittaisdifferentiaaliyhtälöllä
∂φ
∂t +v· ∇φ= 0. (3.2)
Nopeuskentän v komponentti funktion φ tasa-arvopinnan normaalin suuntaan voi- daan kirjoittaa vektorilaskennan perusteella vN = v· |∇φ|∇φ. Sijoittamalla tämä yhtä- löön (3.2) saadaan
∂φ
∂t +vN|∇φ|= 0, (3.3)
missä |·| on euklidinen normi. Yllä oleva yhtälö on välttämätön, jos halutaan seu- rata rajapinnan muutosta tunnetun nopeuskentän vaikutuksessa. Molemmat yhtälöt (3.2) ja (3.3) kuuluvat niin sanottuihin Hamilton-Jacobi -yhtälöihin, joilla on erittäin keskeinen rooli level set -menetelmissä. Hamilton-Jacobi -yhtälöt voidaan yleistää muotoon
∂φ
∂t +H(φx, φy, φz) = 0, (3.4) missä H on Hamiltonin matriisi, joka voi olla sekä paikan että ajan funktio [30, 31]. Viitteissä [30, 31] on esitelty monia erilaisia ja eri tilanteisiin sopivia numeerisia menetelmiä yhtälön (3.4) ratkaisemiseksi.
3.2 Merkki-etäisyysfunktio
Level set -funktiosta φ ei ole vielä sanottu muuta kuin, että sen tulisi olla impli- siittinen ja jossakin määrin sileä. Seuraavaksi esitellään yksi vaihtoehto level set - funktioksi, merkki-etäisyysfunktio (signed distance function). Merkki-etäisyysfunktio on implisiittinen funktio ja sen on osoitettu käyttäytyvän paremmin kuin minkään muun yhtälön (3.1) toteuttavan funktion kun lasketaan numeerisia approksimaatioita Hamilton-Jacobi -yhtälöille. [30]
Level set -funktio φ määritellään merkki-etäisyysfunktion avulla seuraavasti
φ(x) =
−d(x), kunx∈Λ d(x), kunx /∈Λ
0, kunx∈∂Λ
, (3.5)
missä d(x) tarkoittaa etäisyysfunktiota, joka on määritelty
d(x) = min(|x−xB|) ∀xB ∈∂Λ. (3.6) Merkki-etäisyysfunktion käyttöön level set -funktiona liittyy monia muita ominai- suuksia, joita muilla funktiovaihtoehdoilla ei ole [30]. Esimerkiksi
|∇φ|= 1. (3.7) Numeerista laskentaa suoritettaessa level set -funktion φ alkuarvaus ajautuu lopulta erilleen alkuperäisestä arvostaan eikä se välttämättä enää täytä merkki- etäisyysfunktion kriteerejä. Tämä voi aiheuttaa epästabilisuutta ja virhettä lopullises- sa ratkaisussa. Tämän takia on suositeltavaa uudelleen alustaa φ aika ajoin takaisin muistuttamaan merkki-etäisyysfunktiota. Näin ollen varmistutaan level set -funktion riittävästä sileydestä ja ratkaisun stabiilisuudesta. [30, 47, 8]
Artikkelissa [47] on esitetty yksi tapa, jolla φ voidaan uudelleen alustaa muistut- tamaan merkki-etäisyysfunktiota. Uudelleen alustus perustuu yhtälön
∂φ
∂t +S(φ0)(|∇φ| −1) = 0 (3.8)
ratkaisuun numeerisesti. Funktio S(φ0) on merkkifunktio, joka saa arvon 1 alueen Λ ulkopuolella, -1 sisäpuolella ja 0 rajapinnalla Γ.
Luku 4
Materiaalit ja menetelmät
Tämän tutkielman tavoitteena on tutkia level set -menetelmän toimivuutta impe- danssitomografian rekonstruktio-ongelman ratkaisussa. Kaikki simuloinnit ja lasken- ta suoritettiin käyttäen Matlab-ohjelmaa. Tässä luvussa käydään läpi level set - lähestymistapa impedanssitomografiaan sekä simulointien että numeeristen testien yksityiskohdat.
4.1 Level Set-lähestymistapa ja EIT
Tässä tutkielmassa impedanssitomografian käänteisongelma ratkaistaan soveltamal- la level set -menetelmää. Käytettävä lähestymistapa käänteisongelman ratkaisuun on hyvin samankaltainen kuin viitteissä [8, 11, 43, 44, 16] pienillä eroavaisuuksilla.
Pääidea menetelmässä on, että kohteen johtavuus oletetaan paloittain vakioksi ja eri johtavuuksien alueet ja niiden väliset rajapinnat esitetään level set -formuloinnin t.s level set funktioidenφavulla. Erona tässä tutkielmassa on level set -funktionφ valin- ta. Level set funktioksi ei valita varsinaisesti merkki-etäisyysfunktiota vaan level set -funktio määritellään seuraavasti
φ(x) =
−c, kunx∈Λ c, kunx /∈Λ
, (4.1)
missä c ∈ R+ ja se voidaan valita esimerkiksi luvuksi 1. Funktion φ arvojen kehit- tyessä rajapinta Γ alueiden Λ välillä t.s nollataso φ(x) = 0 määritetään kohtaan, jossa funktion positiiviset arvot muuttuvat negatiivisiksi. Level set -funktion φ reini- talisointi tehdään tarvittaessa ja se toteutetaan asettamalla sen hetkiset negatiiviset arvot arvoon −cja positiiviset arvot arvoon c.
22
Seuraavaksi johdetaan level set -lähestymistapa johtavuuden rekonstruktio- ongelman ratkaisemiseksi, kun kohde sisältää kahta eri johtavuuden omaava mate- riaalia sekä tilanne, missä eri johtavuuksia voi olla kolme tai neljä. Ongelman rat- kaisuun kun eri materiaaleja on kaksi riittää yksi level set -funktio, mutta 3-4 eri materiaalin ratkaisuun tarvitaan kaksi level set -funktiota, jotta kaikki mahdolliset eri alueet saataisiin esitettyä. Yleinen sääntö on, ettäK määrällä level set -funktioita voidaan esittää 2K määrä eri alueita, joilla on eri johtavuus tai jokin muu parametrin arvo [8, 16, 53]. Tässä tutkielmassa rajoitutaan tapauksiin, jossa level set -funktioita on yksi tai kaksi.
4.1.1 Kaksi johtavuutta
Oletetaan, että kohde Ω sisältää kahta eri materiaalia Ω1 ja Ω2, joiden johtavuudet ovatσ1 jaσ2, vastaavassa järjestyksessä. Kun level set -funktio on määritelty, voidaan alueet Ω1 ja Ω2 ja niiden välinen rajapinta Γ määrittää sen avulla
Ω1 ={x∈Ω|φ(x)>0} ,
Ω2 ={x∈Ω|φ(x)<0} , (4.2)
Γ ={x∈Ω|φ(x) = 0} .
Kuva 4.1: Kohde Ω, joka sisältää kahta eri materiaalia johtavuuksillaσ1 ja σ2.
Tämä on myös esitetty kuvassa (4.1). Nyt kohteen Ω johtavuusσ(x) voidaan kirjoittaa paloittain vakiona funktiona level set -funktion avulla seuraavasti
σ(x) =
σ1, kunφ≥0 σ2, kunφ <0
=σ1H(φ(x)) +σ2(1−H(φ(x))), (4.3) missä H(x) on askelfunktio ja se on määritelty
H(x) =
1, kunx≥0, 0, kunx <0
. (4.4)
Kun johtavuudelle σ on olemassa esitys niin muokataan luvussa 2 esitetty yleis- tetty Tikhonov regularisointifunktionaalin (2.31) minimointi sisältämään level set - lähestymistavan. Minimoitava funktionaali on muotoa
F(σ) = F(σ(φ, σ1, σ2)) =kL1(V −U(σ))k2+αkL2(σ−σ∗)k2. (4.5) Funktionaali (4.5) on johtavuuden σ funktio, mutta σ riippuu edelleen level set - funktiosta φ ja johtavuuksista σ1 ja σ2 eli σ(φ, σ1, σ2). Jotta iteratiivista algoritmia (2.33) voitaisiin käyttää estimoitavien parametrienφ,σ1 jaσ2päivitykseen täytyy en- sin laskea kuvauksen U(σ) Jacobin matriisi estimoitavien parametrien suhteen. Ko- konaisuudessa Jacobin matriisi on muotoa
Jtot(σ) = (Jφ, Jσ1, Jσ2), (4.6) missä matriisit Jφ ∈ Rm×Nn, Jσ1 ∈ Rm×1 ja Jσ2 ∈ Rm×1 lasketaan seuraavasti deri- voinnin ketjusääntöä hyödyntäen
Jφ= ∂U(σ(φ, σ1, σ2))
∂φ = ∂U
∂σ
∂σ
∂φ = ∂U
∂σ(σ1−σ2)δ(φ), (4.7) Jσ1 = ∂U(σ(φ, σ1, σ2))
∂σ1 =
Z
Ω
∂U
∂σ
∂σ
∂σ1 dx =
Z
Ω
∂U
∂σH(φ)dx , (4.8)
Jσ2 = ∂U(σ(φ, σ1, σ2))
∂σ2 =
Z
Ω
∂U
∂σ
∂σ
∂σ2 dx =
Z
Ω
∂U
∂σ(1−H(φ))dx . (4.9)
Yhtälöiden (4.7)-(4.9) yksityiskohtaisemmat derivoinnit on esitetty liitteessä. Yhtä- lössä (4.7) δ(x) on Diracin deltafunktio, joka on määritelty
δ(x) =
1, kunx= 0, 0, muulloin
. (4.10)
Yhtälöiden (4.7)-(4.9) numeerinen toteutus on seuraavanlainen: Termi ∂U∂σ yhtälöissä (4.7)-(4.9) on perinteisen EIT-ongelman Jacobin matriisi (2.34). Yhtälössä (4.7) del- tafunktio tarkoittaa, että matriisista ∂U∂σ valitaan ainoastaan ne sarakkeet, jotka vas- taavat määritettyjä nollakäyränφ(x) = 0 solmupisteitä laskentahilassa. Nollakäyrä ei välttämättä kulje suoraan solmupisteiden kautta, jolloin nollasolmupisteiksi valitaan nollakäyrän lähellä olevat solmupisteet. Yhtälöissä (4.8)-(4.9) integraalilla ja askel- funktiolla tarkoitetaan matriisin ∂U∂σ niiden sarakkeiden summausta, jotka vastaavat johtavuusalueiden Ω1 ja Ω2 solmupisteitä. Numeerisessa laskennassa deltafunktio ja askelfunktio on järkevää korvata sileämmillä approksimaatioilla, mutta tästä aiheesta lisää myöhemmin tässä luvussa.
4.1.2 Kolme tai neljä johtavuutta
Oletetaan, että kohde Ω voi sisältää enimmillään neljää eri materiaalia, joilla on eri johtavuus σ1, σ2, σ3 ja σ4. Nyt alueiden Ω1,Ω2,Ω3 ja Ω4 esittämiseen tarvitaan kaksi level set -funktiotaφ1 ja φ2. Kirjoitetaan kohteen johtavuus σ(x) neljän johtavuuden tilanteessa nyt paloittain vakiona funktiona level set -formuloinnin avulla
σ(x) =
σ1, kunφ1 ≥0, φ2 ≥0 σ2, kunφ1 <0, φ2 ≥0 σ3 kunφ1 ≥0, φ2 <0 σ4 kunφ1 <0, φ2 <0
=σ1H(φ1)H(φ2) +σ2(1−H(φ1))H(φ2) +...
σ3H(φ1)(1−H(φ2)) +σ4(1−H(φ1))(1−H(φ2)). (4.11) Nyt kohteen Ω johtavuus on level set -funktioiden φ1 ja φ2 sekä eri alueiden johta- vuuksien σi, i = 1,2,3,4, funktio eli σ = σ(φ1, φ2, σ1, σ2, σ3, σ4). Tällöin kuvauksen U(σ) Jacobin matriisi estimoitavien parametrien suhteen edellisen kappaleen esityk- sen nojalla on nyt muotoa
Jtot(σ) = (Jφ1, Jφ2, Jσ1, Jσ2, Jσ3, Jσ4), (4.12) missä matriisit Jφ1, Jφ2, Jσ1, Jσ2, Jσ3 ja Jσ4 voidaan laskea seuraavasti
Jφ1 = ∂U(σ)
∂φ1 = ∂U
∂σ
(σ1−σ2)H(φ2) + (σ3−σ4)(1−H(φ2))
δ(φ1), (4.13) Jφ2 = ∂U(σ)
∂φ2 = ∂U
∂σ
(σ1−σ3)H(φ1) + (σ2−σ4)(1−H(φ1))
δ(φ2), (4.14) Jσ1 = ∂U(σ)
∂σ1 =
Z
Ω
∂U
∂σH(φ1)H(φ2)dx , (4.15)
Jσ2 = ∂U(σ)
∂σ2 =
Z
Ω
∂U
∂σ(1−H(φ1))H(φ2)dx , (4.16)
Jσ3 = ∂U(σ)
∂σ3 =
Z
Ω
∂U
∂σH(φ1)(1−H(φ2))dx , (4.17)
Jσ4 = ∂U(σ)
∂σ4 =
Z
Ω
∂U
∂σ(1−H(φ1))(1−H(φ2))dx . (4.18) Yksityiskohtaisemmat derivoinnit on esitetty liitteessä.
4.1.3 Askelfunktion ja deltafunktion approksimaatiot
Numeerisessa laskennassa johtavuuden σ(x) esityksessä ja Jacobin matriisien yhtä- löissä esiintyvät askelfunktio H ja deltafunktio δ on suotavaa korvata sileämmillä funktioilla. Tässä tutkielmassa käytetään seuraavia funktioita approksimoimaan as- kel -ja deltafunktiota [8, 16, 53]
H1(φ) = 1 π tan−1
φ 1
+1
2 , (4.19)
δ2(φ) = 2
π(φ2+22) , (4.20)
missä parametrit1 ja2 määrittävät approksimaation sileyden. Parametrieni valin- ta vaikuttaa approksimaatioiden tarkkuuteen, mitä pienempii niin sitä terävämpi ja tarkempi approksimaatio on. Kuitenkin numeerisen laskennan kannalta liian pienten parametrien i käyttö ei ole kannattavaa. Etenkin deltafunktion tapauksessa hyvin pienillä parametrin 2 arvoilla approksimaatiosta tulee erittäin terävä ja se voi ai- heuttaa epästabiilisuutta numeerisessa laskennassa. Askelfunktion approksimaatiolle
voidaan ja on suositeltavaakin käyttää kohtuullisen pientä parametrin2 arvoa. Para- metrien i valintaan vaikuttaa myös numeerisessa laskennassa käytettävän hilan ele- menttien keskimääräinen koko. Tämän tutkielman rekonstruktioissa sileysparametrit valittiin pääsääntöisesti seuraavasti hieman riippuen käytetystä hilasta: 1 ∈ [0.01 0.15] ja 2 ∈ [1 5] .
4.1.4 Gauss-Newton algoritmi ja Matlab implementointi
Seuraavaksi esitetään Gauss-Newton algoritmi EIT:n käänteisongelman ratkaisemi- seksi level set -menetelmällä sekä kuinka tutkielman numeeriset testit on toteutettu Matlabilla. Tämä on näytetty vain yhden level set -funktion ja kahden johtavuuden tapaukselle, mutta siirtyminen kahteen level set -funktioon ja useampaan johtavuu- teen on hyvin suoraviivaista.
Minimoitava funktionaali (4.5) on johtavuuden σ funktio, mutta level set - lähestymistavassa johtavuus riippuu vielä level set -funktiosta φ sekä yksittäisistä johtavuuden arvoista σ1 ja σ2. φ, σ1 ja σ2 ovat tuntemattomat parametrit ja ne voi- daan estimoida käyttäen iteratiivistä Gauss-Newton algoritmia. Nyt luvussa 2 esitetty Gauss-Newton algoritmi (2.33) saa muodon
ˆ
xk+1 = ˆxk+hkJkTW1Jk+αW2−1JkTW1(V −U(ˆxk))−αW2(ˆxk−x∗) , (4.21) missä ˆxk = (φk, σ1(k), σ2(k)) ja Jacobin matriisi J on yhtälön (4.6) muotoa. Tässä tut- kielmassa painotusmatriisi L1 ja regularisointimatriisi L2 valitaan yksikkömatriiseik- si, jolloin myös matriisit W1 ja W2 ovat yksikkömatriiseja. Prioriestimaattina käy- tetään nollavektoria x∗ = 0 sekä askelpituus hk pidetään vakiona koko iteraation ajan ja se valitaan väliltä [0.10 0.40] tilanteesta riippuen. Regularisointiparametrin α valinta riippuu useasti eri tekijästä, esimerkiksi mittauksissa olevasta häiriön suuruu- desta, oikean johtavuusjakauman mahdollisesta monimutkaisesta muodosta ja joh- tavuusinkluusioiden sijainnista toisiinsa nähden. Regularisointiparametri valittiin ta- pauskohtaisesti ja parhaan regularisointiparametrin arvon valinta tehtiin tulosten vi- suaalisen tarkastelun perusteella.
Mittausten simulointi ja johtavuuden rekonstruointi Matlabin avulla on esitetty liitteessä pseudokoodi-muodossa kommentein varustettuna.
4.1.5 Hilat
Tässä tutkielmassa level set -menetelmää testataan EIT:n käänteisongelman ratkai- suun niin 2D- kuin 3D-geometriassa. 2D-tilanteessa laskenta-alueena käytetään ym- pyränmuotoista kohdetta, joka sisältää 16 elektrodia reunalla ja 3D-tapauksessa koh- teena on sylinteri, jonka reunalla on yhteensä 64 elektrodia neljässä kerroksessa. Mo- lemmissa tapauksissa numeerinen laskenta toteutetaan käyttäen kolmea eri hilaa. Si- muloidut oikeita mittauksia vastaavat jännitteet lasketaan yhdellä hilalla, suora on- gelma ja Jacobin matriisi lasketaan toisella hilalla ja kolmas hila on käänteisongelman ratkaisua varten. Kuvissa (4.2) ja (4.3) on esitetty 2D- ja 3D-tilanteissa käytetyt hilat sekä niitä vastaavat elementtien ja solmupisteiden määrät.
(a) (b) (c)
Kuva 4.2: 2D hilat: (a) mittausten simulointihila, jossa on 20642 elementtiä ja 10770 solmupistettä, (b) suoran ongelman hila, jossa on 15058 elementtiä ja 7986 solmupis- tettä ja (c) käänteisongelman hila, jossa on 918 elementtiä ja 492 solmupistettä
(a) (b) (c)
Kuva 4.3: 3D hilat: (a) mittausten simulointihila, jossa on 81597 elementtiä ja 17132 solmupistettä, (b) suoran ongelman hila, jossa on 81269 elementtiä ja 17027 solmu- pistettä ja (c) käänteisongelman hila, jossa on 30144 elementtiä ja 5882 solmupistettä
Luku 5
Tulokset
Tässä luvussa esitetään tuloksia level set -menetelmän toimivuudesta sähkönjohta- vuuden rekonstruoimiseksi erilaisissa tilanteissa sekä 2D- että 3D-geometriassa. Tut- kittavia tapauksia ovat esimerkiksi erilaisten johtavuusjakaumien (yksinkertaiset vs.
monimutkaiset) rekonstruoiminen, rekonstruktion lopputuloksen riippuvuus level set -funktioiden alkuarvauksesta ja jännitemittausten häiriötason vaikutus rekonstruk- tioon. Luku on jaettu kahteen osaan 2D- ja 3D-tuloksiin ja tarkemmat yksityiskohdat testitapauksista on kerrottu niitä käsittelevissä kappaleissa.
Simuloituihin jännitemittauksiin ei ole lisätty satunnaiskohinaa kuin ainoastaan niissä tapauksissa missä kohinatason vaikutusta rekonstruktioon on tarkasteltu. Niin sanotun "inverse crimen" välttämiseksi simuloidut mittausjännitteet laskettiin käyt- täen eri hilaa kuin rekonstruktiossa suoran ongelman ratkaisussa.
5.1 2D
2D-geometriassa tarkastellaan tilanteita, joissa kohteessa on kahta tai kolmea eri joh- tavuutta. Kahdelle eri johtavuudelle rekonstruktio voidaan suorittaa yhdellä level set -funktiolla, mutta kolmen johtavuuden tapauksessa tarvitaan kaksi level set - funktiota.
5.1.1 Yksi level set -funktio
Seuraavissa numeerisissa esimerkeissä kuvannettavan kohteen johtavuus σ koostuu kahdesta eri johtavuudesta σ1 ja σ2. Taustan johtavuus σ1 on 0.1 Sm−1 ja epäho- mogeenisuuskohteen tai -kohteiden johtavuuden arvo σ2 on 0.5 Sm−1. Taustan johta- vuus oletetaan melkein kaikissa rekonstruktioissa tunnetuksi, mutta toista johtavuu- den arvoa ei välttämättä tunneta. Yhdessä tapauksessa molemmat johtavuuden arvot
29
oletaan tuntemattomiksi. Rekonstruktioissa tarkastellaan kolmenlaisia tapauksia: 1.
molemmat johtavuuden arvot tiedetään ja estimoidaan pelkkää johtavuusjakauman muotoa 2. taustan johtavuus tiedetään ja estimoidaan samanaikaisesti epähomogee- nisuuden johtavuuden arvoa ja muotoa 3. molempia johtavuuksia estimoidaan level set -funktion rinnalla. Alkuarvauksena mahdollisille tuntemattomille johtavuuden σ1 ja σ2 arvoille käytettiin 0.2 Sm−1 ja 0.25 Sm−1.
Ensimmäisessä esimerkissä todellinen johtavuusjakauma σ sisälsi kaksi epähomo- geenisuuskohdetta. Rekonstruktiossa epähomogeenisuuksien paikat ja muodot saatiin määritettyä täsmällisesti sekä myös johtavuuden σ2 arvo pystyttiin estimoimaan tar- kasti. Tulokset pelkän level set -funktion päivitykselle on nähtävissä kuvassa (5.1) ja kuvassa (5.2) on esitetty rekonstruktio kun on estimoitu myös johtavuutta σ2 level set -funktion lisäksi.
(a) Todellinen johtavuus (b) Alkuarvaus
(c) 10. iteraatio (d) 30. iteraatio
Kuva 5.1: (a) Todellinen johtavuus. (b) Alkuarvaus. (c)-(d) Rekonstruktiot iteraa- tion eri vaiheissa kun on estimoitu pelkkää level set -funktiota.
Seuraavassa numeerisessa esimerkissä on tutkittu tilannetta, jossa kaksi ympy- rän muotoista kohdetta ovat hyvin lähellä toisiaan. Tälläisten kohteiden erottaminen toisistaan perinteisillä pikseli-rekonstruktiomenetelmillä on erittäin hankalaa, mutta
