-14-

Teil B: Grenzschichtbestimmung mittels LIDAR

1. Einleitung Die unterste Schicht unserer Atmosphäre steht in direkter Wechselwirkung mit der Erdoberfläche. Sie ist im Sinne der Strömungsdynamik eine Grenzschicht und wird deshalb als planetare Grenzschicht bezeichnet. Ihre Mächtigkeit ist einem täglichen Rhythmus unterworfen und schwankt zwischen einigen Metern und 5 Kilometer. Die Strömung der Luft ist hier turbulent. Dafür gibt es mechanische und thermische Ursachen, insbesondere Scherströmungen und Konvektion. Da sich unser Leben weitgehend innerhalb der planetaren Grenzschicht abspielt, sind wir mit ihren Eigenschaften gefühlsmäßig vertraut. Fragen nach ihrer vertikalen Struktur und Ausdehnung und danach, worin sie sich genau von anderen Schichten der Atmosphäre unterscheidet, können jedoch in zunehmendem Maße erst durch moderne Fernerkundungsmethoden seit der zweiten Hälfte des 20. Jahrhunderts beantwortet werden [Schwiesow, 1984]. Eine charakteristische Größe der planetaren Grenzschicht ist die Mischungsschichthöhe. Bis hier werden Emissionen vom Boden vertikal gut durchmischt und verteilt. Von der Bestimmung ihrer genauen räumlichen und zeitlichen Struktur hängt nicht nur die Qualität der Vorhersage meteorologischer Größen ab, sondern in ganz entscheidendem Maße auch die Vorhersage der Konzentration von unliebsamen Luftbeimengungen. Für gewöhnlich ist die Mischungsschicht wesentlich stärker verschmutzt als andere Luftschichten. Das Verdünnungspotenzial für bodennah emittierte Spurenstoffe ist von der vertikalen Mächtigkeit des konvektiv durchmischten Raumes am Tage abhängig. Eine niedrige Mischungsschichthöhe führt zu einer Ausprägung bzw. Verschärfung lufthygienisch ungünstiger Bedingungen. Die Mischungsschichthöhe ist eine wesentliche Skalierungsgröße bei der Modellierung des vertikalen turbulenten Austauschkoeffizienten und der zugehörigen turbulenten Vertikalflüsse von Wärme, Feuchtigkeit, Impuls und Stoffkonzentrationen. Nach Sonnenuntergang sinkt die Obergrenze auf wenige hundert Meter ab, und es bleibt eine nächtliche Residualschicht. Sie ist mitbestimmend für die Menge des verfügbaren Ozons, das in der Phase des Übergangs von der nächtlichen stabilen zur tagsüber meist labil geschichteten planetaren Grenzschicht eingemischt werden kann. Betrachtet man mittels LIDAR die vertikalen Profile des Aerosol-Rückstreukoeffizienten, hebt sich die freie Troposphäre deutlich durch ihre Sauberkeit von der bodennahen Grenzschicht ab. Aerosol, das durch Aufwirbelung und chemische Reaktionen der vom Boden stammenden Stoffe gebildet wird, gelangt nur sehr selten über die Barriere zur freien Troposphäre. Die verschieden stark streuenden Schichten machen aber deutlich, dass Luftbeimengungen innerhalb der Mischungsschicht nicht homogen verteilt sind. Farbliche Darstellungen (Vgl. Abb. 1) der zeitlichen Entwicklungen der Höhenprofile geben ein eindrucksvolles Zeugnis der komplizierten Struktur.

1.1 Fragestellung Das Ziel dieses Teils des Praktikums war es, die LIDAR Daten in Bezug zur Grenzschichtmeteorologie zu stellen. Dabei sollte insbesondere folgenden Fragen nachgegangen werden: a. Wie lassen sich aus Ceilometerdaten Grenzschichthöhen bestimmen? Unter welchen meteorologischen Bedingungen ist dies möglich?

Abb. 1: LIDAR-Plot

-15-

b. Wie hängen die gemessenen Schadstoffkonzentrationen von den gemessenen Grenzschichthöhen ab?

2. Methoden und Daten In diesem Kapitel soll zuerst ein kleiner Einblick in die Funktionsweise des LIDAR gegeben werden. Anschliessend soll aufgezeigt werden, mit welchem Ansatz wir aus den Ceilometerdaten versuchten die Grenzschichthöhe zu bestimmen. Zum Schluss wird dann das Matlab-Programm zur Datenauswertung erläutert.

2.1. Bestandteile eines LIDAR Ein Laser-Fernerkundungssystem besteht aus einer Sende- und einer Empfangseinheit. Senderseitig ist der Laser das Herzstück. Hinzu kommen noch Strahldiagnostik und Aufweitungsoptik. Für bestimmte Messungen wird auch eine Vorrichtung benötigt, um die Sendestrahlrichtung ändern zu können. Der Empfänger besteht aus einem Teleskop, den Tageslichtfiltern, den Detektoren und den Signal-Digitalisierern bzw. Photonenzählern. Meist sind noch weitere Optiken zur Signalaufbereitung notwendig. Ein Computer wird für die Steuerung des Laser und der Mechanik, zum Speichern der Daten und zur Überwachung von Systemparametern verwendet. Eine zusätzliche Schwierigkeit besteht darin, dass alle Teile gemeinsam nicht nur im Labor, sondern auch im freien Feld oder auf Expeditionen funktionieren müssen. Dieser Aspekt sollte bei der Auswahl der physikalischen Methoden unbedingt berücksichtigt werden.

2.2. Ansatz Wie in Abbildung 1 ersichtlich, lässt sich die zeitlichen Entwicklungen der Höhenstruktur auf den LIDAR Bilder gut erkennen. Die Bilder lassen erahnen, bis auf welche Höhe sich die Mischungsschicht erstreckt, geben aber keinen objektiven Anhaltspunkt. Ausgehend von der Überlegung, dass sich das LIDAR Signal im Bereich der Oberkante der Grenzschicht am Stärksten ändern müsste, verfolgten wir zur Grenzschichthöhen- bestimmung den Ansatz, nach dem stärksten Gradienten im linearisierten Rückstreusignal zu suchen.

2.2. Vorgehensweise und Programmierung Die Auswertung der Daten erfolgte mittels Matlab. Im Folgenden wird das von uns codierte Programm kurz erklärt. Die ersten fünf Programmmodule (Importieren der Lidarrohdaten, Bildmatrix erstellen, Relative Korrektur mit Response-Function, Echtskaliereung reskalieren, Regenfilter) haben wir vom IAC der ETH Zürich erhalten. Diese beschreiben wir nur kurz, um einen Verständnis für das Programm zu ermöglichen. Die darauf folgenden Module des Programms (Gradientenmatrix erstellen, Der Finder, Das Rauschen, Die Grenzschichthöhe) haben wir erstellt und werden demnach auch ausführlicher erklärt. Die Bilder sollen jeweils den Effekt der einzelnen Programmabschnitte visualisieren.

Importieren der Lidarrohdaten: Die LIDAR Rohdaten für einen Tag werden als Vektor abgespeichert. Dieser wird importiert. worm = textread('21_4.dat','%s'); Filename des Datafiles eingeben a=0; for i=1:1:size(worm,1) if worm{i,1}=='S' a=a+1; sv(a,1)=i;

-16-

end if worm{i,1}=='E' ev(a,1)=i; end end Nun gilt es den Vektor, welcher die Messdaten eines Tages enthält in eine Matrix umzuwandeln. Jede Spalte der Matrix entspricht einer vierminütigen gemittelten Messung über die gesamte Höhe von 4.5 Kilometer. Somit repräsentiert die Bildmatrix die Lidarmessung eines Tages über die Höhe von 4.5 Kilometer. Jeder Pixel im Bild entspricht einem Wert der Matrix und somit einem Messewert des LIDAR. Zwei aufeinander folgende Pixel einer Spalte entsprechen einem Höhenunterschied von 10 Meter.

Bildmatrix erstellen: matrix = NaN(244,size(sv,1)); a=2; for i=1:1:size(sv,1) for j=2:1:(ev(i,1)-sv(i,1)-1) matrix(j-1,i)=str2num(worm{sv(i,1)+j}); end end worm = [];

Relative Korrektur mit Response-Function: Die Korrektur wird benötigt, weil das Gerät eine nicht konstante Empfindlichkeit mit der Entfernung hat. Es wird auf 180 Pixel reduziert. A=0; B=0; for counter=1:length(matrix); A=A+1; B=0; for counter2=1:180 B=B+1; matrix_corr1(B,A)=matrix(B,A)/correction(B,1); end end

Echtskalierungen reskalieren: Um Speicherplatz zu sparen, wurden die ersten 500m mit einer Auflösung von 10m gemessen, die nächsten 1000m mit 20m Auflösung, die nächsten 1500m mit 30m und die restlichen Meter mit 50m Auflösung. Das ganze wird dann wieder auf eine lineare Entfernungsskala zurücktransformiert. skala10(:,:)=matrix_corr1(1:50,:); skala20(:,:)=matrix_corr1(51:100,:);

-17-

skala30(:,:)=matrix_corr1(101:150,:); skala50(:,:)=matrix_corr1(151:180,:); c=0; for count=1:length(matrix_corr1) c=c+1; aa=-1; b=0; d=-2; e=0; f=-4; g=0; for count1=1:50 aa=aa+2; b=b+2; d=d+3; e=e+3; rescale20(aa:b,c)=skala20(count1,c); rescale30(d:e,c)=skala30(count1,c); end for count2=1:30 f=f+5; g=g+5; rescale50(f:g,c)=skala50(count2,c); end end matrix_resc(:,:)=vertcat(skala10,rescale20,rescale30,rescale50); Abbildung 2 zeigt einen typischen Lidarplot, wie er zum Beispiel an Flughäfen benutzt wird, um die Bewölkungssituation zu erkennen.

Regenfilter: Dieses Programm löscht alle Spalten der Matrix, welche einen sehr hohen Wert enthalten. Ein extrem hoher Wert bedeutet eine sehr hohe Reflektivität und lässt somit auf grosse Partikel, sprich Regentropfen schliessen. matrixRegen=NaN(size(matrix_resc)); Dieser Codefetzen bildet eine leere Matrix `matrixRegen` der gleichen Grösse, wie die Matrix `matrix_resc`. NaN wird in alle Zellen der leeren Matrix geschrieben und dient als Platzhalter ohne Wert. In die Matrix `matrixRegen` werden nun an der entsprechenden Stelle alle Spalten der Matrix `matrix_resc` eingefügt, welche keine zu hohen Werte aufweisen. Alle anderen Spalten bleiben leer. Für unsere Anwendung birgt dieser Filter einen zusätzlichen Vorteil. Da wir vor allem Tage betrachten wollen, an denen es nicht geregnet hat, können wir mit diesem Programm schnell eine Vorselektion betreiben. a=0;

Abb. 2: Plot nach Skallierung

-18-

B=0; for counter=1:length(matrix_resc) a=a+1; x=find(matrix_resc(:,a)>100000); y=isempty(x); if y==1 matrixRegen(:,a)=matrix_resc(:,a); else B=B+1; matrixRegen(:,a)=NaN; end x=[]; end

Gradientenmatrix erstellen: Dieses Programm berechnet aus der Matrix `matrixRegen` den Gradienten aller Werte einer Spalte und fügt diese Werte an der entsprechenden Stelle in die neue Matrix `Gradve` ein. Der Gradient: Wir beginnen in der ersten Spalte mit dem untersten Pixel (Messung um Mitternacht auf Höhe 0 km). Wir berechnen den Gradienten, indem wir über die ersten zehn Pixel und den zehn darauf folgenden Pixels mitteln und die Differenz bilden. Dieser Wert wird in der neuen Matrix an der entsprechenden Stelle eingefügt, das heisst in der ersten Spalte an der untersten Stelle. v1=0; for v1=0:1339 v1=v1+1; u1=0; for u1=0:430 u1=u1+1; p=mean(matrixRegen(u1+10:u1+19,v1)); q=mean(matrixRegen(u1:u1+9,v1)); y=q-p; Gradve(u1,v1)= [y]; end end minimum = min(min(Gradve,[],1)); maximum = max(max(Gradve,[],1)); Nun suchen wir jeweils den grössten - und den kleinsten Wert in der neuen Gradientenmatrix. Diese beiden Werte benötigen wir, um die Grenzen für die Farbskala der Bildmatrix zu definieren. imagesc(Gradve, [minimum maximum]); axis xy; Abbildung 3 zeigt die Gradienten an. Deutlich zu erkennen ist eine Schicht mit einem starken Gradienten (durch Pfeil markiert). Darunter sind noch weitere aber bei weitem schwächere

Abb. 3: Plot der Gradienten

-19-

Schichten zu sehen. Weiter erkennt man das sogenannte Rauschen (mit Kreis markiert). Dieses Rauschen ist technische bedingt und kein eigentliches Rückstreusignal.

Der Finder: Die Idee hinter dem Finder-Programm ist, dass betragsmässig grosse Werte (grossen Gradient) einen Übergang von einer Schicht in eine andere Schicht zeigen. Dieses Programm soll nun aus der Gradientenmatrix eben diese Werte suchen, welche auf eine Grenzschicht schliessen lassen. Da die Grenzschicht je nach Tag verschieden stark ausgebildet ist, und somit auch die Grenzwerte variieren, ist es nicht möglich einen absoluten Grenzwert zu definieren. Deshalb haben wir für jeden Tag versucht einen relativen Grenzwert mittels Quantil zu finden. In unserem Fall haben wir nach dem maximalen 25% Quantil pro Spalte gesucht. Allen Werten (und somit Pixels), welche kleiner als dieses Kriterium sind, erhalten in der neuen Matrix `Layervek` an der entsprechenden Stelle den Wert Eins, alle anderen Werte den Wert Null. Grenzwert = prctile(Gradve,[25]); Grenzwert2 = max(max(Grenzwert,[],1)); Layervek=NaN(size(matrixRegen)); w1=0; for w1=0:600 w1=w1+1; z1=0; for z1=0:430 z1=z1+1; x=find([Gradve(z1,w1)]<=Grenzwert2); if x==1 Layervek(z1,w1)= [1]; end end end imagesc(Layervek,[-1 1]); axis xy; Abbildung 4 zeigt nun die maximalen 25% Quntile an. Dieses Bild wurde nicht mehr für den ganzen Tag geplottet, da es für unsere Fragestellung genügt bis 10 Uhr zu plotten.

Das Rauschen: Dieses Programm filtert das Rauschen (siehe Abb.3 und Bemerkung dazu) weg. Oberhalb der Grenzschicht liegen die Grenzschichtwerte (Wert Eins in der Matrix `Layervek`) unzusammenhängend in der Bildmatrix verteilt. Deshalb filtern wir alle Pixelgruppen der Matrix `Layervek` heraus, welche in der Horizontalen pro 20 aufeinander folgenden Pixel weniger als 90 Prozent den Wert Eins haben. Die Werte, welche übrigbleiben, werden neu in der Matrix `Grenzs` an der entsprechenden Stelle eingetragen. w2=0; r=20;

Abb. 4: Plot nach Finder

-20-

for w2=0:(600-(r+1))%von1339 w2=w2+1; z2=0; for z2=0:430 z2=z2+1; e=r*(mean(Layervek(z2,w2:w2+r))); if e>=0.9*r Grenzs(z2,w2)=[1]; else Grenzs(z2,w2)=[0]; end end end imagesc(Grenzs,[-1 1]); axis xy;

Abbildung 5 zeigt nun die Schicht mit dem von uns definierten relativen Grenzwertbereich an.

Die Grenzschichthöhe: Dieses Programm bestimmt zum Schluss noch die mittlere Grenzschichthöhe: Dazu werden die Anzahl Pixel einer Höhe aus der Matrix `Grenzs` (entspricht einer Zeile der Matrix) mit eben dieser Höhe multipliziert. Dies wird für alle Höhen durchgeführt. Summiert man alle Werte auf und teilt man diese Summe durch die gesamte Anzahl Pixel, resultiert daraus die mittlere Grenzschichthöhe. Die Zahlenwerte, welche wir erhalten haben, werden in „3.Resultate“ aufgeführt. h=0; total=0; for h=0:430 h=h+1; zeilensumme=(600-r)*(mean(Grenzs(h,1:(600-r)))); sumvektor(h)=[zeilensumme]; g=h*zeilensumme; total=total+g; end sumsumvektor=sum(sumvektor); GsHoehe=total/sumsumvektor;

3. Resultate Die Bestimmung der Grenzschichthöhe erfolgte für die Periode zwischen dem 1.4.2007 und dem 24.4.2007. Wir haben folgende Grenzschichthöhen erhalten: 01.04.2007: GS nicht bestimmbar 02.04.2007: GS nicht bestimmbar 03.04.2007: Daten nicht vollständig

Abb. 5: Plot ohne Rauschen

-21-

04.04.2007: Daten nicht vollständig 05.04.2007: GS nicht bestimmbar 06.04.2007: GS nicht bestimmbar 07.04.2007: GS nicht bestimmbar 08.04.2007: GS nicht bestimmbar 09.04.2007: GS nicht bestimmbar 10.04.2007: 229 Meter, gut ersichtlich , x 11.04.2007: GS nicht bestimmbar 12.04.2007: 171 Meter, schwache GS, x 13.04.2007: 244 Meter, schwache GS, x 14.04.2007: 250 Meter, schwache GS, x 15.04.2007: 159 Meter, x 16.04.2007: 228 Meter, x 17.04.2007: 89 Meter, gut ersichtlich und durchgehend, x 18.04.2007: GS nicht bestimmbar 19.04.2007: GS nicht bestimmbar 20.04.2007: 96 Meter, durchgehend, x 21.04.2007: 238 Meter, gut ersichtlich, x 22.04.2007: 90 Meter, gut ersichtlich, x 23.04.2007: 225 Meter, durchgehend, x 24.04.2007: 196 Meter, x x: Für den Vergleich mit der NOx Konzentration benutzt worden Je nach Wahl der Selektionskriterien (Quantil, Anzahl zusammenhängender Pixel in der Horizontalen) erhielten wir mehr oder weniger eindeutige Grenzschichten. Das Problem bestand hauptsächlich im Herausfiltern des Rauschens. Wurden die Kriterien zu hart gewählt, so wurde das Rauschen sehr gut eliminiert, unter Umständen aber auch die Grenzschicht. Definierte man sie hingegen zu schwach, so wurde die mittlere Höhe durch das Rauschen stark verfälscht. Beantwortung der Frage a: Lassen sich aus Ceilometerdaten Grenzschichthöhen bestimmen? Unter welchen meteorologischen Bedingungen ist dies möglich? Es scheint so, als wäre es möglich aus Celiometerdaten die Grenzschicht zu bestimmen. Die Ergebnisse, sprich die Struktur und die Höhe der Schicht, sowie deren zeitliche Wandel (Bildung in der Nacht und Auflösen bei Tagesanbruch), deuten stark auf eine Grenzschicht hin. Wir können jedoch nicht mit Gewissheit sagen, ob wir die gesamte Grenzschicht sehen oder ob die Höhe der Grenzschicht wirklich der Realität entspricht, da wir doch einige Filterkriterien formuliert haben. Um das zu überprüfen, müssen unsere Ergebnisse erstens systematisch mit anderen Methoden (z.B. Ballonsondierungen etc.) verglichen und evtl. abgeglichen werden. Zweitens müssten viel mehr Messungen durchgeführt werden. Eine schwäche unsere Methode, bzw. des Programms ist, dass wir nur bei sehr gutem Wetter eine klare Grenzschicht finden. Wenn die Wolken hoch oben sind, scheint es möglich zu sein, gute Messungen durchzuführen. Liege die Wolken aber zu tief, ist es nicht möglich eine Aussage zu machen. Beantwortung der Frage b: Wie hängen die gemessenen Schadstoffkonzentrationen von den gemessenen Grenzschichthöhen ab?

-22-

Es ist zu erwarten, dass bei niedriger Grenzschichthöhe die Schadstoffkonzentrationen höher sind, da das Mischvolumen geringer ist. Dagegen ist bei hoher oder nicht-vorhandener Grenzschicht die Akkumulation der Schadstoffe geringer. Um diesen Sachverhalt zu überprüfen, haben wir die elf verwendbaren Grenzschichthöhen gegen die kumulierte NOx-Konzentration des jeweiligen Tages zwischen Mitternacht und 10:00 aufgetragen. Wir glauben mit den von uns gewonnenen Daten eine Tendenz festzustellen. Die NOx-Konzentration scheint tatsächlich in den Fällen höher zu sein, in denen wir eine Grenzschicht gefunden haben. Auch in diesem Fall ist es unbedingt notwendig weitere Messungen durchzuführen.

Grenzschichthöhe vs. [NOx]

R2 = 0.2496

0

50

100

150

200

250

300

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000

Kummulierte NOx-Konz. [ppm]

Gre

nzsch

ich

thö

he [

m]

In der obigen Graphik (Abb.6) sind nicht nur die absolut wolkenfreie Tage aufgeführt. Wir haben alle Tage einfliessen lassen, bei denen wir glauben eine Grenzschicht bestimmt zu haben (siehe 3.Resultate). Die Korrelation ist mit einem R^2 = 0.2496 sehr tief. Wir haben in einem weiteren Versuch nur die Tage gewählt, welche unserer Ansicht nach perfekt geeignet waren. Dabei erhielten wir eine viel höhere Korrelation zwischen er Grenzschichthöhe und der NOx-Konzentration (R^2 = 0.4576). Da aber die Anzahl verwendeter Messpunkte zu klein war, um eine gute Interpretation zu ermöglichen, haben wir dieses Ergebnis bewusst nicht genauer aufgeführt. Doch kommen wir zurück auf die obige Graphik. Auch hier stellt sich die Frage, wie man dieses Ergebnis interpretieren darf, denn zum Einen bestehen wie oben erwähnt grosse Unsicherheiten in der Grenzschichthöhenbestimmung. Zum Anderen wurden einige wichtige Einflussfaktoren auf die NOx-Konzentration nicht berücksichtigt. Sprich die Windrichtung und Windstärke an den Messtagen, die schwankenden NOx Inputs durch den Verkehr an den Werktagen und an den Wochenenden, die Sonnenstrahlung und die damit verbundene Photochemie und weitere Faktoren. Wir denken, dass diese Faktoren im Vergleich zur Grenzschichthöhe die NOx – Konzentration zu stark beeinflussen, als dass die Anzahl Daten, welche wir generiert haben, genügen, um eine klare Aussage zu machen. Wir wollen aber an dieser Stelle nochmals betonen, dass wir an die Möglichkeit glauben, mit Hilfe des LIDAR eine Grenzschichthöhe bzw. die Grenzschicht selbst bestimmen zu können. Hinsichtlich des Programmierens braucht es noch Verfeinerungen und Kriterien, deren Filtereffekte besser abgeschätzt werden können. Die Bestimmung der Grenzschicht mit einem relativen Grenzwert, den wir für jede vierminütige gemittelte Messung über eine Höhe von

Abb. 6: Grenzschichthöhe vs. NOx-Konzentration

-23-

4.5 Kilometer neu bestimmen, ist sicher eine Stärke der Methode und scheint gut zu funktionieren.

4. Literaturverzeichnis Schwiesow, 1984