Modelldokumentation

E.-W. Reiche & Martin Meyer

(Stand 6.98)

1 Name und Typ des Modells

Modellname: TOPNEWP5

Programm zur Höhenrasterauswertung (Identifizierung von topographische Einzugsgebieten, Hag- und Senkenzuordnung, Bestimmung der Fläche von hangaufwärts gerichteten Einzugsgebieten, Berechnung des Hanglängenfaktors nach USLE, Berechnung des langfristig mittleren Bodenabtrags durch Wasser unter Einbeziehung von Vegetation (C-Faktor), Bodenbeschaffenheit (Kfaktor) und Erosivität der Niederschläge (R-Faktor).

Autor: E.-W. Reiche, Projektzentrum Ökosystemforschung der Universität Kiel, Schauenbur-gerstr. 112, 24118 Kiel; Mitarbeit Martin Meyer.

Das Modell TOPNEWP2 wurde für eine Auswertung von Höhenrastern, wie sie beispielsweise durch die digitalen Höhenmodelle (DHM) des Landesvermessungsamtes Schleswig-Holstein (LVA) angeboten werden, entwickelt. Das Modell bietet u. a. die Möglichkeit, Teileinzugsgebiete einer vorgegebenen Fläche zu berechnen und entsprechende Karten zu erstellen. Es handelt sich um ein statisches Analysewerkzeug, dessen Kern eine auf 8 Himmelsrichtungen ausgerichtete Nachbarschaftsanalyse darstellt. Nach Identifizierung von maximalen Abflußrichtungen - bezogen auf jede Rasterfläche - und lokalen Senkenflächen (= Flächen ohne Abflußrichtung) werden topographische Einzugsgebiete durch eine iterierende Gefälleabfrage definiert. Senken und Hangformen werden auf der Grundlage von Gefällecharakteristika ausgewiesen. Für jede Einzell-Rasterzelle wird ausgezählt, wieviele Rasterzellen dieser in Abflußrichtung vorglagert sind. Dabei finde konvergierende Hangformen Berücksichtigung. Der von AUERSWALD et al beschreibene Hanglängenfaktor zur Berechnung des Bodenabtrags durch Wasser kann auf dieser Grundlage (Teilgefälle, Gesamthanggefälle, Hanglänge) berechnet werden und in Kombination mit C-(Vegetation) und K-Faktoren (Boden) einer Abschätzung der Bodenerosion dienen. Die Berücksichtigung von linienhaften Elemente (Wälle, Dämme, Knicks etc.) in ihrer Funktion als Abflußbarrieren ist bei der Erosionsberechnung möglich.

1. Festlegung der Abflußrichtung jeder Rasterzelle und Berechnung des maximalen Gefälles.

2. Ableitung der Teileinzugsgebiete aus dem DHM (Digitales Höhenmodell).

3. Differenzierung zwischen Senken, Hängen und sonstigen Flächen.

4. Berechnung der Anzahl den Rasterflächen vorgelagerten Teileinzuggebiete.

5. Bestimmung des Hanglängenfaktors mit Hilfe der modifizierten USLE-Gleichung nach Auerswald et al. (1988).

Die räumliche Basis zur Berechnung der o.g. Punkte durch das Programm TOPNEWP2 stellt ein digitale Höhenmodell (DHM) (Quelle: Landesvermessungsamt) dar. Die Maschenweite und Höhengenauigkeit ist grundsätzlich nicht vorgegeben. Tests haben gezeigt, daß eine Genauigkeit des DHM 5 ( 12,5 x 12,5 m) bzw DHM 25 (25 * 25 m) für die genannten Fragestllungen ausreichend ist. Das DHM 50 (50 * 50 m-Raster) erscheint für die differenzierte Reliefanalyse als nicht ausreichend. Bei der Auwertung werden wird rasterorientiertder vorgegangen, wobei jeder Rasterpunkt als Zentrum eine quadratische Fläche repräsentiert (s. Abb. 1):

Abb. 1: Schematische Darstellung des Höhenrasters mit der Abflußrichtung der Rasterzellen

In dem ersten Programmschritt werden die Abflußrichtungen, in der oberen Abbildung durch Abflußrichtungsnummern (entsprechend der Himmelsrichtungen) gekennzeichnet, für jede Zelle bestimmt und das maximale Gefälle berechnet. Die Abflußrichtungskennung wird wie folgt interpretiert: -1 = Senke, 0 = Zentralpunkt, 1-8 = Himmelsrichtungen im Uhrzeigersinn von Nordost bis Nordwest.

Aus dem Höhenraster des DHM erfolgt zunächst die Generierung eines Datensatzes im Rechteckformat, auf dessen Grundlage die Teileinzugsgebiete mit einem von REICHE programmierten Algorithmus berechnet wird. Zur Festlegung der Teileinzugsgebiete werden dann die signifikanten lokalen Abflußminima (lokale Senke) der zu bearbeitenden Karte bzw. Kartenblätter, d.h. Punkte, deren Höhe niedriger als die ihrer acht Nachbarpunkte ist, identifiziert und nach steigender Höhe sortiert.

Unter der Annahme, daß das Oberflächenwasser in Richtung des stärksten Gefälles im Sinne eines Gradienten abfließt, wird für jede lokale Senke geprüft, mit welchen der benachbarten Flächen dieser Rasterpunkte hydrologisch verbunden ist. Ein Gradient liegt vor, wenn die Exposition einer der benachbarten Rasterzellen dem jeweiligen Zentralpunkt zugewendet ist und somit hydrologisch gleichhoch oder höhergelegen, also "upstream" von ihm liegt. Dementsprechend ordnet das Programm den nach der Höhe sortierten Abflußminima , die "upstream" gelegenen Punkte dem nächsten Gradienten zu und legt ihre Teileinzugsgebietsgröße fest. Aufgrund des beschriebenen Zuordnungsverfahrens wird somit allen durch die Rasterpunkte repräsentierten Flächen nur eine Hangneigung zugewiesen, wodurch sie nur zu einem Teileinzugsgebiet gehören können. Alle mit einer lokalen Senke hydrologisch verbundenen Rasterpunkte bilden daher ein Einzugsgebiet und werden anschließend zu Polygonen aggregiert.

Zur Beschleunigung des Algorithmus werden alle höhergelegenen und gleichhohen Nachbarpunkte die einem Zentralpunkt zugeordnet wurden, markiert, so daß bereits "abgearbeitete" Punkte im Laufe der Rechenprozedur nicht weiter berücksichtigt werden.

Ergebnis der Berechnung sind Punktekarten lokaler Abflußminima signifikanter Teileinzugsgebiete mit den Abflußrichtungen der zugeordneten Rasterpunkte. Deren jeweilige Teileinzugsgebiets- und Abbflußrichtungsnummer werden durch das Programm als Variablen abgespeichert.

Das codierte Raster dient als Grundlage zur Erstellung von Teileinzuggebietskarten, wobei aus speicherplatz- und rechenzeittechnischen Gründen auf die Triangulations- und Latticeprozeduren des TIN-Moduls (TIN= Trianguläres Irreguläres Netzwerk) von ARC/INFO zurückgegriffen werden kann. Seit der Version ARC/INFO 6 bietet sich zudem eine Polygon-Generierung, basierend auf der Einzugsgebietsnummer, über das GRID-Modul an.

Die Unterscheidung zwischen Senken, Hängen und sonstigen Flächen, die den erstgenannten nicht zugeordnet werden können, wie beispielsweise ebene Flächen, Kuppen etc. erfolgt nach folgenden Bedingungen.

- Senkenterm:

Die Identifikation der Senken erfolgt durch die Bestimmung eines maximalen Höhenabstandes HD (Höhendifferenz) zwischen der lokalen Senke (S) und den benachbarten Rasterzellen. Benachbart heißt hier, daß sich die Abflußlinien zwischen Senken- und Rasterpunkten, ausgehend von den Abflußminima, ohne Unterbrechung fortsetzen (Gleichungen 1-3).

(1) HS - Hi _ HD

(2) Ri _ | Ri - 1 + 1 |

(3) Hi _ | Hi - 1 + 1 |

mit:

H= Hochwert

R= Rechtswert

HS= Hochwert (lokale Senke)

HD= Höhendifferenz

- Hangterm:

Für die Identifikation von Hängen werden im allgemeinen zwei Kriterien herangezogen Zum einen muß ein Mindestgefälle erreicht sein, zum anderen müssen mehrere Rasterflächen hintereinander gelegen eine Mindesthanglänge erreichen. Letztere kann interaktiv festgelegt werden.

Im einzelnen läuft die Routine folgendermaßen ab:

a. Ermittlung der höchsten Rasterfläche des Einzugsgebietes mit einem Mindestgefälle von 2%,

b. das Minimum an in Hangrichtung zugeordneter Rasterflächen mit Mindestgefälle von 2% entspricht dem festgelegten Wert bzw. mindestens zwei Rasterflächen,

c. werden im Laufe der Kaskadierungsprozedur weitere Zuflüsse aus Nachbarrasterflächen identifiziert, so werden die Rasterflächen, die ihnen vorgelagert und ein Mindestgefälle von 2% überschreiten, dem Hang zugeordnet.

Treffen weder Senken- noch Hangbedingung zu, werden die jeweiligen Rasterzellen als "sonstige Flächen" ausgewiesen. Im allgemeinen lautet die Kennung der unterschiedlichen Relieffunktionen: 2 = Senke, 1 = Hang, 0 = sonstige Fläche. Können in einem Einzugsgebiet keine Hänge identifiziert werden, fallen die Senkenkennungen weg, was den Vorteil hat, daß bei der späteren Verschneidung überflüssige Flächen von vornherein aus der Berechnung ausgeschlossen sind.

Innerhalb des Programms wird zudem die Anzahl der Rasterflächen, die als vorgelagerte Einzugsgebiete jeder einzelnen Rasterzelle zuzuordnen sind, berechnet.

4.5 Bestimmung des Hanglängenfaktors mit Hilfe der modifizierten USLE-Gleichung nach Auerswald et al. (1988)

Zur Bestimmung und mathematischen Beschreibung des Hanglängenfaktors kann im Modell auf die modifizierte USLE-Gleichung nach AUERSWALD et al. (1988) zurückgegriffen werden (Gleichung 4).

(4)

mit:

K = Bodenerodierbarkeit

L = Hanglängenfaktor

S = Hangneigungsfaktor

C = Bedeckungs- und Bewirtschaftungsfaktor

Ki, Ci = K-bzw. C-Faktor des i-ten Teilstückes

Nach dem Aufruf des Programms mit "topnew" muß der Rechts- und Hochpunkt der Anfangskoordinate des zu bearbeitenden Kartenblattes eingegeben werden. Handelt es sich um mehrere Kartenblätter, kann dies in den darauffolgenden Zeilen (Wieviele 5.000er-Blaetter? Wieviele Karten nach rechts? Wieviele Karten hoch?) eingetragen werden, wobei die Summe der angegebenen Kartenblätter als Gesamtkarte angesehen und bearbeitet wird. Mit der Eingabe des Höhendatenfiles (Variablenname "hoehed"), welcher in einfachen Anführungszeichen stehen muß, werden anschließend die Höhendaten-Datei sowie die Koordinaten-Datei (Variablenname "koord") eingelesen. Zu beachten ist, daß der Kartenname höchstens 6 Zeichen umfassen darf. Ein Beispieldatensatz einer Höhendatei ist in der folgenden Tabelle (s. Tab. 1) aufgelistet:

Tab.1: Höhendatei (Rechtswert, Hochwert, Höhe)

Bei Angabe der Mindesthanglänge ist zu beachten, daß das Minimum an in Hangrichtung zugeordneter Rasterflächen mit Mindestgefälle von 2 % dem festgelegten Wert bzw. mindestens zwei Rasterflächen entspricht und somit nur Werte > 2 annehmen kann (vgl. 4.3 b). Nach Eingabe des definierten Senkenmaßes und des Höhendaten-Rundungsfaktores schließt sich die Rechenprozedur an, deren Ergebnis in einer ersten Übersicht Aufschluß über die Gesamtfläche, die Anzahl der Einzugsgebiete sowie den prozentualen Anteil an Nullflächen und Senken gibt.

Das Einlesen des Ausgabefiles erfolgt schließlich mit page oder pg <dateiname> (in diesem Fall: EINZUG) oder über eine Editor (z.B. joe). Die errechneten Parameter sind: Rechtswert, Hochwert, Höhe, Gefälle, Abflußrichtungs- und Teileinzugsgebietsnummer, Relieffunktion und die Anzahl der Entwässerungsraster. Tabelle 2 zeigt ein Beispiel für einen Ausgabefile :

Tab.2: Ausgabefile EINZUG

Wird das Programm mit Erosionsberechnung = 1 gestartet, so werden alle in der Datei aufgeführten Variablen ausgegeben. Anderenfalls beschränkt sich die Ausgabe auf die Variablen-Nummern (Spalte Vab. Nr.) 1 - 9.

Beim Start eines neuen Programmlauf wird eine bereits vorhandene Resultat-Datei (Name: EINZUG) gelöscht und durch eine neue ersetzt. Der Programmaufruf sowie die In- und Outputformate der Parameter bzw. Variablen werden im Anhang beschrieben.

Das Programm TOPNEWP2 wurde in FORTRAN 77 geschrieben, ist rechnerkompilierbar und wird hauptsächlich auf dem UNIX-Betriebssystem eingesetzt. Da das Modell in der Lage ist, mehrere Karten (bis zu 30) gleichzeitig zu bearbeiten, ist anzumerken, daß die Rechenprozeduren generell einen sehr hohen Speicherbedarf innehaben Zur Zeit ist eine maximale Array-Belegung von 1400 *1400 Rasterzellen vorgegeben. Diese kann, wenn erforderlich und wenn die Hard- und Software es zulassen, erweitert werden.

1. Programm-Intro

Aufruf : topnew Stand 6.98

Achtung: bestehende EINZUG-Dateien werden geloescht!

2. Abfrage nach einzulesenden Höhenraster-Dateien

Wieviele 5.000er-Blaetter ?

Wieviele Karten nach rechts ?

1

Wieviele Karten hoch ?

1

File mit Hoehendaten nennen:

rechts: hoch:

'h8298'

3. Abfrage nach Vorgaben für die Hang- und Senkendefinition:

Angabe der Mindesthanglaenge (Faktor *Rasterkantenlaenge): 2

Erklärung: Die Hangidentifizierung erfolgt auf Grundlage der berechneten Gefällewerte. Um die Möglichkeit auszuschließen, viele sehr kurze Hänge, die nur Längen von 1 Raster-Kantenlänge haben, auszuglidern, wird eine Variable "Mindesthanglänge" abgefragt. Diese Angabe bestimmt, ab welcher Länge eine Rastefolge mit entsprechendem Gefälle (> 2%) als Hang ientifiziert wird.

Senkenmaß nennen (vertikale Distanz, z.B. 0.4)

0.4

Erklärung: Wie oben beschrieben, erolgt die Identifizierung in 2 Schritten. Zunächst werden auf Grundlage der Gefälle analyse sogenannte lokale Gefälle-Minima bestimmt, also solche Rasterpunkte für die in alle 8 Hangrichtungen nur ein nach oben gerichtetes Gefälle zu verzeichnen ist. Von diesen Punkten ausgehend werden topographische Einzugsgebiete definiert, daß heißt alle Raster durch eine Kennzahl zusammengefaßt, die entsprechend ihres Gefälles in diesen Punkt münden. Als Senke werden all die Raster gekennzeichnet, deren vertikaler Abstand zum lokalen Minimum das Senkenmaß (Eingabe-Variable) nicht übertrifft.

Rundungsfaktor fuer Hoehendaten nennen

1000= mm, 100= cm, 10= dm

auch Abstimmungen dazw. moeglich

10

Erklärung: Die Höhenangaben werden häufig in einer numerischen Genauigkeit geliefert, die die tatsächliche Meßgenauigkeit beiweitem übersteigt. AUs diesem Grund werden die eingelesenen Höhen auf eine angemessene Genauigkeit gerundet.

Soll mit EROSION-SEDIMENTATION gerechnet werden (0/1)?

0

Nenne allgemeinen K-Faktor (z. B. 0.6):

1

Nenne allgemeinen C-Faktor (z. B. 1):

1

EINGABE ENDE

RECHNUNG

AUSGABE: kurze Übersicht am Bildschirm

z.B.

Gesamtflaeche (ha): 405.0156

Einzugsgebiete : 185

Nullflaechen 20565Proz.: 79.33722

Senken 2225Proz.: 8.583774

File: EINZUG (Format siehe Kap.7)