Bestimmung des Wärmeflusses

Die Bestimmung des Wärmeflusses in den gegebenen Daten wird ausführlich in [Garbe,2001], [Gar-be, C. S. u. a.,2001], [Garbe, C. S. und Jähne, B., 2001] und [Garbe, Christoph S. u. a., 2002] (in Planung: [Brocke,2003]) vorgestellt und soll an dieser Stelle nur kurz zusammengefaßt werden, um die robusten Komponenten des Verfahrens vorzustellen und aufzuzeigen, wie Änderungskarten aus

Abbildung 5.4: Zonen mit Reflektionen und das generelle Wärmemuster der Wasseroberfläche bewegen sich mit unterschiedlichen Geschwindigkeiten (Richtung und Beträge) durch das Bild (vier aufeinanderfolgende Aufnahmen): die rot markierten Zonen der Reflektionen bewegen sich nach oben, während die Flüssigkeit sich eher nach links bewegt.

einer Vorverarbeitung der Daten in dieses Schema eingebaut werden können.

Aus der Kontinuitätsgleichung und den beiden Fick’schen Gesetzen folgt für Transportphänome-ne (Diffusion eiTransportphänome-ner Konzentrationc) allgemein (erste Zeile) und für die WärmemengeV ρcpT bzw.

TemperaturT als spezielle Größe (zweite Zeile)

j=−D∇c _dt^dc=D∆c=−∇j j_h =−D_hρc_p∇T c_pρ_dt^dT =D_h∆T =−∇j_h

Für zwei TemperaturenT₁(z₁)undT₂(z₂)in unterschiedlichen Tiefenz₁undz₂ergibt sich ein Aus-druck, der die Definition einer Transfergeschwindigkeitkh(engl. piston velocity) nahelegt.

jh =−D_hρcp

T₁(z₁)−T₂(z₂) z1−z2

=−k_hρcp∆T mit kh = j_h ρcp∆T

Lösungen für die Transportgleichung der Wärmemenge bzw. Temperatur können unter bestimmten Randbedingungen angegeben werden. Mit den Temperaturen des Wasserkörpers, der Wasseroberflä-che und den RandbedingungenT(z, t= 0) =TbulkundT(z= 0, t) =Tsurferhält man für den Wärme-flußj_h(siehe [Crank,1975])

j_h(0, t) =−D_hρc_pT_surf−T_bulk

√πD_ht

An der Meeresoberfläche gilt eher die Randbedingung eines konstanten Wärmeflusses, was zu

∂(T_surf−T_bulk)

∂z =− j_h

Dhρcp

führt. Mit Analogieschlüssen (siehe [Garbe,2001], vollständige Ableitung bei [Liu, W. T. und Busin-ger, J. A., 1985] [Soloviev, A. V. und Schlüssel, P., 1994]) folgt aus der vorigen Lösung für diese Randbedingung folgenden Zusammenhang ab:

j_h(0, t) = T_surf−T_bulk α√

t mit α= 2

√πD_hc_pρ

Abbildung 5.5: Links: Statistisch verteilte Erneuerungsereignisse bringen die Oberflächentemperatur in unre-gelmäßigen Abständen auf die Temperatur des kälteren Wasserkörpers unterhalb der Grenzschicht. Die mittlere Temperatur ergibt sich aus der aus der Häufigkeit dieser Ereignisse, der TemperaturTtextbulkund dem Aus-kühlungsverhalten an der Oberfläche. Rechts: Schematische Darstellung des Oberflächenerneuerungsmodells

oder umgeformt

Tsurf−Tbulk=αjh

√t−t0 mit: t≥t0

Differenzieren nach der Zeittund erneutes Einsetzen in die Ausgangsgleichung führt schließlich zu j_h =

√2 α

r

(T_surf−T_bulk) d

dt(T_surf−T_bulk)

Als Tracer für den Gasaustausch kann der Wärmefluß an der Oberfläche also bestimmt werden, in-dem∆T:=T_surf−T_bulk und die totale zeitliche Ableitung von∆T bestimmt werden. Beides ist mit Bildfolgen räumlich-zeitlich aufgelöst möglich. Aus dem Wärmeflußjh ergibt sich die Transferge-schwindigkeitk_h, die über die Schmidtzahlen der Medien und den Schmidtzahlexponenten zur Trans-fergeschwindigkeit für CO2-Gas führt. Bevor die Implementierung der Bestimmung von∆Tund der totalen zeitlichen Ableitung davon über die Verarbeitung der IR-Sequenzen und damit auch die da-für genutzte Methode der Ausreißerdetektion vorgestellt wird, soll zuvor das Modell der zeitlichen Verteilung der Oberflächentemperaturen vorgestellt werden.

Für das Zustandekommen der zeitlichen Wärmeverteilung an der Wasseroberfläche gibt es zahl-reiche Modelle (Übersicht in [Jähne,1980]) für den turbulenten Wärmetransport vom Wasserkörper zur Oberfläche, von denen das Oberflächenerneuerungsmodell [Danckwerts,1970] mit statistisch ver-teilten Erneuerungen die möglichen Exponenten der Schmidtzahl korrekt vorhersagt. Dabei wird an-genommen, daß kleine Wasservolumina in der Grenzschicht an der Oberfläche in zufällig verteilten Zeitabständen mit Wasser aus dem tieferen Wasserkörper durchmischt und auf dessen Temperatur gebracht werden. Zusätzlich kommt es zu einer diffusiven Auskühlung der Wasseroberfläche (siehe Abbildung5.5). Dem Modell zugrunde liegt eine Verteilungsfunktionp(t)der Zeitabstände zwischen solchen aufeinanderfolgenden Erneuerungsereignissen für ein Volumenelement. Die Verteilungsfunk-tion wurde von [Soloviev, A. V. und Schlüssel, P.,1994] als ln-Normalverteilung

p(t) = 1

√πσmt/t⁰exp (

−(ln (t/t⁰)−m)² σ_m²

)

vorgeschlagen, wobeimder Mittelwert derlnt/t⁰ ist undσ_m² die zugehörige Varianz (t⁰ist ein Ska-lenfaktor, der für einheitenfreie Größen sorgt). [Haußecker, H. u. a., 2001] belegte diese Verteilung experimentell. Daraus ergibt sich auch eine mittlere Zeitdauer zwischen den Erneuerungen von

hti=

Zu einer festgehaltenen Erneuerungszeittberechnet [Haußecker,1996] eine Verteilung für die Ober-flächentemperatur von

p(T_surf|t) = 2 (T_surf−T_bulk) t α²j²_h

Unter der Annahme, daß die Erneuerungsereignisse ln-normalverteilt sind, ergibt sich daraus eine Oberflächentemperatur-Verteilung von

für den FallTsurf< Tbulk. Kurvenanpassungen dieser Funktion an die Histogramme von einzelnen IR-Bildern sind in Abbildung5.6gezeigt. Die mittlere OberflächentemperaturT¯_surfwird entweder durch eine einfache Mittelung über das einzelne Bild (oder mehrere aufeinanderfolgende, um eine größere Datenbasis zu haben) berechnet oder über die numerische Berechnung des Integrals

T¯_surf=

Das Modell für die Oberflächen-Erneuerung kann anhand der oben hergeleiteten Bestimmungsglei-chung für den Wärmefluß überprüft werden. In

t−t0 =

Tsurf(t−t0)−Tbulk

αj_h

2

können die beiden Temperaturen aus den Ergebnissen der Analyse der Verteilungen (für einen im Bild konstant angenommenen WertT_bulk) und direkt aus den Daten (fürT_surf(t)) bestimmt werden.

Zur Berechnung der zugehörigen Lebensdauerτ :=t−t₀ fehlt nur der Wert des Wärmeflusses j_h. Dieser wird aus der Ableitung der Oberflächentemperatur nach der Zeit

j_h = 2 α

√τ d

dτT_surf(τ)

entnommen, so daß sich die Lebensdauer an jedem Pixel der Sequenz zu τ = 1

2

T_surf(τ)−T_bulk T˙surf(τ)

ergibt. Wieder wird die totale zeitliche Ableitung der Temperatur benötigt, zu deren Berechnung eine robuste Einbettung des erweiterten Strukturtensors dient. Trägt man die Häufigkeiten der in einer Sequenz vorkommenden Lebensdauernτ zur Überprüfung des Oberflächenerneuerungsmodelles auf, so muß sich die vom Modell vorhergesagte ln-Normalverteilung ergeben, was [Garbe, Christoph S.

u. a.,2002] bestätigen konnte.

Die Bestimmung von∆T aus kalibrierten IR-Sequenzen geschieht über die Histogramm-Analyse der Infrarotbilder, bei derT_bulkdurch den Achsenschnittpunkt der Verteilung bildweise als räumlich konstanter Wert bestimmt wird und∆T sich als Differenz daraus zu den Bilddaten als räumlich-zeitlich aufgelöster Wert ergibt. Zur Bestimmung der totalen räumlich-zeitlichen Ableitung von∆T wird eine robuste Einbettung eines Total-Least-Squares-Verfahrens (TLS) des erweiterten Strukturtensors be-nutzt. In der Bestimmungsgleichung (engl. brightness change constraint equation, BCCE) für den optischen Flußf [Horn, B. K. P. und Schunk, B.,1981] lassen sich mit der üblichen Schreibweise für partielle Ableitungen^∂g/∂i:=g_iein Datenvektordund ein Parametervektorpvoneinander trennen.

dg Das Problem ist schlecht gestellt, weil zwei Größenf = (∂x/∂t, ∂y/∂t)^T aus nur einer Gleichung zu bestimmen sind. In der Regel benutzt man daher aus einer Umgebung um den betreffenden Da-tenpunkt insgesamtn Datenpunkte, nimmt den optischen Fluß f in dieser Umgebung als konstant an und löst dann ein Gleichungssystem ausnGleichungen, das sich als Matrix D·p=0 mit einer n×3-DatenmatrixDschreiben läßt. Bein >2wird für das überbestimmte Problem genau der Parame-tervektorp^∗ gesucht, dessen Summe von Residuen-Quadraten (Betragsquadrat des Residuenvektors

r) n minimal ist. Vergleicht man die totale zeitliche Ableitung

d∆T

dt = ∂∆T

∂t + (u∇) ∆T

mit der BCCE, so ist evident, daß∆T nicht die Grauwertintensität ersetzen kann, da die Bedingung

d/dt∆T = 0 nicht gilt. Ganz im Gegenteil wird gerade dieser Wert ungleich Null gesucht. Unter gleichzeitiger Berücksichtigung einer durch das Geschwindigkeitsfeld u gegebenen Bewegung des Wärmemusters an der Wasseroberfläche und einer in den zitierten Arbeiten für diese Gleichungen als konstant angenommenen Quell-Rate von Wärme aus dem Wasserkörper an die Oberfläche ist die totale zeitliche Ableitung von ∆T zu bestimmen. Würde die Wasseroberfläche ruhen, so wäre die totale zeitliche Änderung der Temperatur durch die lokale Änderung ^∂/∂t∆T an einem festen Ort gegeben und aus den Daten durch Differenzbildung im Bildstapel einfach zu ermitteln. Die BCCE muß also erweitert werden, und die Annahme einer konstanten Änderungcführt zu

gt+ (f∇)g−c= (gx, gy, gt,−1)^T ·(f1, f2,1, c) =d^T ·p= 0

Für die Änderung der Temperatur kann der Datenvektor den IR-Sequenzen und deren Ableitungen entnommen werden, im Parametervektor entsprichtc ≡^d/dt∆T der gesuchten Temperaturänderung.

In einer Umgebung mitnPixeln mit konstant angenommenem Flußf≡uwird die Datenmatrix

Die Minimierung (engl. total least squares, TLS) kann so umformuliert werden, daß der erweiterte Strukturtensor und sein Eigenwertproblem erkennbar werden.

|Dp|= Die Randbedingung p^Tp= 1zur Vermeidung der trivialen Lösungp=0wird mit dem Lagrange-Multiplikatorλin das Funktional

L(p, λ) =p^TJp+λ 1−p^Tp

gekoppelt, so daß sich im Lagrange-Formalismus p^∗ = argmin_pL(p, λ) aus den ^∂L/∂pj

= 0! das Gleichungssystem

Jp=λp

als Eigenwertgleichung mit dem kleinsten Eigenwertλschreibt, der auch der minimale Wert im TLS-Verfahren ist:

min|Dp|=p^TJp=p^Tλp=λ

Gesucht ist also der Eigenvektor p^∗ zum kleinsten Eigenwert des erweiterten Strukturtensors. Die totale Ableitung^d/dt∆T≡c = ^p∗4/p^∗₃ ergibt sich aus den letzten beiden Komponenten dieses Eigen-vektors; das Geschwindigkeitsfelduergibt sich aus(^p∗1/p^∗₃,^p∗2/p^∗₃)^T.

Die Reflexe in den Daten führen zu Ausreißern, die mit einer robusten TLS-Methode durch Least Median of Squares Orthogonal Distances (LMSOD) nach [Bab-Hadiashar, Alireza und Suter, David, 1998] unterdrückt werden, so daß nicht alle Werte einer Umgebung in die Bestimmung des Parame-tervektorsp^∗ eingehen. Ausreißer werden dabei anhand der Residuenr=d^Tpermittelt, für deren Berechnung eine möglichst gute Schätzung vonpbekannt sein muß. Generell werden dabei aus den naus einer Umgebung zur Verfügung stehenden Daten zur Bildung der MatrixD= (d₁, . . . ,d_n)^T solange Unterstichproben gebildet, bis eine Lösung fürp mit möglichst kleinen Residuenquadraten bereitsteht. Dann werden für alle Daten der Umgebung die Residuen erneut berechnet, danach sortiert und für die folgende Berechnung des Parametervektors nur die darin nicht als Ausreißer identifizierten Daten benutzt. Im einzelnen heißt das: es wird eine zufällig ausgewählte Unterstichprobe vonk= 4

DatensätzenDˆ = (d1,d2,d3,d4)^T gebildet, für die sich ein Parametervektor pˆ exakt bestimmen läßt. Für alle andereni= 1, . . . , n−kDaten werden die quadrierten Residuenr²_i= d_i^Tpˆ2

und deren MedianMed r²_i

bestimmt. Dies geschieht für

m= ln (1−Π) ln

1−(1−)^k

Unterstichproben, wobei die geschätzte Rate an Ausreißern ist undΠ die sehr nahe 1 angesetzte Wahrscheinlichkeit, daß mindestens eine Unterstichprobe den Parametervektor ausreißerfrei bestimmt hat. Aus allenminitialen Parametervetorenpˆ wird der mit dem zugehörigen kleinsten Medianwert der Residuenquadrate über den Restdatensatzes benutzt, um damit allej= 1, . . . , nResiduenquadrate r²_j= djTpˆ^∗2

zu berechnen und die Daten nach deren Größe zu sortieren. Mit

|r_j| 1,4826

1 +_n−k⁵ q

Medr²_j

≤

> 2,5 (

kein Ausreißer Ausreißer

werden dann Ausreißer identifiziert. Über die nicht als Ausreißer gekennzeichneten n_a Daten (mit k≤na≤n) wird nun abschließend mit der MatrixDa= (d1, . . . ,dna)^T das Eigenwertproblem des Strukturtensors D_a^TD_a berechnet, dessen Elemente nun über n_a statt ngemittelt werden (es wird eine binominale Filtermaske benutzt).

Rechentechnisch findet in der Verarbeitung der IR-Sequenzen an mehreren Stellen eine Mittelung und Abtastung der Daten statt, bevor die robuste Detektion der durch Reflexe (erscheinen fälsch-licherweise als “zu warme” Regionen) eingebrachten Ausreißer diese aus der Parameterschätzung herausnimmt.

Die256×256großen Bilder werden bei [Garbe,2001] auf der ersten Pyramidenstufe verarbei-tet und dazu bildweise geglätverarbei-tet mit einem Filterkernel¹/16(1,4,6,4,1)^T, um die hohen räumlichen Frequenzen so zu unterdrücken, daß auf der folgenden Pyramidenstufe der Größe 128×128 kein Aliasing-Effekt mehr auftritt [Jähne, 2002]. Auf der neu abgetasteten Pyramidenebene findet eine weitere Glättung mit einem Kernel¹/4(1,2,1)^T statt. Alle folgenden Rechnungen finden auf dieser Skalenebene statt. Die Ableitungen in der erweiterten BCCE werden mit5×5×5optimiertem Sobelfil-ter [Scharr,2000] berechnet. Für das TLS Verfahren werden9×9×3Datenpunkte ausgewählt und der Ausreißerdetektion unterzogen. Für das Eigenwertproblem des erweiterten Strukturtensors wird über die verbliebenenna≤nPunkte mit einer binominalen Fensterfunktionwgemittelt. Die Größe der Masken und der Übergang auf eine höhere Pyramidenebene können zur Folge haben, daß eine vor-handene Ausreißerregion mit ihrem additivem Grauwertbeitrag4bis5Pixel von ihrer Kante aus noch Einfluß auf die Daten hat. Schwierigkeiten können dann entstehen, wenn die für die Ausreißeranalyse benutzte Datenprobe sich dadurch nicht mehr in “gute Daten” und Ausreißer einteilen läßt, sondern auch vermeintlich “gute Daten” durch die Nähe zu einer vorhandenen Ausreißerregion bereits gestört sind.

Im Dokument Ruprecht-Karls-Universität Heidelberg (Seite 152-158)