Thomas Dierkes Rekonstruktionsverfahren zur

optischen Tomographie

Für Bene

Angewandte Mathematik

Rekonstruktionsverfahren zur optischen Tomographie

Inaugural-Dissertation zur Erlangung des Doktorgrades

der Naturwissenschaften im Fachbereich

Mathematik und Informatik

der Mathematisch-Naturwissenschaftlichen Fakultät der Westfälischen Wilhelms–Universität Münster

vorgelegt von

Thomas Dierkes aus Münster

Dekan:  Prof. Dr. W. Lange
Erster Gutachter:  Prof. Dr. F. Natterer
Zweiter Gutachter:  Prof. Dr. C.W. Cryer

Tag der mündlichen Prüfungen: Tag der Promotion:

Inhaltsverzeichnis

1 Einleitung 3

2 Physikalische Problemstellung 9 2.1 Diffusionsapproximation . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Quellterm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4 Modellgleichung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.5 Datenmodellierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3 Mathematischer Teil 15

3.1 Direktes Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.1.1 Definitionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.1.2 Existenz und Eindeutigkeit einer Lösung . . . . . . . . . . . . . . 20 3.1.3 Stetigkeit der Dirichlet-zu-Neumann–Abbildung . . . . . . . . . . 23 3.1.4 Differenzierbarkeit der Dirichlet-zu-Neumann–Abbildung . . . . . 25 3.1.5 Greensche Funktion . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.1.6 Born–Näherung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.1.7 Rytov–Näherung . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.1.8 Liouville–Transformation . . . . . . . . . . . . . . . . . . . . . . 32 3.1.9 Transformation der Fréchet–Ableitung . . . . . . . . . . . . . . . 32

INHALTSVERZEICHNIS

3.2 Inverses Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2.1 Eindeutigkeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2.2 Schlecht gestelltes Problem . . . . . . . . . . . . . . . . . . . . . 38 3.2.3 Iterativer Lösungsansatz . . . . . . . . . . . . . . . . . . . . . . . 39 3.2.4 Adjungierte Abbildung . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2.5 Verallgemeinerte Lösung . . . . . . . . . . . . . . . . . . . . . . . 44 3.2.6 Sherman–Morrison–Woodbury–Formel . . . . . . . . . . . . . . . 47 3.2.7 Vorteil der komplexwertigen Rechnung . . . . . . . . . . . . . . . 50 3.2.8 Tikhonov–Regularisierung . . . . . . . . . . . . . . . . . . . . . . 52 3.2.9 Kaczmarz–Newton–Methode . . . . . . . . . . . . . . . . . . . . 53 3.2.10 Konvergenz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4 Numerischer Teil 57

4.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.1.1 Finite Differenzen . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.1.2 Orthogonale Geometrie . . . . . . . . . . . . . . . . . . . . . . . 61 4.1.3 Zylindrische Geometrie . . . . . . . . . . . . . . . . . . . . . . . 62

4.2 Numerische Experimente . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2.1 Rekonstruktionsgebiet . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2.2 Datenerzeugung . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2.3 Rekonstruktionen . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

4.3 Diskussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

5 Zusammenfassung 85

Literaturverzeichnis 87

Kapitel 1

Einleitung

Ausgangspunkt der optischen Tomographie

Ausgangspunkt der optischen Tomographie ist die Verwendung energiearmer Strahlung im nahen Infrarot-Bereich (NIR, 700 - 1000 nm) als Informationsträger, um Kenntnis über die optische Beschaffenheit des Inneren eines für diese Infrarotstrahlung durchlässigen Gegenstandes oder Objektes zu erlangen, ohne diesen Gegenstand aufzuschneiden oder anderweitig zu zerstören. Dazu wird der fragliche Gegenstand von mehreren Seiten mit Hilfe von Glasfasern, die das NIR-Licht leiten, abgetastet und die aus dem Gegenstand austretende Strahlung mittels weiterer Glasfasern aufgezeichnet.

Die Entwicklung und Untersuchung geeigneter numerischer Verfahren zur Rekonstruktion der optischen Beschaffenheit des durchleuchteten Objektes aus diesen zu messenden Daten mit einem Computer ist Gegenstand dieser Arbeit.

Anwendungsbeispiele

Im Bereich der medizinischen Diagnostik ist eine ähnliche Art der Abtastung seit der Erfindung der Computertomographie (CT) in den frühen 60iger Jahren bekannt. Im Gegensatz zur optischen Tomographie wird bei der CT die energiereichere und damit für biologisches Gewebe schädlichere Röntgenstrahlung zur Durchleuchtung benutzt. Neben der erheblichen Kostenreduzierung hätte der Einsatz der optischen Tomographie zu diagnostischen Zwecke auch den Vorteil der geringeren Strahlenbelastung, was insbesondere in bezug auf mammographische Reihenuntersuchungen, dem sogenannten screening, von sehr hohem Interesse ist. In den Industrieländern ist bei Frauen eine der häufigsten Todesursachen Brustkrebs, welcher bei entsprechender Früherkennung mitunter erfolgreicher behandelt werden könnte.

Grundlegend für die Entwicklung der Forschung der optischen Tomographie für medizinische Zwecke kann die Arbeit von F. F. Jöbsis [26] aus dem Jahre 1977 angesehen werden, in der die grundsätzliche Anwendbarkeit der optischen Tomographie zu diagnostischen Zwecken aufgezeigt wurde.

KAPITEL 1. EINLEITUNG

Andere Anwendungen der optischen Tomographie sind darüberhinaus in der Verfahrenstechnik, z.B. zur zerstörungsfreien Prüfung von Werkstücken, und in der chemischen Industrie, z.B. als Kontrollinstrument für chemische Reaktionen in flüssiger Lösung, die für dieses elektomagnetische Spektrum durchlässig sind, zu finden.

Nachteile der optischen Tomographie

Die Verwendung energieärmerer Strahlung bei der optischen Tomographie hat zur Folge, daß wegen der hohen Streuung im Inneren des Objektes die gemessenen Daten erheblich geringeren Informationsgehalt besitzen als etwa CT-Daten.

Darüberhinaus erfährt das Eingangssignal neben der Streuung zusätzlich eine nicht zu vernachlässigende Absorption. Dies macht es erforderlich, daß die Meßgeräte, die die austretene Strahlung aufnehmen, den weiten dynamischen Bereich der Daten abdecken. Es werden also besondere Anforderungen an die verwendete Meßtechnik zu stellen sein, siehe z.B. [24], [4].

Stand der Forschung

Einen recht umfassenden Überblick gibt der Artikel von S. R. Arridge [2] aus dem Jahr 1999.

Grundsätzlich kann die zeitliche Entwicklung der Photonendichte in einem (stark) streuenden Medium mit der linearen Transportgleichung, oder auch bekannt als lineare Boltzmann-Gleichung, beschrieben werden, siehe [28], [15]. Prinzipiell ist diese Gleichung zur Beschreibung der physikalischen Vorgänge und damit für etwaige Rekonstruktionsverfahren zu verwenden, stellt jedoch erhebliche Anforderungen an die vorhandene Rechnerkapazität, da sie numerisch nur sehr aufwendig zu lösen ist.

In etlichen Untersuchungen konnte festgestellt werden, daß die sogenannte Diffusionsnäherung der Boltzmann-Gleichung eine hinreichend genaue Approximation zur Simulation der optischen Tomographiedaten ermöglicht, siehe z.B. [59], die zudem mit viel geringerem numerischen Aufwand zu lösen ist.

Durch die Verwendung der Diffusionsnäherung ergibt sich auch eine einfachere Modellierung des physikalischen Meßprozesses im Frequenzbereich, bei dem im Gegensatz zur temporal aufgelösten Messung die Amplitudendämpfung und die Phasenverschiebung zu einem gegebenen zeitharmonisch modulierten Trägersignal registriert wird. Ein Nachteil dieser Messungen im Frequenzbereich liegt in der aufwendigeren und daher kostspieligeren experimentellen Realisierung begründet. Dennoch sei in dieser Arbeit vornehmlich diese Art der Messung stets angenommen, da hierbei für etwaige Rekonstruktionen erheblich weniger Daten als im zeitaufgelösten Fall zu verarbeiten sind, was erneut einen Geschwindigkeitsvorteil bedeutet. Tatsächlich kann mit den Ergebnissen dieser Arbeit belegt werden, daß Rekonstruktionen mit zeitaufgelösten Daten keine genaueren Ergebnisse liefern als Rekonstruktionen mit den hier vorgeschlagenen Algorithmen, die sämtlich mit Daten aus dem Frequenzbereich operieren.

Eine zusätzliche Schwierigkeit, die in der Literatur bisher kaum beachtet wurde, entsteht bei der hier betrachteten optischen Tomographie dadurch, daß bei der Rekonstruktion aus den gemessen Daten zwei simultan zu bestimmenden Größen auftreten, Streuung und Absorption. Arbeiten, die diesen Umstand nicht vernachlässigen, sind zum Beispiel [69], [58], [34].

Rekonstruktionsmethoden

Zunächst kann man zur Lösung des Problems der optischen Tomographie dieses als inverses Problem einordnen, im Gegensatz zu dem direkten Problem, bei gegebenen optischen Parametern die Meßdaten analytisch oder numerisch zu bestimmen.

In dem Übersichtsartikel [2] werden in Analogie zu bestehenden Verfahren zur (klassischen) Tomographie insgesamt drei Methodenklassen für die optische Tomographie unterschieden:

• SPECT-ähnliche Methoden (\underline{s}ingle \underline{p}hoton \underline{e}mission \underline{c}omputer \underline{t}omography): lineare Transferfunktion, Matrixinversion, siehe z.B. [27, Kap. 4.2] für SPECT.

• CT-ähnliche Methoden (\underline{c}omputer \underline{t}omography):

inverse Radon-Transformation, gefilterte Rückprojektion, siehe [25], [47] für CT.

• EIT-ähnliche Methoden (\underline{e}lectric \underline{i}mpedance \underline{t}omography):

Minimierung eines Kostenfunktionals, Iteration, siehe z.B. [13] für EIT.

Die ersten beiden Methodenklassen basieren jeweils auf einem linearen Modellansatz, wohingegen die dritte Klasse den nichtlinearen Charakter der Problemstellung durch den Ansatz als ein Optimierungsproblem bzw. durch die Verwendung iterativer Verfahren zur Lösung nichtlinearer Probleme berücksichtigt, siehe [35], [29]. In diese dritte Klasse fallen auch die in dieser Arbeit angegebenen Kaczmarz–Newton–Verfahren. Hierbei sei ausdrücklich darauf hingewiesen, daß diese Verfahren keine Gradientenverfahren oder Abstiegsverfahren, wie z.B. die Landweber–Iteration, darstellen, sondern eine Erweiterung der aus der klassischen Computertomographie bekannten linearen ART-Methode (\underline{a}lgebraic \underline{r}econstruction \underline{t}echnique) sind.

Beispiele für Arbeiten über einfache Rückpropagationsmethoden zur optischen Tomographie sind [8], [16].

Eine Vielzahl an Veröffentlichungen benutzt den Ansatz, die nichtlineare Abbildung, die die gesuchten optischen Parameter den zu messenden Daten zuordnet, durch Linearisierung auf eine lineare Abbildung zu reduzieren, um die daraus resultierende lineare Gleichung mit bekannten Methoden direkt oder iterativ zu lösen.

Als nichtlineare iterative Methoden zur optischen Tomographie mit frequenzmodulierten Daten sind zum Beispiel das Newton–Raphson–Verfahren (siehe [54]), das Levenberg–Marquardt–Verfahren (siehe [60]), eine Born–Iterations–Methode (siehe [69]) zu nennen.

Ein spektakulärer Ansatz zur optischen Tomographie mit zeitaufgelösten Daten wird in der Arbeit [31] verfolgt. Ausgehend von der Diffusionsnäherung, einer partiellen Differentialgleichung zweiter Ordnung, bezüglich einer als bekannt vorausgesetzten Startnäherung für die gesuchten Parameter wird ein elliptisches System vierter Ordnung hergeleitet. Die Lösung dieses Systems wird dazu verwendet, die Startnäherung additiv zu korrigieren, um damit als letzten Schritt der Iteration eine neue Diffusionsnäherung zu bestimmen.

Einige Resultate dieser Arbeit

Ein Umstand, der gänzlich bei Rekonstruktionsverfahren zur optischen Tomographie im Frequenzbereich außer Acht gelassen wurde, ist der, daß die gemessenen Daten zwar als komplexwertige Größen gegeben sind, die unbekannten Parameter jedoch als rein reelle Größen zu bestimmen sind. Diese Arbeit zeigt auf, daß die korrekte Berücksichtigung dieses Umstandes in bezug auf das vorgestellte nichtlineare Kaczmarz–Newton– Verfahren nicht in der simplen Einschränkung der Rekonstruktionen auf den Realteil besteht, wie vielleicht zu vermuten wäre.

Eben dies bringt gerade bei der Implementation des hier vorgestellten Verfahrens einen wesentlichen Gewinn in der Qualität der erzielten Rekonstruktionsergebnisse im Vergleich zu Ergebnissen aus bisherigen Veröffentlichungen zur optischen Tomographie. Damit ist auch begründet, daß zeitaufgelöste Daten im Vergleich zu frequenzmodulierten Daten keine zusätzlichen Informationen zur optischen Tomographie enthalten.

Besonders der Effekt des Übersprechens in den Rekonstruktionen der beiden Parameter wird durch die Verwendung des vorgestellten Algorithmus reduziert.

Die Anmerkungen zum Schluß des Artikels [2] hinsichtlich der Problematik der Rekonstruktion reellwertiger Parameter aus komplexwertigen Daten, daß die simple Einschränkung auf den jeweiligen Realteil der zu bestimmenden Größen durchaus korrekt sei, müssen somit revidiert werden. Zwar mag dieser Ansatz physikalisch vollkommen korrekt sein, spielt aber in Abhängigkeit des eingesetzten Rekonstruktionsverfahrens durchaus eine entscheidene Rolle.

Kapitelübersicht

– Kapitel 2 Physik –

In diesem Kapitel werden die physikalischen Hintergründe erläutert, die wesentlich zur mathematischen Modellbildung der optischen Tomographie beitragen.

– Kapitel 3 Mathematik –

Das mathematische Kapitel stellt die wesentlichen Aussagen bereit, die zur Aufstellung der iterativen Rekonstruktionsverfahren benötigt werden, welche in dieser Arbeit für die speziellen Gegebenheiten der optische Tomographie angepaßt werden. Insbesondere wird die Problematik der reellwertigen Rekonstruktion aus komplexwertigen Daten behandelt.

– Kapitel 4 Numerik –

Dieses Kapitel enthält die genaue Beschreibung des implementierten Kaczmarz-Newton-Verfahrens zur Lösung des inversen Problems der optischen Tomographie, samt der Diskretisierungen der elliptischen Randwertprobleme für unterschiedliche Geometrien des Grundgebietes anhand der Methode der finiten Differenzen.

Die konkreten Konstanten für die Berechnungen werden angegeben und numerische Resultate für 2D-Modellrechnungen gezeigt.

Schließlich werden in diesem Kapitel die vorgestellten Ergebnisse erläutert und der Zusammenhang zu bisherigen Resultaten diskutiert.

Danksagung

An dieser Stelle sei Herrn Prof. Dr. F. Natterer ganz herzlich für die Anregung zu dieser Arbeit und für viele hilfreiche Hinweise und Diskussionen gedankt.

Einen nicht unerheblichen Teil am Gelingen dieser Arbeit ist auch der herzlichen und freundlichen Atmosphäre unter den Mitgliedern des Institus für Numerische und instrumentelle Mathematik der Westfälischen Wilhelms-Universität Münster zu verdanken. Ihnen gebührt daher ebenfalls ein großes Dankeschön.

Insbesondere zu nennen sind hier Herr Dipl.-Math. Helmut Sielschott und, als ehemaliges Mitglied des Instituts, Herr Dr. Oliver Dorn, die mir bei der Durchsicht des Manuskripts kritisch und hilfreich zur Seite standen.

KAPITEL 1. EINLEITUNG

8

Physikalische Problemstellung

Vor gut zwei Jahrzehnten hat man festgestellt, daß biologisches Gewebe in einem Wellenbereich von etwa 700 nm bis 1000 nm des elektromagnetischen Spektrums (naher Infrarotbereich) relativ geringere Absorption aufweist, siehe z.B. [26], [53], [43]. Der Gedanke, dieses physikalische Phänomen für diagnostische Zwecke zu nutzen, liegt auf der Hand. Die Idee ist, daß man die Intensität und die Phase des Lichtes aus diesem Spektrum, welches durch ein zu testendes Objekt geschickt wurde, als Daten sammelt. Nun stellt sich die Frage, ob es eine eindeutige Zuordnung dieser Daten zu den optischen Parametern im Inneren des Objektes gibt. Dies konnte bisher nur unter sehr restiktiven Voraussetzungen an die zu bestimmenden optischen Parameter positiv beantwortet werden, [20]. Jedoch sei, gewissermaßen als Ausgangspunkt für diese Arbeit, angenommen, daß die fragliche Injektivität der beschriebenen Datenabbildung stets gegeben sei. Demnach ist als nächstes zu klären, wie eine brauchbare Rekonstruktion der gesuchten optischen Parameter aus gemessenen Daten zu erreichen ist. In dieser Arbeit werden die mathematischen Zusammenhänge betrachtet, die bei solch einer tomographischen Aufgabenstellung auftreten, um daraus numerische Verfahren zur Rekonstruktion zu entwickeln. Zunächst soll das zugrunde gelegte physikalische Modell erläutert werden.

2.1 Diffusionsapproximation

Die zeitliche Dichteverteilung Φ(x, t) eines Stoffes in einem Gebiet Ω ⊂ IR^{n} (n = 2, 3, …) mit einem stark streuenden Medium läßt sich durch die partielle Differentialgleichung

\frac{1}{c} \frac{∂Φ}{∂t} − ∇_{x} · (κ(x)∇_{x}Φ) + μ(x)Φ = q(x, t; y_{0}), y_{0} ∈ ∂Ω × \{0\},  (2.1)

beschreiben. Dabei bedeuten

KAPITEL 2. PHYSIKALISCHE PROBLEMSTELLUNG

c = c_{medium}  die Ausbreitungsgeschwindigkeit des Stoffes

im Medium,

μ(x) = μ_{a}(x)  der mittlere Absorptionsquerschnitt,
κ(x) = \frac{1}{3}[μ′_{s}(x) + μ_{a}(x)]^{−1} der Diffusionskoeffizient und
q(x, t; y_{0}) = q(x − y_{0})δ(t) ein Quellterm abhängig von y_{0} ∈ ∂Ω.

Der Term μ′_{s}(x) = (1 − g)μ_{s}(x) wird der reduzierte Streuquerschnitt genannt, wobei g der mittlere Cosinus einer Streufunktion φ(s, ŝ) ist,

g = \frac{1}{ω_{n}} ∫_{S^{n−1}} φ(s, ŝ)s · ŝdŝ, s ∈ S^{n−1}.  (2.2)

Die Herleitung der Diffusionsgleichung (2.1) basiert auf einem Gleichgewichtsansatz. Ausgehend von der Annahme, daß keine Energie verloren gehen oder aus dem Nichts erzeugt werden kann, gelangt man für die Verteilungsdichte I(x, t, s) eines sich ausbreitenden Stoffes bezüglich einer Richtung s ∈ S^{n−1} zu der Differential-Integralgleichung

\frac{1}{c} \frac{∂I}{∂t} + s · ∇_{x}I + (μ_{s}(x) + μ_{a}(x))I = \frac{μ_{s}(x)}{ω_{n}} ∫_{S^{n−1}} φ(s, ŝ)I(x, t, ŝ) dŝ+ q(x, t, s; y_{0}), (2.3)

die als Boltzmann-Gleichung oder lineare Transportgleichung bezeichnet wird, [15]. Sie stellt eine Bilanzgleichung dar, welche die Änderungsrate der Dichteverteiltung I(x, t, s) zum Zeitpunkt t im Ort x ∈ IR^{n} und bezüglich der Richtung s aufgrund der Abnahme durch Absorption und Streuung einerseits und des Gewinns durch Produktion und Streuung aus anderen Richtungen andererseits beschreibt.

Diese kann unter der Voraussetzung der isotropen Streuung und der Bedingung, daß im zu betrachtenden Gebiet der mittlere Streuquerschnitt sehr viel größer als der Absorptionsquerschnitt ist,

μ_{s}(x) ≫ μ_{a}(x), ∀ x ∈ Ω,  (2.4)

in die geschlossene Form (2.1) gebracht werden, [28], [15].

2.2 Randbedingungen

Es wird ein beschränktes Gebiet Ω ⊂ IR^{n} mit hinreichend glattem Rand ∂Ω angenommen.

Als Anfangs-Randbedingungen für die partielle Differentialgleichung (2.1) kann man zunächst Dirichlet-Bedingungen setzen:

Φ(x, t_{0}) = 0 für x ∈ Ω,  (2.5)
Φ(y, t) = 0 für y ∈ ∂Ω, t_{0} ≤ t.  (2.6)

2.3 QUELLTERM

Physikalisch kommen diese Bedingungen einem total absorbierenden Rand gleich, so daß die Teilchendichte außerhalb Ω Null beträgt.

Falls Ω von einem nichtstreuenden Medium umgeben wird und keine Streuung an der Grenzfläche ∂Ω auftritt, wird dies durch Robin-Bedingungen modelliert:

Φ(x, t_{0}) = 0 für x ∈ Ω,  (2.7)
Φ(y, t) + 2κ(y)\frac{∂Φ}{∂ν_{y}} = 0 für y ∈ ∂Ω, t_{0} ≤ t.  (2.8)

Hierbei soll ν_{y} die äußere Normale zu ∂Ω im Punkt y ∈ ∂Ω bedeuten.

Der Fall, daß die Brechungsindizes des inneren und des äußeren Gebietes nicht übereinstimmen, und also Reflektionseffekte in Richtung des Inneren auftreten können, wird dadurch modelliert, daß die Bedingung (2.8) durch

Φ(y, t) + 2κ(y)\frac{∂Φ}{∂ν_{y}} = K (Φ(y, t) − 2κ(y)\frac{∂Φ}{∂ν_{y}} ) für y ∈ ∂Ω, t_{0} ≤ t,  (2.9)

mit einer Konstanten K ersetzt wird, die die Reflektionseigenschaften des Randes ∂Ω widerspiegelt. Mit A := (1 + K)/(1 − K) läßt sich (2.9) auch schreiben als

Φ(y, t) + 2A κ(y) \frac{∂Φ}{∂ν_{y}} = 0 für y ∈ ∂Ω, t_{0} ≤ t.  (2.10)

Die Konstante K bestimmt man experimentell zu

K ≈ −1.4399n^{−2} + 0.7099n^{−1} + 0.6681 + 0.0636n,  (2.11)

wobei n der Brechungsindex innerhalb Ω bedeutet, [59].

2.3 Quellterm

Zur Berechnung einer Lösung der Gleichung (2.1) ist neben den Anfangs-Randbe-dingungen auch die Angabe eines konkreten Quellterms q(x, t; y_{0}) notwendig. Dies soll jetzt geschehen.

Zu modellieren sei eine zeitabhängige Punktquelle im Ortsraum Ω. Dazu sei zunächst ein Punkt x_{0} ∈ Ω fest. Damit gelte für den Quellterm q(x, t; x_{0}) der Ansatz

∫ q(x, t; x_{0}) dx = f(t), f(t) := q(x_{0}, t; x_{0}).  (2.12)

Mit der Dirac’schen δ-Distribution

⟨δ, φ⟩ = φ(0)

läßt sich der Quellterm somit schreiben als

q(x, t; x_{0}) = f(t)δ(x − x_{0}) = f(t)δ_{x_{0}}(x).  (2.13)

Schließlich erhält man, wenn man x_{0} → y_{0} ∈ ∂Ω gehen läßt,

q(x, t; y_{0}) = f(t)δ_{y_{0}}(x), t_{0} ≤ t.  (2.14)

Für ein mit der Kreisfrequenz ω = 2πf harmonisch moduliertes Signal erhält man sodann

q(x, t; y_{0}) = H(t − t_{0})e^{iω(t−t_{0})}δ_{y_{0}}(x).  (2.15)

Dabei steht H(t) für die Heaviside-Distribution,

H(t) = \{ \begin{matrix}1, & 0 ≤ t, \\ 0, & t < 0. \end{matrix}  (2.16)

2.4 Modellgleichung

Stellt man sich Photonen im Teilchenmodell als einen Stoff vor, der durch biologisches Gewebe diffundiert, so können die Ergebnisse der vorherigen Abschnitte zur Modellbildung herangezogen werden.

Sei q(x, t, y_{0}) = e^{iwt}δ_{y_{0}}(x) und sei daher für die Photonendichte Φ(x, t) = u(x)v(t).

Für den Ortsanteil u(x) ergibt sich als resultierende Gleichung, die die Dichteverteilung der Photonen bezüglich einer zeitlich harmonisch modulierten Punktquelle in einem beschränkten Gebiet mit hinreichend glattem Rand beschreibt,

\begin{aligned}−∇_{x} · (κ(x)∇_{x}u) + (μ(x) + \frac{iω}{c} ) u & = & δ_{y_{0}}(x), x ∈ Ω, \\ u(y) + 2A κ(y) \frac{∂u}{∂ν_{y}} & = & 0, y ∈ ∂Ω. \end{aligned}  (2.17)

Für die theoretischen und praktischen Betrachtungen in den folgenden Kapiteln stellt sich heraus, daß eine leicht abgeänderte Modellgleichung einige Vorzüge in der Darstellung der Theorie aufweist. Sei φ(y) eine (meßbare) Verteilung auf dem Rand ∂Ω. Die modifizierte Modellgleichung lautet dann wie folgt,

\begin{aligned}−∇_{x} · (κ(x)∇_{x}u) + (μ(x) + \frac{iω}{c} ) u & = & 0, x ∈ Ω, \\ u(y) & = & φ(y), y ∈ ∂Ω. \end{aligned}  (2.18)

Diese Gleichung soll nun als Grundlage für die Untersuchungen in den folgenden Kapitel dienen.

2.5 DATENMODELLIERUNG

2.5 Datenmodellierung

Die Photonen, die das Gebiet Ω durch den Rand ∂Ω verlassen, können gemessen werden. Da die Amplitude und die Phase des Eingangssignals bekannt ist, kann somit für die optische Tomographie im Frequenzbereich die Phasen- und Amplitudendifferenz als Daten eines realen Experiments gemessen werden. Für die Modellgleichung (2.17) bzw. (2.18) läßt sich anhand der Lösung u(x) der Differentialgleichung wie folgt berechnen. Sei ν = ν(y) das äußere Normalenfeld an ∂Ω. Die Amplitude r und die Phase ϕ des gemessenen Signals wird dann durch die in Polarkoordinaten ausgedrückte komplexwertige Größe g = r e^{iϕ} der Richtungsableitung

g = (κ ∂_{ν}u)|_{∂Ω}  (2.19)

beschrieben.

KAPITEL 2. PHYSIKALISCHE PROBLEMSTELLUNG 14

Mathematischer Teil

In diesem Kapitel wird zunächst das direkte Problem der optischen Tomographie analysiert. Ein wesentliches Ergebnis wird die Fréchet–Differenzierbarkeit des Datenfunktionals sein, welches die physikalischen Gegebenheiten der Messungen in der optischen Tomographie mathematisch beschreibt.

Daran anschließend werden die Ergebnisse dieser Analyse dazu verwendet, einen iterativen Algorithmus zur Lösung des inversen Problems herzuleiten. Dabei wird insbesondere der Umstand berücksichtigt, daß das inverse Problem sehr schlecht gestellt ist, und also entsprechend eine Regularisierung zur Lösung des inversen Problems vorgenommen werden muß.

Desweiteren ist für das angegebene Verfahren genau darauf zu achten, daß die im Frequenzbereich gemessenen Daten als komplexwertige Größen vorliegen, wohingegen die gesuchten optischen Parameter als Urbild der gegebenen Daten rein reelle Parameter sind. Interessanterweise ergibt sich durch die Ausnützung dieses Umstandes eine sehr effiziente Methode zur Rekonstruktion der beiden optischen Parameter aus der gemessenen Phasen- und Amplitudendifferenz am Rand des zu untersuchenden Gebietes als Daten für das inverse Problem.

3.1 Direktes Problem

Es sei L_{f} : U(Ω) ⟶ Q(Ω), u ↦ L_{f}u ein linearer Differentialoperator abhängig von Parameterfunktionen f = (f_{1}, …, f_{m}) ∈ F_{1}(Ω) × … × F_{m}(Ω). Die Mengen U(Ω), Q(Ω) und F_{j}(Ω), j = 1, …, m, seien geeignete separable Banachräume von Funktionen über einem beschränkten Gebiet Ω ⊂ IR^{n} mit einem (n − 1)-dimensionalen Rand ∂Ω. Schließlich sei φ ∈ R(∂Ω) eine Funktion oder allgemeiner eine Verteilung, die auf dem Rand ∂Ω definiert ist.

Dann versteht man unter einem direkten Randwertproblem die folgende Aufgabe.

KAPITEL 3. MATHEMATISCHER TEIL

Definition 3.1.1 (Direktes Randwertproblem)

Bestimme bei gegebenen Parametern f ∈ F_{1}(Ω) × … × F_{m}(Ω) und φ ∈ R(∂Ω) eine Lösung u_{f} ∈ U(Ω) der Gleichung

L_{f}u = q  in Ω,  u = φ  auf ∂Ω,  (3.1)

wobei q ∈ Q(Ω) ein Term unabhängig von L_{f}, f und u ist.

Es stellt sich natürlich sofort die Frage, ob es überhaupt solche Funktionen u_{f} ∈ U(Ω) gibt, die die Gleichung (3.1) erfüllen. Ist dies der Fall, so kann man weiter fragen, ob es nur genau ein u_{f} ∈ U(Ω) als Lösung für (3.1) gibt.

Diese Fragen nach der Existenz und Eindeutigkeit einer Lösung von (3.1) sollen für das Modell der Lasertomographie beantwortet werden, also für (3.1) mit dem Differentialoperator

L_{f}u ≡ [−∇_{x} · (f_{1}(x)∇_{x}) + (f_{2}(x) + iω_{0})]u,  (3.2)

wobei f_{1}(x) := κ(x), f_{2}(x) := μ_{a}(x), ω_{0} := ω/c und q := 0 gesetzt sind. Für die Beantwortung dieser Fragen kommt es darauf an, geeignete Funktionenräume zu wählen, so daß der Ausdruck L_{f}u überhaupt Sinn macht.

Aus den Beweisen der entsprechenden Sätze über die Existenz und Eindeutigkeit ergeben sich dann auch Aussagen über die Art der Abhängigkeit der Lösungen u_{f} von den Parametern f.

Mit dieser Analyse des direkten Problems erhält man sodann entsprechende Aussagen über die Abbildung

P_{f}(φ) = (f_{1} ∂_{ν}u_{f} )|_{∂Ω},  (3.3)

die einer gegebenen Randverteilung φ auf ∂Ω eine durch f_{1} skalierte Normalenableitung der Lösung u_{f} von (3.1) bezüglich q = q(φ) zuordnet. Mit P_{f} werden also die gemessenen Daten modelliert.

Es stellt sich heraus, daß diese Datenfunktion P_{f} stetig und differenzierbar in den Parametern f ist, wenn sie im schwachen Sinne, d.h. als ein Datenfunktional, definiert wird.

3.1.1 Definitionen

Im folgenden sei stets Ω ⊂ IR^{n} ein beschränktes Gebiet mit glattem Rand oder Ω = IR^{n} der ganze euklidische Raum, versehen mit dem kanonischen Lebesgue–Maß.

Es sei für α = (α_{1}, …, α_{n}) ∈ IN^{n} der Betrag |α| := ∑ α_{j} und die Schreibweise

D^{α}f := \frac{∂^{|α|}f}{∂x_{1}^{α_{1}} ⋯ ∂x_{n}^{α_{n}}}  (3.4)

eingeführt.

Alle Funktionen seien zunächst komplexwertige Funktionen. Die folgenden Definitionen der jeweiligen Funktionenräume seien also komplexwertige Vektorräume.

Sofern reellwertige Funktionen betrachtet werden sollen, wird dies explizit angegeben.

3.1 DIREKTES PROBLEM

Die Räume C_{0}^{∞}(Ω) und D′(Ω).

Für Ω ⊂ IR^{n} sei mit C_{0}^{∞}(Ω), oder auch mit D(Ω), der Raum aller komplexwertigen C^{∞}–Funktionen f(x) auf IR^{n} mit kompaktem Träger supp f ⊂ Ω bezeichnet. Insbesondere verschwinden alle in C_{0}^{∞}(Ω) enthaltenen Funktionen außerhalb einer kompakten Teilmenge von Ω.

Konvergenz in C_{0}^{∞}(Ω) sei wie folgt definiert: die Funktionen f_{j}(x) ∈ C_{0}^{∞}(Ω), j ∈ IN, konvergieren genau dann gegen eine Funktion f(x) ∈ C_{0}^{∞}(Ω), wenn sämtliche Träger der Funktionen f_{j}(x) in einer kompakten Teilmenge Ω′ ⊂ Ω enthalten sind und D^{α}f_{j}(x) für alle α ∈ IN^{n} gleichmäßig gegen D^{α}f(x) konvergiert. Damit sei ein lineares Funktional auf C_{0}^{∞}(Ω), ℓ : C_{0}^{∞}(Ω) → C, f ↦ ⟨ℓ, f⟩ := ℓ(f), genau dann stetig, wenn stets aus der erklärten Konvergenz f_{j} → f in C_{0}^{∞}(Ω) die Konvergenz ⟨ℓ, f_{j}⟩ → ⟨ℓ, f⟩ in C folgt.

Mit D′(Ω) sei dann der Dualraum

D′(Ω) := { ℓ : C_{0}^{∞}(Ω) → C, ℓ linear und stetig in C_{0}^{∞}(Ω) }  (3.5)

bezeichnet.

Die Räume L^{p}(Ω), L_{loc}^{p}(Ω) und L^{∞}(Ω).

Für 1 ≤ p < ∞ sei mit L^{p}(Ω) der Raum aller meßbaren, komplexwertigen Funktionen f(x) auf Ω bezeichnet, für die

∫_{Ω} |f(x)|^{p} dx < ∞  (3.6)

ist. Mit der Norm

∥f∥_{L^{p}} := (∫_{Ω} |f(x)|^{p} dx)^{1/p}  (3.7)

ist L^{p}(Ω) ein Banachraum.

Für 1 ≤ p < ∞ sei entsprechend der lokale Raum

L_{loc}^{p}(Ω) := { f : Ω → C meßbar, ∀ Ω′ ⊂⊂ Ω : ∫_{Ω′} |f(x)|^{p} dx < ∞ }  (3.8)

definiert. Man beachte, daß L_{loc}^{p}(Ω) kein normierter Raum ist.

Mit L^{∞}(Ω) sei der Raum aller fast überall beschränkten, komplexwertigen Funktionen f(x) auf Ω bezeichnet, für die das essentielle Supremum ess sup_{Ω} |f(x)| < ∞ ist.

Zusammen mit der Norm

∥f∥_{L^{∞}} := ess sup_{x∈Ω} |f(x)|  (3.9)

wird L^{∞}(Ω) ein Banachraum.

Die Räume W^{s,p}(Ω) (H^{s}(Ω)), W_{0}^{s,p}(Ω) (H_{0}^{s}(Ω)), und W^{s,p}(∂Ω) (H^{s}(∂Ω)).

Es sollen die Sobolevräume definiert werden, die zur Formulierung der Aussagen in diesem Kapitel benötigt werden. Dazu sei zunächst an den Begriff der schwachen Ableitung erinnert. Mit

W^{k}(Ω) := { f ∈ L_{loc}^{1}(Ω), ∀ |α| ≤ k ∃ g_{α} ∈ L_{loc}^{1}(Ω) :  (3.10)

∫_{Ω} (g_{α} φ) dx = (−1)^{α} ∫_{Ω} (f D^{α}φ) dx für alle φ ∈ C_{0}^{∞}(Ω) }

sei der Raum der k-mal schwach differenzierbaren Funktionen definiert. Die schwache Ableitung D^{α}f von f ∈ W^{k}(Ω) ist dann durch D^{α}f := g_{α} ∈ L_{loc}^{1}(Ω) gegeben.

Für ganzzahlige k ≥ 0 sei mit W^{k,p}(Ω) der Sobolevraum

W^{k,p}(Ω) := { f ∈ W^{k}(Ω), ∀ |α| ≤ k : D^{α}f ∈ L^{p}(Ω) }  (3.11)

bezeichnet. Zusammen mit der Norm

∥f∥_{W^{k,p}} := ( ∑_{|α|≤k} ∫_{Ω} |D^{α}f|^{p} dx)^{1/p}  (3.12)

ist der Raum W^{k,p}(Ω) ein Banachraum. Für nichtganzzahlige s > 0 sei s = k + σ die Zerlegung von s in einen ganzzahligen Anteil k und einen nichtganzzahligen Anteil σ mit 0 < σ < 1. Mit W^{s,p}(Ω) sei dann der Raum

W^{s,p}(Ω) := \{ f ∈ W^{k,p}(Ω), ∀ |α| = k : ∫_{Ω} ∫_{Ω} \frac{|D^{α}f(x) − D^{α}f(y)|^{p}}{|x − y|^{n+σp}} dx dy < ∞ \}

(3.13)

bezeichnet, der durch die Norm

∥f∥_{W^{s,p}} := (∥f∥_{W^{k,p}}^{p} + ∑_{|α|=k} ∫_{Ω} ∫_{Ω} \frac{|D^{α}f(x) − D^{α}f(y)|^{p}}{|x − y|^{n+σp}} dx dy)^{1/p}  (3.14)

zu einem Banachraum wird.

Der Fall p = 2 sei für s > 0 beliebig mit H^{s}(Ω) := W^{s,2}(Ω) bezeichnet, denn mit p = 2 ist W^{s,p}(Ω) ein Hilbertraum bezüglich des skalaren Produktes für ganzzahliges s = k

⟨u | v⟩_{k} := ∑_{|α|≤k} ∫_{Ω} (D^{α}u \overline{D^{α}v}) dx  (3.15)

bezeihungsweise für nichtganzzahliges s = k + σ, 0 < σ < 1,

⟨u | v⟩_{s} := ⟨u | v⟩_{k} + ∑_{|α|=k} ∫_{Ω} ∫_{Ω} \frac{(D^{α}u(x) − D^{α}u(y))(\overline{D^{α}v(x) − D^{α}v(y)})}{|x − y|^{n+2σ}} dx dy. (3.16)

Desweiteren sei für s > 0 mit W_{0}^{s,p}(Ω) := \overline{C_{0}^{∞}(Ω)}^{∥·∥_{W^{s,p}}} der Abschluß von C_{0}^{∞}(Ω) bezüglich der Norm ∥ · ∥_{W^{s,p}} bezeichnet. Dies macht Sinn, da für s > 0 die Menge C_{0}^{∞}(Ω) nicht dicht in W^{s,p}(Ω) liegt. Der Fall p = 2 sei ebenfalls mit H_{0}^{s}(Ω) := W_{0}^{s,2}(Ω) bezeichnet. Zusammen mit dem Skalarprodukt (3.15) bzw. (3.16) ist der Raum H_{0}^{s}(Ω) ein Hilbertraum.

Für s > 0 sei mit W^{−s,p}(Ω) := [W_{0}^{s,q}(Ω)]′, p^{−1} + q^{−1} = 1, der Dualraum von W_{0}^{s,q}(Ω) bezeichnet, d.h. W^{−s,p}(Ω) ist der Raum aller linearen Funktionale ℓ : W_{0}^{s,q}(Ω) → C, die bezüglich der Norm ∥ · ∥_{W^{s,q}} stetig sind. Mit H^{−s}(Ω) := W^{−s,2}(Ω) sei wie oben der Hilbertraumfall p = 2 abgekürzt bezeichnet.

Schließlich sei Γ := ∂Ω der Rand des beschränkten Gebietes Ω ⊂ IR^{n} und es sei Γ_{0} ⊂ Γ eine offene Teilmenge von Γ. Der Rand ∂Ω sei vom Typ C^{k,1}, d.h. es gebe für jeden Punkt x ∈ ∂Ω eine Umgebung U := U(x) ⊂ IR^{n} und ein neues orthogonales Koordinatensystem { y_{1}, …, y_{n} }, so daß

1. die Menge U in den neuen Koordinaten ein Quader ist,

U = { (y_{1}, …, y_{n}), −a_{j} < y_{j} < a_{j} für 1 ≤ j ≤ n };

2. es auf der Menge U′ = { y′ = (y_{1}, …, v_{n−1}), −a_{j} < y_{j} < a_{j} für 1 ≤ j ≤ n } eine C^{k,1}–Funktion ϕ(y′) gibt, für die |ϕ(y′)| ≤ a_{n}/2 für alle y′ ∈ U′ und

Ω ∩ U = { y = (y′, y_{n}) ∈ U, y_{n} < ϕ(y′) },
Γ ∩ U = { y = (y′, y_{n}) ∈ U, y_{n} = ϕ(y′) }

gilt.

Es sei Φ(y) := {y_{1}, …, y_{n−1}, ϕ(y_{1}, …, y_{n−1}) } definiert. Damit ist Γ als eine (n − 1)- dimensionale Untermannigfaltigkeit des IR^{n} vom Typ C^{k,1} erklärt.

Nun bezeichne V die Menge aller Paare (U, ϕ) bezüglich Γ_{0}. Für |s| ≤ k + 1 sei dann mit W^{s,p}(Γ_{0}) der Banachraum

W^{s,p}(Γ_{0}) := { ℓ ∈ D′(Γ_{0}), ∀ (U, ϕ) ∈ V : ℓ ∘ Φ ∈ W^{s,p}(U′ ∩ Φ^{−1}(Γ_{0} ∩ U)) } (3.17)

mit der Norm

∥f∥_{W^{s,p}(Γ_{0})} := (∫_{Γ_{0}} |f|^{p} dσ + ∫_{Γ_{0}}∫_{Γ_{0}} \frac{|f(x) − f(y)|^{p}}{|x − y|^{n−1+sp}} dσ_{x} dσ_{y})^{1/p}  (3.18)

im Fall 0 < s < 1 und der Norm ohne das Doppelintegral analog zu (3.12) im ganzzahligen Fall s = k definiert. Die übrigen Fälle sind gemäß der Bildungsweise (3.14) zu verstehen. Erneut sei für den Fall p = 2 die Bezeichnung H^{s}(Γ_{0}) := W^{s,2}(Γ_{0}) eingeführt.

Eine Funktion f(x) ∈ L_{loc}^{1}(Γ_{0}) sei durch die Injektion f ↪ T_{f} ∈ D′(Γ_{0}),

⟨T_{f}, g⟩ := ∫_{Γ_{0}} (f g) dσ  (3.19)

als Element aus D′(Γ_{0}) erklärt, wobei dσ das Oberflächenmaß auf Γ ist.

3.1.2 Existenz und Eindeutigkeit einer Lösung

Ein wesentliches Hilfsmittel für die Existenzaussage schwacher Lösungen bei elliptischen Differentialgleichungen ist das folgende Lemma.

Lemma 3.1.1 (Lax–Milgram) Es sei H ein Hilbertraum und

a : H × H ⟶ C  (3.20)

eine Sesquilinearform. Weiter gebe es zwei positive Konstanten c und C, so daß

|a[x, y]| ≤ C ∥x∥_{H} ∥y∥_{H}  (3.21)

für alle x, y ∈ H und

Re a[x, x] ≥ c ∥x∥_{H}^{2}  (3.22)

für alle x ∈ H gilt.

Dann existiert für jedes lineare Funktional F ∈ H^{∗} ein eindeutig bestimmtes Element y_{0} ∈ H, so daß die Gleichung

a[x, y_{0}] = F (x)  (3.23)

für alle x ∈ H erfüllt ist. Darüberhinaus gilt für dieses y_{0} ∈ H die Abschätzung

∥y_{0}∥_{H} ≤ \frac{1}{c} ∥F ∥_{H^{∗}}.  (3.24)

Beweis: Zu finden in [70, III.7].

□

Nun sei Ω ⊂ IR^{n} ein Gebiet des Typs C^{1,1}. Der Rand ∂Ω habe das Oberflächenmaß dσ und es sei φ ∈ H^{1/2}(∂Ω). Eine Funktion u ∈ H^{1}(Ω) nehme genau dann die Randwerte φ an, wenn

∫_{∂Ω} (u w) dσ = ∫_{∂Ω} (φ w) dσ  ∀ w ∈ L^{2}(∂Ω)  (3.25)

gilt.

Für reellwertige Funktionen f_{1}, f_{2} ∈ L^{∞}(Ω : IR) gelte inf_{Ω} f_{1}(x) > 0, inf_{Ω} f_{2}(x) > 0. Diese Eigenschaft sei im weiteren mit f_{1} >> 0 und f_{2} >> 0 abgekürzt.

Es folgt, daß für die Funktionen f_{1}, f_{2} ∈ L^{∞}(Ω : IR) mit f_{1} >> 0 und f_{2} >> 0 der reellwertige Differentialoperator

L̃_{f} ≡ [−∇_{x} · (f_{1}(x)∇_{x}) + f_{2}(x)]  (3.26)

in Ω elliptisch ist, und es gilt, daß es dann genau eine schwache Lösung u ∈ H^{1}(Ω) der Gleichung

L̃_{f}u = 0,  u|_{∂Ω} = φ  (3.27)

gibt, siehe [18, Thm. 8.3]. Mit dem Variationsprinzip in H^{1}(Ω),

∫_{Ω} (f_{1}|∇u|^{2} + f_{2}|u|^{2}) dx = inf_{v|_{∂Ω}=φ} ∫_{Ω} (f_{1}|∇v|^{2} + f_{2}|v|^{2}) dx  (3.28)

(siehe z.B. [17, Chap. 2F] mit entsprechenden Modifikationen), und der Eigenschaft, daß es unabhängig von φ ∈ H^{1/2}(∂Ω) und v ∈ H^{1}(Ω) Konstanten C_{∂Ω}, C′_{∂Ω} > 0 gibt, so daß

C′_{∂Ω} ∥φ∥_{H^{1/2}} ≤ inf_{v|_{∂Ω}=φ} ∥v∥_{H^{1}} ≤ C_{∂Ω} ∥φ∥_{H^{1/2}}  (3.29)

gilt (siehe [38, Chap. I, Thm. 9.4]), erhält man aus (3.28) eine Abschätzung für die Lösung u ∈ H^{1}(Ω) der Gleichung (3.27),

∥u∥_{H^{1}} ≤ C_{∂Ω} \frac{max\{∥f_{1}∥_{∞}, ∥f_{2}∥_{∞}\}}{min\{inf_{Ω} f_{1}, inf_{Ω} f_{2}\}} ∥φ∥_{H^{1/2}}.  (3.30)

Analoge Aussagen sollen nun für den mit dem komplexen Anteil i ω_{0} gestörten Differentialoperator L_{f} ≡ L̃_{f} +i ω_{0}, der im Model (2.17) bzw. (2.18) der optischen Tomographie auftritt, nachgewiesen werden. Man hat den

Satz 3.1.1 Es sei L_{f} ≡ L̃_{f} +i ω_{0} mit den genannten Voraussetzungen an f und Ω ⊂ IR^{n} ein beschränktes Gebiet vom Typ C^{1,1}.

Dann gilt:

1. Die Gleichung

L_{f}u = 0 in Ω, u|_{∂Ω} = φ ∈ H^{1/2}(∂Ω),  (3.31)

besitzt eine eindeutige Lösung u ∈ H^{1}(Ω).

2. Mit der Konstanten C_{∂Ω} > 0 aus (3.29) gilt für die Lösung u ∈ H^{1}(Ω) die Abschätzung

∥u∥_{H^{1}} ≤ C_{∂Ω} (\frac{max \{∥f_{1}∥_{∞}, ∥f_{2}∥_{∞}, ω_{0}\}}{min \{ inf_{Ω} f_{1}, inf_{Ω} f_{2}\}} + 1) ∥φ∥_{H^{1/2}}.  (3.32)

Beweis: Zunächst gilt für die zu L_{f} = L̃_{f} + i ω_{0} gehörige Bilinearform

a(u, v) ≡ ∫_{Ω} (f_{1}(∇u · \overline{∇v}) + (f_{2} + i ω_{0}) (u \overline{v})) dx  (3.33)

für alle v ∈ H_{0}^{1}(Ω) die Gårdingsche Ungleichung,

Re a(v, v) = ∫_{Ω} (f_{1}|∇v|^{2} + f_{2}|v|^{2}) dx

≥ ( inf_{Ω} f_{1}) ∫_{Ω} |∇v|^{2} dx + ( inf_{Ω} f_{2}) ∫_{Ω} |v|^{2} dx
≥ ( min { inf_{Ω} f_{1}, inf_{Ω} f_{2}}) ∥v∥_{H^{1}}^{2} − λ_{G}∥v∥_{L^{2}}^{2}, λ_{G} := 0. (3.34)

Da sich mit der Cauchy-Schwarzschen Ungleichung die Abschätzung

|a(u, v)| ≤ ∥f_{1}∥_{∞} ∫_{Ω} |∇u| |∇v| dx + ∥f_{2} + i ω_{0}∥_{∞} ∫_{Ω} |u| |v| dx
≤ ( max {∥f_{1}∥_{∞}, ∥f_{2}∥_{∞}, ω_{0}}) ∥u∥_{H^{1}} ∥v∥_{H^{1}}  (3.35)

ergibt, sind die Voraussetzungen des Lax–Milgram–Lemmas 3.1.1 bezüglich des Hilbertraumes H_{0}^{1}(Ω) und der Bilinearform a(u, v) erfüllt.

Nun sei v ∈ H^{1}(Ω) mit v|_{∂Ω} = φ. Damit sei für ϕ ∈ C_{0}^{∞}(Ω) das lineare Funktional

G_{v}(ϕ) := − ∫_{Ω} (f_{1} (∇v · \overline{∇ϕ}) + (f_{2} + i ω_{0}) (v \overline{ϕ})) dx  (3.36)

definiert. Die stetige Fortsetzung von G_{v} auf H_{0}^{1}(Ω) sei wieder mit G_{v} bezeichnet. Dann folgt anhand des Lax–Milgram–Lemmas 3.1.1, daß die Gleichung

a(w, ṽ) = G_{v}(ṽ), für alle ṽ ∈ H_{0}^{1}(Ω)  (3.37)

eine eindeutige Lösung w ∈ H_{0}^{1}(Ω) besitzt. Für die Lösung w ∈ H_{0}^{1}(Ω) gilt die Abschätzung

∥w∥_{H^{1}} ≤ \frac{1}{min\{inf_{Ω} f_{1}, inf_{Ω} f_{2}\}} ∥G_{v}∥_{H^{−1}} ≤ \frac{max\{∥f_{1}∥_{∞}, ∥f_{2}∥_{∞}, ω_{0}\}}{min\{inf_{Ω} f_{1}, inf_{Ω} f_{2}\}} ∥v∥_{H^{1}}. (3.38)

Daraus folgt, daß die Funktion u := w + v ∈ H^{1}(Ω) die Gleichung (3.31) im schwachen Sinne löst; und es gilt die Ungleichung

∥u∥_{H^{1}} ≤ (\frac{max \{∥f_{1}∥_{∞}, ∥f_{2}∥_{∞}, ω_{0}\}}{min \{ inf_{Ω} f_{1}, inf_{Ω} f_{2}\}} + 1) ∥v∥_{H^{1}}  (3.39)

für alle v ∈ H^{1}(Ω) mit v|_{∂Ω} = φ. Mit der Konstanten C_{∂Ω} > 0 aus (3.29) folgt aus (3.39) und (3.29) die Abschätzung (3.32). Damit ist der Satz bewiesen.

□

Die Abschätzung (3.32) ist für die folgenden Abschnitte wesentlich.

Schließlich ermöglicht die Existenz einer Lösung u ∈ H^{1}(Ω) der Gleichung (3.31) die Definition eines linearen Funktionals ℓ_{φ} : H^{1/2}(∂Ω) ⟶ C,

⟨ℓ_{φ}, ψ⟩ := a(u, \overline{v}_{ψ}) = ∫_{Ω} (f_{1}(∇u · ∇v_{ψ}) + (f_{2} + i ω_{0}) (u v_{ψ})) dx,  (3.40)

wobei v_{ψ} ∈ H^{1}(Ω) die Lösung der Gleichung L_{f}v = 0, v|_{∂Ω} = ψ, ist. Wegen (3.35) und

(3.29) ist ℓ_{φ} auch stetig in H^{1/2}(∂Ω). Folglich ist ℓ_{φ} ∈ H^{−1/2}(∂Ω).

Eine partielle Integration in (3.40), wobei alle auftretenen Ableitungen im schwachen Sinn zu nehmen sind, ergibt nun

⟨ℓ_{φ}, ψ⟩ = ∫_{∂Ω} (ψ f_{1}∂_{ν}u) dσ + ∫_{Ω} (v_{ψ} \underbrace{L_{f}u}_{=0} ) dx,  (3.41)

so daß die sogenannte Dirichlet-zu-Neumann–Abbildung (engl.: DtN-map, siehe z.B.

[32] [46]),

P_{f} : H^{1/2}(∂Ω) ⟶ H^{−1/2}(∂Ω),  P_{f}(φ) = (f_{1} ∂_{ν}u)|_{∂Ω}  (3.42)

im Zusammenhang mit der optischen Tomographie im Frequenzbereich erklärt werden kann. Unter dem Bild der Zuordnung (3.42) sei also wegen (3.41) das stetige, lineare Funktional (3.40) verstanden.

3.1.3 Stetigkeit der Dirichlet-zu-Neumann–Abbildung

Für diesen und den nächsten Abschnitt sei an folgende Definitionen erinnert.

Definition 3.1.2 (stetig, Fréchet–differenzierbar) Es seien X und Y normierte Vektorräume über IR oder C, D ⊂ X eine offene Untermenge, x̃ ∈ D, und T : D ⊂ X ⟶ Y eine (nichtlineare) Abbildung.

1. T heißt stetig in x̃ ∈ D, falls es für jede Zahl ε > 0 eine Zahl δ > 0 gibt, so daß ∥T (x) − T (x̃)∥_{Y} ≤ ε für alle x ∈ D mit ∥x − x̃∥_{X} ≤ δ gilt.

2. T heißt Fréchet–differenzierbar in x̃ ∈ D, falls es einen linearen, beschränkten Operator Ax̃ : X ⟶ Y (abhängig von x̃) gibt, so daß

lim_{h→0} \frac{1}{∥h∥_{X}} ∥T (x̃ + h) − T (x̃) − Ax̃h∥_{Y} = 0  (3.43)

gilt. Man setzt dann T′(x̃) := Ax̃; es gilt somit T′(x̃) ∈ L(X, Y), der Raum der linearen, beschränkten Operatoren X → Y.

3. T stetig Fréchet–differenzierbar in x̃ ∈ D, falls T Fréchet–differenzierbar in einer Umgebung V (x̃) ist und die Abbildung T′ : V ⟶ L(X, Y) stetig in x̃ ist.

Zum Nachweis der Stetigkeit der Abbildung (3.42) in den Parametern f sei für den normierten Raum Y der Definition 3.1.2 der komplexe Raum Y := L(H^{1/2}(Ω), H^{−1/2}(Ω)) mit der Operatornorm ∥ · ∥_{1/2,−1/2} := ∥ · ∥_{H^{1/2},H^{−1/2}} und für den normierten Raum X der reelle Vektorraum X := (L^{∞}(Ω : IR) × L^{∞}(Ω : IR)) mit der Norm ∥f∥_{L^{∞}×L^{∞}} := max{ ∥f_{1}∥_{∞}, ∥f_{2}∥_{∞} } gewählt. Damit läßt sich folgender Satz formulieren.

Satz 3.1.2 Es seien f, f̃ ∈ (L^{∞}(Ω : IR) × L^{∞}(Ω : IR)) zwei reellwertige Parameterpaare, für die es eine Zahl γ > 0 gebe, so daß

0 < \frac{1}{γ} ≤ f_{1,2}(x) ≤ γ und 0 < \frac{1}{γ} ≤ f̃_{1,2}(x) ≤ γ  (3.44)

gilt. Dann gibt es eine Konstante C_{∂Ω,γ,ω_{0}} > 0, so daß gilt

∥P_{f} − P_{f}∥_{1/2,−1/2} ≤ C_{∂Ω,γ,ω_{0}} ∥f − f̃∥_{L^{∞}×L^{∞}}.  (3.45)

Beweis: Für die Differenz (u − ũ) der beiden Lösungen u, ũ ∈ H^{1}(Ω) der Gleichungen

L_{f}u = 0, u|_{∂Ω} = φ, Lf̃ũ = 0, ũ|_{∂Ω} = φ, (3.46)

ergibt sich, daß sie die Gleichung

L_{f}(u − ũ) = −(L_{f} − Lf̃)ũ = −L̃_{(f−f̃)}ũ, (u − ũ)|_{∂Ω} = 0  (3.47)

erfüllt. Definiert man

q(z) := − ∫_{Ω} ((f_{1} − f̃_{1})(∇ũ · ∇z) + (f_{2} − f̃_{2})(ũ z)) dx,  (3.48)

so ist q ∈ H^{−1}(Ω) mit

∥q∥_{H^{−1}} ≤ ( max {∥f_{1} − f̃_{1}∥_{∞}, ∥f_{2} − f̃_{2}∥_{∞}}) ∥ũ∥_{H^{1}}.  (3.49)

Anhand des Lax–Milgram–Lemmas 3.1.1 folgt, daß es genau eine schwache Lösung w ∈ H_{0}^{1}(Ω) der Gleichung

L_{f}w = q, w|_{∂Ω} = 0 (3.50)

gibt, die einer Abschätzung

∥w∥_{H^{1}} ≤ \frac{1}{min \{ inf_{Ω} f_{1}, inf_{Ω} f_{2}\}} ∥q∥_{H^{−1}}  (3.51)

genügt. Somit ist (u − ũ) = w ∈ H_{0}^{1}(Ω) und

∥u − ũ∥_{H^{1}} ≤ \frac{max \{∥f_{1} − f̃_{1}∥_{∞}, ∥f_{2} − f̃_{2}∥_{∞}\}}{min \{ inf_{Ω} f_{1}, inf_{Ω} f_{2}\}} ∥ũ∥_{H^{1}}.  (3.52)

Damit erhält man

sup_{∥ψ∥_{1/2}=1} |⟨(P_{f} − Pf̃)(φ), ψ⟩|

= sup_{∥ψ∥_{1/2}=1} | ∫_{Ω} ((f_{1} − f̃_{1})(∇u · ∇v) + f̃_{1}(∇u · ∇v − ∇ũ · ∇ṽ)

+ (f_{2} − f̃_{2})(u v) + (f̃_{2} + i ω_{0})(u v − ũ ṽ)) dx|
≤ ( max {∥f_{1} − f̃_{1}∥_{∞}, ∥f_{2} − f̃_{2}∥_{∞}} ∥u∥_{H^{1}}) sup_{∥ψ∥_{1/2}=1} ∥v∥_{H^{1}}
+ C_{f̃,ω_{0}} (∥u − ũ∥_{H^{1}} sup_{∥ψ∥_{1/2}=1} ∥v∥_{H^{1}} + ∥ũ∥_{H^{1}} sup_{∥ψ∥_{1/2}=1} ∥v − ṽ∥_{H^{1}})  (3.53)
≤ C_{∂Ω, γ, ω_{0}}( max {∥f_{1} − f̃_{1}∥_{∞}, ∥f_{2} − f̃_{2}∥_{∞}}) ∥φ∥_{H^{1/2}},  (3.54)

wobei die Konstante C_{f̃, ω_{0}} > 0 in (3.53) durch C_{f̃,ω_{0}} = max{ ∥f̃_{1}∥_{∞}, ∥f̃_{2}∥_{∞}, ω_{0} } gegeben ist und die Konstante C_{∂Ω,γ,ω_{0}} > 0 in (3.54) anhand (3.52) und (3.32) aus (3.53) hervorgeht. Aus der Abschätzung (3.54) folgt (3.45), was zu zeigen war.

□

3.1.4 Differenzierbarkeit der Dirichlet-zu-Neumann–Abbildung

Betrachte nun für ein reellwertiges Parameterpaar f ∈ (L^{∞}(Ω : IR) × L^{∞}(Ω : IR)) mit f_{1,2} >> 0 die lineare Abbildung

(A_{f} ·)(φ) : L^{∞}(Ω) × L^{∞}(Ω) ⟶ H^{−1/2}(∂Ω), (A_{f}h)(φ) = (f_{1} ∂_{ν}w)|_{∂Ω}, (3.55)

wobei w ∈ H_{0}^{1}(Ω) die Lösung der Gleichung

L_{f}w = −L̃_{h}u_{f} in Ω, w|_{∂Ω} = 0,  (3.56)

und u_{f} = u ∈ H^{1}(Ω) die Lösung der Gleichung (3.31) bezüglich f ist.

Damit läßt sich mit einer Beweistechnik aus [15] der folgende Satz beweisen.

Satz 3.1.3 Es seien die Voraussetzungen des Satzes 3.1.2 erfüllt.

Dann gibt es eine Konstante C′_{∂Ω,γ,ω_{0}} > 0, so daß die Abschätzung

∥P_{f} − P_{f}− A_{f}h ∥_{1/2,−1/2} ≤ C′_{∂Ω,γ,ω_{0}} ∥h∥_{L^{∞}×L^{∞}}^{2}  (3.57)

für h = f − f̃ gültig ist.

KAPITEL 3. MATHEMATISCHER TEIL

Beweis: Die Differenz (u − ũ − w) erfüllt die Gleichung

Lf̃(u − ũ − w) = −L̃_{h}(u − ũ), (u − ũ − w)|_{∂Ω} = 0, (3.58)

also ist (u − ũ − w) ∈ H_{0}^{1}(Ω) eindeutig bestimmt und es gilt die Abschätzung

∥u − ũ − w∥_{H^{1}} ≤ \frac{max \{∥h_{1}∥_{∞}, ∥h_{2}∥_{∞}\}}{min \{ inf_{Ω} f̃_{1}, inf_{Ω} f̃_{2}\}} ∥u − ũ∥_{H^{1}}
≤ (\frac{max \{∥h_{1}∥_{∞}, ∥h_{2}∥_{∞}\}}{min \{ inf_{Ω} f̃_{1}, inf_{Ω} f̃_{2}\}})^{2} ∥u∥_{H^{1}},  (3.59)

da die Differenz (u − ũ) ∈ H_{0}^{1}(Ω) der Gleichung

Lf̃(u − ũ) = −L̃_{h}u, (u − ũ)|_{∂Ω} = 0, (3.60)

genügt. Damit erhält man

sup_{∥ψ∥_{1/2}=1} |⟨(P_{f} − Pf̃ − Af̃h)(φ), ψ⟩|

= sup_{∥ψ∥_{1/2}=1} | ∫_{Ω} (f_{1}(∇u · ∇v) − f̃_{1}(∇ũ · ∇ṽ)
− h_{1}(∇ũ · ∇ṽ) − f̃_{1}(∇w · ∇ṽ + ∇ũ · ∇w)
+ (f_{2} + i ω_{0})(u v) − (f̃_{2} + i ω_{0})(ũ ṽ)
− (h_{2} + i ω_{0})(ũ ṽ) − (f̃_{2} + i ω_{0})(w ṽ + ũ w)) dx|
= sup_{∥ψ∥_{1/2}=1} | ∫_{Ω} (h_{1}(∇u · ∇v − ∇ũ · ∇ṽ) + h_{2}(u v − ũ ṽ)

+ f̃_{1}[∇(u − ũ − w) · ∇v + ∇w · ∇(v − ṽ) + ∇ũ · ∇(v − ṽ − w)]

+ (f̃_{2} + i ω_{0})[(u − ũ − w) v + w (v − ṽ) + ũ (v − ṽ − w)]) dx|

≤ ( max {∥h_{1}∥_{∞}, ∥h_{2}∥_{∞}}) ×

(∥u − ũ∥_{H^{1}} sup_{∥ψ∥_{1/2}=1} ∥v∥_{H^{1}} + ∥ũ∥_{H^{1}} sup_{∥ψ∥_{1/2}=1} ∥v − ṽ∥_{H^{1}})

+ C_{f̃, ω_{0}}(∥u − ũ − w∥_{H^{1}} sup_{∥ψ∥_{1/2}=1} ∥v∥_{H^{1}} + ∥w∥_{H^{1}} sup_{∥ψ∥_{1/2}=1} ∥v − ṽ∥_{H^{1}}

+ ∥ũ∥_{H^{1}} sup_{∥ψ∥_{1/2}=1} ∥v − ṽ − w∥_{H^{1}})  (3.61)
≤ ( max {∥h_{1}∥_{∞}, ∥h_{2}∥_{∞}}) ×

(\frac{max \{∥h_{1}∥_{∞}, ∥h_{2}∥_{∞}\}}{min \{ inf_{Ω} f̃_{1}, inf_{Ω} f̃_{2}\}} C_{∂Ω, f, ω_{0}}^{2} ∥φ∥_{H^{1/2}}

+ C_{∂Ω, f̃, ω_{0}}^{2} ∥φ∥_{H^{1/2}} \frac{max \{∥h_{1}∥_{∞}, ∥h_{2}∥_{∞}\}}{min \{ inf_{Ω} f_{1}, inf_{Ω} f_{2}\}})

+ C_{f̃, ω_{0}} \{(\frac{max \{∥h_{1}∥_{∞}, ∥h_{2}∥_{∞}\}}{min \{ inf_{Ω} f̃_{1}, inf_{Ω} f̃_{2}\}})^{2} C_{∂Ω, f, ω_{0}}^{2} ∥φ∥_{H^{1/2}}

+ \frac{C_{∂Ω, f, ω_{0}} C_{∂Ω, f̃, ω_{0}} ( max \{∥h_{1}∥_{∞}, ∥h_{2}∥_{∞}\})^{2}}{( min \{ inf_{Ω} f_{1}, inf_{Ω} f_{2}\}) ( min \{ inf_{Ω} f̃_{1}, inf_{Ω} f̃_{2}\})} ∥φ∥_{H^{1/2}}

+ C_{∂Ω, f̃, ω_{0}}^{2} ∥φ∥_{H^{1/2}} (\frac{max \{∥h_{1}∥_{∞}, ∥h_{2}∥_{∞}\}}{min \{ inf_{Ω} f_{1}, inf_{Ω} f_{2}\}})^{2}\}  (3.62)

≤ C′_{∂Ω, γ, ω_{0}}( max {∥h_{1}∥_{∞}, ∥h_{2}∥_{∞}})^{2} ∥φ∥_{H^{1/2}}.  (3.63)

Hierbei berechnet sich in (3.61) mit C_{f̃,ω_{0}} = max{ ∥f̃_{1}∥_{∞}, ∥f̃_{2}∥_{∞}, ω_{0} } die Konstante C_{f̃, ω_{0}} > 0 . Die Konstanten C_{∂Ω, f, ω_{0}}, C_{∂Ω, f̃, ω_{0}} > 0 in (3.62) ergeben sich aus den Abschätzungen (3.52) und (3.32).

Die Konstante C′_{∂Ω, γ, ω_{0}} > 0 in (3.63), die sich aus dem Faktor ergibt, der übrig bleibt,

wenn in (3.62) der Term ( max {∥h_{1}∥_{∞}, ∥h_{2}∥_{∞}})^{2} ∥φ∥_{H^{1/2}} ausgeklammert wird, ist

unabhängig von der Differenz (f − f̃). Daher folgt aus (3.63) die Abschätzung (3.57), was zu zeigen war.

□

3.1.5 Greensche Funktion

Für einfache Gebiete Ω ⊂ IR^{n}, n ≥ 2, läßt sich die Greensche Funktion G_{0}(y; x) (mit Singularität in einem Punkt x ∈ Ω) für den Differentialoperator

L_{f}u ≡ [−∇_{y} · (κ_{0}∇_{y}) + (μ_{0} + i ω_{0})]u  (3.64)

mit konstanten Koeffizienten (κ_{0}, μ_{0}) konkret angeben, [2], [63].

Zum Beispiel lautet die Greensche Funktion bezüglich einer Kugel um den Koordinatenursprung mit Radius R, Ω = B_{R}(0) = {y ∈ IR^{n}, |y| < R}, mit einer Singularität im Punkt x ∈ B_{R}(0)

G_{0}(y; x) = \{ \begin{matrix}γ_{n}(|y − x|) − γ_{n} (\frac{|x|}{R}|y − x^{∗}|) , & x ≠ 0, x^{∗} := R^{2}/|x|^{2} \\ γ_{n}(|y|) − γ_{n}(R), & x = 0, \end{matrix}  (3.65)

wobei die Fundamentallösung γ_{n}(r), r > 0, des Differentialoperators L_{f_{0}} gegeben ist durch

γ_{n}(r) = \frac{c_{n}}{κ_{0}} (\frac{b}{r} )^{(n−2)/2} K_{(n−2)/2}(b r) = \{ \begin{matrix}\frac{1}{2πκ_{0}} K_{0}(b r), & n = 2, \\ \frac{1}{4πκ_{0} r} e^{−br}, & n = 3, \end{matrix}  (3.66)

mit b^{2} := (μ_{0} + i ω_{0})/κ_{0} und c_{n} := 2\sqrt{2^{n−2}/(4π)^{n}}. Dies läßt sich wie folgt herleiten.

Sei durch

W_{t}(x) ≡ W (x, t) = (4πt)^{−n/2}e^{−|x|^{2}/4t} , t > 0, x ∈ IR^{n},  (3.67)

der Gauß-Weierstraß-Kern definiert. Es gilt

W (x, t) = t^{−n/2}W_{1}(t^{−1/2}x) und Ŵ_{t}(ξ) = e^{−4π^{2}|ξ|^{2}t},  (3.68)

wobei die Fourier-Tranformation û(ξ) durch

û(ξ) := ∫_{IR^{n}} u(x)e^{−2πi x·ξ} dx  (3.69)

definiert sei. Da

(a + s)^{−β} = \frac{1}{Γ(β)} ∫_{0}^{∞}τ^{β−1}e^{−(a+s)τ} dτ (Re β > 0, Re a > 0),  (3.70)

ist, folgt mit β = \frac{1}{2}α für

B̃_{α}(x) := \frac{1}{Γ(\frac{1}{2}α)} ∫_{0}^{∞}τ^{α/2−1}e^{−aτ} W (x, τ) dτ

= \frac{1}{Γ(\frac{1}{2}α)(4π)^{n/2}} ∫_{0}^{∞}τ^{(α−n)/2−1}e^{−aτ−|x|^{2}/4τ} dτ

= \frac{a^{(n−α)/2}}{Γ(\frac{1}{2}α)(4π)^{n/2}} ∫_{0}^{∞}τ^{(α−n)/2−1}e^{−τ−a|x|^{2}/4τ} dτ

= b^{n−α}B_{α}(bx),  b^{2} = a,  (3.71)

daß

(aI − Δ)^{α/2}B̃_{α} = δ_{0}  (3.72)

ist. Somit folgt mit den Definitionen a := (μ_{0} + i ω_{0})/κ_{0}, ω_{0} := ω/c, daß γ_{n} = \frac{1}{κ_{0}} B̃_{2} gilt. Die Bezeichnung B_{α}(x) steht dabei für die Bessel-Potentiale, die durch

B_{α}(x) = c_{n,α}|x|^{−(n−α)/2}K_{(n−α)/2}(|x|)  (3.73)

mit gewissen Konstanten c_{n,α} > 0 und den modifizierten Besselfunktionen K_{ν}(r) zweiter Art ν-ter Ordnung definiert werden können. In der Tat ergeben sich aus der Darstellung

K_{ν}(r) = ∫_{0}^{∞}cosh(νt)e^{−r cosh(t)} dt  (r > 0, ν ∈ IR)

= \frac{1}{2} (\frac{r}{2})^{ν} ∫_{0}^{∞}τ^{−ν−1}e^{−τ−r^{2}/4τ} dτ  (τ = \frac{1}{2}re^{−t}).  (3.74)

durch Vergleich mit (3.71) für a = 1 die Werte c_{n,α} zu

c_{n,α} = \frac{2}{Γ(α/2)}\sqrt{\frac{2^{n−α}}{(4π)^{n}}} .  (3.75)

Somit erhält man für b = a^{1/2} und α = 2 das Resultat

γ_{n}(|x|) = \frac{c_{n,2}}{κ_{0}|x|^{(n−2)/2}} K_{(n−2)/2}(b|x|) = \{ \begin{matrix}\frac{1}{2πκ_{0}} K_{0}(b|x|), & falls x ∈ IR^{2} \ \{0\}, \\ \frac{1}{4πκ_{0}|x|}e^{−b|x|}, & falls x ∈ IR^{3} \ \{0\}, \end{matrix} (3.76)

wie behauptet.

Mit der Greenschen Funktion, sofern sie existiert, lassen sich Näherungsformeln für die Lösung des direkten Problems angeben. Dies ist Gegenstand der beiden folgenden Abschnitte.

3.1.6 Born–Näherung

Vorausgesetzt, die Parameter f der Gleichung (3.1) können als kleine Störung h um die konstanten positiven Werte f_{0} = (κ_{0}, μ_{0})^{t}, κ_{0} > 0, μ_{0} > 0, geschrieben werden,

f(x) = f_{0} + h(x).  (3.77)

Zusätzlich gelte supp h ⊂⊂ Ω × Ω und die Gleichung (3.1) sei umgeformt zu

L_{f_{0}}u = q − L̃_{h}u,  u|_{∂Ω} = φ.  (3.78)

Nun existiere G_{0}(x, y) = G_{f_{0}}(x, y) die Greensche Funktion bezüglich des Operators L_{f_{0}} mit homogenen Randwerten φ = 0 (Singularität von G_{0} bei y = x), dann ergibt sich für eine Lösung u_{f} von (3.1) die Integralgleichung

u_{f}(x) = ∫_{∂Ω} φ(y) ∂_{ν_{y}}G_{0}(x, y) dσ_{y} + ∫_{Ω} G_{0}(x, y)(q(y) − (L̃_{h}u_{f})(y)) dy

= v_{φ}(x) + ∫_{Ω} G_{0}(x, y)(∇_{y} · (h_{1}(y)∇_{y}u_{f}) − h_{2}(y)u_{f} ) dy  (3.79)

mit

v_{φ}(x) = ∫_{∂Ω} φ(y) ∂_{ν_{y}}G_{0}(x, y) dσ_{y}.  (3.80)

Die erste Bornsche Näherung besteht hier nun darin, in der Gleichung (3.79) die Funktion u_{f} unter dem Integralzeichen zu ersetzen durch die Funktion v_{φ}, gerechtfertigt durch die Annahme, daß die Störung h klein genug sei. Da supp h_{1} ⊂⊂ Ω vorausgesetzt war, erhält man nach partieller Integration die erste Born–Näherung

u_{B}(x) = v_{φ}(x) − ∫_{Ω} (h_{1}(y) ∇_{y}G_{0}(x, y) · ∇_{y}v_{φ}(y) + h_{2}(y) G_{0}(x, y)v_{φ}(y)) dy. (3.81)

Die weiteren Bornschen Näherungen ergeben sich, indem man u_{B}^{(1)} := u_{B} setzt und rekursiv

u_{B}^{(k+1)}(x) = v_{φ}(x)−∫_{Ω} (h_{1}(y) ∇_{y}G_{0}(x, y)·∇_{y}u_{B}^{(k)}(y)+h_{2}(y) G_{0}(x, y)u_{B}^{(k)}(y)) dy (3.82)

berechnet, siehe [27, Abschn. 6.2.1].

3.1.7 Rytov–Näherung

Es sei mit der Greenschen Funktion G_{0}(x, y) des Differentialoperators L_{f_{0}} bezüglich konstanter Koeffizienten f_{0} = (κ_{0}, μ_{0})^{t} und homogener Randwerte die Funktion

v_{φ}(x) = ∫_{∂Ω} φ(y)∂_{ν_{y}}G_{0}(x, y) dσ_{y}  (3.83)

definiert. Die Rytov–Näherung einer Lösung u des direkten Problems

L_{(f_{0}+h)}u = 0,  u|_{∂Ω} = φ,  (3.84)

wird berechnet, indem man für u den Ansatz

u(x) = v_{φ}(x)e^{w(x)}  (3.85)

mit einer noch zu bestimmenden Funktion w macht. Eingesetzt in die Gleichung (3.84) ergibt das

−κ_{0}(v_{φ} Δw + 2 ∇w · ∇v_{φ} + v_{φ}|∇w|^{2}) =  (3.86)
∇ · (h_{1}∇v_{φ}) − h_{2}v_{φ} + h_{1}(2 ∇w · ∇v_{φ} + v_{φ} Δw + v_{φ}|∇w|^{2}) + v_{φ} ∇h_{1} · ∇w.

Vernächlässigt man hierin den Term |∇w|^{2} bzw. den dritten und vierten Term auf der rechten Seite der Gleichung, so ergibt sich die Rytov–Näherung als

u_{R}(x) = v_{φ}(x)e^{w_{R}(x)}  (3.87)

mit der Lösung w_{R}(x) der Gleichung

−κ_{0} Δ(v_{φ} w) + (μ_{0} + i ω_{0})(v_{φ} w) = −κ_{0}(v_{φ} Δw + 2 ∇w · ∇v_{φ})
= ∇ · (h_{1}∇v_{φ}) − h_{2}v_{φ}.  (3.88)

Somit löst

u_{R}(x) = v_{φ}(x) exp (−\frac{1}{v_{φ}(x)} ∫_{Ω} (G_{0}(x, ξ)v_{φ}(ξ)h_{2}(ξ)  (3.89)

+ (∇_{ξ}G_{0} · ∇_{ξ}v_{φ})h_{1}(ξ)) dξ)

das direkte Problem der optischen Tomographie im Frequenzbereich in der Rytov– Näherung, siehe auch [34], [27, Abschn. 6.2.2].

3.1.8 Liouville–Transformation

Unter der Liouville–Transformierten der Lösung u(x) einer elliptischen Differentialgleichung

[−∇_{x} · (κ(x)∇_{x}) + μ(x) + i (ω/c)]u = 0  in Ω  (3.90)

versteht man die Funktion

ϕ(x) := u(x)\sqrt{κ(x)}.  (3.91)

Dann rechnet man nach, daß diese Funktion ϕ(x) die stationäre Schrödingergleichung

[−Δ_{x} + q(x)]ϕ = 0  in Ω  (3.92)

erfüllt, falls die komplexwertige Potentialfunktion q(x) gesetzt wird zu

q(x) := \frac{1}{κ(x)} (μ(x) + \sqrt{κ(x)}Δ_{x}\sqrt{κ(x)} + i (ω/c)) .  (3.93)

Dabei wird vorausgesetzt, daß das Infimum der reellwertigen Funktionen μ(x) und κ(x) auf Ω strikt positiv bleibt und daß die Koeffizientenfunktionen jeweils hinreichend glatt sind, so daß alle auftretenen Ausdrücke stets existieren.

3.1.9 Transformation der Fréchet–Ableitung

Führt man für die Lösung u(x) in der Abbildung P_{f} = P_{(κ, μ)} die Liouville–Transformation ϕ := u\sqrt{κ} durch, und setzt weiter voraus, daß die Funktion κ(x) in der Nähe von ∂Ω konstant bleibt, so erhält man für die Abbildung

P_{q}^{Sch} : H^{1/2}(∂Ω) ⟶ H^{−1/2}(∂Ω), φ ⟼ (\sqrt{κ} ∂_{ν}ϕ)|_{∂Ω}  (3.94)

die Identität

P_{(κ, μ)} ≡ P_{t(κ, μ)}^{Sch}  für  t(κ, μ) := \frac{μ + \sqrt{κ}Δ\sqrt{κ}}{κ} + i \frac{ω}{cκ}.  (3.95)

Dabei steht der obere Index ^{Sch} als Kennzeichen für die Schrödingergleichung (3.92), die zur Berechnung eines Elementes im Bildbereich zu lösen ist. Die Fréchet–Ableitung (A_{f}h)(φ) an der Stelle f transformiert sich nun gemäß der Kettenregel

(A_{f}h)(φ) = (A_{t(κ,μ)}^{Sch} \frac{∂t}{∂(κ, μ)}h)(φ) = (A_{q}^{Sch}k)(φ)  (3.96)

mit q := t(κ, μ) und k := J_{t}(f)h. Die Jacobi-Matrix J_{t}(f) an der Stelle f = (κ, μ) berechnet sich zu

J_{t}(κ, μ) = \frac{∂t}{∂(κ, μ)} = \frac{1}{κ} ( \begin{matrix}\frac{1}{2}Δ_{x}(Id − \sqrt{κ}) − \frac{μ}{κ} & Id \\ −\frac{ω}{cκ} & 0 \end{matrix} ) .  (3.97)

Damit kann der komplexwertige Vektor k = k^{re} +i k^{im} als zweidimensionale reellwertige Größe k = (k^{re}, k^{im})^{t} identifiziert werden. Explizit ausgeschrieben gilt somit

(A_{f}h)(φ) = (A_{q}^{Sch}k)(φ) = (\sqrt{κ} ∂_{ν}ψ)|_{∂Ω},  (3.98)

wobei ψ(x) ∈ H_{0}^{1}(Ω) die Lösung der Gleichung

[−Δ_{x} + q(x)]ψ = −k(x)ϕ in Ω,

(3.99)

ψ = 0  auf ∂Ω

mit den Koeffizienten q(x) = t(κ(x), μ(x)) und

k(x) = \frac{1}{κ(x)} \{\frac{1}{2}Δ_{x}(h_{κ}(x) − \sqrt{κ(x)}) − \frac{μ(x)}{κ(x)}h_{κ}(x) + h_{μ}(x) − i \frac{ω}{cκ(x)}h_{κ}(x)\} (3.100)

ist.

Zwar wird durch die Liouville–Transformation das Problem der optischen Tomographie in ein rein komplexwertiges Problem gewandelt, jedoch steht wegen (3.100) zu erwarten, daß die Abhängigkeit der Daten von dem Parameter q(x) = q_{r}(x)+i q_{i}(x) in P_{q}^{Sch} stark nichtlinear sein wird, da die Abhängigkeit derselben Daten von den Parametern κ(x) und μ(x) in P_{(κ,μ)} nur schwach nichtlinear ist.

3.2 Inverses Problem

Die Bestimmung der Parameter f eines linearen Differentialoperators L_{f} kann als inverses Problem bezeichnet werden, da im Gegensatz zu einem direkten Problem (3.1) nicht die Lösung u_{f}, sondern der Operator L_{f} anhand gegebener Daten gesucht ist. Dabei seien die Daten

g = M(u_{f})  (3.101)

in Abhängigkeit der Lösung u gegeben. Die Abbildung M : U(Ω) ⟶ G(∂Ω) bilde die Menge der Lösungen linear in eine Menge von (verallgemeinerten) Funktionen auf ∂Ω ab.

Ein inverses Randwertproblem ist dann durch folgende Aufgabe gegeben.

Definition 3.2.1 (Inverses Randwertproblem)

Bestimme bei gegebenen Randwertdaten g ∈ G(∂Ω) und φ ∈ R(∂Ω) die Parameter f ∈ F_{1}(Ω) × … × F_{m}(Ω) des Differentialoperators L_{f}, so daß die Gleichungen

M(u_{f}) = g in ∂Ω,  L_{f}u = q in Ω, u|_{∂Ω} = φ  (3.102)

erfüllt sind. Hierbei ist q ∈ Q(Ω) ein Term unabhängig von L_{f}, f und u.

Im vorliegenden Modell der optischen Tomographie sind die reellwertigen Parameter f = (f_{1}, f_{2}) = (κ, μ_{a}) ∈ L^{∞}(Ω : IR) × L^{∞}(Ω : IR) gesucht. Die Meßdaten M(u_{f}) für diese Aufgabe werden durch die Bildung der Normalenableitung am Rand ∂Ω, also durch die Abbildung P_{f} : H^{1/2}(∂Ω) ⟶ H^{−1/2}(∂Ω),

P_{f}(φ) = (f_{1} ∂_{ν}u_{f} )|_{∂Ω},  (3.103)

modelliert. Für verschiedene φ ∈ H^{1/2}(∂Ω) sei somit M(u_{f, φ}) := P_{f}(φ) gesetzt, wobei u_{f, φ} ∈ H^{1}(Ω) die Lösung des direkten Problems

[−∇_{x} · (κ(x)∇_{x}) + (μ_{a}(x) + i ω_{0})]u = 0 in Ω,  u|_{∂Ω} = φ  (3.104)

bedeute.

Mit R : F_{1}(Ω) × … × F_{m}(Ω) ⟶ G(∂Ω),

R(f) := M(u_{f}) − g  (3.105)

läßt sich formal ein inverses Problem als Bestimmung einer Lösung der Gleichung

R(f) = 0  (3.106)

formulieren. Falls diese Gleichung keine Lösung besitzt, weil zum Beispiel die Meßdaten nur ungenau bekannt sind, beschränkt man sich auf die Bestimmung von

min ∥R(f)∥^{2} = min ∥M(u_{f}) − g∥^{2},  (3.107)

3.2 INVERSES PROBLEM

wobei vorausgesetzt wird, daß G(∂Ω) ein normierter (vollständiger) Raum oder sogar ein Hilbertraum ist. Unter allen möglichen Lösungen von (3.107) soll diejenige mit minimaler Norm ∥f∥^{2} ausgezeichnet werden. Für lineare Abbildungen stellt eine solchermaßen gekennzeichnete Lösung die verallgemeinerte (oder Moore-Penrose) Lösung der Gleichung (3.106) dar. Die Ergebnisse dieser Theorie lassen sich allerdings nicht ohne weiteres auf das zu untersuchende inverse Problem übertragen, da die Zuordnung f ↦ u_{f} im allgemeinen nichtlinear ist, und damit auch die Abbildung R : F(Ω) ⟶ G(∂Ω).

3.2.1 Eindeutigkeit

Dieser Abschnitt beschäftigt sich damit, unter welchen Voraussetzungen es genau eine Lösung f = (κ, μ)^{t} zu (3.107) geben kann. Injektivitätsaussagen dieser Art findet man in [20], [64].

Um den Rahmen dieser Arbeit nicht zu sprengen, soll nur der zweidimensionale Fall, welcher in [20] betrachtet wird, wiedergegeben werden. Wegen der Liouville–Transformation genügt es, die Eindeutigkeitsaussage für die Schrödingergleichung

[Δ_{x} − q(x)]u = 0,  q(x) := \frac{μ(x) + \sqrt{κ(x)} Δ_{x}\sqrt{κ(x)}}{κ(x)} + i \frac{ω}{c κ(x)}  (3.108)

mit der Randbedingung u|_{∂Ω} = −φ\sqrt{κ} zu zeigen. Ein wesentliches Hilfsmittel zum Beweis der Eindeutigkeit ist die folgende Orthogonalitätsrelation.

Satz 3.2.1 Es seien

L_{1,2}w ≡ [∇_{x} · (κ_{1,2}(x)∇_{x}) − (μ_{1,2}(x) + i ω_{0})]w  (3.109)

zwei Differentialoperatoren, die die gleichen Datenfunktionale P_{1}(φ) = P_{2}(φ) erzeugen. Beiden Funktionen κ_{1}(x) und κ_{2}(x) seien auf dem Rand ∂Ω des beschränkten Gebietes Ω ⊂ IR^{2} identisch.

Dann gilt für die entsprechenden Lösungen u_{1,2}(x) ∈ H^{2}(Ω) ∩ C^{1}(Ω) bezüglich der Gleichung (3.108) die Relation

∫_{Ω} (q_{1} − q_{2}) u_{1}u_{2} dx = \frac{1}{2} ∫_{∂Ω} (∂_{ν} log \frac{κ_{1}}{κ_{2}} ) u_{1}u_{2} dσ(x).  (3.110)

Beweis: Anhand der Greenschen Formel und wegen u_{1,2} = w_{1,2}\sqrt{κ_{1,2}} erhält man die Gleichungskette

∫_{Ω} (q_{1} − q_{2}) u_{1}u_{2} dx = ∫_{Ω} (u_{2} Δu_{1} − u_{1} Δu_{2}) dx

= ∫_{∂Ω} (u_{2} ∂_{ν}u_{1} − u_{1} ∂_{ν}u_{2}) dσ  (3.111)

= \frac{1}{2} ∫_{∂Ω} (∂_{ν} log \frac{κ_{1}}{κ_{2}} ) u_{1}u_{2} dσ + ∫_{∂Ω} κ (w_{2} ∂_{ν}w_{1} − w_{1} ∂_{ν}w_{2}) dσ.

Aufgrund der Voraussetzung P_{1}(φ) = P_{2}(φ) verschwindet der letzte Term auf der rechten Seite.

□

Es seien also q_{1}(x), q_{2}(x) ∈ L^{p}(Ω) zwei komplexwertige Verteilungen mit kompakten Trägern in dem beschränkten Gebiet Ω ⊂ IR^{2}, die die gleichen Datenfunktionale

P_{q_{1}}^{Sch}(φ) = P_{q_{2}}^{Sch}(φ) für alle φ ∈ C_{0}^{∞}(∂Ω)  (3.112)

erzeugen. Weiter seien für p > 2 die Konstanten

N_{p} := max\{ ∥q_{1}∥_{L^{p}}, ∥q_{2}∥_{L^{p}} \},  c_{p}(Ω) := d_{Ω} (\frac{2π d_{Ω}^{2}}{2 − (1 − 1/p)^{−1}} )^{\frac{1}{1−1/p}}  (3.113)

mit d_{Ω} := diam(Ω) definiert. Damit ergibt sich die folgende Aussage.

Satz 3.2.2 Für ein 2 < p ≤ 3 gelte die Ungleichung c_{p}(Ω) N_{p} ≤ \frac{1}{10} und es sei (3.112) erfüllt.

Dann gilt q_{1} ≡ q_{2} in Ω.

Beweis: Zum Beweis dieses Satzes sei zunächst eine weitere Transformation der Gleichung (3.108) durch

Φ(x) := u(x) exp(−ζ · x) − 1  (3.114)

mit einem noch näher zu bestimmenden komplexwertigen Vektor ζ ∈ C^{2} vorgenommen. Damit transformiert sich (3.108) zu

[Δ_{x} + 2 (ζ · ∇_{x})]Φ = q(x)(Φ + 1).  (3.115)

Mit der Fundamentallösung

g_{ζ}(x) = (\frac{1}{−4π^{2} |ξ|^{2} + 4πi (ξ · ζ)}) (x)  (3.116)

zu der linken Seite der Gleichung (3.115) läßt sich die Lösung Φ(x) jener Gleichung als Faltungsintegral

Φ(x) = (g_{ζ} ∗ (q Φ + q))(x)  (3.117)

darstellen. Hierbei bedeute

(f(ξ)) (x) = ∫_{IR^{2}} f(ξ) e^{2πi x·ξ} dξ  (3.118)

die inverse Fourier-Transformation im IR^{2}. Nun sei der Operator

[G^{±}Q]f := (g_{ζ^{±}} ∗ (q f))  (3.119)

mit ζ^{±} := πk (−i w ± w^{⊥}), w ∈ IR^{2}, |w|^{2} = 1 und k ≥ 0 definiert. Unter der Voraussetzung p > 2 wird in [20, Prop. 2.3] die Abschätzung

∥G^{±}Q∥ ≤ \frac{5 c_{p}(Ω)}{k diam(Ω)} ∥q∥_{L^{p}}  (3.120)

gezeigt. Mit der Bedingung c_{p}(Ω) ∥q∥_{L^{p}} ≤ \frac{1}{10} folgt daraus für Φ(x) = Φ^{±}(x) die Gleichung

Φ^{±}(x) = ((1 − G^{±}Q)^{−1}(g_{ζ^{±}} ∗ q))(x)  (3.121)

und die Abschätzung

|Φ^{±}(x)| ≤ \frac{10 c_{p}(Ω)}{k diam(Ω)} ∥q∥_{L^{p}},  (3.122)

falls k diam(Ω) ≥ 1 ist. Aus dieser Abschätzung folgt, daß

u^{±}(x) = exp(ζ^{±} · x)(1 + Φ^{±}(x)) ∈ H^{2}(Ω) ∩ C^{1}(Ω)  (3.123)

gilt.

Es sei nun mit Φ_{1}^{±}(x) die Lösungen Φ^{±}(x) ∈ H^{2}(Ω) ∩ C^{1}(Ω) der Gleichung (3.115) bezüglich des Potentials q_{1}(x) bezeichnet, und mit Φ_{2}^{±}(x) entsprechend die Lösungen der Gleichung (3.115) bezüglich q_{2}(x).

Dann ergibt sich für die Differenz r(x) := q_{1}(x) − q_{2}(x) aufgrund der Orthogonalitätsrelation (3.110) mit ξ = kw die Gleichung

r̂(ξ) = I_{∂Ω}(ξ) − ∫_{Ω} r(x) e^{−2πi kw·x} (Φ_{1}^{+} + Φ_{2}^{−} + Φ_{1}^{+} Φ_{2}^{−})(x) dx,  (3.124)

wobei das Randintegral I_{∂Ω} gegeben ist durch

I_{∂Ω}(ξ) = \frac{1}{2} ∫_{∂Ω} (∂_{ν} log \frac{κ_{1}(x)}{κ_{2}(x)}) e^{−2πi kw·x} (1 + Φ_{1}^{+} + Φ_{2}^{−} + Φ_{1}^{+} Φ_{2}^{−})(x) dσ(x). (3.125)

Mit der Bedingung p ≤ 3 lassen sich daraus die Abschätzungen

∥r̂(ξ) − I_{∂Ω}∥_{L^{2}(|ξ|≥ \frac{1}{diam(Ω)})} ≤ 5.43 c_{p}(Ω) N_{p} ∥r̂∥_{L^{2}},  (3.126)
∥r̂(ξ) − J_{∂Ω}∥_{L^{2}(|ξ|≤ \frac{1}{diam(Ω)})} ≤ 0.05 c_{p}(Ω) N_{p} ∥r̂∥_{L^{2}}  (3.127)

herleiten, wobei das Randintegral J_{∂Ω}(ξ) ähnlich wie das Integral I_{∂Ω}(ξ) festgelegt wird, [20, Cor. 2.11, Cor. 2.13].

Da nach Voraussetzung die Randintegrale I_{∂Ω} und J_{∂Ω} verschwinden, folgt aus den letzten beiden Abschätzung die Ungleichung

∥r̂∥_{L^{2}} ≤ 6 c_{p}(Ω) N_{p} ∥r̂∥_{L^{2}},  (3.128)

und damit wegen der Bedingung c_{p}(Ω) N_{p} ≤ \frac{1}{10}, daß r̂ ≡ 0 ist. Daher muß für x ∈ Ω r(x) = q_{1}(x) − q_{2}(x) = 0 sein.

□

3.2.2 Schlecht gestelltes Problem

Nach Hadamard [23] lautet ein Problem, gegeben durch eine (lineare oder nichtlineare) Gleichung Ax = b, gut gestellt, wenn

1. es eine Lösung besitzt,

2. die Lösung eindeutig festgelegt ist und

3. die Lösung stetig von den Daten b abhängt.

Ist eine der drei Bedingungen nicht erfüllt, so spricht man von einem schlecht gestellten Problem.

Schon das linearisierte inverse Problem der optischen Tomographie, formuliert durch die Gleichung

(A_{f_{0}} h)(φ) = g(φ)  (3.129)

mit der Fréchet–Ableitung (A_{f_{0}} ·)(φ) an der Stelle f_{0} := (κ_{0}, μ_{0})^{t} bezüglich der Randverteilungen φ ∈ H^{1/2}(∂Ω), und mit den gegebenen Daten g(φ) ∈ H^{−1/2}(∂Ω), läßt sich als sehr schlecht gestellt erkennen, da der dritte Punkt der Definition eines gut gestellten Problems klar verletzt wird. Dies kann man mit Hilfe der Methode von Calderon [10] wie folgt erkennen.

Es sei u_{0} := u_{f_{0}} ∈ H^{1}(Ω) die eindeutige Lösung der Gleichung L_{f_{0}}u = 0, u|_{∂Ω} = φ_{0}. Die Funktion w ∈ H^{1}(Ω) sei die Lösung von

L_{f_{0}}w = ∇ · (h_{1}(x)∇u_{0}) − h_{2}(x)u_{0},  w|_{∂Ω} = 0.  (3.130)

Aus der zweiten Greenschen Formel

∫_{Ω} (z L_{f_{0}}w − w L_{f_{0}}z) dx = −κ_{0} ∫_{∂Ω} (z ∂_{ν}w − w ∂_{ν}z) dσ_{x}  (3.131)

folgt zunächst mit einer Funktion z ∈ H^{1}(Ω), die L_{f_{0}}z = 0 erfüllt, und einer partiellen Integration die Identität

κ_{0} ∫_{∂Ω} (z ∂_{ν}w) dσ_{x} = ∫_{Ω} (h_{1}(x) ∇u_{0} · ∇z + h_{2}(x) u_{0} z) dx.  (3.132)

Wählt man auf der rechten Seite für z und u_{0} die Ansatzfunktionen

z_{ξ}(x) := e^{x·ζ^{+}},  u_{0, ξ}(x) := e^{x·ζ^{−}},  (3.133)

wobei ζ^{±} ∈ C^{n}, n ≥ 3, die Eigenschaften

ζ^{+} + ζ^{−} = −2πi ξ, ξ ∈ IR^{n},  (3.134)
ζ^{+} · ζ^{−} = −2π^{2} |ξ|^{2} − \frac{μ_{0} + i ω_{0}}{κ_{0}} ,  (3.135)
0 = −4π^{2} κ_{0} ζ · ζ + μ_{0} + i ω_{0} für ζ = ζ^{±}  (3.136)

besitzen sollen. Dies ist für

ζ^{±} := ±(α + i β)η − πi ξ  (3.137)

mit η ∈ ξ^{⊥}, |η| = 1 und

α := \sqrt{\frac{ω_{0}}{2κ_{0}}} (ϵ + \sqrt{1 + ϵ^{2}})^{1/2} , β := \sqrt{\frac{ω_{0}}{2κ_{0}}} (ϵ + \sqrt{1 + ϵ^{2}})^{−1/2} ,

ϵ := \frac{κ_{0}}{ω_{0}} (\frac{μ_{0}}{κ_{0}} + π^{2}|ξ|^{2})  (3.138)

möglich. Die Funktionen z_{ξ} und u_{0, ξ} erfüllen beide die Gleichung L_{f_{0}}u = 0 und u_{0, ξ} besitzt die Randwerte φ_{ξ} := e^{x·ζ^{−}} |_{x∈∂Ω}. Dann erhält man aus (3.132)

G(ξ) := ∫_{∂Ω} (z_{ξ} g(φ_{ξ})) dσ_{x}

= ∫_{Ω} (ζ^{+} · ζ^{−} h_{1}(x) + h_{2}(x))e^{−2πi ξ·x} dx,  (3.139)

also mit der inversen Fourier-Transformation

G̃(x) = ∫_{IR^{n}} G(ξ)e^{2πi x·ξ} dξ  (3.140)

die Gleichung

(\frac{1}{2}Δ_{x} − \frac{μ_{0} + i ω_{0}}{κ_{0}} )h_{1}(x) + h_{2}(x) = G̃(x),  x ∈ Ω.  (3.141)

Daraus lassen sich jetzt die gesuchten Funktionen h_{1}(x) und h_{2}(x) durch Vergleich der Real- und Imaginärteile bestimmen: zuerst die Funktion h_{1}(x) aus den Imaginärteilen, und mit der dann bekannten Funktion h_{1}(x) die zweite Funktion h_{2}(x) aus den Realteilen. Da jedoch für ζ^{+} die Asymptotik ζ^{+} = O(|ξ|) für |ξ| → ∞ gilt, ist

z_{ξ} = O(e^{|x|·x}) für |ξ| → ∞.  (3.142)

Damit ist die inverse Fourier-Transformierte G̃(x) auf keinen Fall stetig von den Daten g(φ_{ξ}) abhängig, und das lineare inverse Problem ist sehr schlecht gestellt.

3.2.3 Iterativer Lösungsansatz

Für die exakte Lösung f = f_{•} einer nichtlinearen Gleichung

R(f) = 0  (3.143)

gibt es kein allgemeingültiges Rezept. Ist die Abbildung R : X ⟶ Y, wobei X und Y Hilberträume seien, (Fréchet–)differenzierbar, so stellt das Newton–Verfahren

f_{k+1} = f_{k} + h,  R′(f_{k})h = −R(f_{k}),  k = 0, 1, 2, …  (3.144)

ein bewährtes Iterationsverfahren dar, um eine Lösung sukzessiv approximieren zu können, falls die vorzugebene Startnäherung f_{0} nahe genug an der Lösung f_{•} liegt.

Bei tomographischen Anwendungen liegen meist unterschiedliche Messungen, z.B. aus unterschiedlichen Richtungen, ein und desselben Objekts als Daten vor. Liegen etwa p Messungen aus p Richtungen vor, so wird das Objekt f durch p Gleichungen

R_{j}(f) = 0,  j = 1, 2, …, p  (3.145)

beschrieben.

Die Durchführung des Netwon–Verfahrens bezüglich des zusammengefaßten Operators R(f) = (R_{1}(f), …, R_{p}(f))^{t} würde die Berechnung und Invertierung der (p × p)-Matrix (R′_{j}(f_{k})^{∗}R′_{m}(f_{k}))_{j,m} mit R′_{j}(f_{k}) der Fréchet–Ableitung an der Stelle f_{k} und R′_{j}(f_{k})^{∗} der zugehörigen Adjungierten in jedem Schritt des Verfahrens erforderlich machen. Dies ist für tomographische Zwecke sehr aufwendig bis nicht praktikabel.

Eine entscheidene Idee zur Vermeidung dieses Aufwandes ergibt sich dadurch, zur iterativen Lösung des Systems (3.145) zu versuchen, jeweils nur einen Schritt des Newton– Verfahrens nacheinander für eine einzelne Gleichung durchzuführen, unter Berücksichtigung der schon gewonnenen Information. Sind alle Gleichungen abgearbeitet, gelangt man so zu der neuen Näherung f_{k+1}, und man beginne die Prozedur erneut, jetzt mit dieser Näherung f_{k+1} als Startwert.

Im Fall der optischen Tomographie treten verschiedene Randverteilungen φ_{j} an die Stelle der unterschiedlichen Richtungen. Die Fréchet–Ableitung ist dann durch

R′_{j}(f_{k}) = (A_{f_{k}}h)(φ_{j})  (3.146)

gegeben, so daß in jedem Einzelschritt dieses Iterationsverfahrens zur Berechnung eines h eine Gleichung der Form

(A_{f_{k}}h)(φ_{j}) = −R_{j}(f_{k}), j ∈ {1, …, p}  (3.147)

durch Invertierung des Operators (A_{f} ·)(φ_{j}) aufzulösen ist.

Diese Art eines sozusagen verschränkten Newton–Verfahrens ist in bezug auf lineare Abbildungen A_{j}(f) = 0, j = 1, …, p, im Bereich der Computertomographie als \underline{a}lgebraische \underline{R}ekonstruktions\underline{t}echnik (ART) oder auch Kaczmarz–Methode bekannt, siehe [25] [47].

Was zur Lösung der Gleichung (3.147) zu beachten ist, wird in den folgenden Abschnitten erläutert. Besonderes Augenmerk sei hierbei auf die Tatsache gerichtet, daß in (3.147) im Fall der optischen Tomographie im Frequenzbereich die Daten auf der rechten Seite zwar komplexwertig, die jeweils gesuchten Lösungen h hingegen reellwertig sind. Es wird sich herausstellen, daß die korrekte Behandlung dieser zusätzlichen Information entscheidenen Einfluß auf das numerische Verhalten der vorgestellten Rekonstruktionsverfahren hat.

3.2.4 Adjungierte Abbildung

Zu jedem linearen, beschränkten Operator A : X ⟶ Y, f ↦ Af existiert eindeutig eine adjungierte Abbildung

A′ : Y′ ⟶ X′, ℓ ⟼ A′ℓ := ℓ ∘ A,  (3.148)

die ebenfalls beschränkt ist. Mit X′ und Y′ seien hierbei die entsprechenden Dualräume bezeichnet.

Ist X ein Hilbertraum, läßt sich mittels der konjungiert linearen Isometrie

J_{X} : X ⟶ X′, x ↦ ℓ(·) := ⟨ · , x⟩_{X}  (3.149)

aus dem Riesz’schen Darstellungssatz der Raum X mit seinem Dual X′ identifizieren. Sind also X und Y Hilberträume, so erhält man mit

A^{∗} := J_{X}^{−1} ∘ A′ ∘ J_{Y}  (3.150)

einen linearen, beschränkten Operator A^{∗} : Y ⟶ X. Dieser Operator ist charakterisiert durch die Beziehung

⟨x, A^{∗}y⟩_{X} = ⟨Ax, y⟩_{Y} ,  ∀ x ∈ X, y ∈ Y.  (3.151)

Diese Charakterisierung macht zumindest dann Sinn, wenn die Räume X und Y Hilberträume über dem gleichen Körper IK sind.

Nun sei X_{C} ein Hilbertraum über dem Körper der komplexen Zahlen C, also ein C-Vektorraum. Durch Einschränkung der Multiplikation auf rein reelle Skalare wird dieser zu einem IR-Vektorraum und sei mit X_{IR} bezeichnet.

Zwischen den Funktionalen ℓ : X_{C} ⟶ C und k : X_{IR} ⟶ IR besteht der folgende Zusammenhang.

Satz 3.2.3 Es gilt (X_{C})′ ≅ (X_{IR})′ im algebraischen Sinne:

1. Es sei ℓ ∈ (X_{C})′.

Dann ist k := Re ℓ ein Element aus (X_{IR})′ und es gilt

ℓ(x) = k(x) − i k( ix),  x ∈ X_{C}.  (3.152)

2. Es sei k ∈ (X_{IR})′.

Dann ist ℓ : x ⟼ k(x) − i k( ix) ein Element aus (X_{C})′.

Beweis: Durch einfaches Verifizieren:

1. Das Funktional k(x) = Re ℓ(x) ist sicherlich IR-linear und es gilt

ℓ(x) = Re ℓ(x) + i Im ℓ(x)
= Re ℓ(x) − i Re i ℓ(x)
= Re ℓ(x) − i Re ℓ( ix)
= k(x) − i k( ix).  (3.153)

2. Umgekehrt gilt für das IR-lineare Funktional ℓ(x) = k(x)−i k( ix) die Gleichungskette

ℓ(ix) = k( ix) − i k(−x)
= i(k(x) − i k( ix))
= i ℓ(x),  (3.154)

d.h. ℓ(x) ist auch C-linear.

□

Der Satz 3.2.3 soll auf die adjungierten Abbildung A_{f}^{∗} angewendet werden, die zur komplexwertigen Fréchet–Ableitung

(A_{f} ·)(φ) : H^{s_{1}}(Ω) × H^{s_{2}}(Ω) ⟶ H^{−1/2}(∂Ω)  (3.155)

gehört. Der nächste Satz bestimmt zunächst diese adjungierte Abbildung A_{f}^{∗}.

Satz 3.2.4 Es seien φ ∈ H^{1/2}(∂Ω) und f = (f_{1}, f_{2})^{t} ∈ L^{∞}(Ω : IR) × L^{∞}(Ω : IR) mit f_{1} >> 0 und f_{2} >> 0. Die Funktion u_{f} ∈ H^{1}(Ω) sei die eindeutige Lösung der Gleichung

L_{f}u = 0,  u|_{∂Ω} = φ.  (3.156)

Es sei h = (h_{1}, h_{2})^{t} ∈ H_{0}^{s_{1}}(Ω) × H_{0}^{s_{2}}(Ω) mit s_{1}, s_{2} ≥ 0, so daß es für jedes h ∈ H_{0}^{s_{1}}(Ω) × H_{0}^{s_{2}}(Ω) eindeutig eine Funktion w ∈ H^{1}(Ω) gibt, welche die Gleichung

L_{f}w = ∇_{x} · (h_{1}(x)∇_{x}u_{f}) − h_{2}(x)u_{f},  w|_{∂Ω} = 0,  (3.157)

erfüllt, und somit die Fréchet–Ableitung an der Stelle f ∈ L^{∞}(Ω : IR) × L^{∞}(Ω : IR)

(A_{f} ·)(φ) : H_{0}^{s_{1}}(Ω) × H_{0}^{s_{2}}(Ω) ⟶ H^{−1/2}(∂Ω), h ⟼ (f_{1}∂_{ν}w)|_{∂Ω}  (3.158)

festgelegt wird.

Dann ist die zu (A_{f} h)(φ) adjungierte Abbildung

(A_{f}^{∗} ·)(φ) : H^{1/2}(∂Ω) ⟶ H^{−s_{1}}(Ω) × H^{−s_{2}}(Ω)  (3.159)

3.2 INVERSES PROBLEM gegeben durch

(A_{f}^{∗} g)(φ) = ( \begin{matrix}\overline{∇_{x}u_{f}} · ∇_{x}z \\ \overline{u_{f}} z \end{matrix} ) ,  (3.160)

wobei z ∈ H^{1}(Ω) die eindeutige Lösung der adjungierten Aufgabe

\overline{L_{f}}z = 0, z|_{∂Ω} = g (3.161)

für g ∈ H^{1/2}(∂Ω) ist. Beweis: Anhand der zweiten Greenschen Formel

∫_{Ω} (z L_{f}w − w L_{f}z) dx = ∫_{∂Ω} (f_{1}w ∂_{ν}z − z f_{1}∂_{ν}w) dσ_{x}  (3.162)

bezüglich des elliptischen Differentialoperators

L_{f}u ≡ [−∇_{x} · (f_{1}(x)∇_{x}) + (f_{2}(x) + iω_{0})]u  (3.163)

folgt für die Lösung w ∈ H^{1}(Ω) der Gleichung (3.157) und die Lösung z ∈ H^{1}(Ω) der Gleichung (3.161) die Gleichungskette

∫_{Ω} (h_{1}(x)(∇_{x}u_{f} · \overline{∇_{x}z}) + h_{2}(x)(u_{f} \overline{z})) dx

= − ∫_{Ω} (∇_{x} · (h_{1}(x)∇_{x}u_{f}) − h_{2}(x)u_{f} )\overline{z} dx
= − ∫_{∂Ω} (f_{1}w \overline{∂_{ν}z} − \overline{z} f_{1}∂_{ν}w) dσ_{x}
= ∫_{∂Ω} (f_{1}∂_{ν}w)\overline{g} dσ_{x}.  (3.164)

Daher gilt mit der Definition (3.160) für φ ∈ H^{1/2}(∂Ω) beliebig die Gleichung

⟨ h, (A_{f}^{∗} g)(φ) ⟩_{X} = ⟨ (A_{f} h)(φ), g ⟩_{Y}  (3.165)

für alle h ∈ H_{0}^{s_{1}}(Ω) × H_{0}^{s_{2}}(Ω) und alle g ∈ H^{1/2}(∂Ω), wobei ⟨·, ·⟩_{X} die duale Paarung bezüglich (H_{0}^{s_{1}}(Ω) × H_{0}^{s_{2}}(Ω), H^{−s_{1}}(Ω) × H^{−s_{2}}(Ω)) und ⟨·, ·⟩_{Y} die duale Paarung bezüglich (H^{1/2}(∂Ω), H^{−1/2}(∂Ω)) bedeutet.

□

Der Satz 3.2.3 liefert nun mit der Definition

(A_{f}^{∗} g)_{r}(φ) := ( \begin{matrix}Re(\overline{∇_{x}u_{f}} · ∇_{x}z) \\ Re(\overline{u_{f}} z) \end{matrix} )  (3.166)

die Zerlegung

(A_{f}^{∗} g)(φ) = (A_{f}^{∗} g)_{r}(φ) − i(A_{f}^{∗} ig)_{r}(φ)

= ( \begin{matrix}Re(\overline{∇_{x}u_{f}} · ∇_{x}z) \\ Re(\overline{u_{f}} z) \end{matrix} ) − i ( \begin{matrix}Re(\overline{∇_{x}u_{f}} · ∇_{x} z̃) \\ Re & ( & \overline{u_{f}} & z̃ & ) \end{matrix} )

= ( \begin{matrix}Re(\overline{∇_{x}u_{f}} · ∇_{x}z) − i Re(\overline{∇_{x}u_{f}} · ∇_{x} z̃) \\ Re(\overline{u_{f}} z) − i Re(\overline{u_{f}} z̃) \end{matrix} ) ,  (3.167)

wobei z̃ ∈ H^{1}(Ω) die Lösung der adjungierten Differentialgleichung

\overline{L_{f}} z̃ = 0,  z̃|_{∂Ω} = i g  (3.168)

ist.

3.2.5 Verallgemeinerte Lösung

Es seien X und Y Hilberträume über einem Körper IK. Falls g ∈ range(A)+range(A)^{⊥} ⊂ Y eines linearen, beschränkten Operators A : X ⟶ Y ist, besitzt die lineare Gleichung

Af = g,  f ∈ X, g ∈ Y  (3.169)

genau ein Element f^{+} ∈ ker(A)^{⊥} = \overline{range(A^{∗})}, welches den Defekt ∥Af − g∥_{Y} minimiert. Dieses Element f^{+} ∈ \overline{range(A^{∗})} ist gekennzeichnet als die Lösung der Gleichung

A^{∗}Af = A^{∗}g,  (3.170)

die minimale Norm ∥f∥_{X} besitzt, und wird verallgemeinerte (oder Moore–Penrose) Lösung genannt. Offensichtlich ist die Zuordnung g → f^{+} linear. Sie sei mit

A^{+} : (range(A) + range(A)^{⊥}) ⟶ X, g ↦ A^{+}g := f^{+}  (3.171)

bezeichnet.

Besteht nun die Aufgabe, eine Lösung der Gleichung Af = g zu finden, die neben (3.169) noch weitere lineare Bedingungen erfüllt, so kann eine verallgemeinerte Lösung unter Nebenbedingungen gefunden werden, wie die nächsten Sätze zeigen. Um technischen Schwierigkeiten aus dem Wege zu gehen, seien sie nur für endlich-dimensionale Hilberträume formuliert.

Satz 3.2.5 Es seien X, Y und Z endlich-dimensionale Hilberträume über IK. Mit einem Element b ∈ Z und einem linearen, beschränkten Operator N : X ⟶ Z sei

N := { f ∈ X, Nf = b }  (3.172)

eine Menge, die die Nebenbedingungen für f ∈ X beschreibt.

Dann ist die Menge der Elemente M ⊂ N, für die das Minimum

min_{f∈N} ∥Af − g∥_{Y}  (3.173)

angenommen wird, gegeben durch

M = { (AP)^{+}(g − AN^{+}b) + N^{+}b + (I − (AP)^{+}A)ξ, ξ ∈ ker(N) },  (3.174)

wobei P := I_{X} − N^{+}N gesetzt ist. Das Minimum selbst berechnet sich zu

min_{f∈N} ∥Af − g∥_{Y} = ∥(I − A(AP)^{+})(AN^{+}b − g)∥_{Y} .  (3.175)

Beweis: Zu finden in [11, Thm. 3.6.2].

□

In der Lösungsmenge M für das Minimierungsproblem mit Nebenbedingungen N läßt sich nun ebenfalls ein Element mit minimaler Norm finden.

Satz 3.2.6 Das Element m^{+} ∈ M mit minimaler Norm ∥m∥_{X} ist gegeben durch

m^{+} = (AP)^{+}(g − AN^{+}b) + N^{+}b.  (3.176)

Beweis: Vergleiche [11, Thm. 3.6.3].

□

Als einfaches Beispiel einer verallgemeinerten Lösung unter Nebenbedingung sei die reellwertige verallgemeinerte Lösung einer komplexwertigen Gleichung bezüglich eines C-linearen, beschränkten Operators zu bestimmen. Dazu seien im folgenden mit X_{C} und Y_{C} komplexe Hilberträume bezeichnet und mit X_{IR} := Re X_{C} und Y_{IR} := Re Y_{C} reelle Hilberträume definiert, die durch Einschränkung der Elemente aus X_{C} bzw. Y_{C} und der Multiplikation hervorgehen. Man beachte den Unterschied zur Definition aus Abschnitt 3.2.4.

Nun sei

B : X_{C} ⟶ Y_{C}, f ↦ Bf := (B_{r} + iB_{i})(f_{r} + if_{i})  (3.177)

ein C-linearer, beschränkter Operator aufgespaltet in Real- und Imaginärteil. Die Gleichung Bf = g werde geschrieben als

( \begin{matrix}B_{r} & −B_{i} \\ B_{i} & B_{r} \end{matrix} ) ( \begin{matrix}f_{r} \\ f_{i} \end{matrix} ) = ( \begin{matrix}g_{r} \\ g_{i} \end{matrix} ) ,  (3.178)

es seien also die Identifizierungen X_{C} ≃ (X_{IR} × X_{IR}) und Y_{C} ≃ (Y_{IR} × Y_{IR}) vorgenommen. Die Einschränkung einer Lösung auf den Realteil wird dann formuliert durch die Gleichung Nf = 0 ∈ X_{IR} mit der Definition

N : (X_{IR} × X_{IR}) ⟶ X_{IR}, (f_{r}, f_{i})^{t} ↦ Nf := ( 0 I ) ( \begin{matrix}f_{r} \\ f_{i} \end{matrix} ) ,  (3.179)

wobei I = I_{IR} die Identität auf X_{IR} bedeuten soll. Die Abbildung N^{+} berechnet sich zu

N^{+} = N^{∗}(NN^{∗})^{−1} = ( \begin{matrix}0 \\ I \end{matrix} ) (( 0 I ) ( \begin{matrix}0 \\ I \end{matrix} ))^{−1} = ( \begin{matrix}0 \\ I \end{matrix} ) ,  (3.180)

und damit

P = I_{X_{IR}×X_{IR}} − N^{+}N = ( \begin{matrix}I & 0 \\ 0 & I \end{matrix} ) − ( \begin{matrix}0 & 0 \\ 0 & I \end{matrix} ) = ( \begin{matrix}I & 0 \\ 0 & 0 \end{matrix} ) .  (3.181)

Somit ist die Lösung m^{+} ∈ (X_{IR} × {0}) der Gleichung (3.178) mit minimaler Norm ∥m∥_{X_{IR}×X_{IR}}, die den Defekt ∥Bm − g∥_{Y_{IR}×Y_{IR}} minimiert, gegeben durch

m^{+} = (BP)^{+}g = ( \begin{matrix}B_{r}^{t} & B_{i}^{t} \\ 0 & 0 \end{matrix} ) ( \begin{matrix}B_{r}B_{r}^{t} & B_{r}B_{i}^{t} \\ B_{i}B_{r}^{t} & B_{i}B_{i}^{t} \end{matrix} )^{−1} ( \begin{matrix}g_{r} \\ g_{i} \end{matrix} ) .  (3.182)

Man beachte, daß im allgemeinen

(BP)^{+} ≠ P^{+}B^{+} = P B^{+}  (3.183)

gilt, so daß der Realteil der verallgemeinerten Lösung f^{+} = B^{+}g nicht die verallgemeinerte Lösung m^{+} ergibt,

m_{r}^{+} ≠ (f^{+})_{r}.  (3.184)

Um dennoch die verallgemeinerte Lösung m ∈ (X_{IR} ×{0}) anhand der komplexwertigen verallgemeinerten Lösung f^{+} = (f_{r}^{+}, f_{i}^{+})^{t} in diesem Beispiel zu berechnen, betrachte man für die zu bestimmende Korrektur h = (h_{r}, 0)^{t} ∈ (X_{IR} × {0}) den Ansatz

( \begin{matrix}B_{r} & −B_{i} \\ B_{i} & B_{r} \end{matrix} ) ( \begin{matrix}f_{r}^{+} + h_{r} \\ 0 \end{matrix} ) = ( \begin{matrix}B_{r}f_{r}^{+} − B_{i}f_{i}^{+} \\ B_{i}f_{r}^{+} + B_{r}f_{i}^{+} \end{matrix} ) ,  (3.185)

also

( \begin{matrix}B_{r} & −B_{i} \\ B_{i} & B_{r} \end{matrix} ) ( \begin{matrix}h_{r} \\ h_{i} \end{matrix} ) = ( \begin{matrix}0 & −I \\ I & 0 \end{matrix} ) ( \begin{matrix}B_{r} & −B_{i} \\ B_{i} & B_{r} \end{matrix} ) ( \begin{matrix}f_{i}^{+} \\ 0 \end{matrix} ) = ( \begin{matrix}−B_{i}f_{i}^{+} \\ B_{r}f_{i}^{+} \end{matrix} ) (3.186)

mit der Nebenbedingung Nh = 0. Der nächste Satz zeigt, wie sich die beiden verallgemeinerten Lösungen m^{+}, h^{+} zu der verallgemeinerten Lösung f^{+} zusammensetzen.

Satz 3.2.7 Es sei f^{+} = (f_{r}^{+}, f_{i}^{+})^{t} ∈ (X_{IR} × X_{IR}) ≃ X_{C} die verallgemeinerte Lösung der Gleichung

Bf = g ∈ (Y_{IR} × Y_{IR}) ≃ Y_{C}  (3.187)

bezüglich eines C-linearen, beschränkten Operators B : X_{C} ⟶ Y_{C}. Weiter seien m^{+}, h^{+} ∈ (X_{IR} × {0}) die verallgemeinerten Lösungen der Gleichungen

Bm = g  bzw.  Bh = (−B_{i}f_{i}^{+}, B_{r}f_{i}^{+})^{t}  (3.188)

unter den Nebenbedingungen Nm = 0 bzw. Nh = 0, wobei N die Abbildung aus (3.179) sein soll.

Dann gilt

f^{+} = ( \begin{matrix}m_{r}^{+} − h_{r}^{+} \\ f_{i}^{+} \end{matrix} ) .  (3.189)

Beweis: Die verallgemeinerte Lösung f^{+} ∈ X_{C} der Gleichung Bf = g ist eindeutig als das Element f ∈ \overline{range(B^{∗})} bestimmt, welches die Gleichung

B^{∗}Bf = B^{∗}g  (3.190)

erfüllt. Da sich die verallgemeinerten Lösungen m^{+}, h^{+} ∈ (X_{IR} × {0}) anhand der Formeln

m^{+} = (BP)^{+}g,  h^{+} = (BP)^{+}(iB)f_{i}^{+}  (3.191)

mit P = (I − N^{+}N) berechnen, folgt für die rechte Seite von (3.189)

B^{∗}B(m_{r}^{+} − h_{r}^{+} + if_{i}^{+}) = B^{∗}B(m_{r}^{+} − h_{r}^{+}) + B^{∗}(iB)f_{i}^{+}
= B^{∗} \underbrace{B(BP)^{+}}_{=I_{C}} [g − (iB)f_{i}^{+}] + B^{∗}(iB)f_{i}^{+}

= B^{∗}g.  (3.192)

Da nach Definition m_{r}^{+}, h_{r}^{+} ∈ \overline{range(P B^{∗})} ist, gilt m_{r}^{+} − h_{r}^{+} + if_{i}^{+} ∈ \overline{range(B^{∗})}.

□

3.2.6 Sherman–Morrison–Woodbury–Formel

Es scheint, das Problem der Umrechnung einer komplexwertigen verallgemeinerten Lösung f^{+} auf die reellwertige m^{+} wird durch den Satz 3.2.7 durch die Einführung eines weiteren Minimierungsproblems mit den gleichen Nebenbedingungen nicht gelöst, sondern nur verlagert: Löse zunächst das komplexwertige Problem, um dann eine Korrektur durch eine verallgemeinerte Lösung h^{+} unter der Nebenbedingung Im(h) = 0 zu bestimmen.

Mit der Sherman–Morrison–Woodbury–Formel, siehe z.B. [19, Abschn. 2.1.3],

(C + UV^{t})^{−1} = C^{−1} − C^{−1} U (I + V^{t}C^{−1}U)^{−1} V^{t} C^{−1}  (3.193)

zur Invertierung einer modifizierten Matrix ergibt sich allerdings eine weitere Möglichkeit, die Korrektur h^{+} näherungsweise zu bestimmen. Dies stellt sich folgendermaßen dar.

Die verallgemeinerte Lösung f^{+} = (f_{r}^{+}, f_{i}^{+})^{t} der Gleichung Bf = g berechnet sich anhand der Formel f^{+} = B^{∗}(BB^{∗})^{−1}g, also mit den Definitionen

W := ( \begin{matrix}B_{r} \\ B_{i} \end{matrix} ) ,  V := ( \begin{matrix}−B_{i} \\ B_{r} \end{matrix} )  (3.194)

durch

( \begin{matrix}f_{r}^{+} \\ f_{i}^{+} \end{matrix} ) = ( \begin{matrix}B_{r}^{t} & B_{i}^{t} \\ −B_{i}^{t} & B_{r}^{t} \end{matrix} ) (( \begin{matrix}B_{r}B_{r}^{t} & B_{r}B_{i}^{t} \\ B_{i}B_{r}^{t} & B_{i}B_{i}^{t} \end{matrix} ) + ( \begin{matrix}B_{i}B_{i}^{t} & −B_{i}B_{r}^{t} \\ −B_{r}B_{i}^{t} & B_{r}B_{r}^{t} \end{matrix} ))^{−1} ( \begin{matrix}g_{r} \\ g_{i} \end{matrix} )

= ( \begin{matrix}B_{r}^{t} & B_{i}^{t} \\ −B_{i}^{t} & B_{r}^{t} \end{matrix} ) (W W^{t} + V V^{t})^{−1} ( \begin{matrix}g_{r} \\ g_{i} \end{matrix} ) .  (3.195)

Falls nun die Inverse (I − V^{t}D^{−1}V)^{−1} existiert, kann durch Anwendung der Formel (3.193) die Inverse für C := W W^{t} zur Bestimmung der verallgemeinerten Lösung m^{+} anhand der inversen Abbildung D^{−1} := (W W^{t} + V V^{t})^{−1} berechnet werden,

C^{−1} = (D − V V^{t})^{−1}
= D^{−1} + D^{−1} V (I − V^{t} D^{−1} V )^{−1} V^{t} D^{−1}.  (3.196)

Die gesuchte Korrektur h^{+} ergibt sich damit zu

h^{+} = ( \begin{matrix}h_{r}^{+} \\ 0 \end{matrix} ) = ( \begin{matrix}B_{r}^{t} & B_{i}^{t} \\ 0 & 0 \end{matrix} ) (W W^{t} + V V^{t})^{−1} ( \begin{matrix}−B_{i} h̃ \\ B_{r} h̃ \end{matrix} ) ,  (3.197)

wobei h̃ die Lösung der linearen Gleichung

(I − V^{t}D^{−1}V )h̃ = k  (3.198)

mit k := V^{t}D^{−1}g ist. Insgesamt erhält man für die verallgemeinerte Lösung m^{+} also

m^{+} = ( \begin{matrix}B_{r}^{t} & B_{i}^{t} \\ 0 & 0 \end{matrix} ) (W W^{t} + V V^{t})^{−1} ( \begin{matrix}g_{r} − B_{i} h̃ \\ g_{i} + B_{r} h̃ \end{matrix} ) .  (3.199)

Vorausgesetzt,es gilt ∥V^{t}D^{−1}V ∥ < 1, dann ist das System (3.198) eindeutig lösbar und die Iteration

h̃_{n} = k + V^{t}D^{−1}V h̃_{n−1}, h̃_{0} = 0  (3.200)

konvergiert für n → ∞ gegen die Lösung h̃. Dann ist mit (3.199) und der Iterationsvorschrift (3.200) ein Verfahren zur näherungsweisen Bestimmung der Korrektur h^{+} anhand der Formel (3.193) gegeben.

Übertragen auf die Abbildung B = (A_{(κ,μ)} ·)(φ) der Fréchet–Ableitung genommen an der Stelle (κ, μ)^{t} lautet die Gleichungen (3.199) zunächst

h̃ − (A_{(κ,μ)}^{∗} ˜ g̃)_{i}(φ) = (A_{(κ,μ)}^{∗} g̃)_{i}(φ)  (3.201)

mit (AA^{∗} ˜ g̃)(φ) = (i A_{(κ,μ)} h̃)(φ) und (AA^{∗} g̃)(φ) = g(φ). Ausführlich geschrieben ergibt sich dies zu

( \begin{matrix}h̃_{1} & − & Im(\overline{∇_{x}u_{φ}} · ∇_{x}z[h̃]) \\ h̃_{2} & − & Im(\overline{u_{φ}} z[h̃]) \end{matrix} ) = ( \begin{matrix}Im(\overline{∇_{x}u_{φ}} · ∇_{x}z[g̃]) \\ Im(\overline{u_{φ}} z[g̃]) \end{matrix} ) ,  (3.202)

wobei zur Vereinfachung z[h̃] := z[(AA^{∗})^{−1}(i A_{(κ,μ)} h̃)(φ)] definiert ist und z[g] dementsprechend die Lösung der adjungierten Gleichung

\overline{L_{(κ,μ)}}z = 0,  z|_{∂Ω} = g  (3.203)

zu den Randwerten g darstellt.

Damit erhält man für die Iteration (3.200) die Vorschrift

( \begin{matrix}h_{1} \\ h_{2} \end{matrix} )_{n} = ( \begin{matrix}Im(\overline{∇_{x}u_{φ}} · ∇_{x}(z[g̃] + z[(h_{1}, h_{2})_{n−1}^{t}])) \\ Im(\overline{u_{φ}} (z[g̃] + z[(h_{1}, h_{2})_{n−1}^{t}])) \end{matrix} ) , ( \begin{matrix}h_{1} \\ h_{2} \end{matrix} )_{0} = 0. (3.204)

Wird diese Iteration nach N Schritten abgebrochen, so gelangt man mit (3.199) zu der angenäherten verallgemeinerten Lösung m^{+} für einen Einzelschritt der noch anzugebenen Kaczmarz–Newton–Methode,

m_{r}^{+} = (A_{(κ,μ)}^{∗} (AA^{∗})^{−1}[g + (i A_{(κ,μ)} h_{N})(φ)])_{r}(φ).  (3.205)

Eine Iteration äquivalent zu (3.200), die nur auf den Daten g operiert, ergibt sich übrigens aus dem Ansatz

(D − V V^{t})^{−1} = D^{−1}(I − V V^{t}D^{−1})^{−1}.  (3.206)

Man erhält die Vorschrift

g_{n} = V V^{t}D^{−1}(g_{n−1} + g), g_{0} = 0,  (3.207)

die unter der Voraussetzung ∥V V^{t}D^{−1}∥ < 1 konvergiert, und gelangt damit zu der zu (3.205) äquivalenten Formel

m_{r}^{+} = (A_{(κ,μ)}^{∗} (AA^{∗})^{−1}[g + g_{N}])_{r}(φ).  (3.208)

Für numerische Berechnungen zur optischen Tomographie ist diese Form der Iteration gegenüber (3.200) vorzuziehen, da diese nur die Randdaten g benötigt.

Numerische Experimente zeigen, daß für die optische Tomographie die Iteration (3.204) schon nach einigen Schritten, etwa N < 20, stationär wird. Eine längere Iteration bringt somit keine wesentliche Verbesserung des Rekonstruktionsergebnisses. Darüberhinaus stellt sich heraus, daß die Anzahl der Iterationsschritte das Übersprechverhalten des Rekonstruktionsverhaltens beeinflußt, wie im numerischen Teil aufgezeigt wird.

3.2.7 Vorteil der komplexwertigen Rechnung

Zur Vereinfachung der Darstellung seien in diesem Unterabschnitt alle betrachteten Räume stets endlich-dimensional. Dies stellt insofern keine wesentliche Einschränkung dar, da zur Durchführung der praktischen Rechnungen auf endlichen Maschinen jedes Problem in ein endlich-dimensionales überführt werden muß.

Die verallgemeinerte Lösung eines linearen, unterbestimmten und komplexwertigen Problems

Ax = b  (3.209)

ist gegeben durch

x^{+} = A^{+}b = A^{∗}(AA^{∗})^{−1}b  (3.210)

oder, getrennt in Real- und Imaginärteil, durch

( \begin{matrix}x_{r}^{+} \\ x_{i}^{+} \end{matrix} ) = ( \begin{matrix}W^{t} \\ V^{t} \end{matrix} ) (W W^{t} + V V^{t})^{−1} ( \begin{matrix}b_{r} \\ b_{i} \end{matrix} ) ,  (3.211)

wobei

W := ( \begin{matrix}A_{r} \\ A_{i} \end{matrix} ) ,  V := ( \begin{matrix}−A_{i} \\ A_{r} \end{matrix} )  (3.212)

definiert ist.

Die verallgemeinerte Lösung m^{+} der Gleichung (3.209) unter der Nebenbedingung, daß der Imaginärteil verschwindet, lautet mit (3.212)

( \begin{matrix}m_{r}^{+} \\ 0 \end{matrix} ) = ( \begin{matrix}W^{t} \\ 0 \end{matrix} ) (W W^{t})^{−1} ( \begin{matrix}b_{r} \\ b_{i} \end{matrix} ) .  (3.213)

Es seien nun QR-Zerlegungen der Matrizen A^{∗} = (W V)^{t} und (W 0)^{t} gegeben, also

(W V)^{t} = ( \begin{matrix}A_{r}^{t} & A_{i}^{t} \\ −A_{i}^{t} & A_{r}^{t} \end{matrix} ) = Q ( \begin{matrix}R \\ 0 \end{matrix} )  (3.214)

und

(W 0)^{t} = ( \begin{matrix}A_{r}^{t} & A_{i}^{t} \\ 0 & 0 \end{matrix} ) = Q̃ ( R̃ ; 0 )  (3.215)

mit den orthogonalen Matrizen Q, Q̃ und rechten, oberen Dreicksmatrizen R, R̃. Dann folgt

R̃^{t}R̃ = W W^{t} = R^{t}R − V V^{t}.  (3.216)

Zur Abschätzung der Singulärwerte der Abbildungen R, R̃ ist der folgende Satz über die Eigenwertverteilung symmetrischer Matrizen hilfreich.

Satz 3.2.8 Für die der Größe nach geordneten Eigenwerte λ_{k}(A) und λ_{k}(A + B), k = 1, …, n, zweier symmetrischer, reeller (n × n)-Matrizen A und A + B gilt

λ_{k}(A) + λ_{n}(B) ≤ λ_{k}(A + B) ≤ λ_{k}(A) + λ_{1}(B), k = 1, …, n.  (3.217)

Beweis:

Zu finden in [68, S. 101f]: Die Aussage wird durch das Courantsche Maximum-Minimum-Prinzip bewiesen,

λ_{k}(A) = max_{dim(S)=k} min_{0≠y∈S} \frac{y^{t}Ay}{y^{t}y} ,  (3.218)

siehe [19, Thm. 8.1.2].

□

Da für die Singulärwerte σ̃_{j} = σ_{j}(R̃) der Abbildung R̃ die Gleichung

σ̃_{j}^{2} = σ_{j}^{2}(R̃) = λ_{j}(R̃^{t}R̃) = λ_{j}(W W^{t})  (3.219)

und für die Singulärwerte σ_{j} = σ_{j}(R) der Abbildung R die Gleichung

σ_{j}^{2} = σ_{j}^{2}(R) = λ_{j}(R^{t}R) = λ_{j}(W W^{t} + V V^{t})  (3.220)

gilt, erhält man mit dem letzten Satz die Ungleichungen

σ_{1}^{2} ≤ σ̃_{1}^{2} + λ_{1}(V V^{t}),  σ̃_{m}^{2} + λ_{m}(V V^{t}) ≤ σ_{m}^{2},  (3.221)

wobei m der Index des kleinsten Singulärwertes ≠ 0 ist. Wegen des Maximum–Minimum– Prinzips (3.218) gilt nun λ_{k}(V V^{t}) = λ_{k}(W W^{t}) = σ̃_{k}^{2}. Es folgt

\frac{σ_{1}}{σ_{m}} ≤ \frac{σ̃_{1}}{σ̃_{m}} .  (3.222)

Daher steht zu erwarten, daß die numerische Bestimmung der komplexwertigen verallgemeinerten Lösung x^{+} besser möglich ist als die Berechnung der reellwertigen Lösung m^{+}, da diese wegen (3.222) ein schlechter gestelltes Problem darstellt.

Als eine nicht ganz ernst gemeinte Randbemerkung zu diesem Ergebnis sei ein Zitat von Kronecker, vgl. [33, S. 84]^{1}, angemerkt, dessen Aussage allerdings in bezug auf die Anwendbarkeit des Residuenkalküls zur Auswertung uneigentlicher Integrale getroffen wurde:

,,Wir brauchen hier übrigens zum Zwecke dieser Beweise das Gebiet der reellen Größen nicht zu verlassen. Der Glaube an die Unwirksamkeit des Imaginären trägt auch hier wie anderweitig gute Früchte.”

1Remmert, R.: Funktionentheorie I, Springer, 1989, Abschn. 14.2.3

3.2.8 Tikhonov–Regularisierung

Für kompakte Operatoren A : X ⟶ Y, X und Y Hilberträume, gilt, daß im allgemeinen die zugehörige inverse Abbildung nicht stetig ist, sofern sie überhaupt existiert. Da mit A auch der Operator A^{∗}A kompakt ist, stellt die Lösung der Normalgleichung ein schlecht gestelltes Problem dar. Der Begriff der Regularisierung bedeutet nun in diesem Zusammenhang, den Mangel an Stetigkeit bei der Inversenbildung bezüglich der Normalgleichung durch die folgende Definition zu beheben.

Unter einer Regularisierung der Abbildung A^{+}, welche die Daten g der Lösung f^{+} ∈ \overline{range(A^{∗})} der Gleichung

A^{∗}Af = A^{∗}g  (3.223)

zuordnet, versteht man allgemein eine Familie (T_{α})_{α>0} linearer, beschränkter Operatoren T_{α} : Y ⟶ X mit der Eigenschaft, daß für g ∈ range(A)+range(A)^{⊥} der Grenzwert

lim_{α→0} T_{α}g = f^{+}  (3.224)

ist.

Die Tikhonov–Phillips–Regularisierung wird nun durch die Familie

T_{α} = (A^{∗}A + αI)^{−1}A^{∗} = A^{∗}(AA^{∗} + αI)^{−1}, α > 0  (3.225)

beschrieben, [47, IV.1.2]. Für jedes α > 0 bildet T_{α} also die Daten g ∈ Y eindeutig auf das Element ab, welches statt (3.223) die Gleichung

(A^{∗}A + αI)f = A^{∗}g,  f ∈ X  (3.226)

erfüllt. Wie man leicht zeigt, ist zu dieser Definition äquivalent die Definition der Abbildung T_{α} durch das Element, welches das Funktional

J_{α}(f) := ∥Af − g∥_{Y}^{2} + α ∥f∥_{X}^{2}  (3.227)

minimiert.

Entsprechend ergibt sich mit β > 0 und der Definition P := I − N^{+}N gemäß (3.181) für die reellwertige verallgemeinerte Lösung m^{+} die regularisierende Gleichung

((BP)^{∗}(BP) + βI)m = (BP)^{∗}g,  m ∈ X_{IR} × {0},  (3.228)

die genau das Element m ∈ X_{IR} × {0} als eindeutige Lösung besitzt, welches das Funktional

J_{β}(m) := ∥(BP)m − g∥_{Y_{IR}×Y_{IR}}^{2} + β ∥m∥_{X_{IR}×X_{IR}}^{2}  (3.229)

in X_{IR} × {0} minimiert.

Mit den Lösungsformeln (3.211) und (3.213) überträgt sich dies wegen

((BP)(BP)^{∗} + βI)^{−1} = (D + βI − V V^{t})^{−1}
= (D + βI)^{−1} +  (3.230)

(D + βI)^{−1}V (I − V^{t}(D + βI)^{−1}V )^{−1}V^{t}(D + βI)^{−1}, wobei D = W W^{t} + V V^{t} = BB^{∗} ist, auf die Iterationsvorschrift (3.200) wie folgt,

h̃_{n} = k + V^{t}(D + βI)^{−1}V h̃_{n−1}, h̃_{0} = 0  (3.231)

und k := V^{t}(D + βI)^{−1}g. Entsprechend ergibt sich

g̃_{n} = V V^{t}(D + βI)^{−1} (g̃_{n−1} + g), g̃_{0} = 0.  (3.232)

Somit wird unter der Voraussetzung, daß ∥V^{t}(D + βI)^{−1}V ∥ < 1 beziehungsweise ∥V V^{t}(D + βI)^{−1}∥ < 1 gilt, durch

m_{β,N} = P^{∗}B^{∗}(D + βI)^{−1}(g + V h̃_{N})
= P^{∗}B^{∗}(D + βI)^{−1}(g + g̃_{N})  (3.233)

eine angenäherte Lösung für die Minimumstelle m ∈ X_{IR} × {0} des Funktionals J_{β}(m)

gegeben, ohne die Inverse ((BP)(BP)^{∗} + βI)^{−1} bilden zu müssen.

3.2.9 Kaczmarz–Newton–Methode

Die relaxierte Kaczmarz–Methode für ein System

( A_{1}f ; A_{2}f ; ⋮ ; A_{p}f ) = ( g_{1} ; g_{2} ; ⋮ ; g_{p} )  (3.234)

mit surjektiven linearen Operatoren A_{j} : X ⟶ Y ist durch

f^{(k+1)} = (P_{p}^{ϱ} ∘ … ∘ P_{1}^{ϱ})f^{(k)},  k = 1, 2, …  (3.235)

und f^{(0)} ∈ X beliebig gegeben, wobei 0 < ϱ < 2 der Relaxationsparameter ist, und P_{j}^{ϱ} = (1 − ϱ) I + ϱ P_{j} die entsprechend relaxierte Form der orthogonalen Projektion P_{j} : X ⟶ U_{j} auf den affinen Unterraum

U_{j} = {x ∈ X, A_{j}x = g_{j}},  j = 1, …, p,  (3.236)

ist. Man rechnet nach, daß sich die orthogonale Projektion P_{j} : X ⟶ U_{j} ergibt zu

x ↦ P_{j}x = x − A_{j}^{∗}[A_{j}A_{j}^{∗}]^{−1}(A_{j}x − g_{j}).  (3.237)

Diese geometrische Interpretation der Kaczmarz–Methode ermöglicht sehr einfach die Erweiterung der Methode auf den Fall eines Systems

( R_{1}(f) ; R_{2}(f) ; ⋮ ; R_{p}(f) ) = 0  (3.238)

mit nichtlinearen, Fréchet–differenzierbaren Abbildungen R_{j} : X ⟶ Y, indem man für die Menge U_{j} den affinen Unterraum

V_{j} = { x ∈ X, R_{j}(f) + R′_{j}(f)(x − f) = 0 }  (3.239)

einsetzt, wobei R′_{j}(f) die Fréchet–Ableitung an der Stelle f ist. Der affine Unterraum V_{j} hängt natürlich von der Stelle f ab, an der die Ableitung genommen wird. Für jedes feste f ∈ X ist nach (3.237) die orthogonale Projektion Q_{j} : X ⟶ V_{j} durch

x ↦ Q_{j}x = x − R′_{j}(x)^{∗}[R′_{j}(x)R′_{j}(x)^{∗}]^{−1}R_{j}x  (3.240)

gegeben. Das Element Q_{j}x genügt somit der linearen Gleichung R′_{j}(f)(x−f) = −R_{j}(f), welche ihrerseits die lineare Näherung der Teilabbildung R_{j}(x) = 0 an der Stelle f darstellt. Dieser Ansatz kommt nach Wissen des Autors erstmalig 1955 in [66] vor.

Mit f_{0} = f^{(k)} und f_{j} = Q_{j}^{ϱ}f_{j−1}, wobei Q_{j}^{ϱ} = (1 − ϱ) I + ϱ Q_{j} ist, läßt sich dann der Iterationsschritt (3.235) übertragen auf das System (3.238) schreiben als

f^{(k+1)} = (Q_{p}^{ϱ} ∘ … ∘ Q_{1}^{ϱ})f^{(k)},  (3.241)

wobei der Startwert f^{(0)} ∈ X beliebig ist und 0 < ϱ < 2 erneut den Relaxationsparameter darstellt. Ausführlich erhält man aus (3.241) mit B_{j} := R′_{j}(f_{j−1}) den Algorithmus

f_{0} = f^{(k)},
f_{j} = f_{j−1} − ϱ B_{j}^{∗}[B_{j}B_{j}^{∗}]^{−1}R_{j}(f_{j−1}),  j = 1, …, p,  (3.242)
f^{(k+1)} = f_{p}

für den Iterationschritt f^{(k)} ⇝ f^{(k+1)}.

Man beachte, daß der Einzelschritt (3.242) auch durch das Lösen der Normalgleichung

B_{j}^{∗}B_{j}h = −B_{j}^{∗}R_{j}(f_{j−1}),  h ∈ \overline{range(B^{∗})}  (3.243)

beschrieben wird. Ist diese Gleichung nicht stetig invertierbar, so erhält man mit α > 0 aus der Gleichung

(B_{j}^{∗}B_{j} + α I)h = −B_{j}^{∗}R_{j}(f_{j−1}),  h ∈ X  (3.244)

eine regularisierte verallgemeinerte Lösung für h. Somit ergibt sich aus (3.241) das regularisierte Verfahren

f_{0} = f^{(k)},
f_{j} = f_{j−1} − ϱ B_{j}^{∗}[B_{j}B_{j}^{∗} + α I]^{−1}R_{j}(f_{j−1}),  j = 1, …, p, (3.245)
f^{(k+1)} = f_{p}

für einen Schritt f^{(k)} ⇝ f^{(k+1)} zur Lösung des nichtlinearen Gleichungssystems (3.238). Unter der Annahme, daß sich der Operator B_{j}B_{j}^{∗} = R′_{j}(f_{j−1}) ∘ R′_{j}(f_{j−1})^{∗} während der Iteration nicht wesenlich verändert, ergibt sich mit C_{j}C_{j}^{∗} := R′_{j}(f^{(0)}) ∘ R′_{j}(f^{(0)})^{∗} schließlich eine deutlich weniger aufwendige Version des Verfahrens (3.245) durch die Rechenvorschrift

f_{0} = f^{(k)},
f_{j} = f_{j−1} − ϱ B_{j}^{∗}[C_{j}C_{j}^{∗} + α I]^{−1}R_{j}(f_{j−1}),  j = 1, …, p, (3.246)
f^{(k+1)} = f_{p},

da die Abbildungen C_{j}C_{j}^{∗} unabhängig von der Iteration vorberechnet werden können.

Tritt nun, wie in der optischen Tomographie, der Fall ein, daß die gegebenen Daten g_{j} = R_{j}(f) zwar komplexwertig, die gesuchten Größen f = (f_{1}, …, f_{m})^{t} jedoch rein reellwertig sind, so sind die Verfahren (3.242) bzw. (3.245) sicherlich unbrauchbar, da offensichtlich in (3.242) bzw. (3.245) nicht die dem Problem entsprechenden orthogonalen Projektionen berechnet werden. Auch die Einschränkung der Formeln (3.242) bzw. (3.245) auf den Realteil alleine würde wegen (3.183) diesen Mißstand nicht beheben.

Wie jedoch am Ende des Abschnitts 3.2.5 aufgezeigt, ist eine näherungsweise Korrektur der Projektion durch die Verwendung der Sherman–Morrison–Woodbury–Formel möglich. Dies wird auch sehr schön durch Ergebnisse numerischer Experimenten, die mit einem derart modifizierten Schema durchgeführt wurden, im nächsten Kapitel bestätigt.

3.2.10 Konvergenz

Für eine Konvergenzaussage über das regularisierte Kaczmarz–Verfahren bezüglich des linearen Systems (3.234) sei zunächst S_{j} ein positiv definiter Operator mit der Eigenschaft

S_{j} ≥ A_{j}A_{j}^{∗},  j = 1, …, p,  (3.247)

d.h. ⟨S_{j}x, y⟩_{Y} ≥ ⟨A_{j}A_{j}^{∗}x, y⟩_{Y} für alle x ∈ X und y ∈ Y. Mit diesem S_{j} stelle man den Algorithmus

f_{0} = f^{(k)},
f_{j} = f_{j−1} − ϱ A_{j}^{∗} S_{j}^{−1}(A_{j}f_{j−1} − g_{j}),  j = 1, …, p,  (3.248)
f^{(k+1)} = f_{p}

auf. Es ergibt sich folgende Aussage.

Lemma 3.2.1 Angenommen, das System (3.234) besitze eine Lösung.

Dann konvergiert die durch (3.248) beschriebene Iterationsfolge (f^{(k)})_{k} für k → ∞ gegen die Lösung des Systems (3.234) mit minimaler Norm, falls 0 < ϱ < 2 und f^{(0)} ∈ range(A^{∗}) = range(A_{1}^{∗}) + … + range(A_{p}^{∗}) gewählt werden.

Beweis: Zu finden in [47, Thm. 3.6] für S_{j} = A_{j}A_{j}^{∗} bzw. [49, Thm. 4.1] für S_{j} ≥ A_{j}A_{j}^{∗}.

□

Über derartig allgemeine, analoge Konvergenzaussagen der Verfahren bezüglich der nichtlinearen Gleichungen (3.238) ist bisher wenig bekannt, siehe z.B. [42], [41], [66].

KAPITEL 3. MATHEMATISCHER TEIL 56

Numerischer Teil

In diesem Kapitel soll die konkrete Umsetzung der bisher hergeleiteten Rekonstruktionsverfahren beschrieben werden und die theoretisch vorhergesagten Resultate anhand einiger Beispiele numerisch bestätigt werden.

4.1 Implementation

Zur iterativen Lösung des inversen, nichtlinearen Problems der optischen Tomographie ist es notwendig, geeignete Löser für die Modellgleichung bereit zu stellen. Dies geschehe durch die Routinen diffuse und back diffuse.

Es seien p Randverteilungen φ_{j}, j = 1, …, p, gegeben, die jeweils in einem Quellpunkt s_{j} ∈ ∂Ω konzentriert sind.

Die Routine diffuse(κ, μ, s_{j}, ω_{0}) stehe für die approximative Lösung u der Diffusionsgleichung

−∇_{x} · (f_{1}(x) ∇_{x}u) + (f_{2}(x) + i ω)u = 0  in Ω,  (4.1)
u = φ  auf ∂Ω

mit den reellwertigen Koeffizienten f_{1} = κ und f_{2} = μ bezüglich der Randverteilung φ = δ_{s_{j}} und des konstanten Imaginärteils ω = ω_{0}.

Der Aufruf back diffuse(κ, μ, r_{j}, ω_{0}) soll entsprechend eine numerische Lösung z der adjungierten Gleichung

−∇_{x} · (f_{1}(x) ∇_{x}z) + (f_{2}(x) − i ω)z = 0  in Ω,  (4.2)
z = g  auf ∂Ω

mit den reellwertigen Koeffizienten f_{1} = κ und f_{2} = μ, dem konstanten Imaginärteil ω = ω_{0} und den Randwerten g = r_{j} liefern.

KAPITEL 4. NUMERISCHER TEIL

Desweiteren sei mit

C_{j} = ( \begin{matrix}B_{r}B_{r}^{t} & B_{r}B_{i}^{t} \\ B_{i}B_{r}^{t} & B_{i}B_{i}^{t} \end{matrix} )  (4.3)

die diskretisierte Fassung der jeweiligen linearen Abbildung (BP)(BP)^{∗} vorhanden, wobei die Abbildung B = B_{r} + i B_{i} = (A_{(κ_{0},μ_{0})} ·)(φ_{j}) die Fréchet-Ableitung bezüglich der Randverteilung φ_{j} genommen an der Stelle der Startwerte (κ_{0}, μ_{0})^{t} ist und die Abbildung P = I − N^{+}N für den Operator aus Abschnitt 3.2.5 steht.

Entsprechend sei die Diskretisierung der C-linearen Abbildung BB^{∗} als komplexwertige Matrix mit

D_{j} = (A_{(κ_{0},μ_{0})}A_{(κ_{0},μ_{0})}^{∗} ·)(φ_{j})  (4.4)

bezeichnet.

Es ist anzumerken, daß sich die Blöcke der Matrix C_{j} anhand der Zerlegung aus Abschnitt 3.2.4 einfach berechnet werden können unter der Verwendung einer Routine frechet(κ_{0}, μ_{0}, h_{κ}, h_{μ}, u_{j}, ω_{0}), die eine angenäherte Lösung w der Gleichung

−∇_{x} · (κ_{0}(x) ∇_{x}w) + (μ_{0}(x) + i ω_{0})w = ∇_{x} · (h_{κ}(x) ∇_{x}u_{j}) − h_{μ}(x)u_{j} in Ω, (4.5)
w = 0  auf ∂Ω

liefert. Dabei wird angenommen, daß u_{j} = diffuse(κ_{0}, μ_{0}, s_{j}, ω_{0}) gilt. Die Ableitung wird dann durch die Verteilung

(A_{(κ_{0},μ_{0})} h)(φ_{j}) = (\frac{∂_{ν}w}{∂_{ν}v_{j}} ) |_{∂Ω}  (4.6)

gegeben.

Schließlich sei bezüglich jeder Randverteilung φ_{j}, konzentriert im Quellpunkt s_{j}, die m-te Komponente des Datenvektors durch

d_{j,m} = d_{j}(y_{m}) = \frac{∂_{ν}u_{j}^{true}(y_{m})}{∂_{ν}v_{j}(y_{m})} − 1,  y_{m} ∈ ∂Ω,  (4.7)

für m = 1, …, M gegeben, wobei v_{j} die Lösung der Gleichung

−∇_{x} · (κ_{0}(x) ∇_{x}u) + (μ_{0}(x) + i ω_{0})v = 0  in Ω,  (4.8)
v = φ_{j} auf ∂Ω

mit den (nicht notwendig konstanten) Startnäherungen κ = κ_{0} und μ = μ_{0} ist.

Damit läßt sich das regularisierte und relaxierte Rekonstruktionsverfahren wie folgt formulieren.

4.1 IMPLEMENTATION

Algorithmus 4.1

1. Lese die Startwerte κ = κ_{0} und μ = μ_{0} ein und setze Werte für ϱ, α, N, ε und ω_{0} fest.

2. Für jede Quellposition s_{j}

(a) lese den Datenvektor d_{j} und die Matrizen C_{j} und D_{j} ein,

(b) v_{j} = diffuse(κ_{0}, μ_{0}, s_{j}, ω_{0}).

3. Für jede Quellposition s_{j} führe aus

(a) u_{j} = diffuse(κ, μ, s_{j}, ω_{0}),

(b) r_{j} = ∂_{ν}(u_{j} − v_{j})/(∂_{ν}v_{j}) − d_{j},

(c) r̃_{j} = 0,

(d) für n = 1 bis N berechne /* N = 0 bedeutet leere Schleife */

i. t = (D_{j} + α I)^{−1}(r_{j} + r̃_{j}), ii. Re(r̃_{j}) ← (C_{j})_{22} Re(t) − (C_{j})_{21} Im(t), iii. Im(r̃_{j}) ← (C_{j})_{11} Im(t) − (C_{j})_{12} Re(t),

(e) ˜ r̃_{j} = (D_{j} + α I)^{−1}(r_{j} + r̃_{j}),

(f) z = back diffuse(κ, μ, ˜ r̃_{j}, ω_{0}),

(g) κ ← κ + ϱ Re(\overline{∇u}_{j} · ∇z),

(h) μ ← μ + ϱ Re(\overline{u}_{j}z),

4. Falls (∥r_{1}∥ + ∥r_{2}∥ + … + ∥r_{p}∥) > εp, gehe zu Schritt 3.

Hierbei ist die Reihenfolge der Quellpositionen s_{j} in Schritt 3 nach Möglichkeit zufällig zu wählen. Dies bestimmt die Reihenfolge der Unterräume, auf die im Verlauf des Kaczmarz–Newton–Verfahren projiziert wird. Eine feste Reihenfolge bedeutet daher, daß die Rekonstruktionen jeweils nach einem Durchlauf des Schritts 3 stets Element des gleichen Unterraumes sind. Eine Analyse des Kaczmarz–Verfahrens angewendet für die Inversion der Radon–Transformation zur Berechnung von CT-Rekonstruktionen zeigt, daß sich genau dieser Umstand ungünstig auf das Konvergenzverhalten des Kaczmarz– Verfahrens auswirkt, siehe [47, Abschn. V.4.3].

Die reine reelle Rekonstruktionsmethode ergibt sich aus Algorithmus 4.1 einfach, indem N = 0 gesetzt wird und der Schritt (e) ersetzt wird durch die beiden Teilschritte

(e′_{1}) Re(˜r̃_{j}) = Re((C_{j} + α I)^{−1}(Re(r_{j}), Im(r_{j}))^{t}),

(e′_{2}) Im(˜r̃_{j}) = Im((C_{j} + α I)^{−1}(Re(r_{j}), Im(r_{j}))^{t}).

Für die rein reelle Methode werden die Matrizen D_{j} nicht benötigt und die Schleife (d) nicht durchlaufen.

Eine Variante dieses Algorithmus ergibt sich, indem diffuse und back diffuse durch

die Routinen diffuse schroedinger und back diffuse schroedinger ersetzt werden

und dementsprechend über das komplexe Potential q(x) = q_{r}(x) + i q_{i}(x), das sich gemäß Abschnitt 3.1.8 aus κ(x), μ(x) und ω_{0} berechnet, iteriert wird. Die Routine diffuse schroedinger(q, s_{j}) liefere demnach die Lösung u der Gleichung

−Δ_{x}u + q(x)u = 0  in Ω,  (4.9)
u = φ auf ∂Ω

mit der Randverteilung φ = δ_{s_{j}}.

Entsprechend berechne die Routine back diffuse schroedinger(q, r_{j}) die Lösung z der zugehörigen adjungierten Gleichung

−Δ_{x}z + \overline{q(x)}z = 0  in Ω,  (4.10)
z = g auf ∂Ω

bezüglich der Randwerte g = r_{j}.

Damit ergibt sich der folgende, verblüffend einfache Algorithmus.

Algorithmus 4.2

1. Lese den Startwert q = q_{0} ein und setze Werte für ϱ, α und ε fest.

2. Für jede Quellposition s_{j}

(a) lese den Datenvektor d_{j} und die Matrix D_{j} ein,

      (b) v_{j} = diffuse schroedinger(q_{0}, s_{j}).

3. Für jede Quellposition s_{j} führe aus

      (a) u_{j} = diffuse schroedinger(q, s_{j}),

(b) r_{j} = ∂_{ν}(u_{j} − v_{j})/(∂_{ν}v_{j}) − d_{j},

(c) r̃_{j} = (D_{j} + α I)^{−1} r_{j},

      (d) z = back diffuse schroedinger(q, r̃_{j}),

(e) Re(q) ← Re(q) + ϱ Re(\overline{u}_{j}z),

(f) Im(q) ← Im(q) + ϱ Im(\overline{u}_{j}z),

4. Falls (∥r_{1}∥ + ∥r_{2}∥ + … + ∥r_{p}∥) > εp, gehe zu Schritt 3.

4.1.1 Finite Differenzen

Man kann zwei prinzipiell voneinander verschiedene Methoden zur Berechung approximativer Lösungen partieller Differentialgleichungen unterscheiden: Bildung finiter Differenzen für die auftretenen Differentialausdrücke und die Methode der Finiten Elemente.

Die Methode der Finiten Elemente besitzt eine größere Flexibilität bei der Anpassung an komplizierte Geometrien der Grundgebiete, auf denen die Differentialgleichung gelöst werden sollen. Da es unter anderem in dieser Arbeit jedoch um die grundsätzliche Realisierbarkeit des vorgestellten iterativen Verfahrens geht, reicht es völlig aus, Beispiele für einfache Geometrien wie Rechtflach oder Kugel zu betrachen. Daher wurde von der Methode der Finiten Elemente abgesehen und Routinen mittels der finiten Differenzen entwickelt. Die Diskretisierung bezüglich kegelartiger Gebiete in IR^{3} sind nicht Standard, und sollen daher hier mit aufgeführt werden, siehe auch [65, Abschn. 4.5].

4.1.2 Orthogonale Geometrie

Es sei ein beschränkter Quader Ω ⊂ IR^{3} gegeben und mit einem äquidistanten Gitter

x_{j} = x_{0} + jh, j = 0, …, J,
y_{k} = y_{0} + kh, k = 0, …, K,  (4.11)
z_{l} = z_{0} + lh, l = 0, …, L,

durchzogen. Dann können die partiellen Ableitungen u_{x}, u_{y}, u_{z} an der Stelle (x_{j}, y_{k}, z_{l}) einer genügend oft differenzierbaren Funktion u(x, y, z) durch

u_{x} = \frac{1}{2h}(u_{j+1,k,l} − u_{j−1,k,l})  (4.12)
u_{y} = \frac{1}{2h}(u_{j,k+1,l} − u_{j,k−1,l})  (4.13)
u_{z} = \frac{1}{2h}(u_{j,k,l+1} − u_{j,k,l−1})  (4.14)

und entsprechend die zweiten Ableitungen u_{xx}, u_{yy}, u_{zz} durch

u_{xx} = \frac{1}{h^{2}} (u_{j+1,k,l} − 2u_{j,k,l} + u_{j−1,k,l})  (4.15)
u_{yy} = \frac{1}{h^{2}} (u_{j,k+1,l} − 2u_{j,k,l} + u_{j,k−1,l})  (4.16)
u_{zz} = \frac{1}{h^{2}} (u_{j,k,l+1} − 2u_{j,k,l} + u_{j,k,l−1})  (4.17)

angenähert werden, wobei u_{j,k,l} := u(x_{j}, y_{k}, z_{l}) zur Abkürzung definiert sei.

Damit ergeben sich für die Lösung der Diffusionsgleichung

− ∇ · (κ(x, y, z)∇u) + (μ(x, y, z) + i ω_{0})u = 0  (4.18)

die linearen Gleichungen

−\frac{1}{4h^{2}} [(κ_{j+1} + 4κ_{j} − κ_{j−1})u_{j+1} + (−κ_{j+1} + 4κ_{j} + κ_{j−1})u_{j−1}

+(κ_{k+1} + 4κ_{k} − κ_{k−1})u_{k+1} + (−κ_{k+1} + 4κ_{k} + κ_{k−1})u_{k−1}  (4.19)
+(κ_{l+1} + 4κ_{l} − κ_{l−1})u_{l+1} + (−κ_{l+1} + 4κ_{l} + κ_{l−1})u_{l−1}
−(4h^{2}(μ_{j,k,l} + i ω_{0}) + 24κ_{j,k,l})u_{j,k,l}] = 0

als Näherung in den Gitterpunkten (x_{j}, y_{k}, z_{l}) für j = 1, …, (J − 1), k = 1, …, (K − 1) und l = 1, …, (L − 1). In (4.19) wurde dabei auf Indizes, die für die jeweiligen zwei übrigen Koordinaten unwesentlich sind, der Übersichtlichkeit halber verzichtet.

Falls in (4.19) für die Indizes j ∈ {1, (J −1)} oder k ∈ {1, (K −1)} oder l ∈ {1, (L−1)} zutrifft, dann sind in (4.19) die Terme mit j ∈ {0, J} oder k ∈ {0, K} oder l ∈ {0, L} durch die Randbedingungen

u_{j,k,l} = φ_{j,k,l}  für j ∈ {0, J} oder k ∈ {0, K} oder l ∈ {0, L}  (4.20)

bekannt. Somit stellt (4.19) ein lineares Gleichungssystem mit (J − 1)(K − 1)(L − 1) Gleichungen für die (J − 1)(K − 1)(L − 1) Unbekannten u_{j,k,l} dar, welches mit den bekannten Methoden entweder direkt (z.B. LR-Zerlegung) oder iterativ (z.B. GMRES) gelöst werden kann, siehe [9], [19].

Die numerische Lösung der adjungierten Gleichung

− ∇ · (κ(x, y, z)∇u) + (μ(x, y, z) − i ω_{0})u = 0  (4.21)

mit den Randwerten

u_{j,k,l} = g_{j,k,l}  für j ∈ {0, J} oder k ∈ {0, K} oder l ∈ {0, L}  (4.22)

durch das finite Differenzenschema (4.19) mit ω_{0} ersetzt durch −ω_{0} kann somit vollkommen analog zu den vorherigen Ausführungen erfolgen.

4.1.3 Zylindrische Geometrie

Hier sei ein kegelzylindrischer beziehungsweise rotationssymmetrischer Bereich Z ⊂ IR^{3} in den Koordinaten Radius r, Winkel φ und Höhe z gegeben. Das diskrete Gitter sei wie folgt aufgebaut,

r_{j}(z) = f(z)jh_{r}, φ_{k} = kh_{φ}, z_{l} = z_{0} + lh_{z},  (4.23)

wobei f(z) eine nullstellenfreie, stetig differenzierbare Funktion bedeute. Zunächst soll jetzt der Divergenzoperator −∇_{x} · (κ∇_{x}) in diese Koordinaten umgerechnet werden. Man betrachte also die Abbildung

Φ : Z′ ⟶ Z, ( \begin{matrix}r \\ φ \\ z \end{matrix} ) ↦ ( \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ) = ( \begin{matrix}f(z) r cos φ \\ f(z) r sin φ \\ z \end{matrix} ) .  (4.24)

Die Jacobi-Matrix dieser Abbildung ergibt sich zu

(∂_{r,φ,z}Φ) = ( \begin{matrix}f(z) cos φ & −f(z)r sin φ & f′(z)r cos φ \\ f(z) sin φ & f(z)r cos φ & f′(z)r sin φ \\ 0 & 0 & 1 \end{matrix} )  (4.25)

mit der Determinaten

\sqrt{g} = | det(∂_{r,φ,z}Φ)| = f(z)^{2}r.  (4.26)

Der zugehörige Maßtensor lautet

(g_{pq}) = (∂_{r,φ,z}Φ)^{t}(∂_{r,φ,z}Φ) = ( \begin{matrix}f(z)^{2} & 0 & f(z)f′(z)r \\ 0 & f(z)^{2}r^{2} & 0 \\ f(z)f′(z)r & 0 & 1 + f′(z)^{2}r^{2} \end{matrix} ) ,  (4.27)

der entsprechende inverse Tensor

(g^{pq}) = \frac{1}{f(z)^{2}} ( \begin{matrix}1 + f′(z)^{2}r^{2} & 0 & −f(z)f′(z)r \\ 0 & 1/r^{2} & 0 \\ −f(z)f′(z)r & 0 & f(z)^{2} \end{matrix} ) .  (4.28)

Damit transformiert sich der Differentialoperator der Diffusionsgleichung für die Koordinaten y = (y_{1}, y_{2}, y_{3}) = (r, φ, z) zu

−∇_{x} · (κ(x)∇_{x}) + (μ(x) + i ω_{0})
= −\frac{1}{\sqrt{g}} ∑_{p,q=1}^{3} ∂_{y_{p}} [(κ ∘ Φ)(r, φ, z)g^{pq}\sqrt{g} ∂_{y_{q}}] + ((μ ∘ Φ)(r, φ, z) + i ω_{0})
= −\frac{1}{f(z)^{2}r} [∂_{r}[(κ ∘ Φ)(r, φ, z)r(1 + f′(z)^{2}r^{2})∂_{r}]
− f(z)f′(z)∂_{r}[(κ ∘ Φ)(r, φ, z)r^{2}∂_{z}]  (4.29)

+ \frac{1}{r}∂_{φ}[(κ ∘ Φ)(r, φ, z)∂_{φ}]

− r^{2}∂_{z}[(κ ∘ Φ)(r, φ, z)f(z)f′(z)∂_{r}]

+ r∂_{z}[(κ ∘ Φ)(r, φ, z)f(z)^{2}∂_{z}]] + ((μ ∘ Φ)(r, φ, z) + i ω_{0}).

Mit den finiten Differenzenapproximationen – die unwesentlichen Indizes sind der Übersichtlichkeit halber fortgelassen –

∂_{r}[D(r, φ, z)∂_{r}u] =  (4.30)

\frac{1}{2h_{r}^{2}} ((D_{j+1} + D_{j})u_{j+1} + (D_{j+1} + 2D_{j} + D_{j−1})u_{j} + (D_{j} + D_{j−1})u_{j−1}),

∂_{φ}[D(r, φ, z)∂_{φ}u] =  (4.31)

\frac{1}{2h_{φ}^{2}} ((D_{k+1} + D_{k})u_{k+1} + (D_{k+1} + 2D_{k} + D_{k−1})u_{k} + (D_{k} + D_{k−1})u_{k−1}),

∂_{z}[D(r, φ, z)∂_{z}u] =  (4.32)

\frac{1}{2h_{z}^{2}} ((D_{l+1} + D_{l})u_{l+1} + (D_{l+1} + 2D_{l} + D_{l−1})u_{j} + (D_{l} + D_{l−1})u_{l−1})

und

∂_{r}[D(r, φ, z)∂_{z}u] =  (4.33)

\frac{1}{4h_{r}h_{z}} (D_{j+1}u_{j+1,l+1} − D_{j+1}u_{j+1,l−1} − D_{j−1}u_{j−1,l+1} + D_{j−1}u_{j−1,l−1}),

∂_{z}[D(r, φ, z)∂_{r}u] =  (4.34)

\frac{1}{4h_{r}h_{z}} (D_{l+1}u_{j+1,l+1} − D_{l+1}u_{j−1,l+1} − D_{l−1}u_{j+1,l−1} + D_{l−1}u_{j−1,l−1})

ergibt sich die Differenzengleichung zu

−\frac{1}{f_{l}^{2}r_{j}} [\frac{1}{2h_{r}^{2}} ((K_{j+1,l} + K_{j,l})u_{j+1} − (K_{j+1,l} + 2K_{j,l} + K_{j−1,l})u_{j}

+(K_{j,l} + K_{j−1,l})u_{j−1})

−\frac{f_{l}f′_{l}}{4h_{r}h_{z}} ((κ_{j+1,l}r_{j+1}^{2})u_{j+1,l+1} − (κ_{j+1,l}r_{j+1}^{2})u_{j+1,l−1}

−(κ_{j−1,l}r_{j−1}^{2})u_{j−1,l+1} + (κ_{j−1,l}r_{j−1}^{2})u_{j−1,l−1})
−\frac{r_{j}^{2}}{4h_{r}h_{z}} ((κ_{j,l+1}f_{l+1}f′_{l+1})u_{j+1,l+1} − (κ_{j,l+1}f_{l+1}f′_{l+1})u_{j−1,l+1}  (4.35)
−(κ_{j,l−1}f_{l−1}f′_{l−1})u_{j+1,l−1} + (κ_{j,l−1}f_{l−1}f′_{l−1})u_{j−1,l−1})

+\frac{1}{2r_{j}h_{φ}^{2}} ((κ_{k+1} + κ_{k})u_{k+1} − (κ_{k+1} + 2κ_{k} + κ_{k−1})u_{k}

+(κ_{k} + κ_{k−1})u_{k−1})

+\frac{r_{j}}{2h_{z}^{2}} ((κ_{l+1}f_{l+1}^{2} + κ_{l}f_{l}^{2})u_{l+1} − (κ_{l+1}f_{l+1}^{2} + 2κ_{l}f_{l}^{2} + κ_{l−1}f_{l−1}^{2})u_{l}

+(κ_{l}f_{l}^{2} + κ_{l−1}f_{l−1}^{2})u_{l−1})]
+(μ_{j,k,l} + i ω_{0})u_{j,k,l} = 0,

wobei zur Abkürzung K(r, φ, z) := κ(r, φ, z)r(1 + f′(z)^{2}r^{2}), also K_{j,k,l} = κ_{j,k,l}r_{j}(1 + (f′_{l})^{2}r_{j}^{2}) definiert sei.

Zur Behandlung der singulären Stelle bei r = 0 soll eine weitere Differenzengleichung hergeleitet werden. Dazu wird angenommen, daß die gesuchte Lösung bei r = 0 unabhängig vom Winkel ist, also daß ∂_{φ}u|_{r=0} = 0 gilt.

Da sich der Gradient ∇_{x} bezüglich kartesischer Koordinaten gemäß der Formel

∇_{x} = (∂_{r}, ∂_{φ}, ∂_{z})(∂_{r,φ,z}Φ)^{−1}

= \frac{1}{f(z)} ( \begin{matrix}cos φ \\ sin φ \\ −f′(z) r \end{matrix} ) \frac{∂}{∂r} + \frac{1}{f(z) r} ( \begin{matrix}− sin φ \\ cos φ \\ 0 \end{matrix} ) \frac{∂}{∂φ} + ( \begin{matrix}0 \\ 0 \\ 1 \end{matrix} ) \frac{∂}{∂z} (4.36)

in die kegelzylindrischen Koordinaten Φ transformiert, lautet die Gaußsche Integralidentität

∫_{Ω} ∇_{x} · (κ(x)∇_{x}u) dx = ∫_{∂Ω} (κ(x)∇_{x}u) · ν dσ_{x} = ∫_{∂Ω} (κ ∂_{ν}u) dσ_{x}  (4.37)

in den Koordinaten Φ für einen kleinen Bereich [0, h_{r}/2]×[0, 2π]×[z_{0} −h_{z}/2, z_{0} +h_{z}/2] um die für die Differentialgleichung singuläre Stelle (r, ·, z) = (0, ·, z_{0})

∫_{z_{0}−h_{z}/2}^{z_{0}+h_{z}/2} ∫_{0}^{2π} ∫_{0}^{h_{r}/2}(\frac{1}{\sqrt{g}} ∑_{p,q=1}^{3} ∂_{y_{p}}((κ ∘ Φ)(r, φ, z) g^{pq}\sqrt{g} ∂_{y_{q}}u)) \sqrt{g} dr dφ dz

= ∫_{z_{0}−h_{z}/2}^{z_{0}+h_{z}/2} ∫_{0}^{2π}((κ ∘ Φ)(h_{r}/2, φ, z)\frac{(1 + f′(z)^{2}h_{r}^{2}/4)∂_{r}u − (f(z)f′(z) h_{r}/2)∂_{z}u}{f(z)\sqrt{1 + f′(z)^{2}h_{r}^{2}/4}}) ×

× (f(z) h_{r}/2 \sqrt{1 + f′(z)^{2}h_{r}^{2}/4}) dφ dz

+ ∫_{0}^{2π} ∫_{0}^{h_{r}/2}(κ ∘ Φ)(r, φ, z_{0} + h_{z}/2)(∂_{z}u − \frac{f′(z_{0} + h_{z}/2) r}{f(z_{0} + h_{z}/2)} ∂_{r}u)f(z_{0} + h_{z}/2) r dr dφ

− ∫_{0}^{2π} ∫_{0}^{h_{r}/2}(κ ∘ Φ)(r, φ, z_{0} − h_{z}/2)(∂_{z}u − \frac{f′(z_{0} − h_{z}/2) r}{f(z_{0} − h_{z}/2)} ∂_{r}u)f(z_{0} − h_{z}/2) r dr dφ

= \frac{h_{r}}{2} ∫_{z_{0}−h_{z}/2}^{z_{0}+h_{z}/2} ∫_{0}^{2π}(κ ∘ Φ)(h_{r}/2, φ, z)((1 + f′(z)^{2}h_{r}^{2}/4) ∂_{r}u − (f(z)f′(z) h_{r}/2) ∂_{z}u) dφ dz

+ ∫_{0}^{2π} ∫_{0}^{h_{r}/2}(κ ∘ Φ)(r, φ, z_{0} + h_{z}/2)(f(z_{0} + h_{z}/2) r ∂_{z}u − f′(z_{0} + h_{z}/2) r^{2} ∂_{r}u) dr dφ

− ∫_{0}^{2π} ∫_{0}^{h_{r}/2}(κ ∘ Φ)(r, φ, z_{0} − h_{z}/2)(f(z_{0} − h_{z}/2) r ∂_{z}u − f′(z_{0} − h_{z}/2) r^{2} ∂_{r}u) dr dφ.

(4.38)

Indem diese Gleichung nun gleich dem Mittelwert des Absorptionsterms der Differentialgleichung über das angegebene Volumen gesetzt wird und die Oberflächenintegrale mittels der Trapezregel beziehungsweise der Mittelpunktregel angenähert werden, erhält man mit der rechten Seite (4.38) die gesuchte finite Differenzengleichung für den Punkt (0, ·, z_{0}).

4.2 Numerische Experimente

Obwohl die Umsetzung der Iterationsalgorithmen bezüglich eines dreidimensionalen Grundgebietes Ω ⊂ IR^{3} prinzipiell keine Schwierigkeiten bereitet, werden hier ausschließlich Ergebnisse bezüglich zweidimensionaler Rechtecksgebiete präsentiert. Der einzige Grund für diese Beschränkung ist die Einsparung an Laufzeit der Programme, da für einen Einzelschritt der Kaczmarz–Verfahren zur Rekonstruktion zwei partielle Differentialgleichungen numerisch gelöst werden müssen und typische Rekonstruktionsläufe etwa 20·p Einzelschritte erfordern, wobei p für die Anzahl der unterschiedlichen Messungen steht.

Die physikalisch relevanten konstanten Werte seien wie folgt gesetzt. Die Ausbreitungsgeschwindigkeit des Lichtes in Ω sei 0.75·c_{0}. Die Modulierungsfrequenz des zeitharmonischen Trägersignals betrage 2π · 110Mhz. Dies ergibt für den konstanten Imaginäranteil ω_{0} den Wert ω_{0} = ω/c = 3.072 · 10^{−2} cm^{−1}. Der konstante Grundwert der Absorption betrage stets μ_{0} = 2.5 · 10^{−2} cm^{−1} und der der Streuung κ_{0} = 2.772 · 10^{−2} cm.

Sämtliche numerische Rechnungen sind in einfacher Genauigkeit (single bzw. float) auf SUN–UltraSparc–Workstations durchgeführt worden.

4.2.1 Rekonstruktionsgebiet

Alle vorgestellten Ergebnisse beziehen sich auf einen festen geometrischen Aufbau. Das zugrunde gelegte Gebiet Ω ⊂ IR^{2} sei ein 6.0 × 4.0cm großes Rechteck. Dieses Gebiet sei mit 49 × 33 Punkten diskretisiert. Dies ergibt eine Maschenweite h = 0.125cm.

Als Randverteilung φ_{j}, 1 ≤ j ≤ p, sind p = 14 Punktquellen gewählt worden, konzentriert jeweils auf ein Pixel. Die Pixelkoordinaten s_{j} := (j, k) dieser Quellen lauten

(12, k), (16, k), (20, k), (24, k), (28, k), (32, k), (36, k),  k ∈ {0, 32},  (4.39)

und sind somit an den langen Seiten des Rechtecks verteilt.

4.2 NUMERISCHE EXPERIMENTE

4.2.2 Datenerzeugung

Zur Datenerzeugung wird hier das finite Differenzenschema (4.19) verwendet, welches auch in den Rekonstruktionsalgorithmen eingesetzt wird. Dieses Vorgehen birgt die Gefahr in sich, daß ein eventueller systematischer Fehler bei der Lösung des direkten Problems unentdeckt bliebe (engl. inverse crime). Demzufolge wurden die numerischen Lösungen aus (4.19) des direkten Problems mit den anderer Methoden, wie z.B. Finiter Elemente [61], oder anderer physikalischer Modelle, wie z.B. Transportgleichung [15], verglichen. Es konnte annähernde Übereinstimmung festgestellt werden.

Für die j-te Punktquelle, 1 ≤ j ≤ p, besteht der Datenvektor d_{j} = (d_{j,1}, …, d_{j,M})^{t} aus den Komponenten

d_{j,m} = ∂_{ν}(u_{j} − v_{j})(y_{m})/(∂_{ν}v_{j}(y_{m})),  m = 1, …, M,  (4.40)

wobei v_{j} die Lösung des direkten Problems bezüglich der konstanten Werte κ_{0} und μ_{0} und u_{j} die Lösung des direkten Problems bezüglich der variablen Koeffizienten κ und μ darstellt. Mit y_{m} sei die Verteilung der Detektorpositionen bezeichnet, an denen die Ableitung in Richtung der äußeren Einheitsnormalen genommen wird. Hier seien alle M = 156 echten Randpunkte des diskretisierten Rechtecks genommen, d.h. die vier Eckpunkte werden nicht berücksichtigt.

Die approximative Bestimmung des Wertes der Richtungsableitung geschieht wiederum anhand finiter Differenzen,

∂_{ν}u(y_{m}) ≃ \frac{1}{h}(φ_{m} − u_{m+1}),  (4.41)

falls das Gebiet Ω im Punkt y_{m} zur Rechten liegt.

4.2.3 Rekonstruktionen

Für die noch freien Parameter der Rekonstruktionsverfahren N und α sind in den hier vorgestellten Beispielen aus experimenteller Erfahrung N = 15 und α = 0.01 gewählt worden. Der Relaxationsparameter ϱ mußte dem Konvergenzverhalten der jeweiligen Testläufe angepaßt werden, und wird in den entsprechenden Bildunterschriften angegeben.

Als Abbruchkriterium der äußeren Schleife wurde jedoch nicht die Fehlerschranke pε benutzt, sondern die Schleife stets nach 20 Durchgängen abgebrochen.

Die aus den Randdaten zu bestimmenden Objekte sind zwei kreisrunde Teilgebiete in Ω mit höherer Absorption bzw. höherer Streuung. Dabei habe die höhere Absorption einen Kontrast von 500% in bezug auf den Hintergrundwert und die Streuung den Kontrast 300%. Man beachte, daß die Abbildungen auf der linken Seite jeweils die Absorption μ_{a} darstellen und auf der rechten Seite den Diffusionskoeffizient κ zeigen.

Der Durchmesser der zu rekonstruierenden Teilgebiete beträgt in diesen Beispielen 10mm bzw. 4mm.

Fall  Hintergrund  Objekte c · χ_{B(x,r)}
μ_{a} [cm^{−1}]  μ′_{s} [cm^{−1}]  μ_{a} [c = 0.1cm^{−1}] μ′_{s} [c = 24.0cm^{−1}]
a)  0.025  12.0  r = 0.5  r = 0.5
x = (2.0, 2.0)  x = (4.0, 1.5)
b) 0.025 + 0.3 sin 0.3t 12.0 + 0.1 sin 0.2t  r = 0.5  r = 0.5
x = (2.0, 2.0)  x = (4.0, 1.5)
c)  0.025  12.0  r = 0.2  r = 0.2
x = (2.0, 1.8)  x = (4.0, 2.0)
d) 0.025 + 0.3 sin 0.3t 12.0 + 0.1 sin 0.2t  r = 0.2  r = 0.2
x = (2.0, 1.8)  x = (4.0, 2.0)

Für die Testläufe werden die Objekte additiv zu zwei verschiedene Hintergründe hinzugefügt: einen Hintergrund mit konstanten Werten und einen in Richtung der Koordinatenachsen sinusartig variierenden Hintergrund (Schachbrettmuster).

Insgesamt ergeben sich so vier Fallbeispiele, die in oben stehender Tabelle zusammengefaßt sind, wobei Radius und Mittelpunkt der Kreise, die die Lage der Objekte beschreiben, in cm bezüglich des 6×4cm großen Rekonstruktionsgebietes angegeben sind.

Ergebnisse des Rekonstruktionsalgorithmus 4.1 für die Fälle a) bis d) zeigen die Bilder 1 bis 4. Der Graph in der letzten Zeile der Rekonstruktionen bildet die normalisierten Residuenwerte

r = \frac{∥r_{j}∥_{ℓ_{2}}^{2}}{∥d_{j}|_{ℓ_{2}}^{2}}  mit  r_{j} = d_{j} − (\frac{1}{∂_{ν}v_{j}} ∂_{ν}(u_{j}^{akt} − v_{j})) |_{∂Ω}  (4.42)

des ersten Durchlaufs (grün) und des zwanzigsten Durchlaufs (blau) ab. Die horizontale Achse entspricht dabei der Reihenfolge der Quellpositionen in den Einzelschritten.

Die Bilder 5 und 6 sollen den Unterschied zwischen den Rekonstruktionsergebnissen der rein reellwertigen Rechnung gemäß der Formel (3.182), der naiven Rechnung, bei der mit N = 0 in Algorithmus 4.1 die innerste Schleife übersprungen wird, und der modifizierten Rechnung mit N = 15 in den Fällen a) und c) verdeutlichen.

Darüberhinaus sind mit den Werten der Fälle a) und c) komplexwertige Potentiale

q(x) = \frac{μ_{a}(x) + \sqrt{κ(x)} Δ_{x}\sqrt{κ(x)}}{κ(x)} + i \frac{ω_{0}}{κ(x)},  κ(x) = \frac{1}{3(μ′_{s}(x) + μ_{a}(x))}  (4.43)

gebildet worden. Die damit erzeugten Randdaten zeigen in der Tat sehr gute Übereinstimmung mit den entsprechenden Datensätze der Fälle a) und c).

Rekonstruktionen mit dem komplexwertigen Kaczmarz–Newton–Algorithmus 4.2 zeigen die Bilder 7 und 8.

Schließlich wurden zum Vergleich Rekonstruktionen für Parameterwerte, die aus dem Artikel [54] entnommen sind, berechnet. Im einzelnen werden dort folgende Werte verwendet: ω = 2π ·300MHz, μ_{a} = 0.1cm^{−1} und μ′_{s} = 10.0cm^{−1} als Hintergrundwerte, und 100% Konstrast bei Objekten in μ_{a} und in μ′_{s}. Der Durchmesser der Objekte in [54] beträgt 10mm bzw. 20mm. Das gesamte Rekonstruktionsgebiet hat dabei die Ausmaße 10 × 10cm und sei mit 33 × 33 Punkten diskretisiert.

Zwar wird in der zitierten Arbeit das direkte Problem anhand eines vollen 3D–Modells erzeugt, die Rekonstruktionen jedoch auf eine Schicht der in z-Richtung konstanten Objekte beschränkt. Insofern ist ein Vergleich mit den hier ausschließlich in 2D durchgeführten Rechnungen gerechtfertigt. Desweiteren sind die geringere Anzahl an Quellpositionen (4) und Detektoren (20) zu beachten. Diese seien gleichmäßig über die vier Ränder des Rekonstruktionsgebietes verteilt.

Die Bilder 9 bis 12 zeigen verschiedenen Rekonstruktionen für Objekte mit diesen Referenzwerten.

Bild 1

Rekonstruktion für Fall a) mit ϱ = 0.4, α = 0.01 und N = 15 nach 20 Durchläufen; links μ_{a}(x), rechts κ(x), von oben nach unten: Farbkeil, Originale, Querschnitte und Verlauf der Residuenwerte im ersten (grün) und im zwanzigsten (blau) Durchlauf.

4.2 NUMERISCHE EXPERIMENTE Bild 2 Rekonstruktion für Fall b) mit ϱ = 0.3, α = 0.01 und N = 15 nach 20 Durchläufen; gestartet mit den konstanten Hintergrundwerten κ_{0}, μ_{0} aus Fall a).

Bild 3

Rekonstruktion für Fall c) mit ϱ = 0.9, α = 0.01 und N = 15 nach 20 Durchläufen; links μ_{a}(x), rechts κ(x), von oben nach unten: Farbkeil, Originale, Querschnitte und Verlauf der Residuenwerte im ersten (grün) und im zwanzigsten (blau) Durchlauf.

4.2 NUMERISCHE EXPERIMENTE Bild 4 Rekonstruktion für Fall d) mit ϱ = 0.9, α = 0.01 und N = 15 nach 20 Durchläufen; gestartet mit den konstanten Hintergrundwerten κ_{0}, μ_{0} aus Fall c).

Bild 5

Vergleich der Rekonstruktionen für Fall a) nach 20 Durchläufen, oben: ϱ = 0.4, α = 0.01 und N = 15, mitte: ϱ = 0.4, α = 0.01 und N = 0, unten: ϱ = 0.0004 und α = 0.02; links μ_{a}(x), rechts κ(x).

Bild 6

Vergleich der Rekonstruktionen für Fall c) nach 20 Durchläufen, oben: ϱ = 0.9, α = 0.01 und N = 15, mitte: ϱ = 0.9, α = 0.01 und N = 0, unten: ϱ = 0.0009 und α = 0.02; links μ_{a}(x), rechts κ(x).

Bild 7

Rekonstruktion für Fall a) mit komplexwertigen Potential q(x) und ϱ = 0.12, α = 0.01 nach 20 Durchläufen; links Re(q(x)), rechts Im(q(x)), von oben nach unten: Farbkeile, Originale, Querschnitte und Verlauf der Residuenwerte im ersten (grün) und im zwanzigsten (blau) Durchlauf.

Bild 8

Rekonstruktion für Fall c) mit komplexwertigen Potential q(x) und ϱ = 0.08, α = 0.01 nach 20 Durchläufen; links Re(q(x)), rechts Im(q(x)), von oben nach unten: Farbkeile, Originale, Querschnitte und Verlauf der Residuenwerte im ersten (grün) und im zwanzigsten (blau) Durchlauf.

Bild 9 Rekonstruktionen mit vier Quellpunkten und gesamtem Rand als Detektor nach 20 Durchläufen, linke Seite μ_{a}(x), rechte Seite κ(x), oben: jeweils Original, reellwertige und korrigierte Rechnung, Farbkeil, mitte: Querschnitte: Original (grün), reellwertige (blau) und korrigierte (rot) Rechnung, unten: Residuenwerte im ersten (grün, blau) bzw. zwanzigsten (rot, pink) Durchlauf.

Bild 10 Rekonstruktionen mit vier Quellpunkten und gesamtem Rand als Detektor nach 20 Durchläufen, linke Seite μ_{a}(x), rechte Seite κ(x), oben: jeweils Original, reellwertige und korrigierte Rechnung, Farbkeil, mitte: Querschnitte: Original (grün), reellwertige (blau) und korrigierte (rot) Rechnung, unten: Residuenwerte im ersten (grün, blau) bzw. zwanzigsten (rot, pink) Durchlauf.

Bild 11 Rekonstruktionen mit vier Quellpunkten und fünf Detektoren je Seite nach 20 Durchläufen, linke Seite μ_{a}(x), rechte Seite κ(x), oben: jeweils Original, reellwertige und korrigierte Rechnung, Farbkeil, mitte: Querschnitte: Original (grün), reellwertige (blau) und korrigierte (rot) Rechnung, unten: Residuenwerte im ersten (grün, blau) bzw. zwanzigsten (rot, pink) Durchlauf.

Bild 12 Rekonstruktionen mit vier Quellpunkten und fünf Detektoren je Seite nach 20 Durchläufen, linke Seite μ_{a}(x), rechte Seite κ(x), oben: jeweils Original, reellwertige und korrigierte Rechnung, Farbkeil, mitte: Querschnitte: Original (grün), reellwertige (blau) und korrigierte (rot) Rechnung, unten: Residuenwerte im ersten (grün, blau) bzw. zwanzigsten (rot, pink) Durchlauf.

4.3 Diskussion

Alle gezeigten Rekonstruktionen sind mit exakten, fehlerlosen Randdaten berechnet worden. Daher kann nichts über die Stabilität der verwendeten Verfahren gesagt werden. Die Rekonstruktionsergebnisse bezüglich des variablen Hintergrundes könnten jedoch als ersten Hinweis derartiger Aussagen gedeutet werden. Da das behandelte inverse Problem äußerst schlecht gestellt ist, siehe Abschnitt 3.2.2, ist eine der Zielrichtungen dieser Arbeit gewesen, zu untersuchen, ob die hier angegebenen Algorithmen zur optischen Tomographie verwendet werden können oder nicht.

Die Fälle a) bis d) sollen nun im einzelnen diskutiert werden.

Fall a)

Wie schon durch V. P. Palamodov [51] bemerkt, zeigt sich bei dieser Rekonstruktion sehr deutlich, daß man unter der gegebenen Parameterwahl nur erwarten kann, höchstens zwei bis drei Fourier–Moden der Objekte aus den gegebenen Randdaten wiederzufinden. Um so erstaunlicher ist die sehr gute übereinstimmung des rekonstruierten Diffusionskoeffizienten κ(x) zum Original. Es ist anzumerken, daß bei einer forgesetzten Iteration über 20 Durchläufe hinaus das Objekt im Diffusionskoeffizient κ(x) mehr und mehr überbewertet wird, so daß schließlich die Elliptizität des zum direkten Problem gehörigen Differentialoperators an dieser Stelle verloren geht. Damit bricht die Iteration zusammen. Der Verlauf der Residuenwerte zeigt jedoch an, daß nach 20 Durchläufen die Daten praktisch erfüllt sind, und also keine weitere Information zur Verarbeitung vorhanden ist.

Fall b)

Es zeigt sich, daß die Rekonstruktion relativ unempfindlich gegenüber der hier gewählten Störung im Hintergrund ist. Dabei sei noch einmal betont, daß der Algorithmus 4.1 mit den konstanten Hintergrundwerten μ′_{s} und μ_{a} aus Fall a) gestartet ist.

Fall c)

Im Vergleich zum Resultat des Falles a) kann man vermuten, daß in diesem Fall die Auflösungsgenauigkeit sowohl in κ(x) als auch μ_{a}(x) erreicht sind. Darüberhinaus zeigen die Residuenwerte an, daß zu erwarten ist, daß die Iteration nach einigen weiteren Durchläufen wegen der großen Rekonstruktionsfehler am Rand des Gebietes divergieren wird. Man kann feststellen, daß dies in der Tat der Fall ist.

Fall d)

In diesem Fall trifft in etwa ähnliches zu wie in Fall c). Der variable Hintergrund scheint jedoch eine glättende Einfluß auf die Rekonstruktion zu haben.

4.3 DISKUSSION

Der Vergleich zwischen der naiven (N = 0) und der korrigierten Fassung (etwa N = 15) des Kaczmarz–Newton–Verfahrens machen eindrucksvoll die Notwendigkeit dieser Korrektur klar deutlich. Auch die Rechnungen bezüglich der Referenzwerte aus [54] legen diesen Schluß nahe.

Die komplexwertigen Rekonstruktionen mit dem Algorithmus 4.2 weisen das in Abschnitt 3.1.9 vermutete Verhalten auf. Allerdings kann bei diesen Rekonstruktionen das sonst zu beobachtene Übersprechverhalten in den gesuchten Koeffizienten überhaupt nicht festgestellt werden. Insofern könnte im Hinblick auf Anwendungen, bei denen es hauptsächlich um die Trennung der beiden optischen Parameter ankommt, dieses Rekonstruktionsverfahren sehr interessant sein.

KAPITEL 4. NUMERISCHER TEIL 84

Zusammenfassung

Diese Arbeit behandelt iterative Verfahren zur Lösung des inversen Problems der optischen Tomographie im Frequenzbereich.

Die Fréchet–Differenzierbarkeit der Dirichlet-zu-Neumann–Abbildung, die die gesuchten reellwertigen optischen Parameter den zu messenden komplexwertigen Daten zuordnet, wird bewiesen.

Aufgrund dieser Eigenschaft wird das Kaczmarz–Newton–Verfahren, eine nichtlineare Erweiterung des linearen ART–Verfahrens, hergeleitet. Es zeigt sich, daß die Besonderheit, reellwertige Lösungen aus komplexwertigen Daten zur Lösung des inversen Problems der optischen Tomographie zu bestimmen, genau zu berücksichtigen ist, um mit diesem Verfahren brauchbare Rekonstruktionen berechnen zu können. Darüberhinaus wird durch die Verwendung der Sherman–Morrison–Woodbury–Formel in dem Kaczmarz–Newton–Verfahren eine Regularisierungsmethode für das schlecht gestellte inverse Problem der optischen Tomographie vorgeschlagen.

Es wird aufgezeigt, daß diese Regularisierungsmethode für das nichtlineare Kaczmarz– Newton–Verfahren bei geringerem numerischen Aufwand vergleichbare oder sogar bessere Resultate liefert als bisher mit gängigen Rekonstruktionsverfahren zur optischen Tomographie erzielten Ergebnisse.

Literaturverzeichnis

KAPITEL 5. ZUSAMMENFASSUNG 86 [1] Alt, H.W.: Lineare Funktionalanalysis, Springer-Verlag, Berlin, 1992.

[2] Arridge, S.R.: Optical tomography in medical imaging, Inverse Problems 15 (1999), R41–R93.

[3] Arridge, S.R.; Schweiger, M.: A gradient-based optimisation scheme for optical tomography, Optics Express 2 (1998), 213–226.

[4] Arridge, S.R.; Hebden, J.C.: Optical imaging in medicine: II. Modelling and reconstruction, Phys. Med. Biol. 42 (1997), 841–853.

[5] Arridge, S.R.: Photon-measurement density functions. Part I: Analytical forms, Appl. Opt. 34 (1995), 7395–7409.

[6] Arridge, S.R.; Schweiger, M.: Photon-measurement density functions. Part II: Finite-Element-Method calculations, Appl. Opt. 34 (1995), 8026–8037.

[7] Arridge, S.R.; Hiraoka, M.; Schweiger, M.; Delpy, D.T.: Performance of an iterative reconstruction algorithm for near-infrared absorption and scatter imaging, in: Chance, B.; Alfano, R.R. (eds.): Photon Migration and Imaging in Random Media and Tissue, Proc. SPIE 1888 (1993), 360–371.

[8] Barber, D.C.; Seagar, A.D.: Fast reconstruction of resistive images, Clin. Phys. P. 8 (1987), 47–57.

[9] Barrett, R.; Berry, M.; Chan, T.F.; Demmel, J.; Donato, J.; Dongarra, J.; Eijkhout, V.; Pozo, R.; Romine, Ch.; van der Vorst, H.: Templates for the solution of linear systems: Building Blocks for iterative methods, SIAM, Philadelphia, 1994.

[10] Calderón, A.P.: On an inverse boundary value problem, Seminar on Numerical Analysis and its Applications to Continuum Physics, Soc. Brasileria de Mathemàtica (1980), Rio de Janeiro, 65–73.

[11] Campbell, S.L.; Meyer, C.D.: Generalized Inverses of Linear Transformations, Pitman, London, 1979.

[12] Chance, B.; Boas, D.A.; Yodh, A.G.: Experimental images of heterogeneous turbid media by frequency-domain diffusing-photon tomography, Opt. Lett. 20 (1995), 426–428.

[13] Cheney, M.; Isaacson, D.; Newell, J.C.: Electrical impedance tomography, SIAM Review 41 (1999), 85–101.

[14] Courant, R.; Hilbert, D.: Methoden der mathematischen Physik, Springer-Verlag, Berlin, 1993.

[15] Dorn, O.: Das inverse Transportproblem in der Lasertomographie, Dissertation, Universität Münster, 1997.

[16] Feng, S.; Zeng, Z.A.; Chance, B.: Photon migration in the presence of a single defect: a perturbation analysis, Appl. Opt. 34 (1995), 3826–3837.

[17] Folland, G.B.: Introduction to Partial Differential Equations, 2nd ed., Princeton Academic Press, Princeton, 1995.

[18] Gilbarg, D.; Trudinger, N.S.: Elliptic Partial Differential Equations of Second Order, Springer-Verlag, Berlin, 1977.

[19] Golub, G.H.; van Loan, Ch.F.: Matrix Computations, The Johns Hopkins University Press, Baltimore, 1996.

[20] Grinberg, N.: Inverse Boundary Problem for 2D Diffusion Equation near a constant Background: Uniqueness, Technical Report 10/98 N, Fachbereich Mathematik der Universität Münster, 1998.

[21] Grisvard, P.: Elliptic Problems in Nonsmooth Domains, Pitman Publishing, London, 1985.

[22] Gruenbaum, F.A.; Patch, S.K.: How many parameters can one solve for in diffuse tomography?, in: Chavent, G. (ed.): Inverse Problems in Wave Propagation, IMA Vol. Math. Appl. 90 (1997), 219–235.

[23] Hadamard, J.: Le problème de Cauchy et les equations aux derivées partielles linéaires hyperboliques, Herman, Paris, 1932.

[24] Hebden, J.C.; Arridge, S.R.; Delpy, D.T.: Optical imaging in medicine: I. Experimental techniques, Phys. Med. Biol. 42 (1997), 825–840.

[25] Herman, G.T.: Image Reconstruction from Projections, Academic Press, New York, 1980.

[26] Jöbsis, F.F.: Noninvasive infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters, Science 198 (1977), 1264–1267.

[27] Kak, A.C.; Slaney, M.: Principles of Computerized Tomographic Imaging, IEEE Press, New York, 1987.

[28] Kaltenbach, J.P.; Kaschke, M.: Frequency and Time-Domain Modelling of Light Transport in Random Media, in: Müller, G. (ed.): Medical Optical Tomography: Functional Imaging and Monitoring, SPIE Optical Engineering Press, Vol. IS11 (1993), 65–86.

[29] Kelley, C.T.: Iterative Methods for Linear and Nonlinear Equations, SIAM, Philadelphia, 1995.

[30] Kirsch, A.: An Introduction to the Mathematical Theory of Inverse Problems, Springer-Verlag, New York, 1996.

[31] Klibanov, M.V.; Lucas, T.R.; Frank, R.M.: A fast and accurate imaging algorithm in optical/diffusion tomography, Inverse Problems 13 (1997), 1341–1361.

[32] Kohn, R.; Vogelius, M.: Determining conductivity by boundary measurements, Commun. Pure Appl. Math. 37 (1984), 289–298.

[33] Kronecker, L: Theorie der einfachen und der vielfachen Integrale, ed. Netto E., Teubner, Leipzig, 1894.

[34] O’Leary, M.A.; Boas, D.A.; Chance, B.; Yodh, A.G.: Simultaneous scattering and absorption images of heterogeneous media using diffusive waves within the Rytov approximation, in: Chance, B.; Alfano, R.R. (eds.): Optical Tomography, Photon Migration, and Spectroscopy of Tissue and Model Media: Theory, Human Studies, and Instrumentation, Proc. SPIE 2389 (1995), 320–327.

[35] Ortega, J.M.; Rheinboldt, W.C.: Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970.

[36] Lin, F.C.; Fiddy, M.A.: On the issue of the Born-Rytov controversy: I Comparing analytical and approximate expressions for the one-dimensional case, J. Opt. Soc. Amer. A 9 (1992), 1102–1110.

[37] Lin, F.C.; Fiddy, M.A.: The Born-Rytov controversy: II Applications to the nonlinear and stochastic scattering problem in a one-dimensional half-space, J. Opt. Soc. Amer. A 10 (1993), 1971–1983.

[38] Lions, J.L.; Magenes, E.: Non-homogeneous Boundary Value Problems and Applications, Vol. I, Springer-Verlag, New York, 1972.

[39] Louis, A.K.: Inverse und schlecht gestellte Probleme, Teubner-Verlag, Stuttgart, 1989.

[40] Maritz, M.F.; Herman, G.T.; Yee, C.: Recovery of the absorption coefficient from diffused reflected light using a discrete diffusive model, SIAM J. Appl. Math. 59 (1999), 58–71.

[41] McCormick, S.F.: An iterative procedure for the solution of constrained nonlinear equations with application to optimization problems, Numer. Math. 23 (1975), 371–385.

[42] Meyn, K.-H.: Solution of Underdetermined Nonlinear Equations by Stationary Iteration Methods, Numer. Math. 42 (1983), 161–172.

[43] Mitic, G.; Kölzer, J.; Otto, J.; Plies, E.; Sölkner, G.; Zinth, W.: Time-gated transillumination of biological tissues and tissuelike phantoms, Appl. Opt. 33 (1994), 6699–6710.

[44] Model, R.; Ortl, M.; Walzel, M.; Hünlich, R.: Reconstruction algorithm for nearinfrared imaging in turbid media by means of time-domain data, Appl. Opt. 14 (1997), 313–324.

[45] Model, R.; Hünlich, R.: Optical imaging of highly scattering media, Z. Angew. Math. Mech. 76 (1996), 483–484.

[46] Nachman, A.I.: Reconstructions from boundary measurements, Ann. Math. 128 (1988), 531–576.

[47] Natterer, F.: The mathematics of computerized tomography, Wiley, New York, 1986.

[48] Natterer, F.; Wübbeling, F.: A propagation-backpropagation method for ultrasound tomography, Inverse Problems 11 (1995), 1225–1232.

[49] Natterer, F.: Numerical Solutions of Bilinear Inverse Problems, Technical Report 19/96 N, Fachbereich Mathematik der Universität Münster, 1996.

[50] Paasschens, J.: On the Transmission of Light through Random Media, Dissertation, Rijksuniversiteit te Leiden, 1997.

[51] Palamodov, V.P.: persönliche Mitteilung, Oberwolfach, 1998.

[52] Patterson, M.S.; Pogue, B.W.; Farrell, T.J.: Forward and inverse calculations for 3D frequency-domain diffuse optical tomography, in: Chance, B.; Alfano, R.R. (eds.): Optical Tomography, Photon Migration, and Spectroscopy of Tissue and Model Media: Theory, Human Studies, and Instrumentation, Proc. SPIE 2389 (1995), 320–327.

[53] Peters, V.G.; Wyman, D.R.; Patterson, M.S.; Frank, G.L.: Optical properties of normal and diseased human breast tissue in the visible and near infrared, Phys. Med. Biol. 35 (1990), 1317–1334.

[54] Pogue, B.W.; Patterson, M.S.; Jiang, H.; Paulsen, K.D.: Initial assessment of a simple system for frequency domain diffuse optical tomography, Phys. Med. Biol. 40 (1995), 1709–1729.

[55] Renardy, M.; Rogers, R.C.: An Introduction to Partial Differential Equations, Springer-Verlag, New York, 1992.

[56] Romanov, V.G.; He, S.: Some uniqueness theorems for mammography-related time-domain inverse problems for the diffusion equation, Inverse Problems 16 (2000), 447–459.

[57] Schomberg, H.: Nonlinear image reconstruction from projections of ultrasonic travel times and electric current densities, in: Herman, G.T.; Natterer, F. (eds.): Mathematical aspects of computerized tomography, Lect. Notes in Med. and Biol. 8, Springer, New York, 1981.

[58] Schweiger, M.; Arridge, S.R.: Comparison of two- and three-dimensional reconstruction methods in optical tomography, Appl. Opt. 37 (1998), 7419–7428.

[59] Schweiger, M.; Arridge, S.R.; Hiraoka, M.; Delpy, D.T.: The finite element method for the propagation of light in scattering media: Boundary and source conditions, Med. Phys. 22 (1995), 1–14.

[60] Schweiger, M.; Arridge, S.R.; Delpy, D.T.: Application of the finite-element method for the forward and inverse models in optical tomography, J. Math. Imag. Vision 3 (1993), 263–283.

[61] Schweiger, M.; Arridge, S.R.: Inverse Methods for Optical Tomography, in: Barrett, H.H., Gmitro, A.F. (eds.): Information Processing in Medicine and Imaging, 13th Int. Conf., Proc. IPMI ’93, Springer, Berlin, 1993, 259–277.

[62] Stakgold, I.: Boundary Value Problems of Mathematical Physics, Vol. I, MacMillan Company, New York, 1967.

[63] Stakgold, I.: Boundary Value Problems of Mathematical Physics, Vol. II, MacMillan Company, New York, 1968.

[64] Sylvester, J.; Uhlmann, G.: A global uniqueness theorem for an inverse boundary value problem, Ann. Math. 125 (1987), 153–169.

[65] Thomas, J.W.: Numerical Partial Differential Equations: Finite Difference Methods, Springer, New York, 1995.

[66] Tompkins, C.: Projection methods in calculation, Proc. Sec. Symp. Lin. Progr., Washington D.C., 1955, 425–448.

[67] Triebel, H.: Höhere Analysis, VEB Deutscher Verlag der Wissenschaften, Berlin, 1972.

[68] Wilkinson, J.H.: The Algebraic Eigenvalue Problem, Clarendon Press, Oxford, 1965.

[69] Yao, Y.; Wang, Y.; Pei, Y.; Zhu, W.; Barbour, R.L.: Frequency-domain optical imaging of absorption and scattering distributions by a Born iterative method, J. Opt. Soc. Am. A 14 (1997), 325–342.

[70] Yosida, K.: Functional Analysis, Springer-Verlag, Berlin, 1980.

Lebenslauf

Persönliche Daten

Name:  Thomas Dierkes
geboren:  6. Januar 1969 in Münster/i.Westf.
Familienstand:  ledig
Name der Eltern:  Ernst Günter Dierkes und

Leonie Katharina Dierkes, geb. Dammann

Schulbildung

1975 – 1979  Nikolaischule Wolbeck
1979 – 1988  Gymnasium Wolbeck
10. Juni 1988  Abitur am Gymnasium Wolbeck in Münster

Studium

Oktober 1988  Beginn des Studiums der Mathematik

(Diplom, Nebenfach Informatik) an der

Westfälischen Wilhelms-Universität (WWU) Münster

April 1991  Vordiplom in Mathematik
Oktober 1993  Beginn der Diplomarbeit am Institut für

Numerische und instrumentelle Mathematik der WWU Münster unter der Betreuung von Prof. Dr. F. Natterer

10. Mai 1995  Diplom im Studienfach Mathematik
August 1995  Beginn der Dissertation am Institut für

Numerische und instrumentelle Mathematik der WWU Münster Betreuer: Prof. Dr. F. Natterer

Tätigkeiten

Jan. 1989 – Dez. 1989 Studentische Hilfskraft am Institut für

Soziologie/Sozialpädagogik (IfSS) der WWU Münster

Okt. 1991 – Mär. 1992

Mai 1992 – Okt. 1994 Studentische Hilfskraft am Universitätsrechenzentrum

der WWU Münster

seit Aug. 1995  Wissenschaftlicher Mitarbeiter am Institut für

Numerische und instrumentelle Mathematik

der WWU Münster