Iterative Bildrekonstruktion für die Ultraschalltomografie

Johanna Ziegeltrum, März 2015

Ultraschalltomografie

Grundlagen

Der Begriff Tomografie setzt sich aus den beiden altgriechischen Wörtern „τóμοσ“ (tomos) für „Schnitt“ und „γράφειν“ (graphein) für „schreiben“ zusammen. ^[1] Mit Hilfe von Ultraschallwellen werden zweidimensionale Schnitte eines Testobjekts rekonstruiert. Die mathematische Theorie geht zurück auf Radon (1917) welcher zeigt, dass die inneren Gegebenheiten eines Objekts mit Hilfe eines Satzes an Projektionen durch das Objekt exakt rekonstruiert werden können. Die Laufzeit eines Ultraschallsignals wird auf verschiedenen Wegen durch das Testobjekt gemessen. Veränderungen im Objekt führen zu einer veränderten Laufzeit und somit kann bei einer ausreichenden Anzahl an Projektionen ein Schnitt durch das Objekt erstellt werden. ^[2] In Abbildung 1 ist das Grundprinzip der Tomografie dargestellt.

Am Ende soll eine Geschwindigkeitsverteilung des Ultraschallsignals über die Schnittfläche entstehen. Hieraus können Rückschlüsse auf die Eigenschaften und Fehlstellen eines Objekts gezogen werden. Um aus den Projektionsdaten auswertbare Bilder zu erhalten gibt es verschiedene Rekonstruktionstechniken. In diesem Artikel sollen iterative Bildrekonstruktionstechniken näher erläutert werden.

Anwendungsbeispiel für die Ultraschalltomografie

Die Ultraschalltomographie wird hauptsächlich im medizinischen Bereich verwendet. Ein Anwendungsbeispiel im Bauingenieurwesen ist das Auffinden von schlecht verpressten Hüllrohren in nachträglich vorgespannten Spannbetonbrücken. ^[2] Viele der Spannbetonbrücken wurden in den 1960er Jahren erbaut. Fehler beim Verpressen der Hüllrohre, steigende Verkehrslasten und das zur Erbauungszeit fehlende Wissen über die korrosive Wirkung von Chloriden auf Stahl haben zu Bedenken hinsichtlich der Standsicherheit von Spannbetonbrücken geführt. Das alkalische Verpressmaterial schützt den Spannstahl normalerweise vor Korrosion. Außerdem kann in einem richtig verpressten Hüllrohr die Vorspannung auch nach dem Reißen eines Spanngliedes in einem bestimmten Abstand wieder vollständig aufgebaut werden. Somit ist die richtige Verpressung eine bedeutende Größe bei der Beurteilung der Standsicherheit von Spannbetontragwerken. Die Hüllrohre, in denen die Spannglieder verlegt und verpresst werden, sind meist aus Metall. Somit eignet sich Radar als Messmethode nicht, da hier die Eindringtiefe stark von der elektrischen Leitfähigkeit des Materials abhängt. Bei einer hohen elektrischen Leitfähigkeit – wie beispielsweise bei Metallhüllrohren – ist die Wellendämpfung sehr groß und die Eindringtiefe somit zu klein, um Hohlstellen hinter der Metallhülle zu erkennen. ^[3] Besser eignet sich für Metallhüllrohre eine Prüfung mittels Impakt-Echo oder Ultraschalltomographie. ^[2]

Iterative Rekonstruktionsmethoden

Grundlegende Überlegungen für iterative Rekonstruktionsmethoden

Bei diesen Verfahren wird ein Schnitt durch das zu untersuchende Objekt als Feld aus unbekannten Größen dargestellt. Im Fall der Ultraschalltomografie ist dies die Dämpfung des Ultraschallsignals auf seinem Weg durch das Testobjekt. Für die Unbekannten wird ein lineares Gleichungssystem aufgestellt, welche als Ergebnis die gemessenen Daten aus den Projektionen besitzen. In Abbildung 2 ist ein quadratisches Gitter dargestellt, das über das unbekannte Bild $//<![CDATA[ \begin{array}{l}f(x,y)\end{array} //]]>$ gelegt wird. Es wird angenommen, dass $//<![CDATA[ \begin{array}{l}f(x,y)\end{array} //]]>$ in jeder Zelle einen konstanten Wert besitzt. Somit beschreibt $//<![CDATA[ \begin{array}{l}f_j\end{array} //]]>$ einen konstanten Wert in Zelle $//<![CDATA[ \begin{array}{l}j\end{array} //]]>$. $//<![CDATA[ \begin{array}{l}N\end{array} //]]>$ beschreibt die gesamte Anzahl an Zellen. Ein Weg des Signals durch das Objekt wird dargestellt durch eine Linie mit der Stärke $//<![CDATA[ \begin{array}{l}τ\end{array} //]]>$. $//<![CDATA[ \begin{array}{l}p_i\end{array} //]]>$ ist die aufsummierte Messgröße entlang eines Pfades, $//<![CDATA[ \begin{array}{l}M\end{array} //]]>$ ist die gesamte Anzahl an Pfaden durch das Objekt. $//<![CDATA[ \begin{array}{l}w_{ij}\end{array} //]]>$ wird Wichtungsfaktor genannt und beschreibt, wie viel Prozent der $//<![CDATA[ \begin{array}{l}j\end{array} //]]>$-ten Zellenfläche der $//<![CDATA[ \begin{array}{l}i\end{array} //]]>$-te Pfad einnimmt. ^[4]

Somit kann folgendes Gleichungssystem aufgestellt werden: ^[4]

$//<![CDATA[ \begin{array}{l}\sum_{j=1}^N w_{ij} * f_j = p_i\end{array} //]]>$

$//<![CDATA[ \begin{array}{l}i = 1, 2, ..., M\end{array} //]]>$

Ausgeschrieben lässt sich das Gleichungssystem folgendermaßen darstellten: ^[4]

$//<![CDATA[ \begin{array}{l}w_{11}f_1 + w_{12}f_2 + w_{13}f_3 + ... + w_{1N}f_N = p_1\end{array} //]]>$

$//<![CDATA[ \begin{array}{l}w_{21}f_1 + w_{22}f_2 + w_{23}f_3 + ... + w_{2N}f_N = p_2\end{array} //]]>$

$//<![CDATA[ \begin{array}{l}...\end{array} //]]>$

$//<![CDATA[ \begin{array}{l}w_{M1}f_1 + w_{M2}f_2 + w_{M3}f_3 + ... + w_{MN}f_N = p_M\end{array} //]]>$

Wäre die Anzahl $//<![CDATA[ \begin{array}{l}M\end{array} //]]>$ und $//<![CDATA[ \begin{array}{l}N\end{array} //]]>$ klein, könnte man herkömmliche Methoden der Gleichungslösung anwenden. In der Realität jedoch ist die Anzahl $//<![CDATA[ \begin{array}{l}N\end{array} //]]>$ meist ein sehr großer Wert, beispielsweise ergibt sich bei einer Auflösung von 256 x 265 Pixel eine Anzahl von 65.536 unbekannten $//<![CDATA[ \begin{array}{l}N\end{array} //]]>$. Auch die Anzahl der Projektionen hat für die meisten Anwendungsfälle eine ähnliche Größenordnung, so dass eine Invertierung der [$//<![CDATA[ \begin{array}{l}w_{ij}\end{array} //]]>$] Matrix zur Lösung des Problems mit herkömmlicher Matrizenrechnung und vertretbarem Rechenaufwand unmöglich wird. Für eine große Anzahl $//<![CDATA[ \begin{array}{l}M\end{array} //]]>$ und $//<![CDATA[ \begin{array}{l}N\end{array} //]]>$ gibt es iterative Lösungsmechanismen, deren Grundlage die Kaczmarz-Methode bildet. ^[5]

Kaczmarz-Methode

Durch ein Netz mit $//<![CDATA[ \begin{array}{l}N\end{array} //]]>$ Zellen erhält das zu rekonstruierende Bild $//<![CDATA[ \begin{array}{l}N\end{array} //]]>$ Freiheitsgrade. In einem $//<![CDATA[ \begin{array}{l}N\end{array} //]]>$-dimensionalen Raum beschreibt jede Reihe des Gleichungssystems eine Hyperebene. Falls eine exakte Lösung des Problems existiert, kann sie über den Schnittpunkt aller Hyperebenen in einem Punkt beschrieben werden. Das Vorgehen bei der Kaczmarz-Methode soll im Folgenden der Einfachheit halber im zweidimensionalen Raum erläutert werden. Zu lösen ist folgendes Gleichungssystem: ^{[1] [4]}

$//<![CDATA[ \begin{array}{l}w_{11}f_1 + w_{12}f_2 = p_1\end{array} //]]>$

$//<![CDATA[ \begin{array}{l}w_{21}f_1 + w_{22}f_2 = p_2\end{array} //]]>$

Die Kaczmarz-Methode - graphisch dargestellt in Abbildung 3 - beginnt mit einer beliebigen Startschätzung. Diese Schätzung wird auf die erste Linie projiziert, die der ersten Gleichung des Gleichungssystems entspricht. Der neu erhaltene Punkt wird wiederum auf die zweite Linie projiziert. Dieses Vorgehen wird so lange wiederholt, bis der Unterschied zwischen den beiden letzten Projektionen einen festgelegten Grenzwert unterschreitet. Gibt es eine eindeutige Lösung, so wird diese nach unendlich vielen Iterationen erreicht.

Die Umsetzung der Methode im $//<![CDATA[ \begin{array}{l}N\end{array} //]]>$-dimensionalen Raum erfolgt analog. Zunächst wird ein Startvektor gewählt:

$//<![CDATA[ \begin{array}{l}\vec{f^{(0)}} = (f_1^{(0)}, f_2^{(0)}, ..., f_N^{(0)})^T\end{array} //]]>$

Dieser enthält in der Praxis meist nur Nulleinträge. Der Startpunkt wird daraufhin auf die erste Hyperebene projeziert, die durch die erste Gleichung des linearen Gleichungssystems repräsentiert wird. So entsteht der Vektor $//<![CDATA[ \begin{array}{l}f^{(1)}\end{array} //]]>$, welcher wiederum auf die zweite Hyperebene projeziert wird, um zu $//<![CDATA[ \begin{array}{l}f^{(2)}\end{array} //]]>$ zu führen. Mathematisch kann dieser Vorgang wie folgt beschrieben werden:

$//<![CDATA[ \begin{array}{l}\vec{f^{(i)}} = \vec{f^{(i - 1)}} - \frac{(\vec{f^{i - 1}} * \vec{w_i} - p_i)}{\vec{w_i} * \vec{w_i}}\end{array} //]]>$

Dieses Vorgehen wird so lange wiederholt, bis $//<![CDATA[ \begin{array}{l}f^{(M)}\end{array} //]]>$ entstanden ist. Dieses wird dann wieder auf die erste Hyperebene und alle folgenden projiziert, bis $//<![CDATA[ \begin{array}{l}f^{(2M)}\end{array} //]]>$ entstanden ist und so weiter. Falls eine eindeutige Lösung $//<![CDATA[ \begin{array}{l}f_s\end{array} //]]>$ existiert, gilt: ^[4]

$//<![CDATA[ \begin{array}{l}\lim_{k\to\infty} \vec{f^{(kM)}} = \vec{f_s}\end{array} //]]>$

Bezüglich der Konvergenz kann festgestellt werden, dass der Winkel zwischen den aufeinanderfolgenden Hyperebenen von großer Bedeutung ist. Würden die beiden Geraden in Abbildung 3 beispielsweise senkrecht aufeinander stehen, so werden nur 2 Iterationen bis zur exakten Lösung benötigt. Schließen die Hyperebenen jedoch nur einen sehr kleinen Winkel ein, so wird eine große Anzahl an Iterationsschritten erforderlich. Dies zeigt, dass der Rechenaufwand durch eine sinnvolle Anordnung der Gleichungen erheblich verringert werden kann. Benachbarte Pfade durch das zu untersuchende Objekt verlaufen nahezu parallel. Deshalb ist es sinnvoller, weit auseinander liegende Pfade im Gleichungssystem aneinander zu reihen. ^[4]

Bei der Bildrekonstruktion sind überbestimmte Gleichungssysteme mit $//<![CDATA[ \begin{array}{l}M > N\end{array} //]]>$ und Rauschen üblich. In diesen Fällen existiert keine eindeutige Lösung. Mit der Kaczmarz-Methode kann jedoch trotzdem eine gute Näherung erreicht werden. In Abbildung 4 wird dies wieder in der zweidimensionalen Ebene gezeigt. ^[4]

Liegen unterbestimmt Systeme mit M<N vor, kann keine eindeutige Lösung gefunden werden, es sind unendlich viele Lösungen möglich. In diesen Fällen beschreibt die Kaczmarz-Methode wie in Abbildung 5 den Punkt, der am nächsten zur Startschätzung liegt. Die iterative Näherung beschreibt somit eine Lösung $//<![CDATA[ \begin{array}{l}f_s\end{array} //]]>$ so dass gilt: ^[4]

$//<![CDATA[ \begin{array}{l}|\vec{f^{(0)}} - \vec{f_s'}|\end{array} //]]>$ ist minimal.

ART (Algebraic Reconstruction Techniques)

Bei der ART Methode handelt es sich im Grunde um die Kaczmarz-Methode. Um die Implementierung zu vereinfachen, wird die Gleichung für die Projektionen leicht abgewandelt:

$//<![CDATA[ \begin{array}{l}f_j^{(i)} = f_j^{(i - 1)} + \frac{p_i - q_i}{\sum_{k=1}^N w_{ik}^2} * w_{ij}\end{array} //]]>$ wobei $//<![CDATA[ \begin{array}{l}q_i = \vec{f^{(i - 1)}} * \vec{w_i} = \sum_{k=1}^N f_{k}^{(i - 1)} * w_{ik}\end{array} //]]>$

In dieser Gleichung ist $//<![CDATA[ \begin{array}{l}p_i\end{array} //]]>$ die gemessene Summe des $//<![CDATA[ \begin{array}{l}i\end{array} //]]>$-ten Pfades, während $//<![CDATA[ \begin{array}{l}q_i\end{array} //]]>$ die berechnete Summe des $//<![CDATA[ \begin{array}{l}i\end{array} //]]>$-ten Pfades mit Hilfe der Daten aus der $//<![CDATA[ \begin{array}{l}(i - 1)\end{array} //]]>$-ten Lösung ist.

$//<![CDATA[ \begin{array}{l}\Delta f_j^{(i)} = f_j^{(i)} - f_j^{(i - 1)} = \frac{p_i - q_i}{\sum_{k=1}^N w_{ik}^2} * w_{ij}\end{array} //]]>$

Das heißt um die $//<![CDATA[ \begin{array}{l}i\end{array} //]]>$-te Projektion zu erhalten wird zu der $//<![CDATA[ \begin{array}{l}(i−1)\end{array} //]]>$-ten Projektion ein Korrekturwert $//<![CDATA[ \begin{array}{l}\Delta f\end{array} //]]>$ hinzugezählt. Dieser Korrekturwert wird berechnet, indem zuerst der Unterschied zwischen gemessener und berechneter Summe normiert und dann mit dem jeweiligen $//<![CDATA[ \begin{array}{l}w_{ij}\end{array} //]]>$ gewichtet wird. ^[4]
Aufgrund der großen Anzahl an zu bestimmenden $//<![CDATA[ \begin{array}{l}w_{ij}\end{array} //]]>$-Werten und dem damit verbundenen enormen Rechenaufwand werden verschiedene Näherungen verwendet. So erhalten die Einträge für $//<![CDATA[ \begin{array}{l}w_{ik}\end{array} //]]>$ in vielen ART Implementierungen entweder eine 1 oder 0, je nachdem ob der Schwerpunkt der $//<![CDATA[ \begin{array}{l}k\end{array} //]]>$-ten Zelle im $//<![CDATA[ \begin{array}{l}i\end{array} //]]>$-ten Pfad liegt oder nicht. Somit ist $//<![CDATA[ \begin{array}{l}\sum_{k=1}^N w_{ik}^2 = N_i\end{array} //]]>$ die Anzahl an Zellen, die im $//<![CDATA[ \begin{array}{l}i\end{array} //]]>$-ten Pfad liegen und der Korrekturwert vereinfacht sich: ^[4]

$//<![CDATA[ \begin{array}{l}\Delta f_j^{(i)} = \frac{p_i - q_i}{N_i}\end{array} //]]>$

Bilder, die mit Hilfe von ART rekonstruiert wurden, weisen oft „salt and pepper“ Rauschen auf. Dabei handelt es sich um weiße und schwarze Pixel, die über das Bild verteilt auftreten, beispielhaft zu sehen in Abbildung 6. ^[6] Entstehen kann dieses Rauschen aufgrund von schlechten Näherungen für die $//<![CDATA[ \begin{array}{l}w_{ik}\end{array} //]]>$-Werte. Der Effekt kann durch Relaxation reduziert werden, indem der Korrekturwert mit $//<![CDATA[ \begin{array}{l}α\end{array} //]]>$ multipliziert wird, wobei $//<![CDATA[ \begin{array}{l}α < 1\end{array} //]]>$. ^[4]

SIRT (Simultaneous Iterative Reconstruction Technique)

SIRT weist eine schlechtere und langsamere Konvergenz als ART auf, produziert jedoch bessere Bilder. ^[5] Der Korrekturwert $//<![CDATA[ \begin{array}{l}\Delta f_j^{(i)}\end{array} //]]>$ wird genauso ermittelt wie bei ART, jedoch wird der Wert der $//<![CDATA[ \begin{array}{l}j\end{array} //]]>$-ten Zelle noch nicht zu dieser Zeit geändert. Zuerst werden alle Gleichungen durchlaufen und am Ende jeder Iteration wird der Zellenwert um einen durchschnittlichen Korrekturwert dieser Zelle verändert. In der zweiten Iteration wird wieder in der ersten Gleichung begonnen. ^{[1] [4]}

SART (Simultaneous Algebraic Reconstruction Technique)

Diese Rekonstruktionstechnik kombiniert die Vorteile der beiden Techniken ART und SIRT. Dies bedeutet, dass Rekonstruktionen mit guter Qualität und numerischer Genauigkeit in nur einer Iteration möglich sind. ^{[1] [4]}

Literatur

• Zabak, F.: Robust Algorithms for Discrete Tomography. Institute of Applied Mathematics, Delft University of Technology. Juni 2012.
• Martin, J.; Broughton, K.; Giannopolous, A.; Hardy, M. und Forde, M.: Ultrasonic tomography of grouted duct post-tensioned reinforced concrete bridge beams . NDT&E International, pp. 107 - 113, 2001.
• Große, Christian: Grundlagen der zerstörungsfreien Prüfung (Kapitel 8: Radar). TU München, SS 2014, pp. 67 - 69.
• Kak, A. C. und Slaney, M.: Principles of Computerized Tomographic Imaging (Algebraic Reconstruction Algorithms). Society of Industrial and Applied Mathematics, 2001, pp. 275 - 296.
• Carrubé, E. C.: Inversion Methods. pp. 63 - 86.
• Wikipedia.org: Salt-and-pepper noise. Zugriff am 05 Februar 2015.
• Begin, G.: Note on Ultrasonic Tomography. Journal of Nondestructive Evaluation, Bd. 2, Nr. 2, pp. 147 - 149, 1981.

Einzelnachweise

1. Zabak, F.: Robust Algorithms for Discrete Tomography. Institute of Applied Mathematics, Delft University of Technology. Juni 2012.
2. Martin, J.; Broughton, K.; Giannopolous, A.; Hardy, M. und Forde, M.: Ultrasonic tomography of grouted duct post-tensioned reinforced concrete bridge beams . NDT&E International, pp. 107 - 113, 2001.
3. Große, Christian: Grundlagen der zerstörungsfreien Prüfung (Kapitel 8: Radar). TU München, SS 2014, pp. 67 - 69.
4. Kak, A. C. und Slaney, M.: Principles of Computerized Tomographic Imaging (Algebraic Reconstruction Algorithms). Society of Industrial and Applied Mathematics, 2001, pp. 275 - 296.
5. Carrubé, E. C.: Inversion Methods. pp. 63 - 86.
6. Wikipedia.org: Salt-and-pepper noise. Zugriff am 05 Februar 2015.