Approksimaatiovirheiden mallinnus fotoakustisessa tomograassa
Teemu Sahlström Pro gradu -tutkielma Sovelletun fysiikan koulutusohjelma Itä-Suomen yliopisto, Sovelletun fysiikan laitos 4. helmikuuta 2019
ITÄ-SUOMEN YLIOPISTO, Luonnontieteiden ja metsätieteiden tiedekunta Sovelletun fysiikan koulutusohjelma, laskennallinen fysiikka
Teemu Sahlström: Approksimaatiovirheiden mallinnus fotoakustisessa tomograassa Pro Gradu -tutkielma, 72 sivua
Tutkielman ohjaajat: Apulaisprofessori Tanja Tarvainen FT, Tutkijatohtori Aki Pulk- kinen FT ja Nuorempi tutkija Jenni Tick FM
Helmikuu 2019
Avainsanat: fotoakustinen tomograa, Bayesilaiset käänteisongelmat, approksimaatio- virhemenetelmä
Tiivistelmä
Fotoakustinen tomograa on viimeisten vuosikymmenten aikana kehitetty biolääketie- teellinen kuvantamismenetelmä, joka perustuu fotoakustiseen ilmiöön, sen mittaamiseen ja mittauksien perusteella tehtävään kuvan rekonstruktioon. Fotoakustisessa kuvantami- sessa yhdistyvät pehmytkudoksen vahva optinen kontrasti ja ultraäänen korkea resoluu- tio. Fotoakustisen kuvan rekonstruktioon on kehitetty useita erilaisia menetelmiä. Malli- pohjaisissa rekonstruktiomenetelmissä kuvan muodostaminen perustuu kuvauskohteessa etenevien ultraääniaaltojen mallintamiseen ja siihen perustuvan käänteisongelman rat- kaisuun. Kuvan laadun varmistamiseksi käänteisongelmassa käytetyn mallin tulee kuvata tilanteeseen liittyvää fysiikkaa riittävällä tarkkuudella, jotta käänteisongelman ratkaisu- na saatava kuva olisi mahdollisimman tarkka. Epätarkan mallin käyttö saattaa johtaa merkittävään virheeseen rekonstruoidussa fotoakustisessa kuvassa.
Tässä pro gradu -tutkielmassa tutkittiin mallinnusvirheiden vaikutusta ja niiden korjaamista fotoakustisessa tomograassa. Fotoakustisen tomograan käänteisongelmaa lähestyttiin Bayesilaisesta näkökulmasta. Mallinnusvirheitä tutkittiin tilanteissa, joissa käänteisongelmassa käytetty malli sisälsi virhettä johtuen joko harvennetusta aikadis- kretisaatiosta tai sensorien sijainnin epävarmuudesta. Mallinnusvirheitä korjattiin Baye- silaisen approksimaatiovirhemenetelmän avulla. Mallinnusvirheiden vaikutusta tutkit- tiin simulaatioilla. Kuvien rekonstruktiot suoritettiin tilanteissa, joissa mallinnusvirhe huomioitiin approksimaatiovirhemenetelmän avulla, sekä tilanteissa, joissa mallinnus- virhettä ei otettu huomioon. Näitä kuvia verrattiin tarkkoihin, mallinnusvirhettä sisäl- tämättömiin kuviin. Rekonstruoitujen kuvien lisäksi tarkasteltiin niiden luotettavuutta marginaalisten posteriorijakaumien ja luottamusvälien avulla.
Aikadiskretisaation harvennuksen huomattiin aiheuttavan rekonstruoituun kuvaan merkittävää virhettä. Approksimaatiovirhemenetelmän todettiin kuitenkin toimivan hy- vin aikadiskretisaation harvennuksesta aiheutuvien mallinnusvirheiden korjaamisessa erityisesti suurien sensorimäärien ja harvan aikadiskretisaation tapauksissa. Näissä ti- lanteissa laadullisesti tarkan kuvan tasoisia kuvia pystyttiin tuottamaan huomattavas- ti nopeammin verrattaessa tarkassa aikahilassa laskettuihin kuviin. Sensoripaikkojen epävarmuuksien tapauksessa mallinnusvirheen aiheuttamien virheiden havaittiin olevan riippuvia käytetystä sensorimäärästä ja geometriasta. Suurilla sensorimäärillä sensori- paikkojen epävarmuuksien aiheuttamat virheet olivat pieniä ja approksimaatiovirhe- mallin käytön todettiin heikentävän kuvan laatua. Pienillä vajaakulmageometriaan ase- tetuilla sensorimäärillä sensoripaikkojen epävarmuuksien aiheuttamat mallinnusvirheet olivat puolestaan merkittäviä, ja niistä aiheutuvia artefakteja pystyttiin korjaamaan te- hokkaasti approksimaatiovirhemenetelmän avulla.
Abstract
Photoacoustic tomography is an emerging biomedical imaging modality, which combi- nes strong optical contrast with high ultrasound resolution. A photoacoustic image is formed by measuring ultrasound waves produced by photoacoustic eect and reconstruc- ting the image. A variety of methods for the reconstruction of the photoacoustic image have been developed. In model based reconstruction techniques, the image is formed by modeling the ultrasonic waves propagating in tissue and solving the resulting inverse problem. In order to produce as good images as possible, the underlying physics has to be modeled accurately. Use of an inaccurate model may result in signicant errors in the reconstructed photoacoustic image.
In this thesis, the eect and compensation of modeling errors in photoacoustic to- mography were studied. Reconstruction of the photoacoustic image was approached in the Bayesian framework. Modeling errors were studied in the cases of reduced time disc- retization and uncertainties in sensor locations. Compensation of modeling errors was carried out by using the Bayesian approximation error approach. The eect of modeling errors were investigated by reconstructing accurate photoacoustic images, images using the inaccurate models and images using the approximation error approach for compen- sating the uncertainties. In addition, several condence intervals and marginal posterior densities were calculated to study the reliability of the reconstructed images.
The modeling errors due to reduced time discretization resulted in signicant arte- facts in the reconstructed images. These artefacts were well compensated by using the approximation error approach, especially in the case of large number of sensors and sparse time discretization. In general, the approximation error approach allowed shorter reconstruction times for producing images comparable to the accurate model recon- structions. In the case of uncertain sensor locations, the errors in the image were found to be dependent on the number of sensors and the sensor geometry. In the case where sensors were set evenly on the domain boundary, the errors due to uncertainty in the sensor locations were small. In addition, the approximation error approach was not able to compensate for these errors. The errors in the reconstructed images were however signicantly larger in the case of limited view sensor geometries. In this situation, the artefacts due to the modeling errors were well compensated using the approximation error approach.
Lyhenteet
AEM Approximation error model, Approksimaatiovirhemalli CEM Conventional error model, Normaali virhemalli
CFL -luku Courant-Friedrichs-Lewy -luku CI Credible interval, Luottamusväli
CM Conditional mean, Ehdollinen odotusarvo CT Computed tomography, Tietokonetomograa EEM Enhanced error model, Parannettu virhemalli FD Finite dierence, Äärellinen dierenssi
FE Finite element, Äärellinen elementti
FFT Fast Fourier transform, Nopea Fourier-muunnos MAP Maximum a posteriori
MRI Magnetic resonance imaging, Magneettikuvantaminen PET Positron emission tomography, Positroniemissiotomograa PML Perfectly matching layer, Akustisen aallon vaimentava kerros TR Time reversal, Ajankääntö
Symbolit
α Kromoforin molaarinen absorptiospektri C Kromoforin konsentraatio
c Äänen nopeus c0 Vakio äänen nopeus cmax Suurin äänen nopeus
d Dimensio
dl Sensorin paikka
dl0 Sensorin permutoimaton paikka δ Diracin delta funktio
∆t Aikahilan hilapisteiden väli
∆x Spatiaalihilan hilapisteiden väli
∆xe Efektiivinen resoluutio
η Odotusarvo
e Mittausvirhe E Vektorijoukko
E Odotusarvo-operaattori Ep0 Suhteellinen virhe ε Approksimaatiovirhe
fˆ Funktion f Fourier-muunnos F Fourier-muunnosoperaattori
fmax Suurin suoran mallin tukema taajuus g Siroamiskulman kosinin keskiarvo γd Dimensiosta riippuva vakio
Γ Kovarianssimatriisi Γˆ Grüneisenin parametri
H Absorboitunut optinen energiatiheys i Imaginaariyksikkö
I Valovirta
κ Valon diuusiokerroin
K Lineaarinen suora malli
K(z) Lineaarinen tarkka suora malli K(z0) Lineaarinen epätarkka suora malli k Aaltovektori
kmax Suurin mahdollinen aaltoluku l Sileysparametri
λ Aallonpituus
λref Referenssiaallonpituus µa Valon absorptiokerroin µs Valon sirontakerroin µ0s Rajoitettu sirontakerroin µs,ref Referenssi-sirontakerroin N Normaalijakauma
Nt Aikapisteiden määrä
Nx,tot Laskenta-alueen pikselien kokonaismäärä Nx,rec Rekonstruktioalueen pikselien kokonaismäärä n Kokonaisvirhe
ν Ulospäin osoittava yksikkönormaali
Ω Alue
∂Ω Alueen reuna Φ Fotonivuon tiheys p Akustinen paine
p0 Fotoakustinen alkupaine
p0,rec Alkupaineen rekonstruoitu estimaatti pt Akustinen paine-aikasarja
π(·) Todennäköisyystiheysfunktio r Paikkavektori
ρ Väliaineen akustinen tiheys ρ0 Väliaineen tiheys lepotilassa σ Keskihajonta
σ2 Varianssi
t Aika
u Akustisen partikkelin nopeus U Tasajakauma
Upix Tasajakauma pikseleinä Uspat Tasajakauma metreinä
U Sensorien permutaatiossa käytetty tasajakauma z Mallin tarkkuuteen vaikuttava tarkka parametri z0 Mallin tarkkuuteen vaikuttava epätarkka parametri
Sisältö
1 Johdanto 1
2 Fotoakustinen kuvantaminen 3
2.1 Fotoakustinen suora ongelma . . . 3
2.2 Fotoakustinen käänteisongelma . . . 7
3 Fotoakustisen tomograan suoran ongelman mallinnus 10 3.1 Pseudospektraali k-avaruusmenetelmä . . . 10
4 Bayesilainen käänteisongelma ja estimointi fotoakustisessa tomograas- sa 16 4.1 Uskottavuusjakauma . . . 17
4.2 Etukäteistieto ja priorijakauma . . . 18
4.3 Posteriorijakauma . . . 19
4.4 Bayesilainen approksimaatiovirhemenetelmä . . . 20
4.5 Piste-estimaateista . . . 24
5 Simulaatiot 26 5.1 Mittausdatan simulointi ja väliaine . . . 26
5.2 Käänteisongelman ratkaisu . . . 26
5.2.1 Spatiaaliset diskretisaatiot . . . 28
5.2.2 Suoran mallin muodostaminen . . . 28
5.2.3 Priorijakauma ja kohinamalli . . . 28
5.3 Aikadiskretisaation harvennus . . . 29
5.3.1 Harvennettujen aikadiskretisaatioiden muodostaminen . . . 29
5.3.2 Approksimaatiovirheen mallinnus . . . 29
5.4 Sensoripaikkojen epävarmuudet . . . 30
5.4.1 Mittausdatan simulointi . . . 30
5.4.2 Käänteisongelman ratkaisu ja approksimaatiovirheen mallinnus . . . 31
6 Tulokset 32 6.1 Aikadiskretisaation harvennus . . . 32
6.2 Sensoripaikkojen epävarmuudet . . . 45
7 Pohdinta 52 7.1 Aikadiskretisaation harvennus . . . 52
7.2 Sensoripaikkojen epävarmuudet . . . 53
8 Yhteenveto 55
Liitteet 56
Viitteet 64
1 Johdanto
Viime vuosisadan aikana kehitetyt lääketieteelliset kuvantamismenetelmät ja niiden so- vellukset ovat arvokkaita työkaluja sekä kliinisessä- että tutkimuskäytössä. Yleisimmät kuvausmenetelmät, kuten magneettikuvantaminen (engl. magnetic resonance imaging, MRI), tietokonetomograa (engl. computed tomography, CT) ja positroniemissiotomogra- a (engl. positron emission tomography PET), ovat mahdollistaneet biologisten kudosten tutkimuksen ennennäkemättömällä tarkkuudella ja nopeudella. Nykypäivänä käytetyt ku- vantamismenetelmät kärsivät kuitenkin kullekin menetelmälle tyypillisistä heikkouksista, kuten metallisten implanttien aiheuttamista rajoitteista (MRI) [15] ja ionisoivan säteilyn aiheuttamista säteilyrasituksesta (CT ja PET) [6, 10] tai korkeista käyttökustannuksista.
Näin ollen uudenlaisille lääketieteellisille ja biolääketieteellisille kuvantamismenetelmille on edelleen tarvetta.
Fotoakustinen kuvantaminen on viimeisten vuosikymmenten aikana kehitetty kuvanta- mismentelmä, joka perustuu fotoakustiseen ilmiöön. Fotoakustisessa ilmiössä kuvauskoh- teeseen kohdistettu valo aiheuttaa kuvauskohteeseen lokaalin lämpölaajenemisen ja sitä seuraavan paineen nousun, joka etenee kuvauskohteessa akustisina ultraääniaaltoina. Fo- toakustinen kuva muodostetaan mittaamalla kuvaukohteessa etenevät ultraääniaallot koh- teen pinnalta ja ratkaisemalla tilanteeseen liittyvä matemaattinen käänteisongelma. Fotoa- kustisessa kuvantamisessa yhdistyvät pehmytkudoksen korkea optinen kontrasti sekä ult- raäänen tarkka resoluutio. Fotoakustisen kuvantamisen sovelluskohteina nähdään erityi- sesti iho- ja rintasyövän kuvantaminen, erilaisten hoitojen seuranta, verisuoniston kuvan- taminen ja pieneläinkuvantaminen biolääketieteellisessä tutkimuksessa. [4, 8, 29, 30, 51]
Lääketieteellisissä kuvantamismenetelmissä ja niiden sovelluksissa kuvan muodostami- nen perustuu usein matemaattiseen käänteisongelmaan, jossa kuvan rekonstruktio laske- taan mallintamalla tilanteeseen liittyvää fysikaalista ilmiötä. Tällaista kuvan rekonstruk- tiomenetelmää kutsutaan mallipohjaiseksi rekonstruktiomenetelmäksi. Mallipohjaiset re- konstruktiomenetelmät ovat erittäin joustavia ja niitä voidaan periaatteessa soveltaa riip- pumatta kuvantamislaitteiston geometriasta [11, 28, 49]. Kuvan laadun ja tarkkuuden takaamiseksi rekonstruktiossa käytetyn mallin tulee kuitenkin mallintaa tilanteeseen fysi- kaalista ilmiötä riittävän tarkasti. Tarkan mallin muodostaminen ja käyttö voi kuitenkin olla joissain tilanteissa laskennallisesti liian raskasta, jolloin käänteisongelman ratkaise- miseksi joudutaan käyttämään niin kutsuttua redusoitua epätarkkaa mallia. Käänteison- gelmassa käytetty redusoitu malli voi tilanteesta riippuen olla esimerkiksi approksimaa- tio tarkasta mallista, tai malli, joka on muodostettu käyttäen harvaa laskentahilan dis- kretisaatiota. Epätarkan mallin käyttö saattaa kuitenkin johtaa merkittäviin virheisiin käänteisongelman ratkaisussa erityisesti huonokuntoisten käänteisongelmien tapauksessa [17, 18].
Tässä tutkielmassa tutkitaan fotoakustisessa kuvantamisessa esiintyviä mallinnusvir- heitä ja niiden korjaamista. Fotoakustisen tomograan käänteisongelmaa lähestytään Baye- silaisesta näkökulmasta, jossa tietoa kohteen tuntemattomista parametreista pyritään määrittämään kuvauskohteesta tehtyjen mittauksien, tilanteeseen liittyvän fysiikan mal- lintamisen ja kuvauskohteesta tunnetun ennakko- eli prioritiedon avulla [17, 31, 41, 43].
Työssä mallinnetaan ultraäänen etenemistä lineaarisella aaltoyhtälöllä, joka ratkaistaan pseudospektraalilla k-avaruusmenetelmällä käyttäen k-Wave toolbox:ia [9, 27, 39, 47].
Mallinnusvirheiden vaikutusta tutkitaan tilanteissa, joissa käänteisongelman ratkaisus- sa käytetty malli sisältää virhettä joko aikadiskretisaation harventamisen tai sensorien paikkojen epävarmuuksien takia. Mallinnusvirheiden aiheuttamaa virhettä pyritään kor- jaamaan Bayesilaisen approksimaatiovirhemenetelmän avulla, jossa mallinnusvirhettä ap- proksimoidaan tilannetta kuvaavien tarkkojen ja epätarkkojen mallien avulla [17, 18, 19, 20]. Mallinnusvirheen vaikutusta rekonstruoituihin fotoakustisiin kuviin tutkitaan muo- dostamalla fotoakustiset kuvat mallinnusvirhettä sisältävän mallin avulla sekä mallin- nusvirheen huomioon ottavan approksimaatiovirhemallin avulla. Näitä kuvia verrataan tarkkaa mallia käyttäen muodostettuihin kuviin.
Tässä tutkielmassa käsitelty fotoakustinen tomograa ja siihen liittyvä suora- ja kään- teisongelma on esitetty kappalessa 2. Ultraäänikentän mallintamiseen käytetty pseudos- pektraali k-avaruus -menetelmä on esitetty kappaleessa 3. Käänteisongelman ratkaisussa käytetty Bayesilainen lähestymistapa ja mallinnusvirheiden korjaamisessa käytetty ap- proksimaatiovirhemenetelmän teoria on esitetty kappaleessa 4. Tutkielmassa suoritetuis- sa simulaatioissa ja rekonstruktioissa käytetyt parametrit on esitetty kappaleessa 5. Tut- kielman tulokset on esitetty kappaleessa 6. Tuloksien pohdinta ja yhteenveto on esitetty kappaleissa 7 ja 8.
2 Fotoakustinen kuvantaminen
Fotoakustinen kuvantaminen on hybridikuvantamismenetelmä, jossa kuvauskohteen omi- naisuuksia tutkitaan valolla tuotetun ultraäänen avulla. Fotoakustisen kuvantamisen toi- mintaperiaate on esitetty kuvassa 1. Menetelmä perustuu fotoakustiseen ilmiöön, jonka toimintaperiaate voidaan jakaa kolmeen osaan: kuvantamisessa käytettyjen fotonien ab- sorptioon, tästä aiheutuvaan lokaaliin lämpötilan ja paineen kasvuun sekä kuvauskohtees- sa eteneviin akustisiin ultraääniaaltoihin. [4, 8, 29, 30, 51]
Fotoakustinen kuvantamisprosessi aloitetaan kuvauskohteeseen kohdistetulla lyhyellä (tyypillisesti nanosekuntien pituisella) valopulssilla. Valopulssin fotonit etenevät kuvaus- kohteessa ja absorboituvat pääosin valoa tehokkaasti absorboiviin kromoforeihin. Tyypil- lisimpiä pehmytkudoksessa esiintyviä kromoforeja ovat hemoglobiini, melaniini, rasvat ja vesi, jotka toimivat pääasiallisena fotoakustisen kuvan kontrastin lähteenä (kuva 2). Täs- tä johtuen fotoakustisessa kuvantamisessa pyritään käyttämään valoa, jolla kromoforien valon aallonpituudesta riippuva absorptio on suurimmillaan ja eri kromoforien absorp- tiospektri voidaan erottaa. Yleisimpiä fotoakustisessa kuvantamisessa käytettyjä valon aallonpituuksia ovat näkyvän punaisen tai lähi-infrapunan aallonpituudet. [4, 29]
Kuvauskohteeseen kohdistetun valon absorptio aiheuttaa kuvauskohteeseen lokaalin lämpötilan kasvun, jonka seurauksena tapahtuu lämpölaajenemista ja mekaanista rasitus- ta erityisesti korkeiden kromoforikonsentraatioiden alueilla. Lämpölaajenemisen aiheutta- man tilavuuden muutoksen seurauksena kuvauskohteeseen syntyy paine-eroja, jotka ete- nevät kohteessa akustisina ultraääniaaltoina. Kohteen kromoforikonsentraatioista tietoa sisältävät ultraääniaallot mitataan kuvauskohteen pinnalta. Lopulta fotoakustinen kuva rekonstruoidaan ultraäänianturien mittaamien paine-aikasarjojen perusteella. [4, 29]
Tyypillinen fotoakustinen mittauslaitteisto koostuu valonlähteestä ja yhdestä tai use- ammasta ultraäänisensorista, joista yleisimpiä ovat pietsosähköiset sensorit. Tomogra- sista mittauslaitteistoista yleisimpiä ovat geometriat, joissa ultraäänisensorit on asetettu pallosymmetrisesti, sylinterisymmetrisesti tai tasogeometriaan kuvauskohteen ympärille.
Valonlähteinä on tyypillisesti käytetty laser-pohjaisia ratkaisuja niiden tehokkaan valaisu- kyvyn takia. Fotoakustisia mittauslaitteistoja on kuitenkin kehitetty myös LED-valaisuun perustuen erityisesti niiden alhaisen hinnan ja helppokäyttöisyyden vuoksi. [4, 26, 30]
2.1 Fotoakustinen suora ongelma
Fotoakustisen kuvantamisen (kuva 1) mallinnusongelmaa voidaan tarkastella kaksiosai- sena. Toisin sanoen mallinnusongelmassa esiintyvät optinen ja akustinen ilmiö voidaan erottaa, koska ultraäänen eteneminen on pehmytkudoksessa valon absorptioon verrattu- na huomattavasti hitaampaa. Fotoakustisen mallinnusongelman ensimmäinen osa on niin kutsuttu optinen ongelma, jossa mallinnetaan kuvauskohteeseen absorboituneen optisen energiatiheyden jakaumaH kohteeseen kohdistetun valon perusteella. Mallinnusongelman toisessa, eli akustisessa, osassa ratkaistaan kuvauskohteesta mitatut fotoakustiset paine- aikasarjat, kun absorboituneen optisen energian aiheuttama alkupaine p0 on tunnettu [8, 36, 41]. Fotoakustisen kuvantamisen suoran ongelman periaate on esitetty kuvassa 3.
Kuva 1: Fotoakustisen kuvantamisen toimintaperiaate. Kuvantamisprosessi aloitetaan ku- vauskohteeseen kohdistetulla lyhyellä valopulssilla (ensimmäinen kuva). Valopulssin foto- nit absorboituvat kudokseen aiheuttaen lokaalin lämpötilan ja paineen nousun (toinen kuva). Lämpötilan nousua seuraavat ultraääniaallot mitataan kuvauskohteen pinnalta (kolmas kuva). Lopulta fotoakustinen kuva voidaan rekonstruoida mitattujen ultraää- niaaltojen tuottaman informaation perusteella (neljäs kuva). (Kuva lainattu muokattuna lähteestä: L. Wang ja J. Yao, A practical guide to photoacoustic tomography in the life sciences, Nature Methods, 13(8):627-638, 2016.)
Kuva 2: Yleisimpien pehmytkudoksessa esiintyvien kromoforien absorptiokertoimia va- lon aallonpituuden funktiona. Hapettunut hemoglobiini (punainen), hapeton hemoglobii- ni (sininen), vesi (musta), rasva (80 %:a kudoksen tilavuudesta ruskealla ja 100 %:a vaa- leanpunaisella), melaniini (musta katkoviiva), kollageeni (vihreä) ja elastiini (keltainen).
(Kuva lainattu lähteestä: P. Beard, Biomedical photoacoustic imaging, Interface Focus, 1:602-631, 2011.)
Kuva 3: Fotoakustisen kuvantamisen suora ongelma. Optinen ongelma (ylempi laatikko) on ratkaista kuvauskohteeseen absorboituneen optisen energiantiheyden jakauma H ja sen perusteella alkupaine p0. Akustinen ongelma (alempi laatikko) on ratkaista kohteen pinnalta mitattava fotoakustinen paine-aikasarja pt, kun alkupaine on annettu. (Kuva lainattu lähteestä: B. Cox, J. Laufer, S. Arridge ja P. Beard. Quantitative spectroscopic photoacoustic imaging: a review. Journal of Biomedical Optics, 17(6):061202, 2012.)
Valon etenemisen mallintamiseksi d-dimensioisessa (2 tai 3 ulotteisessa) alueessaΩja sen reunalla ∂Ω voidaan käyttää esimerkiksi diuusioapproksimaatiota [1, 41]
−∇ ·κ(r, λ)∇Φ(r, λ) +µa(r, λ)Φ(r, λ) = 0, r∈Ω, Φ(r, λ) + 1
2γdκ(r, λ)∂Φ(r, λ)
∂ν =
Is
γd r ∈j 0 r ∈∂Ω\j,
(1)
missä r on paikkavektori, λ on aallonpituus, j on valonlähteen paikka, Is on valovir- ta lähteestä j, γd on dimensiosta riippuva vakio (γ2 = 1/π, γ3 = 1/4), ν on alueesta ulospäin osoittava normaali, Φ(r, λ) on fotonivuon tiheys (engl. uence), µa(r, λ) on va- lon absorptiokerroin. Lisäksi yhtälössä (1) esiintyvä diuusiokerroin κ voidaan esittää muodossa κ(r, λ) = (d(µa(r, λ) +µ0s(r, λ)))−1 sekä sirontakerroin µ0s muodossa µ0s(r, λ) = (1−g)µs(r, λ), missä µs(r, λ) on sirontakerroin jag on siroamiskulman kosinin keskiarvo.
Valon etenemistä kuvaava fotonivuon tiheys riippuu näin ollen epälineaarisesti kuvaus- kohteen optisista absorptio- ja sirontakertoimista µa ja µs.
Absorboitunut optinen energiatiheys H voidaan kirjoittaa diuusioapproksimaation (1) ratkaisuna saatavan fotonivuon tiheydenΦja absorptiokertoimenµaavulla muodossa [4]
H(r;µa(r, λ), µs(r, λ)) =µa(r, λ)Φ(r, λ). (2) Kuvauskohteen absorptiokerroin riippuu siinä olevista kromoforeista ja niiden konsentraa- tioista. Valon absorptiokerroin voidaan esittää kuvauskohteessa olevien K:n kromoforin valon aallonpituudesta riippuvien absorptiokertoimien summana muodossa [4]
µa(r, λ) =
K
X
k=1
Ck(r)αk(λ), (3)
missä Ck(r) on k:nnen kromoforin konsentraatio ja αk on k:nnen kromoforin molaari- nen absorptiokerroinspektri. Tämän lisäksi kuvauskohteen sironnan spektri voidaan Mie sirontateorian perusteella esittää muodossa [37]
µ0s(r, µs,ref, b, λ) = µs,ref(r) λ
λref −b(r)
, (4)
jonka mukaan rajoitettu sirontakerroin µ0s on verrannollinen rajoitettuun referenssisiron- takertoimeenµ0s,ref aallonpituudellaλref ja suhteelliseen aallonpituuteenλ/λref korotettuna potenssiin −b(r).
Kun absorboituneen energiantiheyden jakaumaH on tiedossa, voidaan lämpölaajene- misen seurauksena aiheutuva hetkellinen paineen nousu, eli fotoakustinen alkupaine p0
ajan hetkellä t= 0, kirjoittaa muodossa [4, 8]
p0(r, λ) = ˆΓ(r)H(r, λ), (5) missä Γˆ on Grüneisenin parametri, joka kuvaa hyötysuhdetta absorboituneen energian ja sitä seuraavan paineen nousun välillä. Alkupaine riippuu näin ollen epäsuorasti kaikista optisen energian absorboitumiseen vaikuttavista parametreista. Alkupaineen eteneminen
vaimentamattomassa ja homogeenisessa väliaineessa voidaan mallintaa käyttäen lineaa- rista aaltoyhtälöä ja alkuehtoja [35]
∇2p(r, t)− 1 c2
∂2p(r, t)
∂t2 = 0 p(r, t= 0) =p0(r)
∂
∂tp(r, t= 0) = 0,
(6)
missä con äänen nopeus, t on aika ja pt on akustinen paine ajan hetkellä t. Käytännös- sä fotoakustisen kuvan rekonstruktioissa käytetty paine-aikasarja pt mitataan äärellisissä määriä pisteitä tai pintoja kuvauskohteen ympäriltä, jolloin mittausdatasta on käytettä- vissä paineen p rajoitettu osa.
2.2 Fotoakustinen käänteisongelma
Fotoakustiseen kuvantamiseen liittyvää käänteisongelmaa voidaan suoran ongelman ta- paan tarkastella kaksiosaisena. Fotoakustisen käänteisongelman rakenne on esitetty ku- vassa 4. Käänteisongelman ensimmäisessä, eli akustisessa, osassa tavoitteena on ratkaista kuvauskohteen sisäinen alkupainejakauma p0 kohteesta mitattujen fotoakustisten paine- aikasarjojen pt avulla. Akustisen käänteisongelman ratkaisuna saatava kuva on luonteel- taan kvalitatiivinen, eikä anna tietoa kohteen fysiologisista parametreista. Fotoakustisen käänteisongelman akustista osaa on perinteisesti kutsuttu fotoakustiseksi tomograaksi.
[4, 42, 43]
Käänteisongelman toisessa eli optisessa osassa tavoitteena on ratkaista kuvauskohteen sisäiset kvantitatiivista tietoa antavat optiset sironta- ja absorptiokertoimet µs ja µa tai kromoforien konsentraatiot Ck, kun absorboitunut optinen energiatiheysH ja kohteeseen syötetyn valon määrä on tunnettu. Jos mittaukset tehdään yhdellä aallonpituudella, saa- daan tulokseksi kromoforien absorptio yhdellä aallonpituudella. Toisaalta, jos mittaukset suoritetaan usealla aallonpituudella, voidaan kunkin kromoforin konsentraatio ratkaista käyttämällä tiedossa olevia kromoforien spektrejä. Fotoakustisen käänteisongelman optis- ta osaa kutsutaan kvantitatiiviseksi fotoakustiseksi tomograaksi. [36, 41]
Fotoakustinen käänteisongelma voidaan vaihtoehtoisesti ratkaista yhdessä osassa, jol- loin akustista ja optista käänteisongelmaa ei eroteta. Tässä tilanteessa kuvauskohteen kromoforien konsentraatiot Ck ratkaistaan suoraan fotoakustisten paine-aikasarjojen pt avulla. [36]
Tässä työssä tarkastellaan fotoakustisen tomograan käänteisongelman akustista osaa.
Fotoakustisessa tomograassa alkupaineenp0 rekonstruktioon on kehitetty useita erilaisia menetelmiä. Yksinkertaisimmillaan rekonstruktio voidaan suorittaa niin kutsutuilla ta- kaisinprojektio -algoritmeilla (engl. back projection algorithm) ja parannetuilla suodate- tuilla takaisinprojektio-algoritmeilla, joilla voidaan saavuttaa eksakti rekonstruktio pallo-, sylinteri- ja tasogeometrioihin asetetuille sensorigeometrioille [24, 52]. Takaisinprojektio- algoritmien lisäksi eksakti rekonstruktio voidaan saavuttaa erilaisilla sarjapohjaisilla me- netelmillä, jotka perustuvat mitattujen fotoakustisten signaalien temporaalisiin ja spati- aalisiin hajotelmiin [25, 33]. Rekonstruktio voidaan suorittaa myös niin kutsutuilla ajan- kääntö (engl. time-reversal, TR) -menetelmillä, joissa mitatut fotoakustiset signaalit las- ketaan temporaalisesti takaperin erilaisten akustisten mallien avulla [16, 53].
Kuva 4: Fotoakustisen kuvantamisen käänteisongelma. Käänteisongelman akustisessa osassa (alempi laatikko) tavoitteena on rekonstruoida kuvauskohteessa oleva alkupaine- jakauma p0 kuvauskohteesta mitattujen paine-aikasarjojenpt avulla. Rekonstruoitu alku- painejakauma voidaan edelleen muuttaa optisen käänteisongelman alkudataksi, eli absor- boituneen optisen energiatiheydenH jakaumaksi Grüneisenin parametrin avulla. Optinen käänteisongelman ratkaisuna saadaan kuvauskohteen optisten parametrien µa ja µs sekä kromoforikonsentraatioiden Ck jakaumat. Vaihtoehtoisena lähestymistapana käänteison- gelma voidaan ratkaista yhdessä osassa, jossa akustista ja optista käänteisongelmaa ei eroteta. (Kuva lainattu muokattuna lähteestä: B. Cox, J. Laufer, S. Arridge ja P. Beard.
Quantitative spectroscopic photoacoustic imaging: a review. Journal of Biomedical Optics, 17(6):061202, 2012.)
Fotoakustisen kuvan rekonstruktiota voidaan lähestyä myös niin kutsutuilla mallipoh- jaisilla menetelmillä, joissa kuvan muodostus perustuu kuvantamisessa käytetyn fysikaa- lisen ilmiön mallintamiseen [2, 34, 50, 54]. Olkoon tilanteeseen liittyvä alkupaineen p0 ja mittausdatan pt yhdistävä lineaarinen havaintomalli muotoa
pt=Kp0, (7)
missäK on matriisioperaattori. Mallipohjaisissa käänteisongelmien rekonstruktiomenetel- missä tavoitteena on määrittää kohteen estimoitu alkupainejakauma p0,rec minimoimalla havaintojen pt ja mallin avulla lasketun mittausdatan residuaali
p0,rec = arg min
p0
kpt−Kp0k22. (8) Huonokuntoisten käänteisongelmien tapauksessa yhtälön (8) ratkaiseminen sellaise- naan ei yleensä ole mielekästä. Käänteisongelmien huonokuntoisuutta voidaan pyrkiä hel- pottamaan esimerkiksi Tikhonov regularisaatioon perustuvilla regularisaatiomenetelmil- lä [54]. Huonokuntoisten käänteisongelmien ratkaisua voidaan lähestyä yleisimpien regu- larisaatiomenetelmien lisäksi myös tilastollisesta näkökulmasta, jossa käänteisongelmas- sa esiintyvät parametrit mallinnetaan satunnaismuuttujina. Tässä työssä fotoakustinen käänteisongelmaa lähestytään tilastollisiin käänteisongelmien ratkaisumenetelmiin kuulu- valla Bayesilaisella lähestymistavalla [42, 43]. Fotoakustisen käänteisongelman Bayesilai- nen ratkaisu ja siihen liittyvä teoria on esitetty kappaleessa 4.
3 Fotoakustisen tomograan suoran ongelman mallin- nus
Käänteisongelmien mallipohjaisissa rekonstruktiomenetelmissä käytettyjen mallien muo- dostaminen perustuu tilannetta kuvaavan fysikaalisen ilmiön mallintamiseen. Fotoakusti- sessa tomograassa ultraäänen etenemistä kuvaavan aaltoyhtälön ratkaisemiseen on käy- tettävissä useita eri menetelmiä, joista yleisimpiä ovat dierenssi- (engl. nite dierence, FD) ja äärellisen elementtien menetelmät (engl. nite element, FE). [14]
FD- ja FE-menetelmät saattavat yleisyydestään huolimatta olla laskennallisesti ras- kaita. Yksinkertaisimmillaan aaltoyhtälön gradienttia approksimoidaan FD-menetelmissä lineaarisella interpolaatiolla vierekkäisten pisteiden välillä. Tarkan gradientin arvon saa- vuttamiseksi tarvitaan kuitenkin useampia hilapisteitä, korkeamman kertaluvun dierens- simenetelmiä ja riittävän tiheä laskentahila (kuva 5). Yleisesti FD- ja FE-menetelmissä aaltoyhtälön gradientin tarkan ratkaisun saavuttamiseksi tarvitaan useita laskentahilan pisteitä jokaista simuloitavaa aallonpituutta kohden. Tässä tilanteessa menetelmissä vaa- dittavat laskentaresurssit saattavat kasvaa käytännön kannalta liian suuriksi erityisesti korkeataajuisten akustisten kenttien tapauksissa. [9, 48]
Dierenssimenetelmille ominaista tiheää laskentahilaa voidaan kuitenkin harventaa käyttämällä niin kutsuttua pseudospektraalia menetelmää, jossa gradientin approksimaa- tio muodostetaan Fourier-muunnoksen avulla (kuva 5). Dierenssimenetelmiin verrattuna pseudospektraali menetelmä tarjoaa kaksi merkittävää etua: nopean Fourier-muunnoksen (FFT) laskennallisen tehokkuuden ja sen kantafunktioiden sinusoidaalisuuden, joka mah- dollistaa teoriassa vain kahden hilapisteen käytön yhtä aallonpituutta kohden. [9, 12, 39]
3.1 Pseudospektraali k-avaruusmenetelmä
Yksinkertainen pseudospektraali aaltoyhtälön ratkaisu voidaan muodostaa Fourier-muun- noksen avulla. Jos akustinen väliaine oletetaan homogeeniseksi ja valon aiheuttaman läm- penemisen oletetaan tapahtuvan välittömästi, voidaan akustinen painekenttä mallintaa käyttämällä lineaarista aaltoyhtälöä (6). Aaltoyhtälö voidaan esittää myös kahden ensim- mäisen asteen osittaisdierentiaaliyhtälön ryhmänä [27, 35, 39]
∂u(r, t)
∂t =− 1
ρ0∇p(r, t)
∂ρ(r, t)
∂t =−ρ0∇ ·u(r, t),
(9)
missä p=c20ρ on akustinen paine, c0 on vakio äänen nopeus, ρ on akustinen tiheys, u on akustinen partikkelinopeus ja ρ0 on väliaineen tiheys lepotilassa. Käyttämällä spatiaalis- ta Fourier-muunnosta, voidaan homogeeninen aaltoyhtälö (6) kirjoittaa paikkataajuusa- varuudessa k= (kx, ky, kz)muodossa
Kuva 5: Spatiaalisen gradientin approksimaatio dierenssi- ja pseudospektraalissa me- netelmissä. Kuvassa (a) ensimmäisen asteen etenevä dierenssimenetelmä. Kuvassa (b) neljännen asteen keskitetty dierenssimenetelmä ja kuvassa (c) pseudospektraali mene- telmä. (Kuva lainattu lähteestä: B. Treeby, B. Cox ja J. Jaros, k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave elds - User Manual, 2016.)
∂2p(k, t)ˆ
∂t2 =−c20k2p(k, t),ˆ (10) missäc0on vakio äänen nopeus,pˆon funktionpFourier-muunnos jak2 =k·k=kx2+k2y+kz2 ja k on aaltovektori. 1
Yhtälössä (10) esiintyvä toisen asteen aikaderivaattaa voidaan edelleen approksimoi- da käyttämällä toisen asteen keskitettyä dierenssimenetelmää, jolloin yhtälö saadaan muotoon [27, 39]
ˆ
p(k, t+ ∆t)−2ˆp(k, t) + ˆp(k, t−∆t)
(∆t)2 =−c20k2p(k, t),ˆ (11) missä ∆t on tasavälisen aikahilan pisteiden välinen aika ja c = c0 on vakio äänen no- peus. Yhtälön (11) vasemmalla puolella esiintyvä ajan suhteen muodostettu dierenssi- lauseke kärsii edelleen dierenssimenetelmälle ominaisista ongelmista, kuten liian harvan
1Yhtälössä (10) esitettyyn muotoon päästään käyttämällä rajoitetun funktion derivaatan Fourier- muunnoksen ominaisuutta
F
d
dxf(x)
(k) = ikfˆ(k) missäion imaginääriyksikkö ja kon taajuus.
laskentahilan aiheuttamasta ratkaisun dispersoitumisesta. Dierenssitermiä voidaan kui- tenkin muokata niin kutsutun k-avaruusmenetelmän avulla. Tarkastelemalla spatiaalisen- Fourier avaruuden määrittämiä tasoaaltoja homogeenisessä väliaineessa voidaan yhtälön (11) dierenssiapproksimaatio muokata eksaktiksi suhteessa yhtälön oikeaan puoleen sinc- muotoisella korjauskertoimella [27, 39]
ˆ
p(k, t+ ∆t)−2ˆp(k, t) + ˆp(k, t−∆t)
(∆t)2sinc(c0∆tkkk/2)2 =−c20k2p(k, t).ˆ (12) Aaltoyhtälön ratkaisua yhtälön (12) avulla kutsutaan pseudospektraaliksi k-avaruusme- netelmäksi.
Verrattaessa yhtälössä (11) esiintyvää aikadiskretisaatiota (∆t)2 yhtälössä (12) esiin- tyvään termiin (∆t)2sinc(c0∆tkkk/2)2, huomataan, että pienillä aikaväleillä termit ovat likipitäen samansuuruisia. Suuremmilla aikaväleillä k-avaruus operaattori toimii numee- rista dispersiota korjaavana terminä, joka mahdollistaa harvempien aikahilojen käytön.
Yhtälön (12) mukaisen aikadiskretisaation on osoitettu olevan tarkka homogeenisissa vä- liaineissa. Epähomogeenisen ja vaimentavan väliaineen tapauksessa diskretisaatio ei kui- tenkaan ole enää tarkka, mutta vähentää numeerista dispersiota huomattavasti ja toimii hyvänä approksimaationa. [27, 39]
Pseudospektraalin k-avaruusmenetelmän paikka-avaruuden vastine saadaan järjestele- mällä termejä ja käyttämällä käänteistä Fourier-muunnosta [39, 48]. Yhtälö (12) voidaan tällöin kirjoittaa muodossa
p(r, t+ ∆t)−2p(r, t) +p(r, t−∆t)
(∆t)2 =−c20F−1
k2sinc(c0∆tkkk/2)2F {p(r, t)} , (13) missä F ja F−1 ovat Fourier-muunnos ja käänteinen Fourier-muunnos. Yhtälössä (13) esiintyvää spatiaalista gradienttia kuvaavaa termiä
− F−1
k2sinc(c0∆tkkk/2)2F {p(r, t)} ≡[∇(c0∆t)]2p(r, t) (14) kutsutaan toisen asteen k-avaruusoperaattoriksi [39]. Operaattorissa yläindeksi (c0∆t) viittaa gradientin korvaavan termin riippuvuuteen äänen nopeudesta ja aikahilan pituu- desta. Toisen asteen k-avaruusoperaattori (14) voidaan edelleen jakaa tekijöihin paikka- dimensioiden (esimerkkinä r = (x, y) ∈ R2) suhteen käyttämällä ensimmäisen asteen k-avaruusoperaattoreita muodossa [39]
∂p(r, t)
∂(c∆t)+x =F−1
ikxeikx∆x/2sinc(c0∆tkkk/2)F {p(r, t)}
∂p(r, t)
∂(c∆t)+y =F−1
ikyeiky∆y/2sinc(c0∆tkkk/2)F {p(r, t)}
∂p(r, t)
∂(c∆t)−x =F−1
ikxe−ikx∆x/2sinc(c0∆tkkk/2)F {p(r, t)}
∂p(r, t)
∂(c∆t)−y =F−1
ikye−iky∆y/2sinc(c0∆tkkk/2)F {p(r, t)} ,
(15)
jolloin
∂
∂(c∆t)+x
∂
∂(c∆t)−x+ ∂
∂(c0∆t)+y
∂
∂(c∆t)−y
p(r, t) = [∇(c∆t)]2p(r, t). (16) Sijoittamalla yhtälössä (15) esiintyvät ensimmäisen asteen derivaattaoperaattorit osit- taisdierentiaaliyhtälöryhmään (9), saadaan kaksiulotteiselle pseudospektraalille k-ava- ruusmenetelmälle muoto [27, 39]
ux(r1, t+)−ux(r1, t−)
∆t =− 1
ρ0(r1)
∂p(r, t)
∂(c∆t)+x uy(r2, t+)−uy(r2, t−)
∆t =− 1
ρ0(r2)
∂p(r, t)
∂(c∆t)+y p(r, t+ ∆t)−p(r, t)
∆t =−ρ0(r)c(r)2
∂ux(r1, t+)
∂(c0∆t)−x +∂uy(r2, t+)
∂(c0∆t)−y
,
(17)
missä r1 = (x+ ∆x, y), r2 = (x, y + ∆y), t+ = t + ∆t/2 ja t− = t− ∆t/2. Saadus- ta yhtälöryhmästä huomataan, että akustisen aallon nopeus ratkaistaan niin kutsutussa porrastetussa hilassa. 2
Pseudospektraalin k-avaruusmenetelmän diskreetti muoto voidaan nyt kirjoittaa toi- sen asteenk-avaruusoperaattorin (14) ja yhtälössä (17) esitettyjen yhtälöiden avulla muo- dossa [45]
∂
∂ξpt=F−1
ikξeikξ∆ξ/2sinc(c0∆tkξ/2)2F pt
utξ+ =utξ−− ∆t ρ0
∂
∂ξpt
∂
∂ξutξ+ =F−1
ikξe−ikξ∆ξ/2sinc(c0∆tk/2)2F utξ+
ρt+∆tξ =ρtξ−∆tρ0 ∂
∂ξutξ+,
(18)
missä ξ = x (1D), ξ = x, y (2D) tai ξ = x, y, z (3D), ρ = Pd
i=1ρεi (d on dimensio), ρ=p/c20, ∆ξ on spatiaalihilan pikselien leveys suuntaanξ ja kξ on aaltoluku suuntaan ξ. Yhtälössä (18) esiintyvät diskreetit aaltoluvut on määritelty muodossa [45]
kξ =
−Nξ 2 ,−Nξ
2 + 1, . . . ,Nξ 2 −1
2π
∆ξNξ jos Nξ on pariton
−Nξ−1
2 ,−Nξ−1
2 + 1, . . . ,Nξ−1 2 −1
2π
∆ξNξ jos Nξ on parillinen (19)
missä Nξ on spatiaalihilan pikselien määrä suuntaan ξ.
Vaikka akustinen paine on skalaarisuure, on se jaettu spatiaalisiin komponentteihin, jotta laskenta-alueen reunalle asetettava akustista absorptiota mallintava kerros (engl.
2Porrastettu hila seuraa yhtälöissä (15) esiintyvistä eksponenttitermeistä ja Fourier-muunnoksen siirto-ominaisuudesta
F {f(x±x0)}(k) =e±i2πx0kfˆ(k).
Kuva 6: Pseudospektraalissa k-avaruusmenetelmässä käytetty porrastettu laskentahila ja yhtälössä (18) käytetyn iteraation kulku kahdessa ulottuvuudessa. Iteraation kulku alkaa varsinaisen hilapisteen (mustat pallot) akustisen paineen arvosta, jonka avulla lasketaan akustisen paineen paikkaderivaatat sekä akustisen aallon nopeudet porrastetuissa hila- pisteissä (kolmiot ja ristit). (Kuva on lainattu lähteestä: B. Treeby, B. Cox ja J. Jaros, k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave elds - User Manual, 2016.)
perfectly matched layer, PML) voidaan toteuttaa. PML:n käyttö mahdollistaa laskenta- alueen reunoilta aiheutuvien heijastuksien eliminoinnin. [5, 39] Yhtälön (18) kulku por- rastetussa laskentahilassa on esitetty kuvassa 6.
Pseudospektraalik-avaruusmenetelmä tarjoaa laskennallisesti tehokkaan ja tarkan ta- van aaltoyhtälön ratkaisemiseksi. Numeeristen menetelmien toimivuuden kannalta on kui- tenkin olennaista tarkastella menetelmän stabiiliutta ja siihen liittyviä ehtoja. Menetel- män stabiiliutta voidaan kontrolloida Courant-Friedrichs-Lewy (CFL) -luvun avulla, joka on määritelty muodossa
CFL=c0∆t
∆x. (20)
CFL-luvun avulla ilmaistuna pseudospektraalink-avaruusmenetelmän staabiilisuusehdok- si saadaan [39]
sin
πCFL 2
≤ c0
cmax (21)
missä cmax on suurin mahdollinen äänen nopeus. Homogeenisen väliaineen tapauksessa cmax = c0, jolloin stabiiliusehto on voimassa kaikilla spatiaali- ja aikadiskretisaation ar- voilla. k-avaruusmenetelmä homogeenisessä väliaineessa on näin ollen stabiili riippumat- ta käytetystä aika- tai spatiaalidiskretisaatiosta. Pseudospektraalink-avaruusmenetelmän tapauksessa tyypillinen riittävän tiheän aikadiskretisaation määrittämisessä käytetty CFL-
luku on 0.3, joka tarjoaa hyvän tasapainon laskenta-ajan ja laskennallisen tarkkuuden vä- lillä. [47]
Pseudospektraalissak-avaruusmenetelmässä käytetyt diskretisaatiot ovat menetelmän stabiiliudesta huolimatta rajoitettuja Nyquistin ehdon nojalla, jolloin aika- ja spatiaali- diskretisaatioille on voimassa ehto [27]
∆t≤ 1
2fmax = π
cmaxkmax = ∆x
cmax, (22)
missä fmax on suurin mahdollinen laskentahilan tukema taajuus, cmax on suurin simulaa- tiossa käytetty äänen nopeus ja kmax on suurin mahdollinen aaltoluku. Menetelmä voi näin ollen olla stabiili, vaikka Nyquistin ehto ei täyty. Esimerkiksi liian harvan aikadis- kretisaation tapauksessa menetelmä saattaa olla stabiili, vaikka ratkaisu menettää tiedon akustisessa kentässä etenevistä korkeataajuisista aalloista.
4 Bayesilainen käänteisongelma ja estimointi fotoakus- tisessa tomograassa
Käänteisongelmat ovat matemaattisia ongelmia, joissa pyritään tuottamaan tietoa fysi- kaalisesta kohteesta tai ilmiöstä siitä mitatun datan avulla. Mielenkiinnon kohteena ovat usein suureet, joita ei voida suoraan mitata. Esimerkkejä käänteisongelmista ovat muun muassa tietokonetomograakuvantaminen, jossa kuvauskohteen vaimennuskerroinjakau- ma määritetään kohteesta mitattujen läpivalaisukuvien perusteella ja geotomogranen maanjäristyksien paikantaminen [13, 21]. Fotoakustisen tomograan tapauksessa mielen- kiinnon suureina ovat kuvauskohteen sisäiset alkupaineen p0 arvot, joita pyritään estimoi- maan kohteesta mitattujen fotoakustisten paine-aikasarjojen pt avulla.
Reaalimaailman käänteisongelmat ovat usein huonokuntoisia. Käänteisongelmaa kut- sutaan huonokuntoiseksi, jos se ei täytä hyväkuntoisen käänteisongelman määritelmiä
1. Ratkaisu on olemassa 2. Ratkaisu on yksikäsitteinen
3. Ratkaisu riippuu jatkuvasti mittausdatasta.
Käänteisongelman huonokuntoisuus korostuu usein tilanteessa, jossa mittausdata sisältää merkittäviä mittausvirheitä ja käänteisongelma on epästabiili. Epästabiilin käänteisongel- man tapauksessa pieni muutos mittausdatassa saattaa aiheuttaa merkittävän muutoksen käänteisongelman ratkaisuun. [40]
Huonokuntoisten käänteisongelmien ratkaisemista lähestytään tyypillisesti joko deter- ministisestä tai tilastollisesta näkökulmasta. Deterministisissä menetelmissä käänteison- gelman huonokuntoisuutta pyritään parantamaan käyttämällä erilaisia regularisaatiome- netelmiä, kuten Tikhonov regularisaatiota. Tilastollisessa lähestymistavassa käänteison- gelman ratkaisua tutkitaan mallintamalla kaikki käänteisongelmassa esiintyvät muuttujat satunnaismuuttujina ja muodostamalla ratkaisu todennäköisyysjakaumien avulla. [17, 40]
Deterministisisessä lähestymistavassa käänteisongelman ratkaisuna saadaan yksittäi- siä arvoja, jotka eivät anna tietoa esimerkiksi ratkaistujen parametrien mahdollisesta ja- kaumasta tai todennäköisistä esiintymisväleistä. Käänteisongelmien tilastollisessa lähes- tymistavassa ratkaisu ei puolestaan ole yksittäinen arvo vaan todennäköisyysjakauma, joka tarjoaa deterministisiin menetelmiin verrattuna enemmän tietoa. Ratkaisun tilas- tollinen luonne tuottaa yksittäisten piste-estimaattien lisäksi informaatiota estimoitavien parametrien arvojen todennäköisyyksistä ja varmuusväleistä. [17, 40]
Tässä työssä fotoakustisen tomograan käänteisongelmaa lähestytään Bayesilaisten käänteisongelmien lähtökohdasta, jossa tavoitteena on ratkaista estimoitavien parametrien p0 ja mittauksien pt ehdollinen todennäköisyys- eli posteriorijakauma π(p0|pt). Bayesin kaavan nojalla posteriorijakauman todennäköisyystiheysfunktio voidaan esittää muodossa
π(p0|pt) = π(pt|p0)π(p0)
π(p) ∝π(pt|p0)π(p0), (23)
missä π(p0) on etukäteistietoon perustuva priorijakauma, π(pt|p0) on likelihood- eli us- kottavuusjakauma ja π(pt)on normalisointivakio. Bayesilaisen käänteisongelman ratkaisu voidaan näin ollen jakaa kolmeen osaan: [17]
1. Uskottavuusjakauman muodostaminen suoran mallin ja mittausvirheen perusteella 2. Priorijakauman muodostaminen etukäteistiedon perusteella
3. Posteriorijakauman ja piste-estimaattien määrittäminen.
Bayesilaisen käänteisongelman ratkaisuna saatavan posteriorijakauman määrittämiseksi tarvitaan näin ollen tietoa sekä priori- että uskottavuusjakaumien rakenteesta.
4.1 Uskottavuusjakauma
Uskottavuusjakauma on todennäköisyysjakauma, joka kuvaa todennäköisyyttä saada ti- lanteeseen liittyvä kohinainen mittausdata, kun estimoitavat parametrit on tunnettu. Us- kottavuusjakauma sisältää tietoa käänteisongelmassa käytetystä mallista sekä mittausvir- heestä, jota voidaan mallintaa usealla eri tavalla. Niin kutsutun normaalin virhemallin (engl. conventional error model, CEM) tilanteessa mittausvirhe oletetaan additiiviseksi.
Tässä työssä tarkasteltavassa fotoakustisessa tomograssa havaintomalli on lineaarinen, jolloin se voidaan kirjoittaa normaalin virhemallin avulla diskreetissä muodossa
pt=Kp0+e, (24)
missä K ∈ Rm×n on alkupaineen p0 ∈ Rn ja mittaukset pt ∈ Rm yhdistävä lineaarinen matriisioperaattori.
Tarkastellaan nyt uskottavuusjakauman muodostamista normaalin virhemallin tapauk- sessa, jolloin tilannetta kuvaava lineaarinen havaintomalli on muotoa (24). Kun havain- tomallissa esiintyvät alkupaine p0 ja mittausvirhe e ovat tiedossa, voidaan ehdollinen jakauma π(pt|p0, e) kirjoittaa Diracin deltan δ avulla muodossa [7]
π(pt|p0, e) =δ(pt−Kp0−e). (25) Mittauksien pt, alkupaineen p0 ja mittausvirheen e yhteisjakaumalle voidaan puolestaan kirjoittaa Bayesin kaavan (23) avulla
π(pt, p0, e) =π(pt|p0, e)π(p0, e)
=π(pt|p0, e)π(e|p0)π(p0)
=π(pt, e|p0)π(p0).
(26)
Uskottavuusjakauma π(pt | p0) voidaan tällöin kirjoittaa yhtälön (26) avulla ehdollisena todennäköisyysjakaumana muodossa [17, 20]
π(pt|p0) = Z
π(pt, e|p0) de
= Z
π(pt|p0, e)π(e|p0) de
= Z
δ(pt−Kp0−e)π(e|p0) de
=πe|p0(pt−Kp0|p0),
(27)
missä viimeiseen muotoon päästään Diracin deltan siirto-ominaisuudenR
f(x)δ(x−X)dx= f(X)avulla. Tilanteessa, jossa alkupainep0 ja mittausvirhee oletetaan toisistaan riippu- mattomiksi 3, voidaan uskottavuusjakauma kirjoittaa muodossa [17, 20]
π(pt|p0) =πe(pt−Kp0). (28) Oletetaan nyt, että mittausvirhe on Gaussisesti jakautunutta
e∼ N(ηe,Γe),
missä ηe on mittausvirheen odotusarvo ja Γe on mittausvirheen kovarianssimatriisi. Us- kottavuusjakauma voidaan tällöin kirjoittaa yhtälön (28) nojalla Gaussisessa muodossa [17]
π(pt|p0)∝exp
−1
2(pt−Kp0−ηe)TΓ−1e (pt−Kp0−ηe)
∝exp
−1
2kLe(pt−Kp0−ηe)k22
,
(29)
missä Le on mittausvirheen kovarianssimatriisin käänteismatriisin Γ−1e Choleskyn hajo- telma.
4.2 Etukäteistieto ja priorijakauma
Bayesilaisen käänteisongelman ratkaisussa käytetty etukäteistieto esitetään todennäköi- syysjakauman muodossa. Tilanteeseen sopivan priorijakauman muodostaminen on kään- teisongelman ratkaisemisen kannalta tärkeässä osassa erityisesti huonokuntoisten kääntei- songelmien tapauksissa. Näissä tilanteissa estimoitavia parametreja hyvin kuvaava priori- jakauma ohjaa estimaattia haluttuun suuntaan ja mahdollistaa käänteisongelman ratkai- suna saatavan estimaatin tulkinnan.
Yleisesti priorijakauma pyrkii esittämään alkupaineesta p0 tunnetun etukäteistiedon todennäköisyysfunktioiden avulla niin, että
π(p0)π(p00), kun p0 ∈E, p00 ∈E0,
missä E on odotettujen vektorien joukko ja E0 on odottamattomien vektorien joukko.
Priorijakauman tulisi siis olla keskittynyt niihin arvoihin, joita käänteisongelman ratkai- sun oletetaan sisältävän. Priorijakaumaan implementoitava ennakkotieto voi olla luonteel- taan kvalitatiivista, kuten esimerkiksi rekonstruoidun kuvan arvojen jakaumien muotoi- hin liittyvää tietoa, tai kvantitatiivista, kuten rekonstruoidun kuvan numeerisiin arvoihin liittyvää tietoa. [17]
Yleisiä Bayesilaisessa käänteisongelmien lähestymistavassa käytettyjä prioreita ovat Gaussiset priorijakaumat, jotka voidaan määritellä niiden odotusarvojen ηp0 ja kovarians- simatriisien Γp0 avulla. Esimerkkejä Gaussisista priorijakaumista ovat muun muassa si- leysoletukselliset neliöllisen eksponentin (engl. squared exponential) ja Matérn-tyyppiset
3Tällöin π (e|p ) =π (e).
priorijakaumat sekä estimoitavien parametrien toisistaan riippumattomuuden olettava valkoisen kohinan (engl. white noise) priorijakauma. [38]
Neliöllisen eksponentin priorijakauma on määritetty Gaussisena jakaumana, jossa ko- varianssimatriisi on muotoa [38, 44]
Γp0,ij =σ2p0exp
−kri−rjk2 2l2
, (30)
missä σp20 on estimoitavien parametrien varianssi, ri ja rj ovat estimoitavien parametrien pikselien paikat ja l on priorin karakteristinen pituus, jolla säädetään priorin spatiaalista sileyttä.
Ornstein-Uhlenbeck priori, joka kuuluu Matérn priorien luokkaan, on neliöllisen eks- ponentin tapaan Gaussinen sileyspriori, ja sen kovarianssimatriisi on määritelty muodossa [38, 44]
Γp0,ij =σp20exp
−kri−rjk l
. (31)
Kovarianssimatriisin eksponenttitermi on tässä tapauksessa neliöllisen eksponentin kaltai- nen, lukuun ottamatta pikselien etäisyyksien ja sileysparametrin neliötä. [38]
Gaussinen priori voidaan määrittää myös spatiaalisesti korreloimattomaksi eli valkoi- seksi kohinaksi. Tässä tapauksessa kovarianssimatriisi on muotoa [38, 44]
Γp0 =σ2p0I, (32)
missä I on yksikkömatriisi.
4.3 Posteriorijakauma
Bayesilaisen käänteisongelman posteriorijakauma voidaan muodostaa, kun uskottavuus- ja priorijakaumat ovat tiedossa. Posteriorijakauma yhdistää mittauksista saadun infor- maation kohteesta tunnettuun etukäteistietoon.
Tarkastellaan nyt posteriorijakauman muodostamista tilanteessa, jossa havaintomalli on yhtälössä (24) esitettyä muotoa. Oletetaan, että alkupaine p0 on Gaussisesti jakautu- nutta ja mittausvirhe e sekä alkupaine p0 ovat toisistaan riippumattomia. Posteriorija- kauma voidaan tällöin kirjoittaa Bayesin kaavan (23) nojalla uskottavuusjakauman (28) ja priorijakauman tulona muodossa [17, 20]
π(p0|pt)∝π(pt|p0)π(p0)
∝πe(pt−Kp0)π(p0)
∝exp
−1
2(pt−Kp0−ηe)TΓ−1e (pt−Kp0−ηe)
·exp
−1
2(p0−ηp0)TΓ−1p0 (p0−ηp0)
∝exp
−1
2 kLe(pt−Kp0−ηe)k22+kLp0(p0−ηp0)k22
,
(33)
missä Bayesin kaavassa esiintyvä mittauksien todennäköisyysjakauma π(pt) on siirretty normalisointivakioon jaLp0 on priorijakauman kovarianssimatriisin käänteismatriisin Γ−1p Choleskyn hajotelma. 0
Avaamalla yhtälössä (33) esiintyvät matriisitulot ja siirtämällä alkupaineen arvoja p0
sisältämättömät termit normalisointivakioon, voidaan posteriorijakauma kirjoittaa Gaus- sisesti jakautuneena muodossa (liite 1.1)
p0|pt ∼ N(ηp0|pt,Γp0|pt), (34) missä
(ηp0|pt = (Γ−1p0 +KTΓ−1e K)−1(KTΓ−1e (pt−ηe) + Γ−1p0 ηp0) Γp0|pt = (Γ−1p0 +KTΓ−1e K)−1.
4.4 Bayesilainen approksimaatiovirhemenetelmä
Mallipohjainen käänteisongelmien ratkaisu perustuu vahvasti käytettyyn malliin ja sen tarkkuuteen. Käänteisongelmien ratkaisussa voidaan päätyä tilanteeseen, jossa käytet- ty malli on epätarkka. Tämän kaltaisia tilanteita ovat esimerkiksi harvan spatiaali- tai aikadiskretisaation tai tuntemattomien parametrien, kuten sensorien paikkojen, aiheut- tamat mallinnusvirheet. Epätarkan mallin aiheuttama mallinnusvirhe saattaa vaikuttaa merkittävästi käänteisongelman ratkaisuun. Mallinnusvirheiden aiheuttamia virheitä voi- daan pyrkiä korjaamaan niin kutsutulla Bayesilaisella approksimaatiovirhemenetelmällä.
[3, 17, 22, 31, 32]
Tarkastellaan tilannetta, jossa mallin tarkkuuteen vaikuttavaa parametria, kuten sen- soripaikkoja tai laskentahilan tiheyttä, merkitään muuttujalla z. Tilannetta kuvaava ad- ditiivisen virheen oletuksen sisältävä tarkka lineaarinen havaintomalli on tällöin muotoa [17, 18, 19]
pt=K(z)p0+e, (35)
Käänteisongelmien tapauksessa voidaan kuitenkin päätyä tilanteeseen, jossa tarkkaa mal- lia ei voida käyttää esimerkiksi riittämättömien laskentaresurssien takia. Tällöin käytössä saattaa olla niin kutsuttu redusoitu malli, jossa mallin tarkkuuteen vaikuttava satunnais- muuttuja z oletetaan vakioksi z0 ja havaintomalli on muotoa [17, 18]
pt=K(z0)p0 +e, (36)
Havaintomallin (36) käyttö saattaa johtaa merkittävään mallinnusvirheeseen ja edelleen virheeseen käänteisongelman ratkaisussa.
Mallinnusvirhettä voidaan pyrkiä korjaamaan käyttämällä tarkkaan malliin (35) pe- rustuvaa approksimaatiovirhemallia [17, 18, 19]
pt =K(z)p0+e
=K(z0)p0+ [K(z)p0−K(z0)p0] +e
=K(z0)p0+ε+e,
(37)
missä satunnaismuuttujaa ε = ε(z, p0) kutsutaan approksimaatiovirheeksi. Bayesilaisen käänteisongelman ratkaisu tilanteessa, jossa havaintomalli on muotoa (37), voidaan muo- dostaa normaalin virhemallin tapaan muodostamalla havaintomalliin liittyvät uskottavuus- ja priorijakaumat.
Uskottavuusjakauma
Tarkastellaan approksimaatiovirhemalliin (37) liittyvää uskottavuusjakaumaa. Kun mallit K(z) ja K(z0) on kiinnitetty, ja alkupaine p0 sekä mittausvirhe e ovat tiedossa, voidaan mittauksienptehdollinen todennäköisyysjakauma muiden mallissa esiintyvien muuttujien suhteen kirjoittaa Diracin deltan avulla muodossa [7, 23]
π(pt|p0, z, ε, e) =δ(pt−K(z0)p0−ε−e). (38) Bayesin kaavan nojalla muuttujien pt, p0, z, ε, e yhteisjakaumalle voidaan edelleen kirjoit- taa [23]
π(pt, p0, z, ε, e) = π(pt|p0, z, ε, e)π(p0, z, ε, e)
=π(pt|p0, z, ε, e)π(e, ε|p0, z)π(z|p0)π(p0)
=δ(pt−K(z0)p0−ε−e)π(ε, e|p0, z)π(z |p0)π(p0)
=π(pt, z, ε, e|p0)π(p0).
(39)
Jos approksimaatiovirhe ε oletetaan parametrista z riippumattomaksi satunnaismuuttu- jaksi, voidaan uskottavuusjakauma kirjoittaa yhtälön (39) nojalla integraalimuodossa [23]
π(pt|p0) = Z Z Z
π(pt, z, ε, e|p0) dzdεde
= Z Z
δ(pt−K(z0)p0−ε−e) Z
π(ε, e|p0, z)π(z |p0) dz
dεde
= Z Z
δ(pt−K(z0)p0−ε−e)π(ε, e |p0) dεde
= Z
πe(pt−K(z0)p0−ε−e)πε|p0(ε|p0) dε.
(40)
Merkitään nyt approksimaatiovirhemallissa esiintyvien virheiden jakaumien summaa muut- tujalla
n|p0 =e+ε|p0. (41)
Oletetaan, että mittausvirhee, approksimaatiovirhe ε|p0 ja alkupainep0 ovat Gaussisesti jakautuneita
e∼ N(ηe,Γe) ε|p0 ∼ N(ηε|p0,Γε|p0) p0 ∼ N(ηp0,Γp0),
jolloin myös kokonaisvirhevektori n on Gaussisesti jakautunut. Uskottavuusjakauma (40) voidaan tällöin tulkita kahden Gaussisesti jakautuneen todennäköisyysjakauman konvo- luutiointegraalina, jolle on voimassa [23]
