Tomographie mit wenigen Daten

Bildrekonstruktion aus Daten der Schallpyrometrie

mittels Kollokationsverfahren

Diplomarbeit

vorgelegt von

Helmut Sielschott

Thema gestellt von

Prof. Dr. F. Natterer 8. August 1994

�

Westfälische Wilhelms-Universität Münster

Institut für Numerische und instrumentelle Mathematik

Diese Version der Diplomarbeit entstand am 15. September 1994. Sie unterscheidet sich inhaltlich nicht von der Originalversion vom 8. August 1994. Es wurden lediglich Rechtschreibfehler korrigiert.

Inhaltsverzeichnis

Abbildungsverzeichnis 2

1 Einleitung 3

2 Schallpyrometrie 6 2.1 Beschreibung der Meßmethode . . . . . . . . . . . . . . . . . . . . . 6 2.2 Aufbereitung der Daten . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 A priori Informationen . . . . . . . . . . . . . . . . . . . . . . . . . . 10

3 Rekonstruktionsverfahren 11 3.1 Problembeschreibung . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 Kollokationsverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3 Mehrfachmessungen . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3.3.1 Kanonische Aufstellung des Gleichungssystems . . . . . . . . 15 3.3.2 Reduktion des Gleichungssystems . . . . . . . . . . . . . . . . 15

3.4 Singulärwertzerlegung der diskreten Radontransformation . . . . . . 17 3.5 Wahl geeigneter Ansatzfunktionen . . . . . . . . . . . . . . . . . . . 21

3.5.1 Der Kovarianzenansatz . . . . . . . . . . . . . . . . . . . . . . 21 3.5.2 Andere Ansatzfunktionen . . . . . . . . . . . . . . . . . . . . 23 3.5.3 Singulärfunktionen . . . . . . . . . . . . . . . . . . . . . . . . 26

3.6 Einbringen von a priori Informationen . . . . . . . . . . . . . . . . . 28 3.6.1 Startbilder . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.6.2 Glättung im oberen Temperaturbereich . . . . . . . . . . . . 30

4 Numerische Ergebnisse 32

4.1 Verfahrensvergleiche anhand von Phantombildern . . . . . . . . . . . 32 4.1.1 Die verwendeten Phantome . . . . . . . . . . . . . . . . . . . 32 4.1.2 Kollokation mit φ_{k} = e^{−λ|x−x_{k}|}: Auswirkung der Regularisierung 35 4.1.3 Vergleich der Ergebnisse mit anderen Verfahren . . . . . . . . 40 4.1.4 Einbringen von a priori Informationen . . . . . . . . . . . . . 44

4.2 Rekonstruktionen aus fehlerbehafteten Daten . . . . . . . . . . . . . 48 4.2.1 Die 24-Liniengeometrie aus Neurath . . . . . . . . . . . . . . 49 4.2.2 Rekonstruktionen mit 4 Meßgeräten pro Wand . . . . . . . . 52

4.3 Anwendung auf Meßdaten . . . . . . . . . . . . . . . . . . . . . . . . 56

ABBILDUNGSVERZEICHNIS

4.4 Diskussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5 Zusammenfassung 60

Literaturverzeichnis 62

Abbildungsverzeichnis

2.1 Approximation von Meßdaten durch einen Spline . . . . . . . . . . . 9

3.1 Singulärfunktionen bei einem Kollokationsverfahren . . . . . . . . . 26 3.2 Erste Singulärfunktion bei verschiedenen Werten von λ . . . . . . . . 27 3.3 Abhängigkeit der ersten Singulärfunktion von den Ansatzfunktionen 28 3.4 Vergleich der η_{k} mit den Ansatzfunktionen beim Kovar.-Verfahren . 29

4.1 Verwendete Phantome . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2 Phantome mit 24-Linien Abtastgeometrie . . . . . . . . . . . . . . . 34 4.3 Die zu den Phantomen gehörenden 24 Daten . . . . . . . . . . . . . 35 4.4 Darstellung der verwendeten Filter . . . . . . . . . . . . . . . . . . . 36 4.5 Rekonstruktionen mit 24 Linien und γ = 0.3 . . . . . . . . . . . . . . 37 4.6 Rekonstruktionen mit 24 Linien und γ = 0.05 . . . . . . . . . . . . . 38 4.7 Rekonstruktionen mit 24 Linien ohne Regularisierung . . . . . . . . 39 4.8 Rekonstruktionen mit dem Kovarianzenverfahren . . . . . . . . . . . 40 4.9 Rekonstruktionen mit dem Fourierverfahren . . . . . . . . . . . . . . 42 4.10 Eine hochfrequente Ansatzfunktion beim Fourierverfahren . . . . . . 42 4.11 Rekonstruktionen mit Pixelansatzfunktionen (ähnlich ART) . . . . . 43 4.12 Auswirkung des Blockprofils als Startbild . . . . . . . . . . . . . . . 44 4.13 Rekonstruktionen mit Startbild und Glättung . . . . . . . . . . . . . 45 4.14 Wie in Abb. 4.13, jedoch stärker regularisiert . . . . . . . . . . . . . 46 4.15 Auswirkungen verschiedener Startbilder . . . . . . . . . . . . . . . . 47 4.16 Rekonstruktion von A1 aus 24 fehlerbehafteten Daten . . . . . . . . 50 4.17 Rekonstruktion von AKomplex aus 24 fehlerbehafteten Daten …. 51 4.18 Die zu den Phantomen gehörenden 96 Daten . . . . . . . . . . . . . 52 4.19 Die verwendeten Filter bei 96 Daten . . . . . . . . . . . . . . . . . . 52 4.20 Rekonstruktionen mit 96 Daten und γ = 0.3 . . . . . . . . . . . . . . 53 4.21 Rekonstruktionen mit 96 Daten und γ = 0.005 . . . . . . . . . . . . 54 4.22 Rekonstruktionen mit 96 Daten und γ = 0.005. Datenfehler: 1% … 55 4.23 Rekonstruktionen aus Meßdaten . . . . . . . . . . . . . . . . . . . . 57

Kapitel 1

Einleitung

Tomographie

Unter dem Begriff ”Tomographie“ versteht man in der Mathematik die Lösung des sogenannten ”Radonschen Inversionsproblems“, d.h. die Rekonstruktion von Funktionen aus ihren Linien- oder Flächenintegralen.

Schon im Jahre 1917 hat J. Radon dieses Problem betrachtet und in der Theorie gelöst: Er konnte eine konstruktive Inversionsformel angeben.

Für die Praxis wurde es erst ab 1963 interessant, als neue medizinische Verfahren, die heute als ”Computertomographie“ allgemein bekannt sind, tatsächlich die Linienintegrale einer Funktion als Meßwerte lieferten und eine Rekonstruktion der Funktion forderten.

Jetzt tauchten in der praktischen Anwendung aber Probleme auf. Das Rekonstruktionsproblem ist (im Sinne Hadamards) schlecht gestellt, das heißt, daß sich die unvermeidlichen Meß-, Diskretisierungs- und Rundungsfehler so stark auswirken, daß die Rekonstruktionen unbrauchbar werden können. Daher mussten stabilere Verfahren entwickelt werden.

Der Standardalgorithmus in der Computertomographie ist seit 20 Jahren die Gefilterte Rückprojektion. Im Gegensatz zu den in dieser Arbeit dargestellten Verfahren ist diese aber nicht auf beliebige Abtastgeometrien anwendbar und wurde hier daher nicht untersucht.

Tomographie mit wenigen Daten

Relativ neu ist nun das Auftreten des Radonschen Inversionsproblems in industriellen Anwendungen wie zum Beispiel der Schallpyrometrie. Der große Unterschied zu den medizinischen Anwendungen liegt in der Zahl der gegebenen Meßdaten: Während man in der medizinischen Tomographie weit über 100.000 Daten pro Bild vorliegen hat (häufig mehrere Millionen), gibt man sich in industriellen Anwendungen mit etwa 20–200 Daten pro Bild zufrieden (in dieser Arbeit sind es typischer Weise 24 Daten).

1 EINLEITUNG

Der Kern aller in dieser Arbeit untersuchten Verfahren sind ”Kollokationsverfahren“, d.h. man macht den Ansatz, die zu rekonstruierende Funktion f als Linearkombination vorgegebener Funktionen φ_{k}, k = 1, …, n zu schreiben,

f(x) = ∑_{k=1}^{n} ξ_{k} · φ_{k}(x),

und stellt nun ein lineares Gleichungssystem zur Berechnung der Koeffizienten ξ_{k} auf.

Typische Kollokationsverfahren sind die ”Direkten Algebraischen Verfahren“, die in [12, Kapitel V.4] und [7] untersucht wurden. A. Trächtler hat bereits ein weiteres solches Verfahren auf die in dieser Arbeit untersuchte Geometrie der RWE Energie AG angewandt (vgl. [19]). Auch die bei der RWE Energie AG bereits verwendeten Verfahren, das sogenannte Fourierverfahren und ART, sind Kollokationsverfahren.

Um a priori Wissen auszunutzen, werden diese Kollokationsverfahren in ein anderes Verfahren eingebaut, das sich an die Konvexen Projektionen anlehnt.

Schallpyrometrie

Die Schallpyrometrie ist ein berührungsloses Temperaturmeßverfahren, mit dem die Temperatur von Gasen gemessen werden kann. Bei den hier vorliegenden Anwendungen geht es um die Temperaturbestimmung in einer horizontalen Ebene im Feuerraum eines großen Dampferzeugers. Die Meßgeräte, die jeweils aus einem Lautsprecher und einem Mikrofon bestehen, werden an den Kesselwänden verteilt. Als Meßwerte erhält man (nach einigen Umrechnungen) die Durchschnittstemperaturen auf den gedachten Verbindungslinien zwischen den Meßeinheiten, es liegt also ein klassisches Tomographieproblem vor. Die für diese Arbeit wesentlichen technischen Details dieser Meßmethode finden sich in Kapitel 2.

In Kapitel 3 sind die verwendeten Rekonstruktionsverfahren und die zugrunde liegende Mathematik dargestellt. Insbesondere werden eine Lösung des Problems von Mehrfachmessungen sowie die Singulärwertzerlegung der diskreten Radontransformation angegeben. Die Güte dieser Verfahren wird schließlich in Kapitel 4 anhand von Phantombildern und echten Meßdaten dargestellt und in der Zusammenfassung abschließend bewertet.

Im Rahmen dieser Arbeit wurde für die RWE Energie AG ein Bericht ([16]) erstellt, der zu großen Teilen mit dieser Diplomarbeit identisch ist, jedoch eine wesentlich ausführlichere Darstellung der mathematischen Grundlagen und eine ausführlichere Untersuchung der numerischen Ergebnisse mit zwei, drei und vier Meßgeräten pro Kesselwand enthält. Da die hier verwendeten Verfahren vom Benutzer noch sehr viele Entscheidungen verlangen (Wahl von Ansatzfunktionen, Pixelungen, Startbildern, diversen Parametern...) wurde für die erstellten Programme eine recht komfortable Benutzeroberfläche zusammengestellt. Die Beschreibung findet sich in [17].

Ich möchte mich an dieser Stelle vielmals bedanken, und zwar zunächst bei Herrn Professor Natterer dafür, daß er mir dieses Thema anvertraut hat und bei Herrn Derichs, RWE Energie AG, Abt. KF-F für die sehr motivierende technische Betreuung.

Der größte Teil der Arbeit für dieses Thema war am Rechner zu leisten, und ohne die Programme von Bernd Setzepfandt und Frank Wübbeling wäre sie wahrscheinlich gar nicht zu bewältigen gewesen. Auch ohne die unzähligen Tips und praktischen Hilfen von Frank Wübbeling und Jörg Fraune wären wohl viele Probleme ungelöst geblieben. Vielen Dank!

Schallpyrometrie

2.1 Beschreibung der Meßmethode

Die Schallpyrometrie ist ein berührungsloses Temperaturmeßverfahren, das sich für die Temperaturbestimmung in einer Ebene eines großen Feuerraumes eignet. Die hier beschriebene Meßapparatur war für mehrere Monate im Block E des Braunkohlekraftwerkes Neurath der RWE Energie AG installiert.

Es wurden acht Meßgeräte, die jeweils aus einem Lautsprecher und einem Mikrofon bestehen, in der 54m–Ebene des Kessels eingebaut. Da an jeder Kesselwand zwei Geräte befestigt waren, ergaben sich 24 Verbindungslinien zwischen den Meßgeräten.

Diese Linien werden im folgenden als Meßpfade bezeichnet. Gemessen werden soll die durchschnittliche Schallgeschwindigkeit auf den Meßpfaden. Dafür arbeitet jeweils ein Meßgerät als Sender und erzeugt ein Schallsignal, welches von den anderen 7 Meßeinheiten empfangen wird. Aus der Laufzeit des Schallsignals läßt sich dann die durchschnittliche Schallgeschwindigkeit berechnen.

Mit Hilfe der Gleichung von Laplace kann man nun bei hinreichend genauer Kenntnis der Gaszusammensetzung im Feuerraum die Temperatur (T) aus der Schallgeschwindigkeit (c) berechnen:

T = \frac{M}{κR}c^{2}  (2.1)

2 SCHALLPYROMETRIE

Das Verhältnis der spezifischen Wärmekapazitäten (κ = c_{p}/c_{v}), die allgemeine Gaskonstante (R) und das Molekulargewicht (M) sind hierbei bekannt.

Nach einem Umlauf, d.h. wenn jede Meßeinheit einmal als Sender fungiert hat, liegen für jeden Meßpfad zwei Meßwerte vor, da jeder Pfad in beiden Richtungen durchschallt wurde. Diese Temperaturwerte werden gemittelt, um Schallmitnahmeeffekte durch die horizontale Rotation des brennenden Kohle/Luftgemisches auszugleichen.

Ursachen von Meßfehlern

Für eine Fehlerabschätzung der Meßdaten sind zahlreiche Einzelaspekte wie zum Beispiel die Unsicherheiten in den als konstant angenommenen Größen M und κ in Gleichung (2.1) zu betrachten. Ausführliche Abschätzungen dieser Fehler sind in [3], [5] und [8] zu finden. Eine Fehlerbetrachtung für die in 2.2 beschriebene Vorgehensweise liegt noch nicht vor.

Auf zwei Fehlerquellen sei jedoch gesondert verwiesen, da sie programmiertechnisch behoben werden können:

1. Schallbeugung: Die Schallpyrometrie beruht auf der höheren Schallgeschwindigkeit in heißeren Teilen des Kessels. Dieser Effekt hat aber auch zur Folge, daß die für die Messung entscheidenden Schallsignale nicht gradlinig laufen, sondern eher einen leichten Bogen durch den heißeren Teil des Gases machen (Fermat’sche Beugung). Aus diesem Grund sind z.B. die an den Wandpfaden gemessenen Temperaturen zu hoch und können nicht verwendet werden. Bei einer bekannten Temperaturverteilung könnte man die tatsächlichen Wege der Schallsignale berechnen. Es wäre also denkbar, zunächst mit den Verfahren aus Kapitel 3 eine Temperaturverteilung aus den Daten zu rekonstruieren um dann festzustellen, wie der tatsächliche Laufweg der Schallsignale bei dieser Verteilung gewesen wäre. Mit diesen korrigierten Informationen über die Meßpfade könnte man das Bild dann erneut rekonstruieren. Nach zwei bis drei Iterationen dieser Art sollte man eine Rekonstruktion erhalten, bei der die Fehlerquelle ”Schallbeugung“ ausgeschlossen ist.

2. Für diese Arbeit wurde die Durchschnittstemperatur pro Meßlinie direkt aus der durchschnittlichen Schallaufzeit pro Meßlinie berechnet. Dies ist streng genommen falsch, da die Temperatur proportional zum Quadrat der Schallgeschwindigkeit ist und damit heiße Stellen einen höheren Einfluß auf die Durchschnittstemperatur gewinnen, als kalte. Bei den hiesigen Anwendungen hat sich jedoch gezeigt, daß der Fehler, den man damit in Kauf nimmt, unter einem halben Prozent liegt. Dieses Problem läßt sich beheben, indem zunächst die lokale Schallgeschwindigkeit z.B. mit Hilfe von Kollokationsverfahren für das gesamte Rekonstruktionsgebiet berechnet wird und dann für jedes Pixel des erhaltenen Bildes die Berechnung der Temperatur erfolgt.

2.2 Aufbereitung der Daten

Die Berechnung der Temperatur aus den von den Meßgeräten gelieferten Laufzeitangaben wurde wie in [5] und [8] beschrieben implementiert und soll hier nicht nochmals dargestellt werden.

Ein Umlauf dauert etwa vierzig Sekunden, an einem Tag fallen also Daten von über 2000 Umläufen an. Es gibt verschiedene Möglichkeiten, diese Datenflut zu bewältigen. Bisher wurden bei der RWE Energie AG jeweils Daten aus 10-14 Umläufen durch Mittelwertbildung zusammengefaßt, wobei ”Ausreißer“, also offensichtliche Fehlmessungen, gestrichen wurden. In [8] wird ein Programm vorgestellt, mit dem die Meßdaten in dieser Form aufbereitet und auf Plausibilität überprüft werden können.

Mit den für diese Arbeit erstellten Programmen ist es möglich, für jeden Pfad eine Funktion zu berechnen, die einen realistischen Verlauf der Durchschnittstemperatur auf diesem Meßpfad für einen Tag darstellt. Diese 24 Funktionen lassen sich für jeden Zeitpunkt auswerten, man kann also z.B. Datensätze berechnen, die dem Temperaturprofil im Kessel entsprechen und einen zeitlichen Abstand von je 30 Sekunden haben. Aus diesen lassen sich dann die entsprechenden Rekonstruktionen der Temperaturverteilungen als Bilder berechnen, welche dann wiederum zu einem Film zusammengefaßt werden können. Auf diese Weise erhält man also einen – im Rahmen der Meß- und Rekonstruktionsfehler – realistischen Temperaturverlauf für einen Zeitraum von z.B. einem Tag.

Im Idealfall liegt für jeden Pfad und jede Richtung alle 40 Sekunden ein neuer Meßwert vor. Um einen Film zu erzeugen, der auch zeitlich hochfrequente Details darstellen kann, braucht man aber für jeden Zeitpunkt eine Approximation an die mittlere Temperatur auf dem jeweiligen Pfad. Bei dieser Approximation sollte bereits eine gewisse ”Glättung“ möglich sein, um realistische (zeitlich) mittlere Temperaturen zu erhalten, und um deutliche Fehlmessungen aussondern zu können.

Dafür werden die Meßwerte durch einen kubischen Spline mit einer vorgebbaren Anzahl von äquidistanten Stützstellen approximiert. Die Koeffizienten der einzelnen Polynome werden so berechnet, daß der Spline im quadratischen Mittel eine Bestapproximation an die gegebenen Meßwerte darstellt, d.h. die Summe der Quadrate der Differenzen von Spline–Wert und Meßwert für alle Meßzeitpunkte wird minimiert. Durch die Vorgabe der Anzahl der Stützstellen des Splines läßt sich nun die (zeitliche) ”Glattheit“ der von den Rekonstruktionsverfahren erzeugten Bilder/Filme sehr gut steuern.

Die mathematischen Details dieser Approximation sind in dieser Arbeit nicht dargestellt, da eine fertige Routine aus einer Programmbibliothek (IMSL) verwendet wurde. Nachzulesen sind diese Details in [1]. Alternativ steht ein Algorithmus zur Verfügung, der eine solche Bestapproximation durch Splines mit statistischen Verfahren berechnet.

Dies ist nun auch die richtige Stelle, um das Herausfiltern offensichtlicher Meßfehler zu behandeln. Fehlmessungen tauchen z.B. auf, wenn ein Lautsprecher ganz ausfällt. In diesem Fall wird die Laufzeit 0.0 ms geliefert, diese Fehler sind also bereits von der direkt zum Meßsystem gehörigen Software detektiert worden, und die Werte können einfach ignoriert werden. Ein größeres Problem stellen sogenannte ”Ausreißer“ (engl. outliners) dar, die aufgrund ihrer Temperaturwerte klar als Fehlmessungen gelten müssen, aber noch in den Daten vorhanden sind. Um auch sie herauszufiltern, werden zunächst alle Daten ignoriert, die zu einer Durchschnittstemperatur unter 500 ^{o}C oder über 1600 ^{o}C führen, da solche Werte nicht auftreten können. Anschließend werden alle verfügbaren Daten wie beschrieben durch einen Spline approximiert, um dann jene Meßwerte aus der Liste zu streichen, welche von diesem Spline um mehr als z.B. 50 ^{o}C abweichen.

Durch die nun noch verbleibenden Meßdaten wird dann der endgültige Spline gelegt.

Abbildung 2.1: Approximation von Meßdaten durch einen Spline

Manchmal entstehen hierbei oder durch den zeitweisen Ausfall eines Meßgerätes recht große (zeitliche) ”Lücken“ in den Meßdaten, so daß für einen Pfad z.B. über 10 Minuten kein Meßwert zur Verfügung steht. Dies führt zu Problemen bei der Berechnung des Splines, da dann eventuell zwischen einigen Knoten des Splines gar keine Meßdaten mehr verfügbar sind. Der Spline fällt dann mangels günstiger Alternativen auf Null ab. Das ist natürlich eine katastrophale Mißinterpretation der Situation und führt zu völlig wertlosen Temperaturrekonstruktionen. Um dies zu verhindern, werden die gestrichenen Meßwerte durch linear interpolierte, künstliche Meßwerte ersetzt. Hierfür werden die jeweils letzten fünf Messungen vor- und hinter der Lücke gemittelt, und zwischen diesen Mittelwerten wird linear interpoliert. Abbildung 2.1 zeigt ein Beispiel, bei dem die Streichung von ”Fehlmessungen“ bereits vorgenommen wurde.

Die Skalen bezeichnen die Uhrzeit (9:00 bis 12:30) und die Temperatur in^{o}C. Die schmale Linie zeigt die Meßwerte, die der Deutlichkeit halber gradlinig verbunden wurden, die breitere Linie stellt den Spline dar, der mit dem erwähnten statistischen Verfahren berechnet wurde. Etwa für den Zeitraum 11:10 bis 11:30 Uhr sind keine Meßwerte vorhanden. Diese wurden daher wie beschrieben ersetzt.

2.3 A priori Informationen

Unter a priori Informationen versteht man Zusatzinformationen, die neben den Meßdaten über das zu rekonstruierende Bild vorliegen.

Aufgrund physikalischer Überlegungen gibt es bei der Schallpyrometrie zahlreiche a priori Informationen. Diese konnten aber bisher nur zum Teil in die beschriebenen Verfahren eingebaut werden. Insbesondere sind sie zumeist nur ”weich“ eingebaut, d.h. es werden nicht zwingend Bilder erzeugt, die diesen Vorgaben entsprechen. Das Einbauen von a priori Informationen ist daher noch eine Aufgabe für die weitere Arbeit an dieser Problemstellung:

• Man weiß aus Messungen mit der Absaugepyrometrie (vgl. z.B. [5]), daß die Temperaturverteilung im Innern des Kessels relativ konstant ist und nur am Rand sehr starke Änderungen erfährt. Diese Information wird durch die in 3.6.2 beschriebene Abwandlung der Konvexen Projektionen (sehr) weich eingebaut.

• Speziell bei der Installation der Meßgeräte im Block E des Kraftwerkes Neurath befanden sich einige Meter unter der Meßebene Einblasstutzen für Ausbrandluft (ABL). Diese sorgten in recht gut bekannten Teilen der Meßebene für ein deutliches Absinken der Temperaturen. Da es sich hierbei um räumlich recht kleine Teile der Ebene handelt, werden diese Meßwerte von den Algorithmen etwas ”fehlinterpretiert“, da diese (aufgrund der vorgenommenen Regularisierung) gar nicht erst versuchen, so kleine Details zu rekonstruieren. Durch die Wahl eines entsprechenden Startbildes lassen sich solche Strukturen einbauen.

• Man kann sich Temperaturschranken vorgeben, die von den Rekonstruktionen eingehalten werden sollen. So kann man z.B. sagen, daß die Maximaltemperatur eines Bildes 1250 ^{o}C nicht überschreiten soll, da dies aus physikalischen Überlegungen heraus bekannt ist.

Rekonstruktionsverfahren

3.1 Problembeschreibung

Die Problemstellung ist hier rein reell. Alle verwendeten mathematischen Konzepte gelten aber ebenso im Komplexen, daher sei IK wahlweise IR oder lC. Für Matrizen werden die im Komplexen üblichen Sprechweisen verwendet.

Das zu rekonstruierende Bild (bzw. die Temperaturverteilung) auf dem Rekonstruktionsgebiet Ω wird als Funktion f : Ω → IK, x ↦ f(x) aufgefaßt. Das Gebiet Ω ist hierbei immer ein beschränktes Gebiet im IK^{2}. Im Falle der Schallpyrometrie ist es ein Rechteck mit dem Länge/Breite–Verhältnis des Kessels. Als Daten über f sind m Linienintegrale über gewisse Linien L_{l} ⊂ Ω, l = 1, …, m gegeben:

g_{l} := ∫_{L_{l}} f(x)dx, l = 1, …, m .

Gleichwertig ist natürlich die Kenntnis des Mittelwertes von f auf jeder Linie. Aus diesen Daten soll eine Rekonstruktion f̃ zu f berechnet werden.

3.2 Kollokationsverfahren

Bei allen Kollokationsverfahren macht man den Ansatz, f̃ als Linearkombination von vorgegebenen Funktionen φ_{k} : Ω ↦ IK, k = 1, …, n darzustellen:

f̃(x) = ∑_{k=1}^{n} ξ_{k} · φ_{k}(x) , ξ_{1}, ..., ξ_{n} ∈ IK.

Geeignete Wahlen für die φ_{k} werden in Abschnitt 3.5 angegeben. Die Koeffizienten ξ_{k} aus dem Ansatz für f̃ sollen jetzt so gewählt werden, daß

∫_{L_{l}} f̃(x)dx ≈ ∫_{L_{l}} f(x)dx , l = 1, …, m

möglichst gut erfüllt ist. Wählt man n ≥ m, so kann f̃ die Daten g_{l} zu f im allgemeinen natürlich beliebig genau erfüllen.

⇝ g_{l} = ∫_{L_{l}} f̃(x)dx = ∫_{L_{l}} ∑_{k=1}^{n} ξ_{k} · φ_{k}(x)dx

3 REKONSTRUKTIONSVERFAHREN

= ∑_{k=1}^{n} ξ_{k} · \underbrace{∫_{L_{l}} φ_{k}(x)dx}_{=:a_{lk}} , l = 1, …, m

⇝ g_{l} = ∑_{k=1}^{n} a_{lk}ξ_{k}  , l = 1, …, m

Der Ansatz führt also auf ein lineares Gleichungssystem mit m Gleichungen und n Unbekannten. Dieses soll nun in Matrixschreibweise ausgedrückt werden:

A := (a_{lk})_{{l=1,...,m \atop k=1,...,n}} = ( \begin{matrix}a_{11} & ⋯ & a_{1n} \\ ⋮ & ... & ⋮ \\ a_{m1} & ⋯ & a_{mn} \end{matrix} ) ∈ IK^{m×n}

g := ( g_{1} ; ⋮ ; g_{m} ) ∈ IK^{m},  ξ := ( ξ_{1} ; ⋮ ; ξ_{n} ) ∈ IK^{n}

Das Gleichungssystem hat damit die Form

Aξ = g .  (3.1)

Hierzu wird die (Moore-Penrose) verallgemeinerte Lösung gesucht. Sie ist definiert als der Vektor ξ_{MP}, der ∥Aξ − g∥_{2} minimiert und unter allen Vektoren mit dieser Eigenschaft noch minimale Norm hat. Sie ergibt sich als eindeutige Lösung der Normalengleichungen

A^{∗}Aξ_{MP} = A^{∗}g mit ξ_{MP} ∈ range(A^{∗}).  (3.2)

Um eine Annäherung an diese zu berechnen, wird die Singulärwertzerlegung der Matrix A verwendet, da A bei tomographischen Anwendungen schlecht konditioniert ist.

Satz 3.2.1 (Singulärwertzerlegung einer Matrix) Sei A ∈ IK^{m×n}. Dann gibt es unitäre Matrizen U = (u_{1}, …, u_{m}) ∈ IK^{m×m}, V = (v_{1}, …, v_{n}) ∈ IK^{n×n}, so daß gilt:

U^{∗}AV = Σ,

mit

Σ = diag(σ_{1}, …, σ_{p}) ∈ IK^{m×n}, p = min(m, n)

und σ_{1} ≥ σ_{2} ≥ … ≥ σ_{p} ≥ 0.

Die rechteckige Diagonalmatrix Σ ist dabei so zu verstehen, daß nur die Diagonale, zu der das Element Σ_{1,1} gehört, besetzt ist.

Beweis: [6, Kapitel 2, Abschnitt 2.5] Im weiteren sei p′ die Anzahl der Singulärwerte von A, welche ungleich Null sind. Für k ≤ p gilt:

A^{∗}u_{k} = σ_{k}v_{k} und Av_{k} = σ_{k}u_{k} .  (3.3)

Außerdem ist

A^{∗}g = ∑_{k=1}^{p} σ_{k}(g, u_{k})v_{k} .  (3.4)

Die gesuchte (Moore-Penrose) verallgemeinerte Lösung ergibt sich mit der Singulärwertzerlegung zu

ξ_{MP} = A^{+}g = ∑_{k=1}^{p′} \frac{1}{σ_{k}} (g, u_{k})v_{k} .

Der Beweis kann analog zum Beweis von Theorem 1.2 aus [12, Kapitel IV.1] geführt werden. An dieser Gleichung erkennt man den Vorteil der Singulärwertzerlegung: Sie separiert das Problem der Lösung von Aξ = g in gut- und schlecht lösbare Anteile. Je kleiner der k’te Singulärwert σ_{k} ist, desto problematischer ist es, den Beitrag \frac{1}{σ_{k}} (g, u_{k})v_{k} bei der Lösung zu berücksichtigen, da Meßfehler in g dabei immer größeren Einfluß auf die Lösung ξ_{MP} haben.

Daher führt man an dieser Stelle einen Filter F_{γ} : IR_{0}^{+} → [0, 1] ein, der diese kritischen Anteile dämpfen soll:

ξ^{γ} := ∑_{k=1}^{p′} F_{γ}(σ_{k})\frac{1}{σ_{k}} (g, u_{k})v_{k} .  (3.5)

Der Leser beachte, daß die Schreibweise des Filters hier etwas anders gewählt wurde, als in [12]. Der wohl einfachste sinnvolle Filter ist

F_{γ}(σ) = \{ \begin{matrix}1 & falls σ ≥ γ, \\ 0 & falls σ < γ. \end{matrix} .

Die Lösung des Gleichungssystems mit diesem Filter wird als ”abgeschnittene Singulärwertzerlegung“ bezeichnet.

In der Anwendung auf die Schallpyrometrie hat sich aber ein anderer Filter besser bewährt, der auf der Methode von Tikhonov-Phillips beruht. Anstelle der Moore-Penrose-Lösung sucht man hierbei nicht einen Minimierer von ∥Aξ − g∥_{2}, sondern fordert, daß für ein γ > 0

h(ξ) := ∥Aξ − g∥_{2}^{2} + γ∥ξ∥_{2}^{2}

minimiert wird. Dies läßt sich z.B. durch die Wahl eines geeigneten Filters erreichen.

Satz 3.2.2 Mit dem Filter F_{γ}^{T P} (σ) = \frac{σ^{2}}{σ^{2}+γ} erhält man aus Gleichung (3.5) dasjenige ξ, welches h(ξ) minimiert.

Beweis: Die Funktion h : IK^{n} → IR ist reell differenzierbar, gesucht ist also ein ξ_{0}, für das der Gradient von h verschwindet.

∇h(ξ) = 2(Aξ − g)^{∗}A + 2γξ^{∗}

⇝ ∇h(ξ_{0}) = 0 ⇔ A^{∗}(Aξ_{0} − g) + γξ_{0} = 0
⇔ (A^{∗}A + γI)ξ_{0} = A^{∗}g

Da A^{∗}A positiv semidefinit ist, ist (A^{∗}A+γI) positiv definit und damit invertierbar. Es gilt wegen den Gleichungen (3.3):

(A^{∗}A + γI)v_{k} = (σ_{k}^{2} + γ)v_{k}
⇔ (A^{∗}A + γI)^{−1}v_{k} = \frac{1}{σ_{k}^{2} + γ} v_{k} .

Mit Gleichung (3.4) hat man nun:

ξ_{0} = (A^{∗}A + γI)^{−1} ∑_{k=1}^{p} σ_{k}(g, u_{k})v_{k}
= ∑_{k=1}^{p} σ_{k}(g, u_{k})\frac{1}{σ_{k}^{2} + γ} v_{k}
= ∑_{k=1}^{p} \frac{σ_{k}^{2}}{σ_{k}^{2} + γ} \frac{1}{σ_{k}} (g, u_{k})v_{k}
= ∑_{k=1}^{p} F_{γ}^{T P} (σ_{k})\frac{1}{σ_{k}} (g, u_{k})v_{k} .

Wegen lim_{∥ξ∥→∞} h(ξ) = ∞ ist ξ_{0} Minimumstelle, der Satz ist damit bewiesen.

□

Da zur numerischen Berechnung der Singulärwertzerlegung Unterprogramme aus mathematischen Programmbibliotheken verwendet werden können, gehe ich auf diesen Punkt nicht weiter ein. Die für die Lösung des Gleichungssystems empfehlenswerten Hilfsmittel stehen jetzt also zur Verfügung.

Schließlich sei noch erwähnt, daß Satz 3.2.2 sicherlich nicht zu der — auf einem Computer — schnellst-möglichen Implementierung der Methode von Tikhonov– Phillips führt, er gestattet aber ein sehr gutes Verständnis dieser Methode. Der Zeitfaktor spielt hier keine Rolle, da die Berechnung der Singulärwertzerlegung im Falle m ≤ 100 nur wenige Sekunden dauert und nur einmal nötig ist.

Der folgende Abschnitt behandelt die Verkleinerung des Gleichungssystems bei bestimmten Meßgeometrien.

3.3 Mehrfachmessungen

In der Praxis treten zuweilen Abtastgeometrien auf, in denen einzelne Linien doppelt vorkommen, d.h. während für alle anderen Linien nur ein Meßwert bestimmt wird, werden die Integrale über diese Linien zweifach gemessen.

Die Zusammenfassung solcher Mehrfachmessungen zu einer Messung, bei der das arithmetische Mittel der Meßwerte verwendet wird, würde einen unnötigen Verzicht auf in den Daten enthaltene Informationen bedeuten. Denn wenn z.B. der zu erwartende Meßfehler bei allen Einzelmessungen gleich groß ist, dann liegen für die Linien mit Mehrfachmessungen genauere Daten vor als für die übrigen Linien. Dies würde bei der Rekonstruktion aber nicht beachtet.

Dieser Abschnitt beantwortet die Frage, wie man bei Kollokationsverfahren die mehrfachen Meßwerte eingehen lassen kann, ohne das entstehende Gleichungssystem unnötig zu vergrößern.

3.3.1 Kanonische Aufstellung des Gleichungssystems

Wenn alle Messungen gleich gewichtet werden sollen, dann ist dies zu erreichen, indem man das Gleichungssystem wie in 3.2 erläutert aufstellt und die Mehrfachmessungen jeweils wie die Einzelmessungen – also unabhängig voneinander – einbringt. Linien, die hierbei mehrfach vorkommen, führen bei der Systemmatrix A zu Zeilen, die jeweils identisch sind. Die Meßwerte selbst werden (unbehandelt) im Vektor b gesammelt. Gesucht ist dann wieder eine verallgemeinerte Lösung des Gleichungssystems

Ax = b.

Es wird die durch die Normalengleichungen

A^{∗}Ax_{MP} = A^{∗}b und x_{MP} ∈ range(A^{∗})  (3.6)

eindeutig bestimmte Moore-Penrose-Lösung verwendet. Im folgenden wird die Matrix A auf eine kleinere Matrix Ã reduziert, in der keine Zeilen mehrfach vorkommen, die jedoch auf dieselbe Lösung x_{MP} führt wie die Gleichung (3.6).

3.3.2 Reduktion des Gleichungssystems

Der folgende Satz beinhaltet für Kollokationsverfahren die Lösung des Problems der Mehrfachmessungen. Da sich die Moore-Penrose-Lösung beim Vertauschen von Gleichungen nicht ändert, kann ohne Beschränkung der Allgemeinheit davon ausgegangen werden, daß die letzten k Zeilen von A übereinstimmen. Es gelte A ∈ IK^{m×n} und r = m − (k − 1). Es zeigt sich, daß man dann die letzten (k −1) Zeilen streichen und den Wert b_{r} durch \frac{1}{k} ∑_{i=r}^{m} b_{i} ersetzen kann. Wird nun noch die r-te Gleichung mit dem Faktor \sqrt{k} multipliziert, dann erhält man mit diesem kleineren Gleichungssystem die alte Moore-Penrose-Lösung. Dies besagt der folgende Satz.

Satz 3.3.1 Seien A = (a_{ij}) ∈ IK^{m×n}, b ∈ IK^{m} wobei die letzten k Zeilen von A identisch sind, r := m − (k − 1). Es sei x_{MP} definiert durch

A^{∗}Ax_{MP} = A^{∗}b und x_{MP} ∈ range(A^{∗}).

Nun wähle man Ã = (ã_{ij}) ∈ IK^{r×n} mit

ã_{ij} = a_{ij} und ã_{rj} = \sqrt{k} a_{rj}

für 1 ≤ i < r, 1 ≤ j ≤ n sowie b̃ = (b̃_{1}, ..., b̃_{r})^{t} ∈ IK^{r} mit

b̃_{i} = b_{i} und b̃_{r} = \sqrt{k} \frac{1}{k} ∑_{j=r}^{m} b_{j} = \frac{1}{\sqrt{k}} ∑_{j=r}^{m} b_{j}

für 1 ≤ i < r. Dann gilt

Ã^{∗}Ãx_{MP} = Ã^{∗}b̃ und x_{MP} ∈ range(Ã^{∗}).

Beweis: Das Bild einer Matrix ist gleich dem von ihren Spalten aufgespannten Vektorraum. Es gilt also range(A^{∗}) = range(Ã^{∗}), da Ã^{∗} aus A^{∗} gewonnen wurde, indem identischen Spalten gestrichen wurden und eine Spalte mit dem Faktor \sqrt{k} ≠ 0 multipliziert wurde. Daher genügt es, A^{∗}A = Ã^{∗}Ã und A^{∗}b = Ã^{∗}b̃ zu zeigen. Seien hierfür i, j ∈ {1, ..., n}:

(A^{∗}A)_{ij} = ∑_{l=1}^{m} a_{li}a_{lj}
= ∑_{l=1}^{r−1} a_{li}a_{lj} + \underbrace{k · a_{ri}a_{rj}}_{=(\sqrt{k}a_{ri})·(\sqrt{k}a_{rj})}

= ∑_{l=1}^{r−1} ã_{li}ã_{lj} + ã_{ri}ã_{rj}
= ∑_{l=1}^{r} ã_{li}ã_{lj}
= (Ã^{∗}Ã)_{ij}

(A^{∗}b)_{j} = ∑_{l=1}^{m} a_{lj}b_{l}
= ∑_{l=1}^{r−1} a_{lj}b_{l} + ∑_{l=r}^{m} a_{lj}b_{l}
= ∑_{l=1}^{r−1} ã_{lj}b̃_{l} + \underbrace{a_{rj}\sqrt{k}}_{=ã_{rj}} · \underbrace{\frac{1}{\sqrt{k}} ∑_{l=r}^{m} b_{l}}_{=b̃_{r}}

= ∑_{l=1}^{r} ã_{lj}b̃_{l}
= (Ã^{∗}b̃)_{j}
□

Die Multiplikation einer Gleichung mit einer reellen Zahl hat also Auswirkungen auf die verallgemeinerte Lösung des Gleichungssystems. Dies wirft die Frage nach einer geeigneten Normierung der Koeffizienten auf. Dies wurde aber bei den hier betrachteten Gleichungssystemen nicht weiter ”künstlich“ beeinflußt, da in allen betrachteten Fällen die Ansatzfunktionen jeweils vergleichbare Größenordnungen haben und damit die Zeilen der Systemmatrix A als sinnvoll normiert angenommen werden können. Sollte bekannt sein, daß manche Meßwerte mit größeren Fehlern behaftet sind, als andere, so ließe sich dies aber an dieser Stelle berücksichtigen, indem die entsprechenden Gleichungen schwächer gewichtet werden.

3.4 Singulärwertzerlegung der diskreten Radontransformation

Wesentlicher Gegenstand der Untersuchungen in diesem Abschnitt sind diskrete Abwandlungen der Radontransformation R. Ist Ω ⊂ IR^{n} beschränkt, so ordnet diese einer L_{2}(Ω)–Funktion ihre Integrale über sämtliche Hyperebenen des IR^{n} zu:

R : L_{2}(Ω, W^{−1}) ⟶ L_{2}(S^{n−1} × IR, w^{−1})
Rf(θ, s) = ∫_{x·θ=s} f(x)dx = ∫_{θ^{⊥}} f(sθ + y)dy

Typische Wahlen der Gewichtsfunktionen für die L_{2}–Räume sind W = idIR^{n} und w(s) = \sqrt{1 − s^{2}}. Diese werden im folgenden jedoch keine Rolle spielen.

Abgesehen von der Schreibweise ist die Radontransformation im zweidimensionalen identisch mit der zweidimensionalen X-Ray-Transformation

P f(θ, x) = ∫_{−∞}^{∞} f(x + tθ)dt,

die jeder Funktion sämtliche Linienintegrale über Geraden im IR^{n} zuordnet.

Angesichts der tomographischen Anwendungen wäre es naheliegender, mit der X-Ray-Transformation anstelle der Radontransformation zu arbeiten, da die Schallpyrometrie ebenso wie die Meßgeräte der medizinischen Tomographie bei dreidimensionalen Problemstellungen immer noch Linienintegrale messen und keine Flächenintegrale. Da in dieser Arbeit jedoch nur zweidimensionale Probleme behandelt werden und in der Literatur üblicherweise mit der Radontransformation gearbeitet wird, soll dies auch hier geschehen. Im folgenden wird daher nur noch die zweidimensionale Radontransformation verwendet, ohne daß dies in der Schreibweise nochmals betont wird.

In praktischen Anwendungen wird niemals die ganze Radontransformation einer Funktion gemessen, sondern die Meßgeräte liefern immer nur endlich viele Linienintegrale. Daher ist die Betrachtung der diskreten Radontransformation sinnvoll. Jede Gerade im IR^{2} läßt sich beschreiben durch ein Tupel (θ, s) ∈ S^{1} × IR. Ist L = {L_{l} := (θ_{l}, s_{l}) ∈ S^{1} × IR : l = 1, ..., m} ein Satz von m Linien und H ein endlichdimensionaler Unterraum von L_{2} mit Skalarprodukt < ·, · >_{H}, so bezeichne ich die Abbildung

\begin{matrix}R_{d} : H & ⟶ & IR^{m} \\ (R_{d}f)_{l} & = & ∫_{L_{l}} f(x)dx \end{matrix}  (3.7)

als diskrete Radontransformation. Die Bedingung dim(H) < ∞ wird hier verlangt, da im folgenden nur für diesen Fall die Singulärwertzerlegung von R_{d} angegeben wird.

Die Problemstellung der Tomographie lautet also nun: Gegeben ist eine diskrete Radontransformation R_{d} einer Funktion f. Rekonstruiere f!

Dieses Problem ist in zweierlei Hinsicht schlecht gestellt: Erstens ist R_{d} im (interessierenden) Fall dim(H) > m nicht injektiv, es gibt also keine eindeutige Lösung, zweitens hängt die Lösung des Inversionsproblems der (kontinuierlichen) Radontransformation nicht stetig von den Daten ab. Dieses Problem überträgt sich auf R_{d}, es ist also eine Regularisierung nötig. Bei der Wahl von H bieten sich zwei Wege an:

1. Falls aus der Anwendung heraus ein hinreichend kleindimensionaler Raum von Lösungsfunktionen bekannt ist, sollte man diesen natürlich verwenden.

2. Gibt die Anwendung keinen geeigneten Raum H vor, so muß ein Raum H ⊂ L_{2} gefunden werden, in dem das Problem eindeutig und mit möglichst guter Kondition gelöst werden kann. Darüber hinaus sollten die Lösungen möglichst ”glatt“ sein.

Im folgenden wird ein Unterraum H ⊂ L_{2} gesucht, der den Anforderungen aus 2. genügt. Hierfür benötigen wir die Singulärwertzerlegung eines linearen Operators.

Definition 3.4.1 Es sei A ein linearer Operator zwischen separablen Hilberträumen (H_{1}, <, >_{H_{1}}), (H_{2}, <, >_{H_{2}}):

A : H_{1} ⟶ H_{2}

Das Tupel {(v_{n}), (u_{n}), (σ_{n})}_{n≥0} heißt Singulärwertzerlegung (kurz: SVD) von A, falls

1. (v_{n})_{n≥0} ist ein vollständiges Orthonormalsystem von H_{1},

2. (u_{n})_{n≥0} ist ein Orthonormalsystem in H_{2},

3. (σ_{n})_{n≥0} ist eine monoton fallende Folge nicht-negativer reeller Zahlen,

4. ∀n ∈ IN : Av_{n} = σ_{n}u_{n} und A^{∗}u_{n} = σ_{n}v_{n}.

Falls nicht beide Hilberträume unendlich–dimensional sind, so sind die entsprechenden Folgen endlich. Damit Punkt 4 noch sinnvoll ist, müssen die Folgen dann eventuell durch Folgenglieder σ_{n} = 0 und u_{n} = 0 bzw. v_{n} = 0 ergänzt werden.

Die Singulärwertzerlegung der Radontransformation ist bekannt, vergleiche z.B. [11] und [12]. In der Tomographie mit wenigen Daten erscheint es aber riskant, solche asymptotischen Ergebnisse auf Abtastgeometrien mit weniger als 100 Linien zu übertragen. Im folgenden wird daher die Singulärwertzerlegung für die diskrete Radontransformation hergeleitet. Dies stellt sich als relativ leicht heraus, da die diskrete Radontransformation eine lineare Abbildung zwischen endlichdimensionalen Vektorräumen darstellt und damit auf die bekannte Singulärwertzerlegung von Matrizen zurückgegriffen werden kann. Diese wurde im Satz 3.2.1 bereitgestellt. Die Spaltenvektoren von V und U bilden dabei die in Definition 3.4.1 geforderten Orthonormalsysteme und erfüllen die Bedingung 4.

Ist nun die diskrete Radontransformation R_{d}f = g gegeben, so wird eine Funktion f̃ ∈ H gesucht, für welche g = R_{d}f̃ gilt und die im Raum aller Lösungen eine minimale Norm hat. Sei dim(H) = n und H = span(φ_{1}, ..., φ_{n}). Dann gibt es ein ξ = (ξ_{1}, …, ξ_{n})^{t} ∈ IR^{n} mit

f̃ = ∑_{k=1}^{n} ξ_{k}φ_{k}

Man erhält nun wie in 3.2 beschrieben ein Gleichungssystem und eine entsprechende Systemmatrix. Mit Hilfe dieser Matrix kann man nun die Singulärwertzerlegung der diskreten Radontransformation gewinnen, falls zusätzlich vorausgesetzt wird, daß eine Orthonormalbasis von H gegeben ist. Dies liefert der folgende

Satz 3.4.1 (Singulärwertzerlegung von R_{d}) Sei χ_{1}, ..., χ_{n} eine Orthonormalbasis von H, A = (a_{lk}) ∈ IK^{m×n} mit

a_{lk} = (R_{d}χ_{k})_{l} für l = 1, ..., m, k = 1, ..., n.

Sei weiter A = UΣV^{∗} die SVD von A gemäß Satz 3.2.1, V = (v_{kj}) und schließlich

ψ_{j} = ∑_{k=1}^{n} v_{kj}χ_{k} , j = 1, ..., n.

Dann ist

{(ψ_{j})_{j=1,...,n}, (u_{k})_{k=1,...,m}, (σ_{k})_{k=1,...,min(m,n)}}

eine Singulärwertzerlegung von R_{d}.

Beweis: Zu Zeigen:

1. {ψ_{j}} ist eine vollständiges orthonormales System von H

2. {u_{j}} ist ein orthonormales System in IR^{m}

3. σ_{1} ≥ σ_{2} ≥ ... ≥ σ_{n} ≥ 0

4. (a) R_{d}ψ_{j} = σ_{j} · u_{j} (b) R_{d}^{∗}u_{j} = σ_{j} · ψ_{j}

Zu 1: Die Orthonormalität der ψ_{j} erhält man aus der Orthonormalität der χ_{k} und der Orthogonalität der Matrix V: (seien i, j ∈ {1, ..., n})

< ψ_{j}, ψ_{i} >_{H} = < ∑_{k=1}^{n} v_{kj}χ_{k}, ∑_{l=1}^{n} v_{li}χ_{l} >_{H}
= ∑_{k,l=1}^{n} v_{kj} · v_{li} \underbrace{< χ_{k}, χ_{l} >_{H}}_{=δ_{kl}}

= ∑_{k=1}^{n} v_{kj} · v_{ki}
= δ_{ji} ,  da V · V^{T} = I_{n}.

Die Vollständigkeit ergibt sich aus dim(H) = #{ψ_{j}, j = 1, ..., n}. Zu 2, 3: Da A = UΣV^{∗} eine Singulärwertzerlegung von A ist, ist dies bereits bekannt. Zu 4a):

(R_{d}ψ_{j})_{l} = ∫_{L_{l}} ψ_{j}(x)dx = ∫_{L_{l}} ∑_{k=1}^{n} v_{kj}χ_{k}(x)dx
= ∑_{k=1}^{n} v_{kj} · ∫_{L_{l}} χ_{k}(x)dx
= ∑_{k=1}^{n} v_{kj} · a_{lk} = ∑_{k=1}^{n} a_{lk} · v_{kj}
= (A · V)_{lj} = (U · Σ)_{lj}
= (σ_{j} · u_{j})_{l} .

Für 4b) wird gezeigt, daß die Abbildung F : IR^{m} → H, ∑ α_{j}u_{j} ↦ ∑ α_{j}σ_{j}ψ_{j} mit R_{d}^{∗} identisch ist. Es sei < ·, · >_{m} das Skalarprodukt im IR^{m}. Z.z.: ∀h ∈ H, u ∈ IR^{m} : < R_{d}h, u >_{m}=< h, F u >_{H} Sei h = ∑_{i=1}^{n} β_{i}ψ_{i}, u = ∑_{j=1}^{m} α_{j}u_{j}. Wegen

< R_{d}h, u >_{m} = < R_{d} ∑_{i=1}^{n} β_{i}ψ_{i}, ∑_{j=1}^{m} α_{j}u_{j} >_{m}
= ∑_{i=1}^{n} ∑_{j=1}^{m} β_{i}α_{j} < R_{d}ψ_{i}, u_{j} >_{m}
< h, F u >_{H} = < ∑_{i=1}^{n} β_{i}ψ_{i}, F ∑_{j=1}^{m} α_{j}u_{j} >_{H}
= ∑_{i=1}^{n} ∑_{j=1}^{m} β_{i}α_{j} < ψ_{i}, F u_{j} >_{H}

genügt es, zu zeigen, daß < R_{d}ψ_{i}, u_{j} >_{m}=< ψ_{i}, F u_{j} >_{H} gilt:

< R_{d}ψ_{i}, u_{j} >_{m} =^{4a)}< σ_{i}u_{i}, u_{j} >_{m}= σ_{i}δ_{ij} = σ_{j}δ_{ij} =< ψ_{i}, σ_{j}ψ_{j} >_{H} =^{def}< ψ_{i}, F u_{j} >_{H} .

□

Die Singulärfunktionen (ψ_{j})_{j} erfüllen nun zwei Orthogonalitätsrelationen. Zum einen gilt < ψ_{i}, ψ_{j} >_{H}= δ_{ij}, d.h. sie bilden eine Orthonormalbasis von H. Definiert man sich weiterhin eine Bilinearform auf H mit Hilfe der diskreten Radontransformation durch < φ, ψ >_{R_{d}}:=< R_{d}φ, R_{d}ψ >_{m}, so gilt auch < ψ_{i}, ψ_{j} >_{R_{d}}= δ_{ij} aufgrund der Orthogonalität der Matrix V.

Wählt man im Satz anstelle der Orthogonalbasis χ_{1}, ..., χ_{n} eine beliebige Basis φ_{1}, ..., φ_{n} von H, so sind die ψ_{j} orthonormal bezüglich des Skalarproduktes, bezüglich dem die Basis φ_{1}, ..., φ_{n} orthonormal ist.

Die Berechnung einer Rekonstruktion mittels Gleichung (3.5) ist damit gleichwertig zur direkten Verwendung der Lösungsformel mit der Singulärwertzerlegung der diskreten Radontransformation:

f̃(x) = ∑_{k=1}^{p′} F_{γ}(σ_{k}) · \frac{1}{σ_{k}} · (g, u_{k}) · ψ_{k}(x)

Wird also z.B. die Singulärwertzerlegung nach dem zehnten Singulärwert abgeschnitten, dann werden genau die ersten 10 Singulärfunktionen (vgl. Abb. 3.1) zur Berechnung der Rekonstruktion verwendet.

3.5 Wahl geeigneter Ansatzfunktionen

3.5.1 Der Kovarianzenansatz

Bei diesem Ansatz wird die Bildrekonstruktion aus einer völlig anderen Richtung angegangen, und zwar als statistisches Problem. Dies wird in [12] dargestellt. Das Bild f und die Daten g werden als Familien von normalverteilten Zufallsvariablen angesehen, wobei f von dem Parameter x abhängt. Zwischen ihnen besteht der Zusammenhang

g = Rf + n .

Hierbei repräsentiert n eine normalverteilte Familie von Zufallsvariablen mit Mittelwert 0, welche von f und g unabhängig ist und die Meßfehler darstellt. R ist ein linearer Operator. Der Mittelwert von f sei f̄ und die Kovarianz von f ist

F (x, y) = E(f(x) − f̄(x))(f(y) − f̄(y))

Die Kovarianz von n werde mit Σ bezeichnet. Ein gängiges Modell in der Bildverarbeitung ist die Annahme von

F (x, y) = σe^{−λ|x−y|} .

Der im Sinne einer Kleinste-Quadrate Approximation beste lineare Schätzer ergibt sich zu

f_{B} = f̄+ F R^{∗}ξ mit ξ = (RF R^{∗} + Σ)^{−1}(g − Rf̄) ∈ IR^{m} .  (3.8)

Sei nun L̃_{l} ⊂ Ω ein ”schmaler“ Streifen im Rekonstruktionsgebiet Ω, der die l’te Linie enthält und χ_{l} die charakteristische Funktion dieses Streifens. Für diesen Abschnitt sei der Operator R : L_{2} → IR^{m} definiert durch

(Rf)_{l} = ∫L̃_{l} f(x)dx = ∫_{Ω} χ_{l}(x)f(x)dx .

Für sehr schmale Streifen L̃_{l} und hinreichend glatte Funktionen f kann dies als gleichwertig zur Kenntnis von ∫_{L_{l}} f(x)dx angesehen werden. Diese Ungenauigkeit erspart an dieser Stelle die Einführung der Dirac’schen δ-Funktion und damit den umfangreichen Theorieapparat der Distributionen.

Mit f ∈ L_{2} und g ∈ IR^{m} ist dann R̃g := ∑_{i=1}^{m} g_{i}χ_{i}(x) wegen

< f, R̃g > = ∫_{Ω} f(x) · R̃g(x)dx
= ∑_{i=1}^{m} ∫_{Ω} f(x)χ_{i}(x)dx · g_{i}
= < Rf, g >

der zu R adjungierte Operator, also R̃ = R^{∗}. Um auf Funktionen angewandt werden zu können, wird F in Gleichung (3.8) durch den Integraloperator

F f(x) = ∫_{Ω} f(y)e^{−λ|x−y|}dy

ersetzt. Hierbei wird σ = 1 gesetzt, da dieser Koeffizient die Lösung des Gleichungssystems nicht beeinflußt. Für n wird weißes Rauschen angenommen, d.h. die Meßfehler sind unabhängig voneinander. Damit ist Σ = γI, γ > 0. Der in (3.8) auftauchende Operator F R^{∗} ergibt sich also zu:

F R^{∗}ξ = F ∑_{i=1}^{m} ξ_{i}χ_{i}(x)
= ∫_{Ω} ∑_{i=1}^{m} ξ_{i}χ_{i}(x)e^{−λ|x−y|}dy
= ∑_{i=1}^{m} ξ_{i} \underbrace{∫L̃_{l} e^{−λ|x−y|}dy}_{=:φ̃_{i}(x)}

= ∑_{i=1}^{m} ξ_{i}φ̃_{i}(x) .

Damit ist

(RF R^{∗}ξ)_{j} = (R ∑_{i=1}^{m} ξ_{i}φ̃_{i}(x))_{j} = ∑_{i=1}^{m} ξ_{i} \underbrace{∫L̃_{j} φ̃_{i}(x)dx}_{=a_{ji}} = ∑_{i=1}^{m} a_{ji}ξ_{i} = (Aξ)_{j}

RF R^{∗} kann also mit der Matrix A aus Gleichung (3.1) identifiziert werden. Abgesehen von der Verschiebung der Meßdaten um Rf̄ wird in (3.8) also folgendes Problem gelöst:

f_{B} = ∑_{i=1}^{m} ξ_{i}φ̃_{i}(x) mit ξ = (A + γI)^{−1}g

Diese Methode ist für symmetrische, positiv definite Matrizen als ”Lavrent’ev“– Regularisierung bekannt. Die Matrix A ist hier symmetrisch, da

a_{ij} = ∫_{L_{i}} ∫_{L_{j}} e^{−λ|x−y|}dydx .

Für großes λ ist die Matrix diagonal dominant mit positiven Diagonalelementen und damit positiv definit. Für symmetrische Matrizen ist die Singulärwertzerlegung identisch mit der Schurschen Normalform, es gilt also U = V. Damit kann man ähnlich wie im Beweis von Satz 3.2.2 nachrechnen, daß sich die Lavrent’ev–Regularisierung mit Gleichung (3.5) und dem Filter

F_{γ}(σ) = \frac{σ}{σ + γ}

berechnen läßt. Für die Berechnung von Rekonstruktionen wurde in dieser Arbeit generell die Methode von Tikhonov–Phillips verwendet, in Abschnitt 4.1.3 wird aber auch auf die Regularisierung nach Lavrent’ev eingegangen.

Aus Effizienzgründen werden die Streifen L̃_{l} für die Berechnung der Ansatzfunktionen wieder auf die Linien L_{l} reduziert.

Diese Sichtweise des Problems bringt also die Empfehlung, Ansatzfunktionen der Form

φ_{k}(x) = ∫_{L_{k}} e^{−λ|x−y|}dy

zu verwenden.

3.5.2 Andere Ansatzfunktionen

Bei der Untersuchung verschiedener Ansatzfunktionen haben sich die folgenden Kriterien für gute Rekonstruktionsverfahren als wünschenswert und erfüllbar erwiesen:

1. Es sollten nicht grundsätzlich problemunabhängige Strukturen in den Rekonstruktionen sichtbar sein.

2. Die Rekonstruktionen sollten möglichst ”glatt“, d.h. niederfrequent sein.

3. Eine einzelne Rekonstruktion sollte schnell berechnet werden können.

ART

ART ist ein Klassiker unter den Tomographieverfahren. Hierbei wird das Rekonstruktionsgebiet in p × q Pixel eingeteilt, die von 1 bis pq durchnummeriert werden. Die Ansatzfunktionen φ_{k} haben dann im k’ten Pixel den Wert 1 und sonst den Wert 0, und das entstehende Gleichungssystem wird mit dem Verfahren von Kaczmarz gelöst.

Diese Wahl der Ansatzfunktionen ist jedoch keineswegs optimal. Falls ein Pixel von keiner Linie getroffen wird, dann liegen für dieses Pixel keine Informationen vor und als Wert in der Rekonstruktion wird 0 gewählt. Daher muß die Pixelung so grob vorgenommen werden, daß jedes Pixel von mindestens einer Linie getroffen wird. Bei der in Kapitel 2 dargestellten 24-Linien-Geometrie kann man daher kaum mehr als 5 × 5 Pixel wählen.

Diese grobe Pixelung stellt aber bereits eine unnötige, problemunabhängige Einschränkung dar und führt zu unrealistischen Bildern.

Fourierverfahren

Beim Fourierverfahren — wie es bisher bei der Schallpyrometrie eingesetzt wird — wird die Rekonstruktion f̃ in eine zweidimensionale Kosinusreihe entwickelt und diese wird z.B. nach den ersten 4 · 4 = 16 Koeffizienten abgeschnitten, so daß ein überbestimmtes Gleichungssystem zu lösen ist. Die Ansatzfunktionen haben also die Form

φ_{k_{1},k_{2}} = cos(\frac{π}{k_{1}} x) · cos(\frac{π}{k_{2}} y), k_{1}, k_{2} = 1, 2, 3, ...

falls auf [0, 1] × [0, 1] rekonstruiert wird.

Das Problem bei diesem Verfahren sind die zu schnell wachsenden Frequenzen der Ansatzfunktionen. Die Ansatzfunktionen weisen sehr bald Detailgrößen auf, die deutlich unterhalb der zu erwartenden Auflösungsgrenze liegen und führen damit zu z.T. katastrophalen Artefakten. Die Kosinusreihe muß deshalb sehr früh abgebrochen werden und kann dann den Informationsgehalt der Daten nicht mehr ausschöpfen, da überbestimmte Gleichungssysteme gelöst werden müssen.

Die Ansatzfunktionen φ_{k}^{(1)} = e^{−λ|x−x_{k}|} und φ_{k}^{(2)} = e^{−λ|x−x_{k}|^{2}}

Mit diesen Funktionen wurden — neben dem Kovarianzenansatz — die besten Ergebnisse erzielt. Das Rekonstruktionsgebiet wird dabei in p × q Pixel eingeteilt und die x_{k} sind die Mittelpunkte dieser Pixel. Die Pixelung sollte dabei so fein gewählt werden, daß deutlich mehr Ansatzfunktionen als Meßlinien vorliegen (pq ≥ 4m hat sich als günstig erwiesen). Es wird also ein stark unterbestimmtes Gleichungssystem gelöst. Da hierfür die Moore-Penrose Lösung — also eine Minimum-Norm Lösung — bzw. die Methode von Tikhonov–Phillips verwendet wird, erfüllen die Rekonstruktionen automatisch die gewünschten Glattheitsbedingungen und weisen keine fatal kleine Detailgröße auf, falls λ nicht zu groß gewählt wird.

Bei hinreichend feiner Pixelung ergeben sich bei den beiden Funktionentypen annähernd identische Singulärfunktionen. Die Funktionen φ_{k}^{(2)} sind jedoch numerisch schneller berechenbar, weil keine Wurzel berechnet werden muß.

Der Nachteil dieser Verfahren ist jedoch der hohe Rechenaufwand, da ein großes Gleichungssystem zu lösen ist und sehr viele Funktionsauswertungen vorgenommen werden müssen, um f̃(x) = ∑_{k=1}^{n} ξ_{k} · φ_{k}(x) zu berechnen.

Um diese Funktionsauswertungen zu beschleunigen, können die Funktionen einfach einmalig in der gewünschten Pixelung in einer (dann sehr großen) Datei abgespeichert werden. Nun ist es aber erst recht wünschenswert, die Zahl der Ansatzfunktionen möglichst klein zu halten.

Eine Lösung besteht darin, die Singulärfunktionen abzuspeichern. Dies ist bei etwa bis zu 1000 Meßlinien mit vertretbarem Rechenaufwand möglich, stellt sich aber als unnötig aufwendig heraus, da hierfür zunächst die Singulärwertzerlegung einer relativ großen Matrix berechnet werden muß. Die folgenden Überlegungen liefern eine geeignetere Basis des relevanten Unterraums von span(φ_{1}, .., φ_{n}).

Die Lösung f̃ wird als f̃(x) = ∑_{k=1}^{n} ξ_{k} · φ_{k}(x) berechnet, wobei ξ die Moore-Penrose Lösung von Aξ = g ist. ξ erfüllt also

A^{∗}Aξ = A^{∗}g und ξ ∈ range(A^{∗}).

Sind a_{1}^{t}, ..., a_{m}^{t} die Zeilen von A, so gilt range(A^{∗}) = span(a_{1}, ..., a_{m}). Daher existieren α_{1}, ..., α_{m} ∈ IR derart, daß ξ = ∑_{i=1}^{m} α_{i}a_{i} also ξ_{k} = ∑_{i=1}^{m} α_{i}a_{ik}. Damit gilt

f̃(x) = ∑_{k=1}^{n} ξ_{k} · φ_{k}(x) = ∑_{k=1}^{n} ∑_{i=1}^{m} α_{i}a_{ik} · φ_{k}(x) = ∑_{i=1}^{m} α_{i} \underbrace{∑_{k=1}^{n} a_{ik} · φ_{k}(x)}_{=: η_{i}(x)} .

Die Lösungen f̃ liegen also in span(η_{1}, ..., η_{m}). Nach der Berechnung der Matrix A können die η_{i} berechnet und abgespeichert werden, um dann als Ansatzfunktionen zu dienen. Die resultierende Systemmatrix ist aus IR^{m×m} und hat damit minimale Größe.

Die η_{k} haben sehr große Ähnlichkeit mit den φ_{k} = ∫_{L_{l}} e^{−λ|x−y|}dy aus dem Kovarianzenansatz (vgl. Abb. 3.4) und führen daher zu ähnlichen Rekonstruktionen. Die Ansätze können also als gleichwertig gelten.

Über zwei Parameter wurde noch fast gar nichts gesagt: Die Feinheit der Pixelung und die Wahl von λ. Beide haben sich als wenig einflußreich erwiesen. Die Pixelung sollte so fein gewählt werden, daß die im folgenden Abschnitt besprochenen hochfrequenten Singulärfunktionen noch aus den vorgegebenen Ansatzfunktionen zusammengesetzt werden können. Es hat sich hierfür bewährt, n ≈ 4m zu wählen, also etwa viermal soviele Ansatzfunktionen zur Verfügung zu stellen, wie Meßlinien vorhanden sind. Eine weitere Verfeinerung der Pixelung sorgt bei deutlich steigendem Rechenaufwand für quasi identische Ergebnisse. Überraschender Weise besitzt auch der Parameter λ kaum Einfluß auf Rekonstruktionen, falls er nicht extrem zu klein oder zu groß gewählt wird. Ist er zu klein, dann unterscheiden sich die Ansatzfunktionen kaum und die Kondition der Systemmatrix wird sehr schlecht. Ist er zu groß, dann werden die Singulärfunktionen zu reinen Abbildern der Abtastgeometrie. Dies ist in Abbildung 3.2 auf Seite 27 dargestellt. Man sollte λ also so wählen, daß sich die Ansatzfunktionen noch deutlich voneinander unterscheiden, andererseits aber die Strukturen der Ansatzfunktionen groß sind und deutlich aus ihren Tomographiedaten hervorgehen. Innerhalb der hierdurch gegebenen Grenzen zeigt der Parameter λ keinen Einfluß auf die Rekonstruktionen.

3.5.3 Singulärfunktionen

Die Singulärfunktionen bei einer festen Wahl der Ansatzfunktionen zeigen jeweils deutliche Strukturen. Interessant ist nun, daß diese Strukturen relativ unabhängig

Abbildung 3.1: Singulärfunktionen bei einem Kollokationsverfahren

von der Wahl der Ansatzfunktionen sind, falls n >> m gilt, d.h. falls wesentlich mehr Ansatzfunktionen vorgegeben werden, als Meßlinien vorhanden sind. Abbildung 3.1 zeigt die Singulärfunktionen bei Wahl von φ_{k}(x) = e^{−4|x−x_{k}|} und einem 10 × 10 Gitter von Punkten x_{k}. Die Singulärfunktionen sind in üblicher Lese-Anordnung nach kleiner werdenden Singulärwerten sortiert.

Es zeigt sich deutlich, welche Strukturen leicht zu rekonstruieren sind, und wo es kritisch ist. Im wesentlichen sind es die hochfrequenten Anteile, die problematisch sind, wobei die Frequenzanteile bei den Singulärfunktionen einer Radontransformation wesentlich langsamer steigen als etwa bei den Ansatzfunktionen einer zweidimensionalen Kosinusreihe.

Mit Ausnahme der Ersten haben alle Singulärfunktionen genau oder annähernd den Mittelwert 0. Durch eine ”Verschiebung“ der Daten mittels

R̃^{+}g := f̄+ R^{+}(g − Rf̄),

wobei f̄ ein konstantes Bild darstellt, wird daher vor allem der Anteil dieser ersten Singulärfunktion gesteuert. Da bei der Schallpyrometrie glatte Bilder erwartet werden, wurde der Wert von f̄ als das Mittel der Meßdaten gewählt, so daß die erste Singulärfunktion nur schwach eingeht.

Abbildung 3.2: Erste Singulärfunktion bei verschiedenen Werten von λ

Abbildung 3.2 zeigt die erste Singulärfunktion für verschiedene Werte des Parameters λ. Der Parameter zeigt sich wie bereits oben gesagt, als wenig einflußreich, falls er nicht deutlich zu groß gewählt wird. In diesem Fall entarten die Singulärfunktionen zu Abbildern der Abtastgeometrie und führen zu wertlosen Rekonstruktionen.

Mit Abbildung 3.3 soll die ”Unabhängigkeit“ der Singulärfunktionen von den gewählten Ansatzfunktionen exemplarisch dargestellt werden. In der ersten Zeile ist jeweils eine Ansatzfunktion dreidimensional abgebildet, in der zweiten Zeile findet sich dieselbe Ansatzfunktion dann zweidimensional, und darunter sieht man die jeweils erste Singulärfunktion bei dieser Wahl von Ansatzfunktionen. Folgende Ansatzfunktionen sind dargestellt:

Erste Spalte: φ_{47}(x) = e^{−4.0|x−x_{47}|}, wobei die x_{k} ein 10 × 10–Gitter bilden,

Zweite Spalte: φ_{3}(x) = ∫_{L_{3}} e^{−4·|x−y|}dy (Kovarianzenansatz)

Abbildung 3.3: Abhängigkeit der ersten Singulärfunktion von den Ansatzfunktionen

Dritte Spalte: Pixelansatzfunktionen, Breite der Pixel: 0.4·(Seitenlänge von Ω), 10 × 10 Pixel.

Die Ähnlichkeit der Singulärfunktionen ist klar erkennbar. Bei den Ansatzfunktionen in der linken und mittleren Spalte sind alle 24 Singulärfunktionen ”optisch fast identisch“. Dies führt dazu, daß beide Verfahren praktisch identische Ergebnisse liefern. Dies wird auch noch einmal durch die Abbildung 3.4 unterstrichen. Sie zeigt in der ersten Zeile die Funktionen η_{1}, ..., η_{4}, die am Ende von Abschnitt 3.5.2 eingeführt wurden, und zwar für φ_{k}(x) = e^{−4|x−x_{k}|}. Die zweite Zeile zeigt die ersten vier Ansatzfunktionen beim Kovarianzenverfahren mit λ = 4.0 und in der untersten Zeile sieht man die Funktionen η_{1}, ..., η_{4}, die zum Kovarianzenansatz gehören.

Die deutlichen Ähnlichkeiten begründen die ähnlichen Resultate. Beim Kovarianzenansatz lohnt sich die Berechnung der η_{k} natürlich nicht, da hier ja bereits eine minimale Anzahl von Ansatzfunktionen vorgegeben wird.

3.6 Einbringen von a priori Informationen

Da die realen Temperaturverteilungen in einem Feuerraum leider im allgemeinen nicht in span(ψ_{1}, ..., ψ_{m}) liegen, müssen weitere Informationen in die Rekonstruktionen eingebracht werden.

Zwei Konzepte finden im folgenden Anwendung: Zum einen die Idee, dem Rechner ein ”Startbild“ zu geben, von dem man weiß, daß es der echten Situation im Abbildung 3.4: Vergleich der η_{k} mit den Ansatzfunktionen beim Kovar.-Verfahren

Kessel bereits recht nahe kommt, zum anderen die nachträgliche ”Glättung“ einer berechneten Temperaturverteilung im oberen Temperaturbereich.

Die Idee, Startbilder zu verwenden, wurde aus [14] übernommen. Man kann sie ohne weiteres bei jedem tomographischen Verfahren verwenden.

Als Abwandlung dieses Konzeptes wäre es möglich, als ersten Schritt eine Rekonstruktion mit einem Kollokationsverfahren durchzuführen, wobei man sich sehr wenige Ansatzfunktionen vorgibt, von denen man dafür aber sicher ist, daß sie sehr realitätsnah sind. So hätte man an dieser Stelle mehr Freiheiten, um auf die Daten zu reagieren, aber auch eine sehr klare Einbringung von a priori Informationen. Dieser Ansatz wurde aber nicht weiter verfolgt.

Die ”Glättung im oberen Temperaturbereich“ ist eine einfache Anpassung der Methode der Konvexen Projektionen an die spezielle Situation in der Schallpyrometrie^{1}. Da über die Konvexen Projektionen vor allem Informationen über die Fouriertransformation des Bildes eingebracht werden können, solche Informationen hier aber nicht bekannt sind, wurden die Konvexen Projektionen in dieser Arbeit nicht weiter verfolgt. Eine ausführliche Beschreibung ist in [15] und [20] zu finden.

Es sei noch darauf hingewiesen, daß die Rekonstruktionen mit den hier dargestellten Methoden unabhängig von den Daten beeinflußt werden können. Dies ist natürlich wünschenswert, da häufig in den Tomographiedaten nicht genügend Informationen vorhanden sind, andererseits kann die fahrlässige Anwendung dieser Methoden zu Rekonstruktionen führen, die künstlich eingebaute (also auch plausible) Details enthalten, obwohl diese Details im Original gar nicht vorhanden waren.

1Es wurde lediglich die Methode übernommen, Näherungen zu korrigieren, denn bei der beschriebenen Glättung liegt gar keine Projektion vor.

3.6.1 Startbilder

Im Falle der Schallpyrometrie weiß man bereits grob im vorhinein, wie die zu berechnenden Bilder aussehen müssen. Dieses Wissen entspringt einerseits anderen Meßmethoden, wie der Absaugepyrometrie, andererseits physikalischen und chemischen Überlegungen. Die Idee ist nun, dem Rechner eine solche Temperaturverteilung in Form eines Bildes vorzugeben, um ihn dann prüfen zu lassen, wie gut dieses Bild bereits den Daten entspricht, und wie stark es dann noch abgewandelt werden muß, damit das Ergebnis wirklich mit den Daten übereinstimmt. Ist h : Ω → IR eine für realistisch gehaltene Temperaturverteilung und R_{d}^{+} der verallgemeinerte inverse Operator zu R_{d}, so wird h mittels

R̃_{d}^{+}g = ζ(κ_{0}h) + R_{d}^{+} (g − ζR_{d}(κ_{0}h))

in die Rekonstruktion eingebaut. Formal muß R_{d} hierfür in kanonischer Weise auf L_{2} fortgesetzt werden.

Die reelle Zahl ζ ∈ (0, 1] dient als Parameter und steuert die Gewichtung, mit der das Startbild eingeht.

Der Faktor κ_{0} wird zur geeigneten Normierung von h eingebaut. Er wird so berechnet, daß ∥g−κ_{0}R_{d}h∥_{2}^{2} minimal ist, d.h. (κ_{0}·h) soll die Daten g bereits möglichst gut erfüllen. Zur Berechnung von κ_{0} sei nun R_{d}h = (h_{1}, ..., h_{m})^{t} ∈ IR^{m}. Dann soll κ_{0} die Funktion r(κ) = ∑_{i=1}^{m}(g_{i} − κh_{i})^{2} minimieren:

⇝ r′(κ_{0}) = \frac{dr}{dκ}(κ_{0}) = 0
⇔ \frac{d}{dκ} ∑_{i=1}^{m}(g_{i} − κh_{i})^{2}|_{κ=κ_{0}} = 0
⇔ ∑_{i=1}^{m}(−2)h_{i}(g_{i} − κ_{0}h_{i}) = 0
⇔ − ∑_{i=1}^{m} h_{i}g_{i} + κ_{0} ∑_{i=1}^{m} h_{i}^{2} = 0
⇔ κ_{0} = \frac{∑_{i=1}^{m} h_{i}g_{i}}{∑_{i=1}^{m} h_{i}^{2}} .

Wegen lim_{|κ|→∞} r(κ) = ∞ ist dies wirklich eine Minimumstelle.

3.6.2 Glättung im oberen Temperaturbereich

Aus Messungen mit Absaugepyrometern ist bekannt, daß die Temperaturverteilung im inneren (und damit im heißen) Bereich eines großen Feuerraumes relativ konstant ist. Die Kollokationsverfahren neigen jedoch dazu, das Gebiet mit hohen Temperaturen als schmalen Bereich mit einem relativ hohen Maximalwert zu ”gestalten“, da dies zu einer kleineren Norm des Bildes führt. Mit einer ähnlichen Korrektur einer gegebenen Näherung an die endgültige Rekonstruktion wie in 3.6.1 kann man dafür sorgen, daß der Bereich hoher Temperatur tendenziell eher als breites Gebiet mit einem kleineren Maximalwert dargestellt wird.

Es sei bereits eine Rekonstruktion h(x) zu den Meßdaten berechnet. Um diese zu glätten, wird sie nun einfach z.B. 50^{o}C unterhalb ihres Maximums abgeschnitten! Wenn also

T = max_{x} h(x) − 50

ist, dann wird das Bild h an allen Stellen, wo es den Wert T überschreitet, auf T gesetzt. Das so geänderte Bild bezeichne ich als h̄. Man erhält damit ein Bild, welches im oberen Temperaturbereich wunderbar glatt ist, aber natürlich stimmt dieses Bild nun nicht mehr mit den Meßdaten überein. Es muß daher wieder mit der Methode aus 3.6.1 korrigiert werden, um die Meßdaten zu erfüllen. Damit erhält man die Lösung zu R_{d}f = g als

f = h̄ + R_{d}^{+}(g − R_{d}h̄)

Es zeigt sich, daß dieses korrigierte Bild in der Tat einen breiteren Bereich hoher Temperatur aufweißt, wobei die Maximaltemperatur kleiner ist.

Dieses Verfahren kann man nun beliebig oft wiederholen, um den gewünschten Effekt zu verstärken. So wird es zu einem iterativen Verfahren. Es hat sich aber gezeigt, daß es nicht sinnvoll ist, die vorhandene Näherung häufiger als 2-3 Mal abzuschneiden, da dann quasi nicht-regularisierte Rekonstruktionen erzeugt werden.

Es sei noch angemerkt, daß mit dieser Methode natürlich viele Arten von a priori Informationen eingebracht werden können: Man rekonstruiere ein Bild, ändere dies entsprechend der Vorkenntnisse und korrigiere das geänderte Bild wieder so, daß es den Tomographiedaten entspricht.

Numerische Ergebnisse

4.1 Verfahrensvergleiche anhand von Phantombildern

Im Kapitel 3 wurden Kollokationsverfahren mit verschiedenen Ansatzfunktionen und das Einbringen von a priori Informationen über Startbilder und Glättung/Korrektur dargestellt.

Diese Verfahren werden nun hier anhand von verschiedenen Testbildern, sogenannten Phantomen, verglichen. Es wird dargestellt, wo in der Praxis die Auflösungsgrenze bei Rekonstruktionen mit der in Kapitel 2 beschriebenen Abtastgeometrie liegt.

Die Vergleiche der Verfahren werden nur mit exakten Daten und einer einzigen Abtastgeometrie (24 Meßpfade) durchgeführt, also ohne die Simulation von Meßfehlern. Das Verfahren, das sich hierbei am besten bewährt, wird dann auf seine Stabilität gegenüber Datenfehlern untersucht.

4.1.1 Die verwendeten Phantome

Unter ”Phantomen“ versteht man in der Tomographie künstlich erzeugte Bilder, die verwendet werden, um die Güte von Rekonstruktionsverfahren zu untersuchen. Hierfür werden aus den Phantomen die Tomographiedaten, also die jeweiligen Linienintegrale, berechnet, um dann daraus wieder Rekonstruktionen zu erstellen. So kann man Original und ”Fälschung“ vergleichen und sich ein Bild davon machen, was mit den verwendeten Verfahren wirklich zu erreichen ist.

Für diese Arbeit wurden sieben Phantome ausgewählt. Sie sind in Abbildung 4.1 dargestellt.

Die Bilder in der oberen Reihe dienen der Auflösungsuntersuchung, also der Frage: Wie große Details können rekonstruiert werden? Die übrigen Bilder sind einigermaßen realistische Temperaturverteilungen, wie sie bei der Schallpyrometrie auftauchen könnten. Oben rechts und unten ist die verwendete Farb- bzw. Graustufenskala dargestellt. Schwarz entspricht dem Wert 850, weiß dem Wert 1250. Dieser Wertebereich der Bilder wird zwar im folgenden immer einheitenlos betrachtet, er gestattet aber die Interpretation der Werte als Temperaturangaben in ^{o}C. Nun zu

4 NUMERISCHE ERGEBNISSE

Abbildung 4.1: Verwendete Phantome

den einzelnen Bildern:

Die Seitenlänge aller Bilder sei 1. In das erste Bild wurden Kreise mit den Durchmessern 0.2, 0.3, 0.4 und 0.5 gezeichnet. Der Hintergrund hat den Wert 900, alle Kreise liegen bei 1200. Das zweite Bild zeigt ebenfalls einen Kreis mit dem Durchmesser 0.5, allerdings an einer für die Abtastgeometrie günstigeren Position. Der Kreis im dritten Bild hat den Durchmesser 0.7. Ich nenne diese Bilder im folgenden A1, A2 und A3.

Die Bilder in der zweiten Reihe werde ich (von links nach rechts) als ”Blockprofil“ (kurz: Block), ”asymmetrisches Blockprofil“ (kurz: ABlock), ”komplexes Profil“ (kurz: Komplex) und ”asymmetrisches komplexes Profil“ (kurz: AKomplex) bezeichnen.

Das Blockprofil entspricht ziemlich genau den Wünschen an eine tatsächliche Temperaturverteilung im Kessel. In der Praxis wird diese Verteilung aber kaum auftreten, da von den 8 Brennern im Feuerraum eines tangential gefeuerten Braunkohlekessels immer mindestens einer aus Wartungsgründen außer Betrieb ist. In der Realität kann man also keine so symmetrischen Verteilungen erwarten. Das asymmetrische Blockprofil wurde durch eine Abwandlung des Blockprofils berechnet. Mit seiner Hilfe und dem Blockprofil kann man feststellen, wie die verschiedenen Verfahren auf symmetrische Daten reagieren. Das komplexe Profil ist ebenfalls realistisch. Bei der Installation der Meßgeräte im Kraftwerk Neurath befanden sich mehrere Meter unterhalb der Ebene mit den Meßgeräten Einblasstutzen für Ausbrandluft, also kältere Luft, die hier eingeblasen wird um die Verbrennungstemperatur zu senken. Da dieses Bild wieder viele Symmetrien aufweist, wurde auch eine asymmetrische Abwandlung davon aufgenommen: Das asymmetrische komplexe Profil. Es unterscheidet sich zum einen durch das Fehlen einer Stelle mit kalter Luft und zum anderen durch ein schwaches Abfallen der Temperatur im heißen Bereich von links nach rechts.

Bildwerte  Datenwerte
Name  Min Mittel Max Min arithm. Mittel gew. Mittel Max
A1  900 1027 1200 900  1038  1041  1200
A2  900 959 1200 900  950  958  1055
A3  900 1015 1200 900  998  1015  1111
Block  900 1114 1230 1027  1106  1120  1175
ABlock  885 1056 1184 961  1051  1060  1109
Komplex 862 1125 1200 1041  1121  1134  1193
AKomplex 862 1115 1199 1019  1111  1124  1189

Interessante Daten über die Bilder sind ihr Minimum, der Mittelwert und das Maximum. In der obigen Tabelle sind außerdem die wichtigsten Informationen über die Daten, also ebenfalls Minimum, Mittelwert und Maximum angegeben. Beim Mittelwert der Daten erscheint es sinnvoll, die Einzeltemperaturen der Pfade entsprechend der Länge der Pfade zu gewichteten. Es wurden daher sowohl das arithmetische Mittel als auch das gewichtete Mittel aufgelistet. Bei den stetigen Ansatzfunktionen (Block, ABlock, Komplex und AKomplex) zeigt sich, daß der Durchschnittswert des Bildes zwischen diesen beiden Mittelwerten liegt. Beide Angaben sind also sinnvoll.

Abbildung 4.2: Phantome mit 24-Linien Abtastgeometrie

Für die Verfahrensvergleiche wird die Abtastgeometrie benutzt, die im Kraftwerk Neurath tatsächlich verwendet wurde. Um eine Vorstellung davon zu geben, aus welchen Daten die Bilder rekonstruiert werden, zeigt die Abbildung 4.2 die Phantome mit der eingezeichneten Meßgeometrie.

Anhand dieser Darstellung kann man sich bereits überlegen, welche Details bestenfalls rekonstruiert werden können. Ein Beispiel: Die rechte untere Ecke bei A1 zeigt den Kreisrand und ein kleines Gebiet mit dem Hintergrundwert 900. Dies ist für die Rekonstruktionsalgorithmen unsichtbar, da keine Linie dieses Gebiet trifft! Jede Rekonstruktion ohne a priori Informationen muß daher in der ganzen Ecke einen konstanten Wert zeigen, und zwar ungefähr den Wert der kurzen Ecklinie. Ähnliches gilt für die Ecken des komplexen Profils: Die dunkle Stelle wird nur von einer Linie berührt, genauere Strukturen sind also auch hier nicht zu erwarten.

Die Daten selbst werden graphisch dargestellt, indem die Meßpfade in einer ihrem Meßwert entsprechenden Farbe gezeichnet werden (Abb. 4.3). Hierbei werden die Mittelwerte auf den Linien dargestellt, nicht die Linienintegrale.

Abbildung 4.3: Die zu den Phantomen gehörenden 24 Daten

4.1.2 Kollokation mit φ_{k} = e^{−λ|x−x_{k}|}: Auswirkung der Regularisierung

Bei dieser Wahl der Ansatzfunktionen gibt es drei Parameter, an denen man ”drehen“ kann:

• λ, also die Detailgröße in den Ansatzfunktionen,

• die Feinheit des Gitters der x_{k} und

• γ, also die Stärke der Regularisierung.

Bei den beiden ersten Parametern hat es sich gezeigt, daß ihre Veränderung wenig Auswirkungen hat. Wählt man das Gitter der x_{k} so fein, daß man mindestens vier mal soviele Ansatzfunktionen wie Meßlinien hat, dann werden die Strukturen der einzelnen Ansatzfunktionen in den Rekonstruktionen bereits nahezu unauffindbar. Eine weitere Verfeinerung des Gitters hat lediglich eine leichte Glättung des Bildes zur Folge. Die Variation von λ wurde für diese Arbeit untersucht, zeigte aber keine bemerkenswerten Ergebnisse. Auf Bilder zur Demonstration wurde daher verzichtet. Man sollte lediglich gemäß Kapitel 3 darauf achten, λ nicht so groß zu wählen, daß die Singulärfunktionen zu bloßen Abbildern der Abtastgeometrie werden. Dies wird vermieden, indem λ so groß gewählt wird, daß alle Strukturen in den Ansatzfunktionen aus den Tomographiedaten deutlich hervorgehen.

Bei der Regularisierung mit der Methode von Tikhonov–Phillips gibt es noch ein praktisches Problem: Wie gibt man an, wie stark gefiltert wurde? Der Parameter hierfür ist natürlich die Zahl γ. Für diesen gibt es aber keine gute Normierung. Wie schon in Kapitel 3 erläutert, findet man mit dem Filter Fγ^{TP}(σ) = \frac{σ^{2}}{σ^{2}+γ} die Lösung ξ des Minimierungsproblems

{min \atop ξ ∈ IR^{n}}(∥Aξ − g∥^{2} + γ∥ξ∥^{2}) .

Die Größenordnung von ξ hängt wegen dem Ausdruck ∥Aξ − g∥^{2} von der Matrix A, also den gewählten Ansatzfunktionen ab. Multipliziert man z.B. alle Ansatzfunktionen φ_{k} mit dem Faktor 2, dann müssen die Komponenten von ξ entsprechend kleiner sein, um ∥Aξ − g∥^{2} zu minimieren. Bei gleichem γ verliert der Term γ∥ξ∥ damit an Bedeutung. Um dieselbe Regularisierung zu erzielen, müsste γ also ver-

Abbildung 4.4: Darstellung der verwendeten Filter

größert werden! Da die Stärke der Regularisierung zudem auch noch von der Größe der Singulärwerte selbst abhängt, ist die Angabe eines einzelnen Zahlenwertes mangels einer passenden Normierung also nicht sinnvoll. In Abbildung 4.4 werden daher die Gewichtungsfaktoren Fγ^{TP}(σ_{k}), k = 1 … n graphisch gegen k dargestellt (durchgezogene Linien). Bei Wahl eines 20 × 20 Gitters und λ = 4.0 sind z.B. die Werte γ = 0.05 und γ = 0.3 eine sinnvolle Wahl.

Mit diesen beiden Werten von γ wurden bei einem 20 × 20 Gitter und λ = 4.0 Rekonstruktionen berechnet. Abbildung 4.5 zeigt die Rekonstruktionen mit γ = 0.3.

Rekonstruktionen  Daten
Name  Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  918  1042  1192  111  25
A2  866  959  1223  33  9
A3  855  1019  1263  79  17
Block  1050  1123  1227  35  16
ABlock  979  1062  1143  30  12
Komplex  1073  1137  1217  45  25
AKomplex  1044  1126  1223  73  20

Abbildung 4.5: Rekonstruktionen mit 24 Linien und γ = 0.3

Hier wurde relativ stark regularisiert, auf die Darstellung von kleinen Details wird also verzichtet, um eine größere Stabilität gegenüber Meßfehlern zu erreichen. Der kleine Kreis in A1 ist höchstens noch andeutungsweise vorhanden. Er liegt offenbar deutlich unter der Auflösungsgrenze. Dies wird auch bei einer schwächeren Regularisierung so bleiben, ohne a priori Informationen lassen sich Details in dieser Größe nicht rekonstruieren, selbst wenn sie an einer günstigen Stelle der Abtastgeometrie liegen. Auch der zweitkleinste Kreis von A1 ist nicht mit seinem korrekten Funktionswert rekonstruiert worden. Dies gilt erst ab dem dritten Kreis. Der Kreis mit dem Durchmesser 0.5 schließlich wird mit seinem Funktionswert richtig erkannt.

Da die Rekonstruktionen bei diesem Verfahren und diesen Ansatzfunktionen immer stetig sind, kann man nicht erwarten, unstetige Objekte in der korrekten Größe mit dem richtigen Funktionswert zu erkennen, da der Abfall auf ein anderes Niveau immer stetig, also ”glatt“ sein muß. Wie erwartet füllt die Rekonstruktion des größten Kreises die ganze untere rechte Ecke aus, da in den Tomographiedaten keine besseren Informationen über die Ecke vorhanden sind. A2 und A3 werden einigermaßen erkannt. Es ist aber nicht mehr zu erkennen, daß der Kreis in A3 im Innern einen konstanten Wert hat.

Die Originale zu den Bildern in der zweiten Reihe besitzen zumeist Details, die bei dieser Regularisierung weit unter der Auflösungsgrenze liegen. Die Rekonstruktionen scheinen daher sehr weit von ihren Phantomen entfernt zu sein. Vor allem das Verhalten am Rand wird sehr schlecht wiedergegeben, da dieser Bereich sehr schmal ist.

Neben den Daten über das Bild selbst sind noch die Größen der Rekonstruktionsfehler angegeben. Hierfür wurden die Tomographiedaten aus den rekonstruierten Bildern neu berechnet und mit den vorgegebenen Daten verglichen. ”max. Fehler“ gibt den Betrag der maximalen Abweichung an, ”ø–Fehler“ das arithmetische Mittel der Beträge der Fehler. Je stärker regularisiert wird, desto größer sind diese Fehler.

Rekonstruktionen  Daten
Name  Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  889  1038  1253  64  12
A2  859  958  1268  9  4
A3  825  1017  1245  36  8
Block  1018  1119  1239  16  6
ABlock  950  1060  1176  14  4
Komplex  1026  1133  1263  24  13
AKomplex  998  1122  1250  41  11

Abbildung 4.6: Rekonstruktionen mit 24 Linien und γ = 0.05

Abbildung 4.6 zeigt Rekonstruktionen mit γ = 0.05. Bei diesen Bildern wurde also nur schwach regularisiert. Die Tomographiedaten werden daher bereits recht gut erfüllt. Der kleine Kreis in A1 erscheint mit etwa dem halben Funktionswert (etwa 1050 vor einem Hintergrund von 900, wobei das Phantom dort den Wert 1200 besitzt). Im komplexen Profil erkennt man gut die negativen Auswirkungen der schwachen Regularisierung: Um die Daten gut zu erfüllen müssen Details eingebaut werden, die unterhalb der zuverlässigen Auflösung liegen. In diesem Fall sind das vier heiße Stellen, die für einen Maximalwert sorgen, der im Vergleich zum komplexen Profil um 63 zu groß ist^{1}. Zu Vergleichszwecken werden in Abbildung 4.7 auch die nicht regularisierten Rekonstruktionen gezeigt. Die noch vorhandenen Fehler in den Daten sind jetzt reine Diskretisierungs– und Rundungsfehler:

Rekonstruktionen  Daten
Name  Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  764  1032  1328  8  1
A2  856  956  1291  7  1
A3  804  1015  1250  8  1
Block  987  1116  1262  9  1
ABlock  931  1057  1202  8  1
Komplex  952  1128  1320  9  1
AKomplex  930  1118  1303  9  1

Abbildung 4.7: Rekonstruktionen mit 24 Linien ohne Regularisierung

Bei diesen Rekonstruktionen (Abbildung 4.7) wurde nun gar nicht ”glättend“ eingegriffen. Sie zeigen relativ kleine Details, doch diese sind zum großen Teil Artefakte. Die exakten Lösungen neigen dazu, ”Peaks“ zu besitzen, also kleine Gebiete mit viel zu hohen oder niedrigen Temperaturen. Hier kommen häufig Abweichungen von mehr als 100 vor.

1In der Bildverarbeitung ist dieser Effekt unter dem Namen ”Aliasing“ bekannt. Damit sind großflächige Auswirkungen von kleinen Details gemeint, die selbst nicht rekonstruierbar sind, jedoch durchaus die jeweiligen Daten beeinflussen. Hier sind das die ”kleinen“ dunklen Stellen am Rand, die bei den komplexen Profilen nicht rekonstruiert werden können.

Die großen Probleme beim Verzicht auf die Regularisierung werden aber erst sichtbar, wenn man fehlerbehaftete Daten verwendet.

4.1.3 Vergleich der Ergebnisse mit anderen Verfahren

Der Kovarianzenansatz

Rekonstruktionen  Daten
Name  Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  906  1041  1231  97  17
A2  858  959  1244  26  6
A3  844  1018  1247  68  12
Block  1040  1122  1221  29  9
ABlock  980  1062  1161  24  7
Komplex  1061  1136  1243  39  17
AKomplex  1019  1126  1238  64  14

Abbildung 4.8: Rekonstruktionen mit dem Kovarianzenverfahren

Beim Kovarianzenansatz wird mit den Ansatzfunktionen

φ_{k}(x) = ∫_{L_{k}} e^{−4.0|x−y|}dy

rekonstruiert. Für die Regularisierung bieten sich hier sowohl die Methode von Tikhonov–Phillips, als auch die Regularisierung nach Lavrent’ev an. In Abbildung 4.4 sind die entsprechenden Filter (vgl. Abschnitt 3.5.1) für geeignete Wahlen von γ dargestellt (gestrichelte Linien^{2}). Sie zeigen ein leicht unterschiedliches Verhalten: Die Regularisierung nach Lavrent’ev dämpft die Singulärfunktionen zu großen Singulärwerten stärker als Tikhonov–Phillips, die Singulärfunktionen zu kleinen Singulärwerten werden aber schwächer gedämpft. Damit liegt die Methode von Tikhonov–Phillips zwischen der Regularisierung nach Lavrent’ev und der abgeschnittenen Singulärwertzerlegung. Sie erscheint daher am günstigsten, da sie die einigermaßen sicheren Anteile bei der Rekonstruktion kaum dämpft, die bei Meßfehlern sehr unsicheren Anteile aber fast wegschneidet. Gleichzeitig hat sie nicht den Nachteil der abgeschnittenen Singulärwertzerlegung, die letzten noch mitgenommenen Singulärfunktionen überzubewerten.

2Die absoluten Werte von γ sind hierbei nicht vergleichbar!

Zur Berechnung von Rekonstruktionen wurde daher die Methode von Tikhonov-Phillips verwendet. Die Ergebnisse sind nun annähernd identisch zu den Ergebnissen bei der Wahl von φ_{k} = e^{−4|x−x_{k}|}. Es wird daher nur ein Beispiel (Abbildung 4.8) angegeben. Die oben gemachten Aussagen über die Regularisierung und die Auflösungsgrenze gelten bei diesem Verfahren unverändert. Der Parameter γ, der die Stärke der Regularisierung angibt, liegt hier aber in einer anderen Größenordnung. Gewählt wurde γ = 0.001 (vgl. Abb. 4.4). Die Regularisierung liegt also zwischen γ = 0.3 und γ = 0.05. Die Rekonstruktionen in Abbildung 4.8 sind vergleichbar mit denen in Abbildung 4.6 auf Seite 38.

Fourierverfahren

Auch hier wird nur ein Beispiel vorgestellt, diesmal allerdings, weil die Ergebnisse deutlich schlechter sind als in 4.1.2. In der Form, wie das Fourierverfahren bei der Schallpyrometrie eingesetzt wird, löst man hier nicht ein unterbestimmtes, sondern ein überbestimmtes Gleichungssystem. Es werden also weniger Ansatzfunktionen bereitgestellt, als Meßlinien zur Verfügung stehen. Auf eine Regularisierung wird verzichtet. Beschränkt man sich in x- und y-Richtung jeweils auf die 4 höchsten Frequenzen der Fourierreihenentwicklung, dann hat man 16 Ansatzfunktionen, also ein überbestimmtes Problem. Abbildung 4.9 zeigt die entsprechenden Rekonstruktionen.

Da relativ wenige Ansatzfunktionen vorgegeben werden, kann das Fourierverfahren auch nur sehr schlecht auf Details in den Bildern reagieren. Dies wird bei der Rekonstruktion von A1 in Abbildung 4.9 besonders deutlich. Die links in Abbildung 4.10 dargestellte Ansatzfunktion des Fourierverfahrens zeigt einige Strukturen von A1. Diese Ansatzfunktion ist daher in der Rekonstruktion zu A1 deutlich wiederzuerkennen. Für eine bessere Annäherung fehlt dem Rekonstruktionsverfahren der Spielraum. Es gelang auch nicht, die unteren beiden Kreise dort darzustellen, wo sie liegen sollten. Statt dessen wurde der ”heiße“ Bereich in die Ecken verschoben und zeigt ein Maximum, daß um knapp 200 zu hoch ist. A2 und A3 werden im Vergleich zu A1 gut erkannt. Der gute Wert des Maximums bei der Rekonstruktion von A2 muß als Zufall gewertet werden. Offenbar liegt der Kreis in A2 noch unterhalb der Auflösungsgrenze und wird daher relativ klein dargestellt^{3}, denn der deutlich

3Es ist ein Erfahrungswert aus der medizinischen Tomographie, daß auch Objekte mit etwa dem halben Durchmesser der Auflösungsgrenze in den Rekonstruktionen noch zu sehen sind, allerdings unterscheiden sie sich dann deutlich weniger von ihrer Umgebung, als beim Phantom.

Rekonstruktionen  Daten
Name  Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  879  1050  1396  105  34
A2  783  957  1218  23  9
A3  752  1015  1301  31  12
Block  1063  1124  1287  43  20
ABlock  950  1064  1177  39  16
Komplex  1074  1138  1269  55  30
AKomplex  1023  1126  1242  62  24

Abbildung 4.9: Rekonstruktionen mit dem Fourierverfahren

größere Kreis in A3 zeigt ein viel zu hohes Maximum. Selbst dieser Kreis ist also noch ein zu kleines Detail. Die ”Ecken“ des Blockprofils und die Bumerang-Struktur im asymmetrischen Blockprofil sind nicht darstellbar.

Abbildung 4.10: Eine hochfrequente Ansatzfunktion beim Fourierverfahren ART

Im engeren Sinne versteht man unter ART ein Kollokationsverfahren mit Pixelansatzfunktionen, bei dem das auftretende Gleichungssystem mit dem Verfahren von Kaczmarz gelöst wird. Da hier nur kleine Gleichungssysteme auftauchen, wurde das Gleichungssystem hier mit Hilfe der SVD gelöst. Ich gebe wiederum nur ein Beispiel an, da sich die Auflösung bei Pixelansatzfunktionen mit disjunkten Pixeln als zu grob erweist. Für Abbildung 4.11 wurden 5 × 5 Pixel und eine relativ schwache Regularisierung (Tikhonov Phillips) gewählt.

Rekonstruktionen  Daten
Name  Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  815  1035  1323  26  12
A2  837  953  1324  10  5
A3  821  1009  1268  21  7
Block  1012  1114  1246  13  5
ABlock  955  1056  1183  13  4
Komplex  1019  1126  1280  35  15
AKomplex  977  1115  1259  30  12

Abbildung 4.11: Rekonstruktionen mit Pixelansatzfunktionen (ähnlich ART)

Die Darstellung ist aber ganz offensichtlich ungeeignet, um Bilder wie die hier gewählten Phantome zu rekonstruieren. Als Verbesserung würde man die Anzahl der Pixel erhöhen wollen, und da irgendwann nicht mehr alle Pixel von mindestens einer Meßlinie getroffen werden, würde man als nächste Verbesserung zu Pixeln übergehen, die sich überschneiden. Spätestens an dieser Stelle ist klar, daß stetige Ansatzfunktionen bei der vorliegenden Problemstellung die natürlichere und bessere Wahl sind.

4.1.4 Einbringen von a priori Informationen

Auswirkung eines Startbildes

Rekonstruktionen  Daten
Name  Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  840  1034  1264  67  14
A2  800  952  1268  19  7
A3  789  1011  1247  38  11
Block  925  1114  1228  4  2
ABlock  865  1054  1179  5  1
Komplex  941  1128  1253  23  13
AKomplex  924  1117  1256  35  9

Abbildung 4.12: Auswirkung des Blockprofils als Startbild

Als Startbild wurde in Abbildung 4.12 eine leicht geglättete Abwandlung des Blockprofils gewählt (vgl. den Abschnitt ”Auswahl von Startbildern“). Es ist oben rechts anstelle der Farbskala dargestellt, die natürlich nicht verändert wurde. Dieses Bild ist günstig, da es zumindest am Rand den meisten Phantomen ähnelt. Die Rekonstruktion erfolgte mit φ_{k} = e^{−4|x−x_{k}|} wie in 4.1.2 bei γ = 0.1. Da es explizit vorgegeben wurde, stimmt das Aussehen der Rekonstruktionen am Rand nun sehr gut mit den Phantomen überein. Das Blockprofil selbst wird natürlich perfekt ”rekonstruiert“, da das Startbild bereits mit den Daten übereinstimmte. Die Ecke unten rechts in A1 ist ein weiteres Beispiel für die Auswirkung des Startbildes: Das Abfallen des Bildwertes zur Ecke hin ist mit den Tomographiedaten allein nicht zu erreichen.

Natürlich kann man mit einem ungünstigen Startbild extrem ”schlechte“ Rekonstruktionen erzeugen. Dies wird im Abschnitt ”Auswahl von Startbildern“ gezeigt.

Die weißen Stellen z.B. im komplexen Profil sind nach wie vor Artefakte. Auch hier ist aus den Tomographiedaten nicht zu ersehen, daß das komplexe Profil im

Rekonstruktionen  Daten
Name  Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  800  1033  1251  46  9
A2  798  953  1250  8  3
A3  783  1011  1232  24  5
Block  923  1113  1212  5  3
ABlock  864  1054  1150  5  2
Komplex  933  1126  1241  18  11
AKomplex  920  1116  1236  27  7

Abbildung 4.13: Rekonstruktionen mit Startbild und Glättung

oberen Temperaturbereich flach ist. Dies muß nun ebenfalls als a priori Information eingebracht werden.

Glättung im oberen Temperaturbereich

Die Bilder in Abbildung 4.13 wurden zunächst ebenso berechnet, wie im Abschnitt ”Startbilder“ (Abb. 4.12). Anschließend wurde jedes Pixel, das einen Wert von mehr als

T = Maximaltemperatur − 50

hat, auf T gesetzt. Zu diesem Bild wurden dann wie in 3.6.2 beschrieben die Fehlerdaten berechnet und das daraus rekonstruierte Bild hinzuaddiert. Regularisierung: γ = 0.1. Da jedoch die Fehlerdaten sehr viel kleinere Werte haben als die Ausgangsdaten, ist der hier übrigbleibende Fehler bei derselben Regularisierung relativ klein. Insgesamt wurde also nur relativ schwach regularisiert.

Diese Glättung hat u.a. den Einfluß des Startbildes gemindert: Bei A1 füllt der große Kreis nun die Ecke wieder fast aus.

Man könnte dieses Verfahren häufiger iterieren, würde dann aber sehr schnell quasi fehlerfreie — also nicht regularisierte — Rekonstruktionen erhalten. Da Meßdaten aber immer fehlerbehaftet sind, würde dies zu deutlichen Artefakten führen.

Bei der Untersuchung von fehlerbehafteten Datensätzen werden sich diese Bilder als eher zu wenig regularisiert herausstellen. Daher wurden die Bilder im folgenden noch einmal mit einer stärkeren Regularisierung rekonstruiert. Hierfür wurde das Startbild mit der Gewichtung 0.5 eingebracht (also mit dem Faktor 0.5 multipliziert, um ein geringeres Gewicht zu haben) und es wurde γ = 0.3 gewählt.

Rekonstruktionen  Daten
Name  Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  876  1037  1208  85  17
A2  837  955  1222  15  5
A3  816  1014  1234  43  11
Block  986  1117  1208  14  7
ABlock  918  1057  1136  13  5
Komplex  1005  1131  1213  28  18
AKomplex  986  1120  1213  51  13

Abbildung 4.14: Rekonstruktion mit Startbild, Glättung und stärkerer Regularisierung

Auswahl von Startbildern

Die Wahl eines vernünftigen Startbildes ist ein mindestens ebenso gewichtiges Problem wie die Wahl eines annehmbaren Regularisierungsparameters. Ein Startbild soll

1. im Rahmen der a priori Informationen realistisch sein, 2. kleine Details einbringen, die mit Sicherheit vorhanden, aber nicht aus den

Tomographiedaten ersichtlich sind,

3. bei unsicheren Details korrigierbar sein, d.h. Details, die nach den Tomographiedaten nicht vorhanden sein dürfen, müssen wieder ausgeglichen werden können.

Abbildung 4.15: Rekonstruktionen der in der obersten Reihe dargestellten Phantome mit Hilfe der Startbilder in der linken Spalte

Leider stellt sich heraus, daß Details in der Größe der dunklen Stellen im komplexen Profil durch Startbilder zu fast 100% beeinflußt werden können. Um dies zu demonstrieren, werden das komplexe Profil und glattere Abwandlungen davon als Startbild verwendet, um das Blockprofil, das komplexe Profil selbst und das asymmetrische komplexe Profil zu rekonstruieren. Die leider recht eindrucksvollen Resultate sind in Abbildung 4.15 zu sehen.

Wenn man das komplexe Profil verwendet, um das Blockprofil zu rekonstruieren, dann bleiben alle kalten Stellen gegenüber dem Startbild so gut wie unverändert, lediglich im heißen Bereich gibt es leichte Veränderungen. Für das Rekonstruktionsverfahren sehen also das Blockprofil und das komplexe Profil fast gleich aus. Da aufgrund der Regularisierung keine exakte Lösung des Problems verlangt wird, sondern eine möglichst glatte, werden die Unterschiede zum großen Teil ignoriert oder großflächiger ausgeglichen. Auf keinen Fall erfüllt das komplexe Profil als Startbild die dritte der geforderten Bedingungen.

Bei der Rekonstruktion des asymmetrischen Blockprofils interessiert die linke Wand: Was passiert bei der im Phantom nicht vorhandenen kalten Stelle? Das Rekonstruktionsverfahren muß hier ein vorhandenes Loch ”auffüllen“. Da dieser dunkle Bereich aber unter der Auflösungsgrenze liegt, wird in einem größeren Bereich um das ”Loch“ herum ”aufgefüllt“. Dies führt zu einem hohen, im Bild weißen Ring um die auszugleichende Stelle herum.

Die beiden anderen Startbilder wurden durch unterschiedlich starke Glättung des komplexen Profils gewonnen. Beim dritten Startbild kann bei der Rekonstruktion des asymmetrischen komplexen Profils der künstliche dunkle Bereich an der linken Wand fast wieder ausgeglichen werden. Auch das Blockprofil wird wieder einigermaßen erkannt. Der Unterschied zur Rekonstruktion des komplexen Profils ist aber nach wie vor gering.

Startbilder sollten also keine Details besitzen, die kleiner sind als die Details des zweiten Startbildes (mittlere Glättung). Diese Details sind dann aber auf jeden Fall in der Rekonstruktion wiederzufinden. Ihr Einbau sorgt aber dafür, daß diese Details nicht fehlinterpretiert werden, da ohne dieses Startbild gröbere Details rekonstruiert würden, die dann weiter von der realen Situation entfernt wären.

Unsichere Details müssen beim Startbild noch stärker ”weggeglättet“ werden, damit sie vom Programm wieder entfernt werden können, falls sie den Tomographiedaten entsprechen.

Diese Rekonstruktionen mit verschiedenen Startbildern sind recht ernüchternd, was die Aussagekraft der Tomographiedaten betrifft. In Abbildung 4.15 fällt auf, daß die Rekonstruktionen, die jeweils in einer Zeile stehen, kaum unterschiedlich sind. Dagegen unterscheiden sich die Rekonstruktionen in je einer Spalte sehr stark und ähneln jeweils eher dem jeweiligen Startbild als dem Phantom.

Man würde sich dieses Verhalten gerade andersherum wünschen, doch leider hat das gewählte Startbild bei diesen Detailgrößen offenbar einen größeren Einfluß auf die Rekonstruktionen, als die gemessenen Daten!

4.2 Rekonstruktionen aus fehlerbehafteten Daten

Im folgenden wurden alle Rekonstruktionen mit φ_{k}(x) = e^{−4|x−x_{k}|} und einem 20×20 Gitter von Punkten x_{k} erstellt. Das Problem in diesem Abschnitt ist die Flut von Bildern, die dargestellt werden muß. Es werden zwei der sieben Phantome mit verschieden großen Fehlern, zwei verschiedenen Regularisierungen und verschiedenen Startbildern rekonstruiert. Da zu den Daten ein zufälliger Fehler addiert wird, wurde dies jeweils vierfach wiederholt, um einen Eindruck davon zu geben, was alles passieren kann.

Der Fehler bei Meßdaten der Schallpyrometrie setzt sich folgendermaßen zusammen: (angegeben werden jeweils die geschätzten Standardabweichungen)

• reine Meßfehler: 1 %

• Fehler aufgrund von stochastischen hochfrequenten Temperaturschwankungen im Feuerraum: 2.3 %

• Fehler aufgrund von niederfrequenten Temperaturschwankungen: 3 %

Der dritte Punkt ist vor allem für Bilder wichtig, die aus 10–Minuten-Mitteln gebildet wurden. Bei der Aufbereitung der Daten für einen Film kann man davon ausgehen, daß langsame Schwankungen einigermaßen korrekt erfaßt werden. Addition der Fehler nach der Gaußschen Fehlerquadratmethode liefert für die beiden ersten Fehler eine Gesamtstandardabweichung von etwa 2.5 %, für alle drei Fehler knapp 4 %.

Zu den exakten Daten wurden daher normalverteilte Fehler mit einer Standardabweichung von 1 %, 2.5 % und 4 % addiert.

4.2.1 Die 24-Liniengeometrie aus Neurath

Der dargestellte Wertebereich bei den Abbildungen 4.16 und 4.17 liegt nach wie vor zwischen 850 und 1250, entspricht also Temperaturangaben in ^{o}C. Bei jedem Datensatz wurden zwei verschiedene Regularisierungen verwendet. Für die Bilder im linken Drittel wurde zu den Daten jeweils ein normalverteilter Fehler mit einer Standardabweichung in einer Größenordnung von einem Prozent addiert. Dies geschah viermal, um einen Eindruck davon zu vermitteln, ”was alles so passieren kann“. Der genaue mittlere Fehler schwankt natürlich um den Wert von einem Prozent, da 24 Zahlen eine sehr kleine Stichprobe sind. In der ganz linken Spalte wurde relativ exakt rekonstruiert (entsprechend Abb. 4.13), also schwach regularisiert und in der Spalte rechts daneben wurde stärker regularisiert (Abb. 4.14). Dasselbe wurde mit Standardabweichungen in der Größe von 2.5 % und 4 % wiederholt.

Nun aber zu den Bildern in Abbildung 4.16: Beginnen wir mit einem groben Überblick und gehen dann zu den Details über. Offenbar wirkt sich der Fehler von einem Prozent noch fast gar nicht aus, er wird einfach ”weggeglättet“. Fehler bis 2.5 % sind bereits sichtbar, verdecken aber noch nicht die wesentlichen Strukturen. Bei einem Fehler mit einer Standardabweichung von rund 4 % dagegen kommen sehr große Veränderungen vor.

Der Einfluß der Regularisierung wird an diesen Bildern gut deutlich: Auf die Darstellung kleiner Details wird verzichtet, da man diese bei fehlerhaften Daten ohnehin kaum noch von Artefakten unterscheiden kann, die aufgrund der Meßfehler

Abbildung 4.16: Rekonstruktion von A1 aus 24 fehlerbehafteten Daten

eingebaut werden. Dies läßt sich z.B. am kleinen Kreis in A1 beobachten. Die Regularisierung führt zu glatteren Bildern die weniger die Tendenz haben, ”Peaks“ zu besitzen. Somit liegt der Maximalwert der stärker regularisierten Rekonstruktionen generell niedriger als bei den Rekonstruktionen, die die Daten genauer erfüllen. Bei einem Datenfehler von einem Prozent liegen die Maximalwerte der Bilder in der linken Spalte von Abb. 4.16 zwischen 1255 und 1265, rechts davon liegen sie zwischen 1198 und 1221.

Wenn größere Meßfehler erwartet werden müssen, dann werden die Ergebnisse natürlich schlechter. Ein normalverteilter Fehler von 4 % stellt sich aber als vernichtend heraus. Die Bilder in der zweiten und dritten Zeile von Abb. 4.16 zeigen, daß Details in der Größe des kleinen Kreises von A1 zu reinen Ergebnissen des Zufalls werden: In der zweiten Zeile fehlt der Kreis ganz, in der dritten ist er ebenso groß, wie der Kreis oben rechts! Der Bereich des Maximums liegt nun bei schwacher Regularisierung (links) zwischen 1279 und 1353, bei stärkerer Regularisierung immer noch bei 1232 bis 1284.

Die Abbildung 4.17 zeigt die entsprechenden Rekonstruktion aus fehlerbehafteten Daten des asymmetrischen komplexen Profils. Als Startbild wurde das mittelstark geglättete komplexe Profil verwendet. Hier ist wiederum interessant, ob die falsche Vorgabe des Startbildes an der linken Seite kompensiert werden konnte. Es zeigt sich jedoch, daß die falsche Vorgabe bereits bei einem Fehler von 2.5 % von den Artefakten aufgrund von Meßfehlern deutlich übertroffen werden.

Abbildung 4.17: Rekonstruktion des Phantoms AKomplex aus 24 fehlerbehafteten Daten. Startbild: komplexes Profil.

Die Tendenz ist auch hier folgende: Ein Datenfehler von einem Prozent führt zu durchaus annehmbaren Ergebnissen, in denen die groben Details noch mit dem jeweiligen Phantom übereinstimmen. Ein Fehler von 2.5 % kann wohl als die Grenze des Annehmbaren angesehen werden. Hier zeigen sich noch die groben Strukturen es Bildes, man darf aber nicht der Versuchung verfallen, Einzelheiten interpretieren zu wollen. Bei einem Fehler von 4 % ist das Ergebnis fast völlig vom Zufall bestimmt. Der Betrachter hat nun keine Chance mehr, aufgrund der Rekonstruktion zu erraten, wie das Phantom wohl aussah.

4.2.2 Rekonstruktionen mit 4 Meßgeräten pro Wand

Abbildung 4.18: Die zu den Phantomen gehörenden 96 Daten

Anhand einer weiteren Abtastgeometrie soll gezeigt werden, daß bei einer Erhöhung der Datenanzahl mit diesen Verfahren durchaus gute Rekonstruktionen auch der komplexen Profile zu erreichen sind. Analog zur Geometrie aus Neurath werden vier Meßgeräte pro Wand verwendet. Damit ergeben sich pro Datensatz 96 Meßwerte und damit deutlich bessere Auflösungsmöglichkeiten.

Die Daten zu den Phantomen sind in Abbildung 4.18 dargestellt. Grobe Details lassen sich nun bereits mit dem bloßen Auge in den Daten wiederfinden. Die dunklen Stellen bei den komplexen Profilen fallen hier aber noch nicht auf.

Abbildung 4.19: Die verwendeten Filter bei 96 Daten

Rekonstruiert wird mit einem 20 × 20 Gitter von Punkten x_{k} und φ_{k} = e^{−4|x−x_{k}|}. Es werden eine starke Regularisierung (γ = 0.3, Abb. 4.20) und eine recht schwache (γ = 0.005, Abb. 4.21) verwendet. Die Filter Fγ^{TP}(σ_{k}) sind in Abbildung 4.19 dargestellt.

Rekonstruktionen  Daten
Name  Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  813  1027  1266  187  18
A2  827  959  1246  65  10
A3  801  1015  1230  70  10
Block  917  1114  1208  9  2
ABlock  873  1056  1151  24  2
Komplex  959  1125  1227  43  15
AKomplex  941  1115  1207  60  13

Abbildung 4.20: Rekonstruktionen mit 96 Daten und γ = 0.3

Problematisch sind nach wie vor die Maximalwerte der Bilder, vor allem bei den Rekonstruktionen von A1. Jedoch liegt der kleine Kreis jetzt offenbar dicht an der Auflösungsgrenze und wird sicher identifiziert. Dies wird sich auch bei den Rekonstruktionen aus fehlerbehafteten Daten zeigen.

Man kann hier wiederum viele Beobachtungen machen, die bereits bei weniger Daten auffielen. So wird z.B. das asymmetrische Blockprofil auch hier durch die Glättung gemäß 3.6.2 eher verfälscht als verbessert. Demgegenüber profitieren die komplexen Profile aber deutlich von diesem Verfahren. Selbst bei fast 100 Daten zeigt das Einbringen von a priori Informationen also noch merkliche Auswirkungen. Man muß daraus den Schluß ziehen, daß z.B. die Information über die Maximaltemperatur in den Tomographiedaten einfach nicht vorhanden ist. Diese läßt sich bei der Annahme stetiger Originale allenfalls bis auf eine Differenz von ±30 bestimmen.

Bei den Rekonstruktionen von unstetigen Phantomen fällt das Gibbs–Phänomen auf. Besonders bei der Rekonstruktion von A3 sieht man einen schwarzen – also tiefen – Rand um den weißen Kreis und der Rand des Kreises liegt deutlich zu hoch.

Rekonstruktionen  Daten
Name  Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  805  1026  1255  96  6
A2  842  959  1219  25  3
A3  854  1015  1230  37  4
Block  910  1114  1218  2  0
ABlock  892  1056  1166  5  1
Komplex  923  1127  1224  52  7
AKomplex  924  1117  1219  49  6

Abbildung 4.21: Rekonstruktionen mit 96 Daten und γ = 0.005

Sind die Datenfehler klein, so kann man durchaus mit γ = 0.005 regularisieren, wie dies in Abbildung 4.22 geschehen ist. Trotz der Addition eines Fehlers mit einer Standardabweichung von einem Prozent werden noch alle kalten Stellen von AKomplex wiedergefunden. Im heißen Teil, der eigentlich sehr glatt verlaufen sollte, zeigen sich schwache hochfrequente Artefakte. Man erkennt diese als unzuverlässige Details, da sie sich kaum vom ”Hintergrund“ unterscheiden. Die kalten Stellen am Rand sind zwar von ihrer Ausdehnung her nicht größer, jedoch aufgrund ihres Funktionswertes deutlich hervorragender und damit auch sicherer.

Es sei noch darauf verwiesen, daß noch kleinere Details in den Phantomen möglicherweise zu schlechteren Rekonstruktionen geführt hätten. Die hier verwendeten Phantome zeigen keine Details mehr, die nicht erkannt werden können, und sind somit vielleicht für die 96-Linien Geometrie etwas optimistisch.

Rekonstruktionen  Daten
Name  Spalte Minimum Mittelwert Maximum max. Fehler ø–Fehler
A1  1  757  1028  1261  63  6
A1  2  797  1026  1264  97  8
A1  3  783  1025  1290  114  8
A1  4  735  1025  1280  67  7
AKomplex  1  912  1115  1227  59  7
AKomplex  2  855  1115  1230  73  7
AKomplex  3  824  1116  1241  91  8
AKomplex  4  882  1117  1251  111  8

Abbildung 4.22: Rekonstruktionen mit 96 Daten und γ = 0.005. Datenfehler: 1%

4.3 Anwendung auf Meßdaten

In diesem Abschnitt werden Rekonstruktionen zu Meßdaten dargestellt, die bei Messungen im Kraftwerk Neurath gewonnen wurden. Es ist nicht möglich, im Rahmen dieser Arbeit eine detailierte Auswertung dieser Meßergebnisse anzufertigen. Hier soll es lediglich darum gehen, dem Leser einen Eindruck von den bei realen Messungen auftretenden Bildern zu vermitteln.

Die Daten wurden nicht aus den in Kapitel 2 erwähnten Filmen entnommen, sondern aus der zu diesem Thema bereits existierenden Literatur. Qualitativ ähneln diese Rekonstruktionen den Filmen natürlich.

Es wurde aber dieselbe Darstellung wie in den Filmen verwendet, um auch hier die wichtigsten Daten auf einen Blick zur Verfügung zu haben.

Die ”Glattheit“ der Rekonstruktionen deuten darauf hin, daß der normalverteilte Fehler bei den Daten kaum über einem Prozent liegt. Es genügt daher, relativ schwach zu regularisieren, und zwar entsprechend 4.1.4 (Startbild: Blockprofil, einmalige Glättung, γ = 0.1).

Um den Speicherbedarf der Filme möglichst niedrig zu halten, wurde die Beschriftung auf ein Minimum reduziert. Die Angaben haben folgende Bedeutung:

D-Fehler: (Durchschnittsfehler) Aus den Rekonstruktionen wurden wiederum die Tomographiedaten berechnet. Als Durchschnittsfehler ist der Mittelwert der Abweichung von den Tomographiedaten der Rekonstruktion und den gegebenen Meßdaten (in Grad Celsius) angegeben. Ebenso wie der maximale Fehler läßt dies einen Rückschluß auf die vorgenommene Regularisierung zu.

max. Fehler: Der Betrag der maximalen Abweichung der Tomographiedaten der Rekonstruktion von den Meßdaten.

Daten: Minimum, gewichtetes Mittel und Maximum der Daten. Bei der Mittelbildung wurden die Daten entsprechend der Länge der jeweiligen Pfade gewichtet.

Temp: Minimum, Mittelwert und Maximalwert des Bildes.

unten links: Angabe der Uhrzeit. (entfällt hier)

unten, Mitte: Angabe des Datums. (entfällt hier ebenfalls, der Platz dient der Quellenangabe)

unten rechts: Angabe des Meßortes. (auch dies entfällt hier)

Die Daten zu den folgenden Rekonstruktionen sind aus [4] entnommen. Die jeweilige Nummer der Abbildung ist im unteren Teil der Bilder angegeben. In [4] wird der Zusammenhang zwischen verschiedenen Kohletypen, verschiedenen Kesselfahrweisen, der Temperaturverteilung sowie der Schadstoffproduktion (NO_{X}) untersucht. Die Bilder haben einen etwas anderen Wertebereich als die Phantome, daher wurde die Farb- bzw. Graustufenskala jetzt auf der Bereich 950 ^{o}C bis 1350 ^{o}C verschoben.

4 NUMERISCHE ERGEBNISSE

Abbildung 4.23: Rekonstruktionen mit Daten aus [4]

4.4 Diskussion

Es wurden verschiedene Rekonstruktionsverfahren anhand von exakten Tomographiedaten verglichen. Hierbei stellten sich eindeutig die Verfahren als besser heraus, die mit einer großen Zahl von glatten Ansatzfunktionen arbeiten, da die Singulärfunktionen bei diesen Ansatzfunktionen gut an die Abtastgeometrie angepaßt sind.

Wenn Ansatzfunktionen verwendet werden, die selbst gar nicht an die Geometrie angepaßt sind, dann müssen relativ viele Funktionen verwendet werden (größenordnungsmäßig die vierfache Anzahl der Meßlinien). Verwendet man angepaßte Funktionen, wie beim Kovarianzansatz, so genügen bereits soviele Ansatzfunktionen, wie Meßlinien vorhanden sind.

Weiterhin sollten alle Ansatzfunktionen lediglich Details besitzen, die oberhalb der Auflösungsgrenze mit exakten Daten liegen. Dies ist nicht überraschend, denn wenn die Funktionen Details enthalten, die für die Rekonstruktionsalgorithmen unsichtbar sind, dann werden diese Details in die rekonstruierten Bilder eingebaut, ohne daß dieses Problem im Programm auffiele. Dies geschieht z.B. beim Fourierverfahren.

Die Auflösungsgrenze in Zahlen anzugeben, ist sehr schwierig. Bei Verwendung exakter Daten und der Geometrie aus Neurath liegt sie ohne Verwendung von a priori Informationen beim 20m × 20m Querschnitt des Kessels bei etwa 8 Metern, d.h. Details, die diese Größe haben, werden mit ihrem Funktionswert wiedergegeben. Dieser ist aber keineswegs konstant innerhalb des Details.

Interessante Größen sind das Minimum, der Mittelwert und das Maximum der jeweiligen Temperaturverteilung. Das Minimum läßt sich aus den Meßdaten der Schallpyrometrie nicht wiedergeben, da die Minimaltemperatur in einem sehr schmalen Bereich am Kesselrand angenommen wird. Hier kommt es zu Abweichungen von mehreren hundert Grad. Viel interessanter wäre auch die Maximaltemperatur, doch auch hier kommt es noch zu Abweichungen, die auch bei vielen Meßdaten und dem Einbringen von a priori Informationen bei 30 ^{o}C liegen. Hierbei wird benutzt, daß die Temperatur im heißen Teil des Kessels einigermaßen konstant ist. Wenn dies in Wirklichkeit aus irgendwelchen Gründen nicht der Fall sein sollte, dann kann die Abweichung wesentlich höher werden.

Bei der Rekonstruktion der Phantome wurde der Mittelwert der Bilder immer bis auf einen Fehler von maximal 10 ^{o}C erkannt. Da dies auch die Größenordnung der reinen Meßfehler ist, ist der Mittelwert also hinreichend gut rekonstruierbar. Um diesen herauszufinden genügt aber auch bereits die einfache Mittelung der Eingangsdaten.

Durch das Einbringen von a priori Informationen kann die Auflösungsgrenze verbessert werden. Hier ist aber Vorsicht geboten, da falsche a priori Informationen häufig nicht aufgrund der Tomographiedaten ausgeglichen werden können.

Bei der Rekonstruktion aus fehlerbehafteten Daten zeigte es sich deutlich, daß die Schlechtgestelltheit der Problemstellung auch bei verbesserten Algorithmen ihre Auswirkungen behält. Ein zufälliger Fehler in den Daten mit einer Standardabweichung von 4 % zerstört jede Möglichkeit, halbwegs detailreiche Rekonstruktionen zu gewinnen. Eine Standardabweichung von 2.5 % führt nur bei sehr vielen Meßdaten zu noch annehmbaren Verfälschungen, bei einem Fehler von rund 1 % gehen die Meßfehler aber in den Regularisierungsfehlern unter, die ja absichtlich eingebaut werden.

Speziell auf die Problemstellung der Schallpyrometrie bezogen zeigt sich, daß 24 Messungen nicht ausreichen, um bei Rekonstruktionen der Temperaturverteilung aus diesen Daten die nötige Auflösung zu erreichen. Wenn es nicht nur um die Feststellung von Schieflagen, sondern z.B. um die Auswirkung von Ausbrandluft geht, dann sind sicherlich drei oder gar vier Meßgeräte pro Wand notwendig. Leider haben kleine Details wie die kalten Stellen bei den komplexen Profilen durchaus Auswirkungen auf die Daten, auch wenn sie bei zu wenigen Daten nicht aufgelöst werden können. Dies führt zu Mißinterpretationen der Daten und etwa bei den komplexen Profilen zu überhöhten Werten in der Mitte der Bilder, die nur durch eine relativ starke Regularisierung ausgeglichen werden können. Dies zeigt, daß die Abtastgeometrie unter Umständen an die erwartete Detailfülle der realen Situation angepaßt werden muß, um realistische Bilder zu erhalten.

Zusammenfassung

Die Schallpyrometrie ist ein Verfahren zur Temperaturmessung, das als Daten die Durchschnittstemperaturen auf gewissen Linien innerhalb der Meßebene liefert. Aus diesen Werten soll die ursprüngliche Temperaturverteilung zurückgerechnet werden. Bei der RWE Energie AG wurden hierfür bereits Programme entwickelt bzw. angewandt, welche auf ART bzw. dem Fourierverfahren basieren.

In dieser Arbeit wurden Verfahren vorgestellt, die man als direkte Weiterentwicklung dieser Verfahren ansehen kann, obwohl der Werdegang ein anderer war. Die Verbesserungen wurden durch

• die Verwendung einer großen Zahl von detailarmen, stetigen Ansatzfunktionen, und

• die Benutzung der Methode von Tikhonov–Phillips zur Lösung des Gleichungssystems

erreicht. Hiermit erhält man Kollokationsverfahren, die einigermaßen stabil gegen Fehler in den Daten sind und weniger als ART und das Fourierverfahren dazu neigen, Artefakte zu produzieren.

In dieser Arbeit wurden Kollokationsverfahren allgemeiner untersucht. In Kapitel 3 wurde hierzu die Singulärwertzerlegung der diskreten Radontransformation beschrieben. Man kann die verwendeten Lösungsmethoden damit als direkte Anwendung von bekannten Lösungsformeln mit Hilfe der Singulärwertzerlegung eines linearen Operators ansehen. Weiterhin wurde eine Lösung des Problems von Mehrfachmessungen bei Verwendung der (Moore-Penrose) verallgemeinerten Lösung dargestellt. Beides läßt sich mit elementaren Methoden der Linearen Algebra behandeln, dies ist jedoch in der mir bekannten Literatur noch nicht geschehen.

Schließlich wurden in Kapitel 3 verschiedene Wahlen von Ansatzfunktionen vorgestellt und bewertet.

Da realistische Temperaturverteilungen im allgemeinen nicht in dem 24-dimen-sionalen Unterraum von L_{2} liegen, aus dem die Rekonstruktionen jeweils ausgewählt werden, mußten auch Möglichkeiten zum Einbringen von a priori Informationen eingeführt werden. Im Text werden sie als Verwendung von Startbildern und Glättung

5 ZUSAMMENFASSUNG

im oberen Temperaturbereich bezeichnet. Diese zeigen deutliche Einflüsse auf die Rekonstruktionen und geben dem Anwender damit die Möglichkeit, die erzeugten Bilder unabhängig von den Tomographiedaten zu regulieren.

In Kapitel 4 wurden die dargestellten Verfahren anhand von Phantombildern verglichen. Als empfehlenswerte Ansatzfunktionen haben sich φ_{k}(x) = e^{−λ|x−x_{k}|} und φ_{k}(x) = ∫_{L_{k}} e^{−λ|x−y|}dy herausgestellt. Die auflösbare Detailgröße stellte sich bei nur 24 Meßdaten als denkbar gering heraus und ist schlechter, als die erwartete Detailgröße der realen Temperaturverteilung bei der Installation der Meßgeräte im Kraftwerk Neurath, wo unter anderem der Einfluß von Ausbrandluft festgestellt werden sollte. Um die Rekonstruktionen zu verbessern, kann man die erwähnten Methoden zur Einbringung von a priori Informationen verwenden. Hierbei müssten jedoch Details in der Größenordnung der Abkühlbereiche durch Ausbrandluft durch entsprechende Startbilder zu 100 % vorgegeben werden, da solche Objekte in den Bildern zu Fehlinterpretationen der Daten führen. So wurde bei den hier verwendeten Phantomen häufig die Temperatur im Innern des Kessels viel zu hoch wiedergegeben, weil im Randbereich Strukturen vorlagen, die zu klein waren.

Eine Erhöhung der Anzahl der Meßlinien schafft hier bessere Abhilfe, mit 16 Meßgeräten und damit resultierenden 96 Meßpfaden konnten auch bei fehlerbehafteten Daten gute Rekonstruktionen erreicht werden.

Literaturverzeichnis

[1] Carl de Boor Splinefunktionen Birkhäuser-Verlag, 1990

[2] W. Derichs, F. Heß et.al.

Abschlußbericht des Arbeitskreises ”Schallpyrometrie“ RWE Energie AG, Dezember 1992

[3] W. Derichs, J. König Die Schallpyrometrie — ein Meßverfahren zur Bestimmung der Temperaturverteilung in Kesselfeuerungen DVV-Kolloquium, 1990

[4] W. Derichs, F. Heß, K. Menzel, E. Reinartz Acoustic Pyrometry: A correlation between temperature distributions and the operating conditions in a brown coal fired boiler 2^{nd} International Conference On Combustion Technologies For A Clean Environment. Lissabon, Juli 1993

[5] Vitus Dewenter Bestimmung eines Temperaturfeldes in einem Kessel-Feuerraum mittels der Schallpyrometrie und Erstellung eines Rechnerprogramms zur Temperatursimulation Diplomarbeit. Aachen, März 1992

[6] G.H. Golub, Ch. van Loan Matrix Computations Second Edition. North-Holland Verlag, 1989

[7] Willehad Klaverkamp Tomographische Bildrekonstruktion mit direkten algebraischen Verfahren Dissertation. Münster, 1991

[8] Ute Kretschmann PC-gestützte Fehler-, und Plausibilitätsabfrage zur Bewertung schallpyrometrisch ermittelter Temperaturdaten im Strahlungsraum eines Braunkohlekraftwerkskessels Diplomarbeit. Gelsenkirchen, 1993 [9] M. M. Lavrent’ev, V. G. Romanov, S. P. Shishat-skii

Ill-posed Problems of Mathematical Physics and Analysis American Mathematical Society; Volume 64, 1980

[10] Thorsten Loehr Rekonstruktion, synthetische Erzeugung und Darstellung rechteckiger Temperaturverteilungen Studienarbeit. Braunschweig, Juli 1992

[11] Peter Maaß Singular Value Decompositions for Radon Transforms in: Lecture Notes in Mathematics, Vol. 1497 Herman, Louis, Natterer (Eds.), Mathematical Methods in Tomography Proceedings. Oberwolfach, 1990

[12] Frank Natterer The Mathematics of Computerized Tomography B. G. Teubner, 1986

[13] F. Natterer, B. Setzepfandt, F. Wübbeling Tomography of flames in combustion chambers in: ECAPT: Process Tomography – A Strategy for Industrial Exploitation 1993

[14] Bernd Setzepfandt ESNM: Ein rauschunterdrückendes EM-Verfahren für die Emissionstomographie Dissertation. Münster, Juli 1992

[15] M. I. Sezan, H. Stark Image Restoration by the Method of Convex Projections: Part 2 — Applications and Numerical Results IEEE Transactions on Medical Imaging, Vol. MI-1, No. 2, Oktober 1982

[16] Helmut Sielschott Tomographie mit wenigen Daten: Anwendung von Kollokationsverfahren in der Schallpyrometrie Abschlußbericht. Münster, Juni 1994

[17] Helmut Sielschott Das Programmpaket ”collocation” Bedienungsanleitung. Münster, Juni 1994

[18] Dieter Strackenbrock Entwicklung eines Programmes zur Berechnung der Temperaturverteilung in Brennkammern und Rauchgaskanälen aufgrund von Laufzeitmessungen mit verschiedenen auf den Umfang angeordneten Schallquellen und Schallempfängern Diplomarbeit. Braunschweig, 1987 [19] A. Trächtler

Use of Symmetries in Tomography in: ECAPT: Process Tomography – A Strategy for Industrial Exploitation 1993

[20] D.C. Youla, H. Webb Image Restoration by the Method of Convex Projections: Part I — Theory IEEE Transactions on Medical Imaging, Vol. MI-1, No. 2, Oktober 1982