Approksimaatiovirheiden mallintaminen raudoitteita sisältävän kohteen resistanssitomografiassa

Approksimaatiovirheiden mallintaminen raudoitteita sisältävän kohteen resistanssitomograassa

Mikko Räsänen

Pro gradu -tutkielma

Sovelletun fysiikan koulutusohjelma

Itä-Suomen yliopisto, Sovelletun fysiikan laitos

16. helmikuuta 2019

ITÄ-SUOMEN YLIOPISTO, Luonnontieteiden ja metsätieteiden tiedekunta Sovelletun fysiikan koulutusohjelma, laskennallinen fysiikka

Mikko Räsänen: Approksimaatiovirheiden mallintaminen raudoitteita sisältävän kohteen resis- tanssitomograassa

Pro Gradu-tutkielma, 64 sivua

Tutkielman ohjaajat: FT Aku Seppänen (pääohjaaja) ja FT Tuomo Savolainen Helmikuu 2019

Avainsanat: resistanssitomograa, käänteisongelma, approksimaatiovirhemenetelmä, raudoite

Tiivistelmä

Resistanssitomograa (ERT) on sähköinen kuvantamismenetelmä, jossa kuvannettavan kohteen ulkopintaan kiinnitettyjen elektrodien välillä syötetään vaihtovirtaa, ja elekt- rodien välisten jännitteiden amplitudit mitataan. Mittausten perusteella lasketaan esti- maatti kohteen sisäiselle sähkönjohtavuusjakaumalle. Sähkönjohtavuusjakauman estimoin- ti on epälineaarinen ja huonosti asetettu käänteisongelma, jonka ratkaisu on hyvin herkkä mittaus- ja mallinnusvirheille. Käänteisongelman ratkaisemiseen voidaan käyttää tilastol- lista eli Bayesilaista lähestymistapaa, jossa mittausvirheen lisäksi voidaan huomioida mal- linnusvirheitä niin sanotulla approksimaatiovirhemenetelmällä.

Resistanssitomograa soveltuu esimerkiksi betonirakenteiden kuvantamiseen, ja sillä voidaan havaita halkeamia ja kosteuden etenemistä rakenteessa. Betonirakenteet sisältä- vät usein raudoitteita, joiden sijainnit ja muodot täytyy tuntea resistanssitomograan mittausten tarkkaa mallintamista varten. Jos raudoitteiden sijainteja ja muotoja ei tunne- ta tarkasti, mallissa esiintyvä virhe aiheuttaa virheitä kohteen sähkönjohtavuusjakauman estimaattiin.

Tässä Pro Gradu - tutkielmassa käytettiin approksimaatiovirhemenetelmää resistans- sitomograan käänteisongelman ratkaisemiseen tilanteessa, jossa mittauksia kuvaava malli on virheellinen johtuen betonirakenteen sisäisen raudoitteen huomiotta jättämisestä. Ap- proksimaatiovirhemenetelmää testattiin simulaatioiden avulla sylinterigeometriassa, jossa elektrodit sijaitsivat sylinterin muotoisen kohteen ulkoreunalla, ja laattageometriassa, jossa elektrodit sijaitsivat laatan muotoisen kohteen yläpinnalla. Sylinterigeometriassa raudoite oli yksinkertainen suora tanko, ja laattageometriassa useammasta tangosta muodostuva verkko.

Tulosten perusteella kohteen sähkönjohtavuusjakaumalle saadaan tarkempia estimaat- teja mallintamalla raudoitteen huomiotta jättämisestä aiheutuvaa approksimaatiovirhet- tä. Etenkin laattageometriassa sähkönjohtavuusjakauman estimaatit ovat virheellisiä, jos approksimaatiovirhettä ei mallinneta. Huomioimalla approksimaatiovirhe saadaan hyvin kohteen sähkönjohtavuusjakaumaa vastaava estimaatti. Työn tulokset osoittavat, että ap- proksimaatiovirhemenetelmä soveltuu kohteen sisäisen raudoitteen epätarkasta mallinta- misesta aiheutuvan mallinnusvirheen korjaamiseen.

Lyhenteet

EIT Impedanssitomograa (Electrical Impedance Tomography) ERT Resistanssitomograa (Electrical Resistance Tomography) FEM Äärellisten elementtien menetelmä (Finite Element Method) MAP Maximum a Posteriori-estimaatti

CEM Tavallinen virhemalli (Conventional Error Model) EEM Paranneltu virhemalli (Enhanced Error Model)

AEM Approksimaatiovirhemalli (Approximation Error Model)

Merkinnät

Ω Kuvannettava kohde

∂Ω Kohteen Ωreuna

⃗r pisteen paikkavektori σ(⃗r) sähkönjohtavuusjakauma

e_l elektrodi, jonka indeksi on l∈N u(⃗r) sähköinen potentiaali

zl elektrodin el ja kohteen välinen kontakti-impedanssi Ul elektrodin el sähköinen potentiaali

V mitattu jännitevektori

ν pinnan ulospäin suunnattu yksikkönormaalivektori I_l sähkövirta elektrodin e_l läpi

j(⃗r) sähkövirran tiheys B(·,·) bilineaarimuoto H¹(Ω) Sobolev-avaruus

π_x(·) satunnaismuuttujan xtiheysfunktio x∗ satunnaismuuttujan xodotusarvo

Γx Satunnaismuuttujan x kovarianssimatriisi π_x,y(·) satunnaismuuttujien x ja y yhteistiheysfunktio π_x|y(·) ehdollisen satunnaismuuttujan x|y tiheysfunktio x∗|y ehdollisen satunnaismuuttujan x|y odotusarvo

⟨a, b⟩ Vektorien a ja b pistetulo

||a|| Vektorin a Euklidinen normi trace(A) matriisin A jälki

Sisältö

1 Johdanto 5

2 Resistanssitomograa 7

2.1 Resistanssitomograan suora malli . . . 8

2.2 Suoran mallin ratkaisun FEM-approksimaatio . . . 9

2.3 Käänteisongelma . . . 12

3 Approksimaatiovirhemenetelmä 18 3.1 Approksimaatiovirheen mallintaminen . . . 19

3.2 Approksimaatiovirheen aiheuttavien parametrien estimointi . . . 23

3.3 Raudoitteisiin liittyvä approksimaatiovirhe . . . 27

4 Simulaatiot 29 4.1 Simulaatiot sylinterigeometriassa . . . 29

4.1.1 Priorit ja approksimaatiovirhemallien muodostaminen . . . 30

4.1.2 Simuloidut kohteet ja estimaattien laskenta . . . 35

4.1.3 Tulokset ja pohdinta . . . 36

4.2 Simulaatiot laattageometriassa . . . 43

4.2.1 Simuloidut kohteet ja estimaattien laskenta . . . 46

4.2.2 Tulokset ja pohdinta . . . 48

5 Johtopäätökset 58

1 Johdanto

Resistanssitomograa on ainetta rikkomaton sähköinen kuvantamismenetelmä, jossa estimoi- daan kuvannettavan kohteen sähkönjohtavuusjakaumaa. Menetelmässä kohteen pintaan kiinni- tetään elektrodeja, joiden välillä syötetään vaihtovirtaa. Virransyöttöjen aikana mitataan elekt- rodien välisten jännitteiden amplitudit. Kohteen sähkönjohtavuusjakauma estimoidaan syötet- tyjen sähkövirtojen ja vastaavien mitattujen jännitteiden perusteella. Estimoidusta sähkönjoh- tavuusjakaumasta saadaan tietoa kohteen sisäisistä rakenteista, joiden sähkönjohtavuudet eroa- vat toisistaan. Resistanssitomograalla on sovelluksia lääketieteellisessä kuvantamisessa, mis- sä sitä voidaan hyödyntää hengityksen seurantaan, pehmytkudosten kuvantamiseen ja aivojen kuvantamiseen [13]. Teollisuudessa resistanssitomograaa voidaan käyttää sekoitustankkien kuvantamiseen ja putkien sisäisten virtausten kuvantamiseen [47].

Resistanssitmograalla voidaan kuvantaa sähköä johtavia kiinteitä materiaaleja [814], jo- ten menetelmä soveltuu esimerkiksi betonirakenteiden kuvantamiseen ja rakenteiden kunnon arviointiin. Resistanssitomograalla voidaan havaita esimerkiksi betonin halkeamia [8] ja veden etenemistä betonissa [13]. Tutkimuksissa menetelmää on testattu simulaatioiden avulla sekä pienikokoisilla muurauslaastista tai betonista valmistetuilla kohteilla. Resistanssitomograa on mahdollisesti käyttökelpoinen menetelmä myös oikeiden, suurikokoisten betonirakenteiden ku- vantamiseen. Rakenteiden kuvantamista varten tutkitaan myös sensing skin - tekniikkaa, jossa ehjän kohteen pintaan muodostetaan sähköä johtava kerros, jota kuvannetaan resistanssitomo- graalla [15]. Kohteen pintaan asti ulottuvat halkeamat halkaisevat tällöin myös pinnoitteen, ja halkeamat näkyvät pinnoitteen estimoidussa sähkönjohtavuusjakaumassa.

Betonirakenteissa käytetään usein raudoitteita, joilla rakenne saadaan kestämään venyttäviä voimia. Raudoite on resistanssitomograan mittauksissa tasapotentiaalissa, joten se vaikuttaa sähköiseen potentiaaliin kohteen sisällä ja edelleen elektrodien väleillä mitattaviin jännitteisiin.

Raudoitteen vaikutus mittausdataan mahdollistaa esimerkiksi raudoitteen paikantamisen resis- tanssitomograan avulla [8]. Toisaalta betonin kuvantamisessa kiinnostuksen kohde on yleensä raudoitteiden sijaan jokin sähkönjohtavuusjakaumaan vaikuttava ominaisuus, kuten betonin kosteusjakauma. Tällöin raudoitteet muodostavat mittauksia häiritsevän tekijän, joka voidaan ottaa huomioon mallintamalla raudoitteiden vaikutusta mittausdataan. Raudoitteiden mallin- tamista varten niiden sijainnit ja muodot täytyy joko tuntea ennakkoon tai estimoida saman- aikaisesti sähkönjohtavuusjakauman kanssa. Raudoitteiden sijaintien ja muotojen estimointi on laskennallisesti raskasta, ja se tuottaa ylimääräistä laskentatyötä tilanteessa, jossa tarvitaan tietoa vain kohteen sähkönjohtavuusjakaumasta. Estimoinnin sijaan raudoitteita voidaan mal- lintaa ennakkotiedon perusteella, mutta epätarkka ennakkotieto johtaa virheelliseen malliin.

Virheellisen mallin käyttö aiheuttaa edelleen virheitä sähkönjohtavuusjakauman estimaattiin.

Sähkönjohtavuusjakauman estimointi on huonosti asetettu käänteisongelma, jonka ratkaise- miseen voidaan soveltaa Bayesilaista eli tilastollista lähestymistapaa. Bayesilainen lähestymis-

tapa sallii mittausvirheen lisäksi myös mallinnusvirheiden huomioimisen estimoinnissa, jolloin saadaan tarkempia estimaatteja kuin huomioimalla vain mittausvirhe [16]. Mallinnusvirheitä kutsutaan usein myös approksimaatiovirheiksi. Resistanssitomograassa ja muissa käänteison- gelmissa on tutkittu usean eri tyyppisen approksimaatiovirheen huomiointia.

Tässä työssä testataan approksimaatiovirheen huomiointia tilanteessa, jossa resistanssito- mograan mittauksia kuvaavan mallin approksimaatiovirhe johtuu kohteen sisäisen raudoitteen huomiotta jättämisestä. Työn tavoitteena on tutkia, saadaanko kohteen sähkönjohtavuusjakau- man estimaatteja tarkennettua ottamalla approksimaatiovirhe huomioon estimoinnissa.

Kappaleessa 2 esitellään resistanssitomograan suora malli, suoran mallin ratkaisun ap- proksimaatio äärellisten elementtien menetelmällä ja työssä käytetty käänteisongelman ratkai- sumenetelmä. Kappaleessa 3 esitellään approksimaatiovirheiden mallintamisen periaate ja so- velletaan sitä raudoitteista johtuviin approksimaatiovirheisiin. Kappaleessa 4 esitellään työssä tehtyjen simulaatioiden tulokset ja pohdintaa saaduista tuloksista. Lopuksi kappaleessa 5 on lyhyt yhteenveto työstä.

2 Resistanssitomograa

Resistanssitomograa (Electrical Resistance Tomography, ERT) on sähköinen kuvantamisme- netelmä, jossa kuva muodostuu kohteen estimoidusta sähkönjohtavuusjakaumasta. Sähkönjoh- tavuusjakauma estimoidaan kohteen pinnalta tehtyjen sähköisten mittausten perusteella. Mit- taukset suoritetaan kiinnittämällä kohteen pintaan elektrodeja, joiden kautta kohteeseen syö- tetään vaihtovirtaa. Virtaa syötetään tavallisesti kahden elektrodin välillä kerrallaan, ja kun- kin virransyötön aikana mitataan kaikkien elektrodien välisten jännitteiden amplitudit. Myös useampaa kuin kahta elektrodia voidaan käyttää samanaikaisesti virransyötöissä [17, 18]. Ku- vassa 1 on esitelty ERT:n mittausasetelma.

Jos käytetyllä mittauslaitteistolla on mahdollista mitata myös jännitteiden vaihe-eroja, voi- daan mitatun amplitudi- ja vaihedatan avulla estimoida kohteen kompleksista admittiivisuus- jakaumaa, jonka reaaliosa sähkönjohtavuus on. Tällöin menetelmää kutsutaan impedanssito- mograaksi (Electrical Impedance Tomography, EIT). Lyhenteellä EIT viitataan usein myös resistanssitomograaan. Resistanssitomograalla on käytännön sovelluksia lääketieteessä [13], prosessikuvantamisessa [47] ja geofysiikassa [19]. Lisäksi se on käyttökelpoinen menetelmä ai- netta rikkomattomassa testauksessa [814].

Ω σ(~r)

−I +I

V₁ V₂

V₃

e₁ e₂ e₃

e₄

Kuva 1: Resistanssitomograan mittausasetelma. Ω on kuvannettava kohde, σ(⃗r) on kohteen sähkönjohtavuusjakauma, I on syötetyn sähkövirran amplitudi ja V_i on elektrodien e_i ja e_i+1 välillä mitatun jännitteen amplitudi.

2.1 Resistanssitomograan suora malli

Resistanssitomograan suora malli yhdistää kohteen sähkönjohtavuusjakaumanσ(⃗r)ja määrät- tyä virransyöttöä vastaavat elektrodien väliset jännitteetV. Malli koostuu sähkönjohtavuudelle σ(⃗r) ja sähköiselle potentiaalille u(⃗r) muodostetusta osittaisdierentiaaliyhtälöstä ja reunaeh- doista. ERT:n suora malli on [20]

∇ ·(σ(⃗r)∇u(⃗r)) = 0, ⃗r∈Ω (2.1)

u(⃗r) +z_lσ(⃗r)∂u

∂ν(⃗r) =U_l, ⃗r ∈e_l, l= 1, ..., L (2.2)

∫︂

el

σ(⃗r)∂u

∂ν(⃗r)dS =I_l, l = 1, ..., L (2.3) σ(⃗r)∂u

∂ν(⃗r) = 0, ⃗r∈∂Ω\

L

⋃︂

l=1

e_l, (2.4)

missä σ(⃗r) on kohteen sähkönjohtavuusjakauma, u(⃗r) on kohteen sisäinen sähköinen potenti- aali,⃗r on pisteen paikkavektori,z_l on elektrodin e_l ja kohteen välinen kontakti-impedanssi, U_l on elektrodin e_l potentiaali ja I_l on elektrodin e_l kautta syötetty sähkövirta. Lisäksi ν on ulos- päin suunnattu yksikkönormaali ja L on elektrodien lukumäärä. Yhtälöiden (2.1)-(2.4) lisäksi sähkövarauksen tulee säilyä, joten virransyötölle asetetaan ehto

L

∑︂

l=1

I_l = 0. (2.5)

Lisäksi sähköisen potentiaalin referenssitaso tulee kiinnittää, jotta suoran mallin ratkaisu on yksikäsitteinen. Referenssitaso voidaan kiinnittää esimerkiksi asettamalla elektrodien potenti- aaleilleU_l ehto

L

∑︂

l=1

U_l = 0. (2.6)

Kun sähkönjohtavuus σ(⃗r), kontakti-impedanssit zl ja virransyöttö on määrätty, yhtälöiden (2.1)-(2.6) ratkaisuna saadaan kohteen sisäinen potentiaali u(⃗r) ja elektrodien potentiaalit Ul.

Suoran mallin yhtälö (2.1) voidaan johtaa Maxwellin yhtälöistä. Kohteen oletetaan koostu- van lineaarisesta ja isotrooppisesta väliaineesta, ja syötetyn sähkövirran oletetaan olevan aika- harmoninen funktio. Lisäksi sähkömagneettisen induktion vaikutus sähkökenttään jätetään huo- miotta, mikä on ERT:ssä käytetyillä matalilla taajuuksilla hyvä approksimaatio. Tällöin sähkö- kenttäE(⃗⃗ r)approksimoidaan pyörteettömäksi, ja se voidaan esittää muodossaE(⃗⃗ r) = −∇u(⃗r). Myös kapasitiiviset ilmiöt jätetään huomiotta, eli sähkövirran approksimoidaan muodostuvan vain varausten liikkeestä eikä lisäksi sähkövuon tiheyden ajallisesta muutoksesta. [20]

sa ∇ ·j(⃗r) = 0, mikä tarkoittaa että kohteen sisällä ei ole sähkövirran lähteitä. Ensimmäinen reunaehto (2.2) tarkoittaa, että elektrodin el potentiaali Ul poikkeaa kohteen potentiaalista u(⃗r)elektrodin alla. Potentiaalin epäjatkuvuus johtuu elektrodin ja kohteen välisestä kontakti- impedanssista z_l. Kontakti-impedanssi aiheutuu käytännössä siitä, että elektrodien ja kohteen rajapinnalla sähkövirta muuttuu elektronien virtauksesta ionien virtaukseksi [5]. Toinen reu- naehto (2.3) tarkoittaa, että virrantiheyden normaalikomponentinj_ν(⃗r) =σ(⃗r)^∂u_∂ν(⃗r) integraali elektrodin e_l pinnan yli on kyseisen elektrodin läpi syötetyn sähkövirran suuruinen. Viimeinen reunaehto (2.4) tarkoittaa, että virrantiheyden normaalikomponenttij_ν(⃗r) = 0kohteen reunalla muualla kuin elektrodien kohdalla.

Yllä esitetyssä mallissa reunaehdot (2.2)-(2.4) muodostavat niin sanotun täydellisen elekt- rodimallin [21], jossa otetaan huomioon elektrodien ja kohteen väliset kontakti-impedanssit.

Täydellisen elektrodimallin on havaittu mallintavan mittauksia mittaustarkkuuden rajoissa, toisin kuin yksinkertaisempien reunaehtojen [22]. Yksinkertaisempia reunaehtoja ovat esimer- kiksi jatkuva malli, jossa sähkövirran tiheyden normaalikomponentti on kohteen reunalla jat- kuva funktio, ja oikosulkumalli, jossa ei oteta huomioon elektrodien ja kohteen välisiä kontakti- impedansseja.

2.2 Suoran mallin ratkaisun FEM-approksimaatio

ERT:n suoralle mallille pystytään muodostamaan analyyttinen ratkaisu vain tapauksissa, jois- sa kohteen geometria ja sähkönjohtavuusjakauma ovat riittävän yksinkertaisia [23, 24]. Ana- lyyttisen ratkaisun sijasta mallin ratkaisua voidaan approksimoida numeerisilla menetelmil- lä, kuten elementtimenetelmällä (Finite Element Method, FEM). ERT:n suoran mallin FEM- approksimaatiota on käsitelty esimerkiksi lähteissä [5, 8, 20, 25, 26]. Tässä kappaleessa FEM- approksimaation periaate esitetään lyhyesti.

Elementtimenetelmässä käytetään suoran mallin variationaalimuotoa, jota sanotaan myös mallin heikoksi muodoksi. Variationaalimuodossa suora malli esitetään integraaliyhtälönä, jossa esiintyy sähköisen potentiaalin u(⃗r) ja elektrodien potentiaalien U_l lisäksi testifunktiot v ja V_l, joille pätee samat reunaehdot kuin u(⃗r):lle ja U_l:lle. Variationaalimuodolle muodostetaan diskreetti approksimaatio, joka kirjoitetaan usean eri testifunktioparin (v, V_l) avulla. Tällöin päädytään yhtälöryhmään, jonka ratkaisuna saadaan approksimaatiotu(⃗r):lle ja U_l:lle.

Variationaalimuodon johtaminen aloitetaan kertomalla yhtälöä (2.1) puolittain testifunk- tiolla v ja integroimalla saatu yhtälö Ω:n yli. Vastaavasti yhtälö (2.2) kerrotaan puolittain testifunktiolla V_l, ja saatu yhtälö integroidaan elektrodine_l yli. Näin saadut kaksi yhtälöä voi- daan sieventää käyttämällä reunaehtoja (2.3)-(2.4) sekä Greenin kaavaa. Yhdistämällä tulokset

saadaan suoran mallin (2.1)-(2.4) variationaalimuoto, joka on

B((u, U),(v, V)) =

L

∑︂

l=1

I_lV_l, (2.7)

missä (v, V) ∈ H ovat testifunktioita, H = H¹(Ω)×R^L ja H¹(Ω) on Sobolev-avaruus. Varia- tionaalimuodossa esiintyvä bilineaarimuoto B :H×H →R on

B((u, U),(v, V))) =

∫︂

Ω

σ∇u· ∇vdx+

L

∑︂

l=1

1 z_l

∫︂

e_l

(u−U_l)(v−V_l)dS. (2.8) On osoitettu, että variationaalimuodolle on olemassa yksikäsitteinen ratkaisu ja että element- timenetelmä on numeerisesti stabiili [21].

Elementtimenetelmässä variationaalimuodon (2.7) ratkaisulle haetaan approksimaatiota ää- rellisulotteisesta ratkaisuavaruudesta, joka määrätään valitsemalla avaruuden kantafunktiot.

Kannan valinnan lisäksi Ω jaetaan äärellisen kokoisiin elementteihin, jotka ovat tässä työssä tetraedrin muotoisia. Sähköinen potentiaali esitetään kantafunktioiden lineaarikombinaationa

u(⃗r)≈u^h(⃗r) =

Nu

∑︂

i=1

β_iφ_i(⃗r), (2.9)

missä N_u on kantafunktioiden lukumäärä, β_i ∈ R ja funktiot φ_i ovat kantafunktioita. Tässä työssä kanta valitaan siten, että kantafunktioita φ_i on yhtä monta kuin elementtihilassa on solmupisteitä, ja kantafunktioille pätee ehto

φ_i(⃗r_j) =

⎧

⎨

⎩

1, kuni=j 0, kuni̸=j

. (2.10)

Tällöin solmupisteessä⃗r_j potentiaalin arvo on β_j. Kantafunktiot voivat olla esimerkiksi lineaa- risia funktioita, jolloinφ_i:n arvo laskee lineaarisesti nollaan siirryttäessä solmupisteestä⃗r_i kohti mitä vain viereistä solmupistettä ⃗r_j. Tässä työssä käytetään 2. asteen kantafunktioita, joilla potentiaalia voidaan approksimoida tarkemmin kuin lineaarisilla kantafunktioilla elementtien määrän ollessa vakio [26]. Suoran mallin numeerista ratkaisua varten myös sähkönjohtavuusja- kauma esitetään muodossa

σ(⃗r)≈σ^h(⃗r) =

Nσ

∑︂

i=1

σ_iψ_i(⃗r), (2.11)

missä N_σ on kantafunktioiden lukumäärä. Potentiaali u(⃗r) ja sähkönjohtavuus σ(⃗r) voidaan esittää eri kannoissa. Tässä työssä kantafunktiot ψ_i ovat lineaarisia. Elektrodien potentiaalit

esitetään muodossa

U ≈U^h =

L−1

∑︂

j=1

γjnj, (2.12)

missäγ_j ∈Rjan₁ = (−1,1,0, ...,0)^T,n₂ = (−1,0,1,0, ...,0)^T ∈R^L×1 jne. Tällä kantavektorien nj valinnalla varmistetaan, että ehto (2.6) pätee. Tällä valinnalla elektrodien potentiaalit ovat

U₁ =

L−1

∑︂

l=1

γ_l

U₂ =−γ₁ U₃ =−γ₂

...

U_L=−γL−1.

Äärellisulotteisen ratkaisuavaruuden kantafunktiot ovat (φ_i,0), ..,(φ_Nu,0),(0, n₁), ..,(0, nL−1), ja näitä funktioita käytetään myös testifunktioina variationaalimuodossa. Variationaalimuotoon (2.7) sijoitetaan sähkönjohtavuuden σ, elektrodien potentiaalien U ja sähköisen potentiaalin u kantafunktioesitykset. Variationaalimuodon täytyy päteä erikseen jokaiselle testifunktiolle, joten siitä saadaan lineaarinen yhtälöryhmä. Yhtälöryhmä voidaan esittää matriisimuodossa

Aθ =I, (2.13)

missä θ = (β, γ)^T, β = (β₁, ..., β_Nu) ja γ = (γ₁, ..., γL−1). Lisäksi I = (0^1×Nu,I)̂^T, Î= (I₁ − I₂, I₁−I₃, ..., I₁−I_L). Matriisi A on muotoa

A=

(︄B C C^T D

)︄

, (2.14)

missä lohkomatriisit B, C ja D ovat

B(i, j) =

∫︂

Ω

(︄_Nσ

∑︂

k=1

σkψk(⃗r) )︄

∇φi· ∇φjdx+

L

∑︂

l=1

1 z_l

∫︂

el

φiφjdS, i, j = 1, ..., Nu (2.15) C(i, j) =−

(︄1 z₁

∫︂

e1

φ_idS− 1 z_j+1

∫︂

ej+1

φ_idS )︄

, i= 1, ..., N_u, j = 1, ..., L−1 (2.16)

D(i, j) =

L

∑︂

l=1

1 z_l

∫︂

e_l

(n_i)_l(n_j)_ldS=

⎧

⎨

⎩

|e1|

z1 , i̸=j

|e1|

z1 +^|e_z^j+1|

j+1 , i=j

i, j = 1, ..., L−1, (2.17) missä|e_i|on elektrodinipinta-ala. Yhtälöissä (2.15)-(2.17) esiintyvien integraalien laskemiseen käytetään numeerisia menetelmiä [25].

Diskretoidun variationaalimuodon ratkaisun θ = A⁻¹I laskemisen jälkeen potentiaali u(⃗r) ja elektrodien potentiaalit Ul voidaan esittää yhtälöiden (2.9) ja (2.12) avulla. Elektrodien po- tentiaaleista U_l voidaan laskea mittauksia vastaavat elektrodien väliset jännitteet V yhtälöllä V =M U, missä M on mittauksia vastaava dierenssimatriisi. Esimerkiksi jos jännitteet mita- taan vierekkäisten elektrodien väliltä,M on muotoa

M =

⎡

⎢

⎢

⎢

⎢

⎣

1 −1 0 · · · 0 0 1 −1 0 · · · 0

...

0. . . 0 1 −1

⎤

⎥

⎥

⎥

⎥

⎦ .

Suoran mallin FEM-approksimaatiota käytetään käänteisongelman ratkaisumenetelmissä. Kään- teisongelman ratkaisuna saadaan sähkönjohtavuusjakauman σ(⃗r) kantafunktioesityksen (2.11) kertoimet σ_i. Jatkossa näille kertoimille käytetään merkintääσ = (σ₁, ..., σ_Nσ)^T.

2.3 Käänteisongelma

Resistanssitomograan käänteisongelma on ratkaista kohteen sisäinen sähkönjohtavuusjakauma σ(⃗r)ja elektrodien kontakti-impedanssitzl, kun virransyöttöjä vastaavat jännitemittaukset on tehty. Kontakti-impedanssien estimointi samanaikaisesti sähkönjohtavuusjakauman kanssa on tärkeää käytettäessä oikeaa mittausdataa, koska väärin asetetut kontakti-impedanssit aiheutta- vat virheitä sähkönjohtavuusjakauman estimaattiin [5]. Toisaalta estimoinnin sijasta kontakti- impedansseille voitaisiin asettaa approksimatiiviset arvot, ja ottaa tähän liittyvä approksimaa- tiovirhe huomioon estimoinnissa [27]. Tässä työssä käytetään simuloitua mittausdataa, ja ole- tetaan yksinkertaisuuden vuoksi kontakti-impedanssit z_l tunnetuiksi. Jatkossa ERT:n suora malli kirjoitetaan lyhyesti muodossaV =U(σ(⃗r);z), missä puolipisteellä erotetaan tuntematon sähkönjohtavuusjakauma σ(⃗r) ja tunnettut kontakti-impedanssitz.

ERT:n käänteisongelma on huonosti asetettu, mikä tarkoittaa että usealla hyvin erilaisel- la sähkönjohtavuusjakaumalla σ(⃗r) simuloidut jännitemittaukset sopivat mitattuun jännite- dataan V mittaustarkkuuden rajoissa [16]. Tämän vuoksi käänteisongelmaa ei voida ratkais- ta esimerkiksi pienimmän neliösumman menetelmällä, jossa etsitään johtavuusjakauma jolle

||V−U(σ(⃗r);z)||² minimoituu. Huonosti asetettujen käänteisongelmien tapauksessa pienimmän neliösumman ratkaisu on hyvin herkkä mittaus- ja mallinnusvirheille, joten ratkaisu ei vastaa kohteen oikeaa sähkönjohtavuusjakaumaa. Huonosti asetettujen käänteisongelmien ratkaisemi- seen voidaan käyttää regularisointimenetelmiä, joissa ongelmaa muokataan siten, että saadaan hyvin asetettu ongelma. ERT:ssä regularisointia on käytetty esimerkiksi lähteissä [17,20,25,26].

Regularisoinnin sijaan käänteisongelmien ratkaisemiseen voidaan käyttää myös Bayesilaista eli tilastollista lähestymistapaa, jossa tuntemattomia muuttujia käsitellään satunnaismuuttu- jina. Bayesilaisessa lähestymistavassa käänteisongelman ratkaisu on suoran mallin ,

mittausvirheen statistiikan ja prioritiedon perusteella muodostettu tuntematonta muuttujaaσ kuvaava todennäköisyystiheysfunktio, jota sanotaan posterioritiheydeksi. Usein lasketaan niin sanottu Maximum a Posteriori-estimaatti (MAP), joka on posterioritiheyden maksimikohta.

Koko posterioritiheyden tarkastelu sallii myös estimaatin luotettavuuden arvioinnin, mutta tä- mä tarkastelu vaatii posterioritiheyden näytteistämistä Markovin ketjuihin perustuvilla mene- telmillä. Pelkän MAP-estimaatin laskeminen on huomattavasti helpompi toteuttaa kuin poste- rioritieyden näytteistäminen. ERT:n käänteisongelmaan liittyvää posterioritiheyden näytteistä- mistä on tehty esimerkiksi lähteessä [28]. Tässä kappaleessa esitetään MAP-estimaatin lasken- nan periaate.

Bayesilaisessa lähestymistavassa käänteisongelman ratkaisemiseen käytetään mittausdatan lisäksi saatavilla olevaa ennakkotietoa kohteen sähkönjohtavuusjakaumasta. Ennakkotietoa ku- vataan sähkönjohtavuuden prioritodennäköisyystiheydellä π_σ(σ). Ennakkoon voidaan esimer- kiksi tietää, että kohteen sähkönjohtavuus on sileä, tai että sähkönjohtavuusjakauma on epäjat- kuva [29]. Lisäksi saatavilla voi olla kvantitatiivista ennakkotietoa, kuten tietoa sähkönjohta- vuusjakauman mahdollisista arvoista. Tässä työssä sähkönjohtavuudelle käytetään katkaistuja Gaussisia prioreja, jotka ovat muotoa

π_σ(σ)∝Ξ(σ) exp {︃

−1

2(σ−σ∗)^TΓ⁻¹_σ (σ−σ∗) }︃

, (2.18)

missäΓ_σ ja σ_∗ ovat Gaussisen satunnaismuuttujanσ kovarianssimatriisi ja odotusarvo. Lisäksi Ξ(σ) on positiivisuuspriori, jonka määritelmä on

Ξ(σ) =

⎧

⎨

⎩

1, kun min(σ)≥0 0, kun min(σ)<0

. (2.19)

Priorin valintaa tässä työssä käsitellään tarkemmin Kappaleessa 4.

Bayesilaisessa lähestymistavassa mittausvirhettä käsitellään satunnaismuuttujana, jolla on todennäköisyystiheys π_e(e). Kirjoitetaan ERT:n havaintomalli lyhyesti muodossa

V =U(σ;z) +e, (2.20)

missä V on mitatut jännitteet, U(σ;z) on ERT:n suoran mallin FEM-approksimaatio, ja e on additiiviseksi oletettu mittausvirhe. Tässä työssä mittausvirheen oletetaan olevan Gaussinen, joten sen tiheysfunktio on muotoa

π_e(e)∝exp {︃

−1

2(e−e∗)^TΓ⁻¹_e (e−e∗) }︃

, (2.21)

missä Γ_e on virheen kovarianssimatriisi ja e∗ on virheen odotusarvo. Mittausvirheen e olete- taan olevan riippumaton sähkönjohtavuudesta σ. Yhteistiheys π_V,σ,e(V, σ, e) voidaan kirjoittaa

ehdollisen tiheyden määritelmän avulla muodossa

π_V,σ,e(V, σ, e) = π_V|σ,e(V|σ, e)πe|σ(e|σ)π_σ(σ)

=δ(V −U(σ;z)−e)πe|σ(e|σ)π_σ(σ), (2.22) missä δ on Diracin deltafunktio ja π_V|σ,e(V|σ, e) = δ(V −U(σ;z)−e) havaintomallin (2.20) perusteella, kunσjaeon kiinnitetty. Toisaalta ehdollisen tiheyden määritelmän avulla voidaan myös kirjoittaa

π_V,σ,e(V, σ, e) = πV,e|σ(V, e|σ)π_σ(σ). (2.23)

Mittausten V ehdollinen tiheys on yhtälöiden (2.22) ja (2.23) nojalla π_V|σ(V|σ) =

∫︂

π_V,e|σ(V, e|σ)de

=

∫︂

δ(V −U(σ;z)−e)π_e|σ(e|σ)de

=πe|σ(V −U(σ;z)|σ)

=π_e(V −U(σ;z)). (2.24)

Integraali on tässä ratkaistu muodollisesti. Kyseessä on Lebesgue-integraali, jossa Diracin delta- funktiota käytetään mittana jonka suhteen integraali lasketaan. Viimeisessä välivaiheessa käy- tettiin oletusta, että e jaσ ovat keskenään riippumattomia. Mittausvirheene tiheysfunktio on muotoa (2.21), joten

π_V|σ(V|σ) = π_e(V −U(σ;z))

=∝exp {︃

−1

2(V −U(σ;z)−e∗)^TΓ⁻¹_e (V −U(σ;z)−e∗) }︃

, (2.25)

missä M on jännitemittausten lukumäärä.

Tiheysfunktiotaπ(V|σ)kutsutaan likelihood-funktioksi [16]. Tämän funktion avulla voidaan laskea ns. maximum likelihood-estimaatti (ML-estimaatti), joka on

̂

σ_ML = arg max

σ

π_V|σ(V|σ). (2.26)

Maximum likelihood-estimaatin laskeminen vastaa painotettua pienimmän neliösumman mene- telmää, jossa mittausvirheen kovarianssimatriisiaΓ_e käytetään eri mittausten suhteelliseen pai- nottamiseen. Resistanssitomograan käänteisongelma on huonosti asetettu, joten ML-estimaatti ei käytännössä vastaa kohteen oikeaa sähkönjohtavuusjakaumaa [30]. Likelihood-funktion lisäksi estimoinnissa täytyy käyttää sähkönjohtavuuden prioritiheyttäπ_σ(σ), jolloin saadaan muodos-

tettua posterioritiheys.

Posterioritiheys saadaan kirjoittamalla sähkönjohtavuuden ja mitattujen jännitteiden yh- teistiheys kahdessa muodossa ehdollisen tiheyden määritelmän avulla.

πσ,V(σ, V) = πV|σ(V|σ)πσ(σ) = πσ|V(σ|V)πV(V). (2.27) Tästä saadaan Bayesin kaava

πσ|V(σ|V) = π_σ(σ)π_V|σ(V|σ)

π_V(V) ∝π_σ(σ)π_V|σ(V|σ), (2.28) missä π_V(V)on normalisointivakio, joka voidaan jättää huomiotta. Tiheysfunktiota πσ|V(σ|V) kutsutaan posterioritiheydeksi [16]. Yhtälöistä (2.18) ja (2.25) saadaan posterioritiheydelle muoto

π_σ|V(σ|V)∝Ξ(σ) exp {︃

−1

2(V −U(σ;z)−e_∗)^TΓ⁻¹_e (V −U(σ;z)−e_∗)

−1

2(σ−σ_∗)^TΓ⁻¹_σ (σ−σ_∗) }︃

.

(2.29)

Posterioritiheydestä voidaan laskea MAP-estimaatti, joka on posterioritiheyden maksimikohta.

MAP-estimaatti siis minimoi posterioritiheyden (2.29) negatiivisen eksponentin. Huomioimalla sähkönjohtavuuden σ positiivisuusrajoite, saadaan

̂

σ_MAP= arg max

σ

π(σ|V)

= arg min

σ≥0

{(V −U(σ;z)−e∗)^TΓ⁻¹_e (V −U(σ;z)−e∗) + (σ−σ∗)^TΓ⁻¹_σ (σ−σ∗)}

= arg min

σ≥0

{(V −U(σ;z)−e∗)^TL^T_eL_e(V −U(σ;z)−e∗) + (σ−σ∗)^TL^T_σL_σ(σ−σ∗)}

= arg min

σ≥0

{(L_e(V −U(σ;z)−e∗)^TL_e(V −U(σ;z)−e∗) + (L_σ(σ−σ∗))^TL_σ(σ−σ∗)}

= arg min

σ≥0

{||L_e(V −U(σ;z)−e_∗)||²+||L_σ(σ−σ_∗)||²}, (2.30) missä Le ja Lσ ovat matriisien Γ⁻¹_e ja Γ⁻¹_σ Cholesky-tekijät. MAP-estimaatin laskeminen on siis optimointiongelma, joka voidaan ratkaista esimerkiksi Gauss-Newton-algoritmilla. Ilman positiivisuusrajoitetta Gauss-Newton-iteraatio olisi muotoa

̂

σk+1 = ̂σk+dk(J^TΓ⁻¹_e J+L^T_σLσ)⁻¹(J^TΓ⁻¹_e (V −U(̂σk;z)−e∗)−Γ⁻¹_σ (̂σk−σ∗))

=: ̂σ_k+d_kp_k, (2.31)

missä σ̂_k on estimaatti iteraatiokierroksella k, d_k ∈ R on askelpituus ja J on mallin U(σ;z) Jacobin matriisi laskettuna pisteessä σ̂_k. Iterointi aloitetaan alkuarvauksesta σ̂₀, ja estimaatti siirtyy jokaisella iterointikierroksella suuntaanp_k askelpituudend_k verran. Tässä työssä alkuar-

vauksena käytetään priorin odotusarvoa. Lisäksi suuntavektori pk lasketaan ensin, jonka jäl- keen askelpituus dk valitaan viivahakumenetelmällä. Jacobin matriisin J laskenta σ:n suhteen on esitetty esimerkiksi lähteessä [5].

Positiivisuusrajoite toteutetaan tässä työssä menetelmällä, jossa minimoitavaan funktionaa- liin (2.30) lisätään sakkofunktioq(σ), jolloin minimointiongelma on muotoa

̂

σ_MAP = arg min

σ≥0

{||L_e(V −U(σ;z)−e_∗)||²+||L_σ(σ−σ_∗)||² +q(σ)}. (2.32) Sakkofunktio q(σ)on tässä työssä muotoa

q(σ) = ∑︂

σ_min,2≤σ_i≤σ_min,1

aσ_i²+bσ_i+c, (2.33)

missä σ_min,2 ja σ_min,1 ovat pieniä positiivisia lukuja siten että σ_min,2 < σ_min,1. Lisäksi a, b ja c ovat vakioita, joiden arvot voidaan laskea kun kiinnitetään arvot σ_min,2, σ_min,1, q(σ_min,2)ja q:n nollakohta q(σ_min,1) = aσ²_min,1 +bσ_min,1+c= 0. Sakkofunktioq(σ) siis kasvattaa minimoitavan funktionaalin (2.32) arvoa josσ̂_k sisältää arvoja, jotka ovat pienempiä kuinσ_min,1. Iteraatiokier- roksellak vektorinσ̂_k negatiiviset komponentit projisoidaan arvoon σ_min,2, jolloin sakkofunktio pakottaa näitä komponentteja kasvamaan iteraatiokierroksellak+ 1. Positiivisuusrajoite toteu- tetaan siis sisäpistemenetelmällä (Interior Point Method), mikä tarkoittaa että σ:n estimaat- ti pysyy jokaisella iteraatiokierroksella ei-negatiivisena. Sisäpistemenetelmää käytetään, koska suoran mallin ratkaiseminen vaatii ei-negatiivisen sähkönjohtavuusjakauman.

Sakkofunktion q(σ)gradientin komponentit ovat

[∇q]i =

⎧

⎨

⎩

0, kunσ_i > σ_min,1

2aσi+b, kunσi ≤σ_min,1

, (2.34)

ja q(σ):n Hessen matriisi on diagonaalimatriisi, jonka diagonaalialkiot ovat muotoa

[Hq]i,i =

⎧

⎨

⎩

0, kunσ_i > σ_min,1 2a, kunσi ≤σ_min,1

. (2.35)

Optimointiongelman (2.32) ratkaisun Gauss-Newton iteraatio on muotoa

̂

σk+1 = ̂σk+dk(J^TΓ⁻¹_e J +L^T_σLσ +Hq)⁻¹(J^TΓ⁻¹_e (V −U(̂σk;z)−e∗)−Γ⁻¹_σ (̂σk−σ∗)− ∇q)

=: ̂σ_k+d_kp_k. (2.36)

Iteraatio on siis samanlainen kuin ilman positiivisuusrajoitteen huomiointia (2.31), mutta nyt iteraatioon on lisätty sakkofunktion q(σ)gradientti ja Hessen matriisi.

ERT:n suoran ongelman epälineaarisuudesta sekä ei-negatiivisuudesta johtuen posterioriti- heys ei ole Gaussinen, joten MAP-estimaatin luotettavuuden tarkka arviointi vaatii posterio- ritiheyden näytteistämistä. Näytteistämisen sijaan MAP-estimaatille voidaan laskea approksi- matiiviset luottamusvälit approksimoimalla posterioritiheyttä Gaussisena. Tällöin posterioriti- heyden kovarianssille saadaan approksimaatio [16]

Γ_σ|V ≈(Γ⁻¹_σ +J^TΓ⁻¹_e J)⁻¹, (2.37) missä Jacobin matriisiJ on sama kuin yhtälössä (2.31), mutta laskettu pisteessäσ̂_MAP. Approk- simatiivisesta posterioritiheyden kovarianssista voidaan laskea kullekin MAP-estimaatin kom- ponentilleσ̂_MAP,k esimerkiksi kolmen keskihajonnan luottamusväli, joka on [σ̂_MAP,k−3√︁

Γ_σ|V,kk,

̂

σ_MAP,k+ 3√︁

Γσ|V,kk].

3 Approksimaatiovirhemenetelmä

Kun käänteisongelmien ratkaisemiseen käytetään Bayesilaista lähestymistapaa, malleissa esiin- tyviä virheitä voidaan ottaa huomioon approksimaatiovirhemenetelmällä [16,31]. Osittaisdie- rentiaaliyhtälöihin perustuvissa malleissa virhettä aiheutuu esimerkiksi mallin diskretoinnista elementtimenetelmällä ratkaisemista varten. Mallin reunaehdot voivat olla huonosti tunnettuja, etenkin jos kohteen reunan muoto on myös huonosti tunnettu. Reunan muodon lisäksi kuvan- tamisongelmissa sensorien paikat voivat olla epätarkasti tunnettuja. Malli voi sisältää huonosti tunnettuja parametreja, joille voidaan ennakktiedon perusteella asettaa vain approksimatiiviset arvot. Virheellisen suoran mallin käyttö käänteisongelman ratkaisemisessa voi johtaa huonoon ratkaisuun, joka ei vastaa kohdetta.

Resistanssitomograassa approksimaatiovirhemenetelmää on sovellettu usean erityyppisen mallinnusvirheen huomiointiin. Näitä virheitä ovat diskretointivirhe [19,27,32], approksimatii- visen reunan muodon aiheuttama virhe [1,19,27,3234], tuntemattomien reunaehtojen aiheut- tama virhe [19] ja tuntemattomien kontakti-impedanssien aiheuttama virhe [27]. Menetelmää on myös käytetty tilanteessa, jossa suoraa mallia on yksinkertaistettu approksimoimalla sähkön- johtavuuden olevan muuttumaton jonkin paikkakoordinaatin suhteen [35]. Lisäksi menetelmää voidaan käyttää tilanteessa, jossa suoran mallin FEM-approksimaatiota yksinkertaistetaan esit- tämällä sähkönjohtavuus ja sähköinen potentiaali pienellä määrällä kantafunktioita [36]. Ap- proksimaatiovirhemenetelmällä voidaan käsitellä useaa erityyppistä mallinnusvirhettä saman- aikaisesti [27].

Approksimaatiovirhemenetelmää on sovellettu ERT:n lisäksi useisiin muihin käänteisongel- miin. Approksimaatiovirheitä on käsitelty ERT:n lisäksi esimerkiksi diuusissa optisessa to- mograassa, jossa kohteen absorptio- ja sirontakerroinjakaumat estimoidaan kohteen reunalla tehtävien optisten mittausten perusteella. Tässä tapauksessa approksimaatiovirhemenetelmää on sovellettu esimerkiksi diskretointivirheen [37], tuntemattomaan reunan muodon [3739], ja säteilynsiirtoyhtälön ja diuusioapproksimaation välisen mallinnusvirheen huomiointiin [40].

Fluoresenssitomograassa approksimaatiovirhemenetelmää on sovellettu approksimatiivisesti mallinnettujen väliaineen optisten ominaisuuksien aiheuttaman virheen huomiointiin [41].

Approksimaatiovirhemenetelmää on käytetty myös tilanteissa, joissa jokin suorassa mallissa esiintyvä jakauma approksimoidaan homogeeniseksi, kun estimoidaan tämän jakauman sijas- ta jotain muuta mallissa esiintyvää suuretta. Lähteessä [42] käytettiin approksimaatiovirhe- menetelmää diuusissa optisessa tomograassa tilanteeseen, jossa kohteen sirontakerroinjakau- maa approksimoitiin homogeenisena ja estimoitiin vain absorptiokerroinjakauma. Lähteessä [43]

approksimoitiin väliaineen permittiivisyysjakauma homogeeniseksi, kun estimoitiin sähkömag- neettisia aaltoja sirottavan kohteen sijaintia.

Approksimaatiovirheitä on käsitelty myös epästationäärisissä käänteisongelmissa, joissa ajas- ta riippuvaa kohteen tilaa estimoidaan eri ajanhetkillä tehtyjen mittausten perusteella [4,6,44,

45].

Puutteellisesta tiedosta johtuvien mallinnusvirheiden lisäksi malli voi sisältää tietoisesti tehtyjä mallinnusvirhetä. Käyttämällä harvaa diskretointia suorassa mallissa voidaan vähen- tää käänteisongelman ratkaisuun tarvittavaa laskentatehoa- ja aikaa. Harvan diskretoinnin li- säksi laskentaa voidaan nopeuttaa asettamalla mallin ei-kiinnostaville parametreille approk- simatiiviset arvot, ja estimoimalla ainoastaan kiinnostavat tuntemattomat parametrit. Jos ei- kiinnostaville parametreille asetetut arvot poikkeavat niiden oikeista arvoista, käänteisongelman ratkaisu kiinnostaville parametreille voi poiketa hyvin paljon parametrien oikeista arvoista. Ap- proksimaatiovirhemenetelmässä otetaan huomioon suorassa mallissa esiintyvät virheet, jolloin kiinnostaville parametreille saadaan tarkempia estimaatteja.

3.1 Approksimaatiovirheen mallintaminen

Suora malli voi sisältää käänteisongelman ratkaisun kannalta ei-kiinnostavia tuntemattomia parametreja, joita merkitään jatkossa vektorilla ξ. Käytetäänσ:n jaξ:n yhteiselle posterioriti- heydelle merkintää π_σ,ξ|V(σ, ξ|V). Tämä posterioritiheys voidaan marginalisoida integroimalla ei-kiinnostavien parametrien ξ suhteen

πσ|V(σ|V) =

∫︂

πσ,ξ|V(σ, ξ|V)dξ, (3.1)

jolloin saadaan posterioritiheys sähkönjohtavuudelle σ. Tämä tapa σ:n estimointiin vaatii ti- heysfunktionπ_σ,ξ|V näytteistämistä, mikä on laskennallisesti raskasta. Approksimaatiovirheme- netelmässä posterioritiheysπσ,ξ|V marginalisoidaan approksimatiivisesti parametrienξ suhteen etukäteen [42]. Tällöin parametreille ξ voidaan asettaa approksimatiiviset arvot, eikä niitä tar- vitse estimoida samanaikaisesti σ:n MAP-estimaatin laskennan kanssa. Tässä kappaleessa esi- tellään approksimaatiovirhemenetelmän käyttö.

Oletetaan, että U_Ω,δ(σ, ξ)on riittävän tarkka laskennallinen malli, jonka approksimaatiovir- he on mittausvirheeseen verrattuna pieni. Tässä mallissa kohteen Ω geometria on mallinnettu tarkasti. Alaindeksi δ viittaa FEM-hilan diskretointiin; δ voi olla esimerkiksi suurin sallittu elementin tilavuus, joka on tarkalle mallille pieni. Lisäksi ξ on parametrivektori, joka voi si- sältää mitä vain mallin parametreja, joita σ:n lisäksi tarvitaan suoran mallin ratkaisemiseen.

Tässä työssä kontakti-impedansseja z pidetään tunnettuina parametreina, joten ne muodos- tavat osan parametrivektorista ξ. Näiden lisäksi ξ voi sisältää esimerkiksi kohteen geometriaa kuvaavia parametreja. Olkoon UΩ,h̃ (σ, ξ₀) approksimatiivinen malli, jossa kohteen geometria voi olla mallinnettu virheellisesti (merkitään Ω̃), diskretointi voi olla harvempi kuin tarkassa mallissa (h≥ δ), ja parametrivektorille ξ on asetettu approksimatiivinen arvo ξ₀. Näistä mal- linnusvirheistä johtuen approksimatiivisella mallillaUΩ,h̃ (σ, ξ₀)lasketut mittaukset ovat malliin UΩ,δ(σ, ξ) verrattuna virheellisiä.

Tarkkaan laskennalliseen malliin liittyvä havaintomalli (2.20) voidaan kirjoittaa muodossa V =U_Ω,δ(σ, ξ) +e

=UΩ,h̃ (σ, ξ₀) + [U_Ω,δ(σ, ξ)−UΩ,h̃ (σ, ξ₀)] +e

=UΩ,h̃ (σ, ξ0) +ε(σ, ξ) +e

=UΩ,h̃ (σ, ξ₀) +η(σ, ξ), (3.2)

missä ε(σ, ξ) =U_Ω,δ(σ, ξ)−UΩ,h̃ (σ, ξ₀) on approksimaatiovirhe ja η(σ, ξ) =ε(σ, ξ) +e. Havain- tomalli (3.2) on siis samaa muotoa kuin havaintomalli (2.20), mutta nyt myös mallinnusvirhe on huomioitu additiivisena virheenä. Approksimaatiovirhemenetelmän periaatteena on käyttää havaintomallia (3.2), ja huomioida mallinnusvirheenεstatistiikka MAP-estimaatin laskennassa.

Johdetaan seuraavaksi posterioritiheys käyttäen havaintomallia (3.2). Mittausvirhettäe kä- sitellään riippumattomana muuttujista (σ, ξ), kuten kappaleessa 2.3. Lisäksi mittausvirheen e ja mallinnusvirheenεoletetaan olevan keskenään riippumattomia. Johdetaan aluksi likelihood- funktio. Ehdollisen todennäköisyyden avulla yhteistiheysπ(V, σ, ξ, e, ε)voidaan kirjoittaa muo- dossa

π_V,σ,ξ,e,ε(V, σ, ξ, e, ε) = π_V|σ,ξ,e,ε(V|σ, ξ, e, ε)π_σ,ξ,e,ε(σ, ξ, e, ε)

=δ(V −UΩ,h̃ (σ, ξ₀)−e−ε)πe,ε,ξ|σ(e, ε, ξ|σ)π_σ(σ), (3.3) missä π_V|σ,ξ,e,ε(V|σ, ξ, e, ε) = δ(V −UΩ,h̃ (σ, ξ₀)−e−ε) havaintomallin (3.2) perusteella, kun (σ, ξ, e, ε)on kiinnitetty (ξ:n kiinnitettyä arvoa merkitäänξ₀:lla). Toisaalta ehdollisen tiheyden määritelmän avulla voidaan myös kirjoittaa

π_V,σ,ξ,e,ε(V, σ, ξ, e, ε) = πV,ξ,e,ε|σ(V, ξ, e, ε|σ)π_σ(σ). (3.4) Yhtälöiden (3.3) ja (3.4) avulla saadaan ehdollinen tiheysfunktio π_V|σ(V|σ) vastaavasti kuin kappaleessa 2.3 (yhtälö (2.24))

π_V|σ(V|σ) =

∫︂ ∫︂ ∫︂

π_V,ξ,e,ε|σ(V, ξ, e, ε|σ)dξdedε

=

∫︂ ∫︂ ∫︂

δ(V −UΩ,h̃ (σ, ξ₀)−e−ε)π_e,ε,ξ|σ(e, ε, ξ|σ)dξdedε

=

∫︂ ∫︂

δ(V −UΩ,h̃ (σ, ξ₀)−e−ε)πe,ε|σ(e, ε|σ)dedε

=

∫︂ ∫︂

δ(V −UΩ,h̃ (σ, ξ₀)−e−ε)π_e(e)πε|σ(ε|σ)dedε

=

∫︂

π_e(V −UΩ,h̃ (σ, ξ₀)−ε)πε|σ(ε|σ)dε (3.5)

Likelihood-funktio π_V|σ(V|σ) saadaan siis kahden tiheysfunktion konvoluutiointegraalina ε:n suhteen. Näistä kahdesta tiheysfunktioista ensimmäinen onπe(V −U(σ, ξ0)), jonka huomataan olevan satunnaismuuttujan e+U(σ, ξ₀) tiheysfunktio evaluoituna pisteessä V. Tässä σ ja ξ₀ olivat kiinnitettyjä, joten U(σ, ξ₀) on vakio, joka on lisätty satunnaismuuttujaane. Konvoluu- tiointegraalin toinen tiheysfunktio on yksinkertaisesti satunnaismuuttujan ε|σ tiheysfunktio.

Kuten edellä, pidetään mittausvirhettä Gaussisena, eliπ_e(e) =N(e∗,Γ_e). Integraalissa (3.5) esiintyvälle tiheysfunktiolleπε|σvoidaan muodostaa Gaussinen approksimaatio approksimoimal- la ensin yhteistiheyttä π_ε,σ(ε, σ) Gaussisena

π_ε,σ(ε, σ)∝exp

⎧

⎨

⎩

−1 2

(︄ε−ε∗

σ−σ∗

)︄T(︄

Γ_ε Γ_εσ Γ_σε Γ_σ

)︄−1(︄

ε−ε∗

σ−σ∗

)︄⎫

⎬

⎭

. (3.6)

Tällöin saadaan, että π_ε|σ on muotoa π_ε|σ(ε|σ) =N(ε_∗|σ,Γ_ε|σ), missä [16]

ε∗|σ =ε∗+ ΓεσΓ⁻¹_σ (σ−σ∗) (3.7)

Γε|σ = Γ_ε−Γ_εσΓ⁻¹_σ Γ_σε. (3.8)

Kokonaisvirheelleηpäteeη|σ =e+ε|σ, joten sen tiheysfunktioksi saadaanπη|σ =N(η∗|σ,Γη|σ), missä

η∗|σ =e∗+ε∗+ Γ_εσΓ⁻¹_σ (σ−σ∗) (3.9) Γ_η|σ = Γ_e+ Γ_ε−Γ_εσΓ⁻¹_σ Γ_σε. (3.10) Tiheysfunktiot πe ja πε|σ ovat Gaussisia ja e ja ε ovat keskenään riippumattomia satun- naismuuttujia. Edelleen myös e+ U(σ, ξ₀) ja ε|σ ovat keskenään riippumattmia. Konvoluu- tiointegraali (3.5) on siis kahden keskenään riippumattoman ja Gaussisen satunnaismuuttujan tiheysfunktioiden konvoluutio. Tällöin tuloksena saadaan tunnetusti Gaussinen tiheysfunktio, jonka odotusarvo on satunnaismuuttujien e+U(σ, ξ₀) ja ε|σ odotusarvojen summa, ja jonka kovarianssimatriisi on vastaavasti e+U(σ, ξ₀):n ja ε|σ:n kovarianssimatriisien summa. Siten likelihood-funktioπ_V|σ(V|σ) on

π(V|σ) = N(e∗+UΩ,h̃ (σ, ξ₀) +ε∗|σ,Γ_e+ Γε|σ)

=N(UΩ,h̃ (σ, ξ₀) +η_∗|σ,Γ_η|σ)

∝exp {︃

−1

2(V −UΩ,h̃ (σ, ξ0)−η∗|σ)^TΓ⁻¹_η|σ(V −UΩ,h̃ (σ, ξ0)−η∗|σ) }︃

. (3.11)

Kun käytetään katkaistua Gaussista prioria kuten kappaleessa 2.3, posterioritiheys on muotoa π(σ|V)∝Ξ(σ) exp

{︃

−1

2(V −UΩ,h̃ (σ, ξ₀)−η∗|σ)^TΓ⁻¹_η|σ(V −UΩ,h̃ (σ, ξ₀)−η∗|σ)

−1

2(σ−σ∗)^TΓ⁻¹_σ (σ−σ∗) }︃

.

(3.12)

Huomioimallaσ:n positiivisuusrajoite kuten kappaleessa 2.3, MAP-estimaatti on

̂

σ_MAP = arg min

σ≥0

{||Lη|σ(V −UΩ,h̃ (σ, ξ₀)−η∗|σ)||²+||L_σ(σ−σ∗)||²+q(σ)}, (3.13) missä L^T_η|σL_η|σ = Γ⁻¹_η|σ. Estimaatti voidaan laskea Gauss-Newton-algoritmilla, joka on nyt muo- toa

̂

σ_k+1 = ̂σ_k+d_k(J^TΓ⁻¹_η|σJ +L^T_σL_σ+H_q)⁻¹(J^TΓ⁻¹_η|σ(V −UΩ,h̃ (̂σ_k, ξ₀)−η_∗|σ)−Γ⁻¹_σ (̂σ_k−σ_∗)− ∇q)

=: ̂σ_k+d_kp_k. (3.14)

Approksimaatiovirheen sisältävää virhemallia kutsutaan approksimaatiovirhemalliksi (Approxi- mation Error Model, AEM). Usein approksimoidaan lisäksi, että σ ja ε ovat keskenään riippu- mattomia, jolloinΓ_σε= 0 jaΓ_εσ = 0. Tällöin kokonaisvirheenη odotusarvo (3.9) ja kovarianssi (3.10) ovat muotoa

η∗|σ =e∗+ε∗ (3.15)

Γη|σ = Γ_e+ Γ_ε. (3.16)

Tällöin virhemallia kutsutaan parannelluksi virhemalliksi (Enhanced Error Model, EEM). σ:n jaε:n approksimointi riippumattomiksi satunnaismuuttujiksi on karkea approksimaatio, mutta suurimmassa osassa approksimaatiovirhemenetelmään liittyvistä tutkimuksista on käytetty tä- tä virhemallia. Tämän virhemallin käytön on havaittu parantavan saatuja MAP-estimaatteja verrattuna vain mittausvirheen sisältävän virhemalliin (2.20) käyttöön [19,33,35,43].

Käytännössä approksimaatiovirheen odotusarvo ε∗, kovarianssi Γ_ε sekä ristikovarianssit Γ_σϵ jaΓϵσ lasketaan simuloitujen approksimaatiovirheen näytteiden perusteella. Tämä simulaatio ja statistiikan laskenta voidaan suorittaa ennen kuin kuvannettavasta kohteesta tehdään mittauk- sia. Statistiikan laskeminen vaatii suoran ongelman ratkaisemista sekä tarkalla mallillaU_Ω,δ(σ, ξ) että approksimatiivisella mallilla UΩ,h̃ (σ, ξ₀). Statistiikan laskemisen jälkeen MAP-estimointiin tarvitaan vain approksimatiivista mallia UΩ,h̃ (σ, ξ₀).

Statistiikan laskenta aloitetaan generoimalla näytteitä parametrivektorille (σ, ξ)^T priori- jakaumasta. Olkoon näytteiden määrä N_s ja merkitään indeksiä l ∈ N vastaavaa näytettä

(σ^(l), ξ^(l))^T. Tätä näytettä vastaava approksimaatiovirhe on

ε^(l) =U_Ω,δ(σ^(l)_δ , ξ^(l))−UΩ,h̃ (σ_h^(l), ξ₀). (3.17) Approksimaatiovirheen odotusarvo ε∗, kovarianssi Γ_ε ja ristikovarianssit Γ_εσ ja Γ_σε lasketaan näytteiden perusteella yhtälöillä

ε_∗ = 1 N_s

Ns

∑︂

l=1

ε^(l) (3.18)

Γ_ε= 1 N_s−1

Ns

∑︂

l=1

(ε^(l)−ε∗)(ε^(l)−ε∗)^T (3.19)

Γ_εσ = 1 N_s−1

Ns

∑︂

l=1

(ε^(l)−ε∗)(σ^(l)−σ∗)^T (3.20)

Γ_σε= Γ^T_εσ. (3.21)

Jos käytetään approksimaatiota, jossa εja σ ovat keskenään riippumattomia, ei tarvitse laskea ristikovarianssimatriiseja Γ_εσ ja Γ_σε.

3.2 Approksimaatiovirheen aiheuttavien parametrien estimointi

Approksimaatiovirhemenetelmässä tuntemattomalle parametrivektorille ξ asetetaan approksi- matiivinen arvo ξ0, ja tähän approksimaatioon ja muihin virhelähteisiin liittyvä mallinnusvir- he otetaan huomioon estimoitaessa vektoriaσ. Approksimaatiovirhemenetelmän käytön sijasta tuntemattomat parametrit ξ voitaisiin estimoida samanaikaisesti σ:n kanssa, mutta erityises- ti kohteen geometriaa kuvaavien parametrien MAP-estimointi on vaikeaa, koska tällöin suoraa mallia täytyy derivoida näiden parametrien suhteen [46,47]. Geometriaa kuvaavien parametrien estimointi iteratiivisilla menetelmillä on laskennallisesti raskasta, koska parametrien estimaatin muuttuessa kohteelle täytyy muodostaa uusi FEM-hila.

Käytettäessä approksimaatiovirhemenetelmää parametreille ξ voidaan muodostaa MAP- estimaatti σ:n estimoinnin jälkeen, jolloin parametreja ξ ei tarvitse estimoida samanaikaisesti σ:n kanssa [1, 34, 48]. Menetelmä perustuu approksimaatiovirheen ε realisaation estimointiin samanaikaisesti kiinnostavan parametrivektorin σ kanssa, jonka jälkeen estimoidun approksi- maatiovirheen realisaation perusteella saadaan estimaatti parametreilleξ. Tämä menetelmäξ:n estimointiin on laskennallisesti kevyempi kuin ξ:n suora estimointi samanaikaisesti σ:n kans- sa. Menetelmää on käytetty resistanssitomograassa kohteen ulkoreunan muodon estimointiin jälkikäteen, kunσ:n estimaattia laskettaessa käytetään ulkoreunan muodolle yksinkertaista ap- proksimaatiota [1,34,48]. Tässä kappaleessa esitelläänξ:n estimoinnin periaate approksimaatio- virhemenetelmän avulla. Kappaleessa käytetään approksimaatiota, jossa approksimaatiovirhe

ε ja sähkönjohtavuus σ ovat keskenään riippumattomia.

Parametrivektorinξestimoinnissa approksimaatiovirhemenetelmällä hyödynnetään approk- simaatiovirheen kovarianssimatriisin Γ_ε ominaisarvohajotelmaa

Γε =W DW⁻¹, (3.22)

missä matriisin W sarakkeet ovat Γ_ε:n ominaisvektorit ja diagonaalimatriisin D päädiagonaa- lilla ovat Γε:n ominaisarvot. Kovarianssimatriisit ovat symmetrisiä, joten matriisi W on orto- gonaalinen eli W⁻¹ =W^T. Tällöin Γ_ε:n ominaisarvohajotelma voidaan kirjoittaa summana

Γ_ε =

m

∑︂

k=1

λ_kw_kw^T_k, (3.23)

missä λ_k, k = 1, ..., m ovat ominaisarvot järjestyksessä suurimmasta pienimpään, ja pystyvek- toritw_k ovat vastaavat ominaisvektorit. Ominaisvektorit muodostavat kannan avaruudelleR^m, joten approksimaatiovirhe ε voidaan esittää yksikäsitteisesti muodossa

ε=ε∗+

p

∑︂

k=1

α_kw_k+

m

∑︂

j=p+1

β_jw_j =:ε∗+ε₁+ε₂, (3.24) missä kertoimet α_k, β_j saadaan skalaariprojektioina, α_k =⟨ε−ε_∗, w_k⟩,β_j =⟨ε−ε_∗, w_j⟩.

Josε₁riittää kuvailemaan approksimaatiovirheen vaihtelua, approksimaatiovirheen realisaa- tiolle saadaan pienidimensioinen estimaatti, kun kertoimet α =: (α₁, α₂, ..., α_p)^T estimoidaan samanaikaisestiσ:n kanssa. Estimointia varten kertoimille α muodostetaan Gaussinen prioriti- heys, jonka odotusarvo ja kovarianssi saadaan ominaisarvohajotelman ominaisuuksien avulla.

Kertoimen α_k odotusarvo on

E{α_k}=⟨E{ε} −ε∗, w_k⟩

=⟨0, w_k⟩

= 0.

Kertoimienα kovarianssimatriisi on siten muotoa

Γ_α =

⎡

⎢

⎢

⎢

⎢

⎣

E{α₁α₁} E{α₁α₂}. . . E{α₁α_p} E{α₂α₁} E{α₂α₂}. . . .

... ... ...

E{α_pα₁} . . . E{α_pα_p}

⎤

⎥

⎥

⎥

⎥

⎦ ,

missä diagonaalialkiot ovat

E{α_kα_k}=E{⟨ε−ε∗, w_k⟩⟨ε−ε∗, w_k⟩}

=E{⟨w_k, ε−ε∗⟩⟨ε−ε∗, w_k⟩}

=E{︁

w_k^T(ε−ε∗)(ε−ε∗)^Twk

}︁

=w^T_kE{︁

(ε−ε∗)(ε−ε∗)^T}︁

w_k

=w^T_kΓ_εw_k

=w^T_kλkwk

=λ_k⟨w_k, w_k⟩

=λ_k, ja diagonaalilta poikkeavat alkiot ovat

E{α_kα_j}=E{⟨ε−ε∗, w_k⟩⟨ε−ε∗, w_j⟩}

=w^T_kΓ_εw_j

=w^T_kλ_jw_j

=λ_j⟨w_k, w_j⟩

= 0.

Kertoimienα prioritiheys on siis N(0,Γ_α), missä Γ_α =diag(λ₁, λ₂, ..., λ_p).

Käyttämällä kovarianssimatriisin Γ_ε ominaisarvohajotelmaa, approksimaatiovirheeen sisäl- tävä havaintomalli (3.2) voidaan kirjoittaa muodossa

V =UΩ,h̃ (σ, ξ₀) +ε+e

=UΩ,h̃ (σ, ξ₀) +

p

∑︂

k=1

α_kw_k+ε_∗+ε₂+e

=UΩ,h̃ (σ, ξ₀) +W_pα+ε∗+ε₂+e, (3.25) missä α= (α₁, α₂, ..., α_p)^T ja W_p = [w₁, w₂, . . . , w_p]. Vastaavasti kuin Kappaleessa 3, vektorille (σ^T, α^T)^T saadaan posterioritiheys

π_σ,α|V(σ, α|V)∝Ξ(σ) exp {︃

−1

2(V −UΩ,h̃ (σ, ξ₀)−W_pα−e_∗−ε_∗)^TΓ⁻¹_e+ε

2· (V −UΩ,h̃ (σ, ξ₀)−W_pα−e_∗−ε_∗)− 1

2(σ−σ_∗)^TΓ⁻¹_σ (σ−σ_∗) }︃

,

(3.26)