Pestizidverteilung im Weinbau
Diese vorliegende Seite fasst ausgewählte Inhalte meiner Masterarbeit zum Thema „Raum-zeitliche mathematische Modellbildung der Stoffverteilung für ausgewählte Agrochemikalien im Pfälzer Weinbau“ zusammen. Die einzelnen Modellierungsschritte werden hier frei zur Verfügung gestellt, damit diese Arbeit weiter entwickelt, optimiert oder an ein anderes betrachtetes Gebiet oder Pestizide verwendet werden kann.
Einleitung
[Bearbeiten]In dieser Arbeit wurde mit Hilfe von GNU Octave drei Modellierungszyklen erstellt. Mit Hilfe von RLP 3D wurde ein Weinbaugebiet in Landau in der Pfalz ausgesucht und vermessen. Die Verteilung des Pflanzenschutzmittels wurde mit Hilfe der Finiten-Differenzen-Methode simuliert. Dabei wurden das explizite Eulerverfahren und die Diffusionsgleichung verwendet.
In dieser Programmierung wurden wesentliche Parameter, beispielsweise die Größe des betrachteten Gebietes oder das Verwendete Pflanzenschutzmittel, variabel gehalten, wodurch sie leicht auszutauschen und anzupassen sind. Dies hatte den Vorteil, dass simulierte Effekte für kleinere Gebiete besser erkennbar waren. Zudem wurde die Bezeichnung der Variablen so verwendet, dass sie die Funktion erkennen lässt. Der ganze Code wurde außerdem mit Erläuterungen in Form von Kommentaren (mit % markiert) ergänzt, damit es leichter nachvollziehbar ist.
Weitere Optimierungsmöglichkeiten sind im Kapitel Ausblick zusammengefasst.
Erster Modellierungszyklus
[Bearbeiten]In diesem Modellierungszyklus wird das Verspritzen eines Pflanzenschutzmittels in der Mitte des Weinbaugebietes in einem 14-tägigen Rhythmus simuliert und die räumlich-zeitliche Ausbreitung untersucht. Es wird ein rechteckiges Weinbaugebiet mit einer Breite von 140 m und einer Länge von 70 m gewählt. Damit ein äquidistantes Gitter für die Finite-Differenzen-Methode entsteht, wird eine konstante Schrittweite von 2 m gewählt. Jeder betrachtete Gitterpunkt steht hierbei für eine Weinrebe. Anschließend wird aus den so konstruierten Punktefolgen (B=1, 3,…, 140) und (L=1, 3,…,70) ein Gitter mit dem Befehl meshgrid erstellt. Als letzter Schritt wird der zentrale Punkt in der Mitte des Weinbaugebietes mit den Koordinaten (xcenter, ycenter)=(35,18) berechnet. Durch den Befehl round, wird die Division auf die nächst größerer ganzer Zahl gerundet.
%Anfangsbedingungen
%Gitter aufstellen, Anzahl der Punkte = Anzahl der Weinreben im Gebiet
Breite=140; %Breite des Gitters in Meter
Laenge=70; %Länge des Gitters in Meter
h1=2; %Schrittweite in x-Richtung, Abstand Reben 2m
h2=2; %Schrittweite in y-Richtung, Abstand Reihen 2m
B=[1:h1:Breite]; %Gitterindizes in x-Richtung sind es 70 Punkte, die einen Abstand von 2m haben
L=[1:h2:Laenge]; %Gitterindizes in y-Richtung sind es 35 Punkte, die einen Abstand von 2m haben
[x,y]=meshgrid(B,L); %Aufbau eines Gitters
Bx=length(B); %Anzahl der Gitterpunkte der Breite in x-Richtung
Ly=length(L); %Anzahl der Gitterpunkte der Länge in y-Richtung
xcenter=round(Bx/2); %x- und y-Koordinate des zentralen Punktes in der Mitte des Gitters
ycenter=round(Ly/2);
Der Diffusionskoeffizient gibt an wie gut sich ein Stoff verbreitet, der von Kupfer beträgt 0,733*10-9 m2/s. Zudem werden in dieser Modellierung 120 Tage betrachtet, da in einem Zeitraum von 4 Monaten von Mai bis August die Pestizide auf die Weinreben aufgetragen wird.
%konstanter Diffusionskoeffizient
D=0.733*(10^(-9)); %von Kupfer m^2/s
%D=0.731*(10^(-9)); %von Schwefel m^2/s
%Definition Zeit
T=120; %betrachtetes Zeitintervall in Tage
deltat=2; %alle zwei Tage betrachten
Für die Finiten-Differenzen-Methode wird eine Systemmatrix A benötigt, die Bx*Ly Zeilen und Spalten hat. Auf der Hauptdiagonalen sind -4 und deren direkten nebendiagonalen sind Einsen zu finden. Zudem sind auf den Nebendiagonalen noch zusätzlich Einheitsmatrizen eingefügt. Die Struktur dieser Systemmatrix kommt vom zentralen Differenzenquotienten. Um die erforderliche Blockmatrix Struktur zu erhalten, werden einige Einträge der Matrix noch Null gesetzt. Danach wird die Systemmatrix mit dem Diffusionskoeffizienten multipliziert und durch die beiden Schrittweiten dividiert. Im letzten Schritt wird die resultierende Systemmatrix A geplottet.
%Systemmatrix A
V=ones(Bx*Ly,1); %erzeugt einen Spaltenvektor
A=diag(-4.*V,0)+diag(ones(Bx*Ly-1,1),1)+ diag(ones(Bx*Ly-1,1),-1)+ diag(ones(Ly*(Bx-1),1),Ly)+ diag(ones(Ly*(Bx-1),1),-Ly);
%erste drei Terme geben die Matrix B an, letzten beiden die Einheitsmatrizen I
%Differenzenquotienten als Annäherung für u mit FDM (Sternchenschema in 2D)
%Größe Bx^2*Ly^2 mit -4 auf der Diagonale und 1 auf der Nebendiagonale und 1 auf den Bx und -Bx-ten Nebendiagonale
%Korrektur der Matrix (Blockdiagonalität)
for i=1:(Ly-1)
A(i*Bx+1,i*Bx)=0; %Nuller einfügen auf der unteren Nebendiagonale von B
A(i*Bx,i*Bx+1)=0; %Nuller einfügen auf der oberen Nebendiagonale von B
endfor
%Matrix mit D/h^2 multiplizieren
A=A*D/(h1*h2);
figure(1)
spy(A)
xlabel('Spalten')
ylabel('Zeilen')
zlabel('MatrixA')
title('Systemmatrix A')
Zur besseren Übersicht der Struktur wird ein kleineres Gebiet betrachtet. Zudem wird die betrachtete Matrix mit konkreten Einträgen ausgegeben.
Im nächsten Schritt werden zwei Funktionen programmiert. Die erste Funktion (Quellfunktion) prüft mittels einer if-Abfrage die Bedingungen (Euklid-Norm der Differenz des Punktes (px, py) zum zentralen Punkt (xcenter, ycenter) ist kleiner-gleich eins, Zeitpunkt ist Vielfaches von 14 über modulo 14) und bestimmt damit, ob das Pflanzenschutzmittel mit einem Wert von 0,3675 aufgetragen wird. Dieser Wert wurde über die Maximalkonzentration von Kupfer bestimmt, der bei 3 kg pro ha pro Jahr liegt.
%Quellfunktion gibt an zu welchem Zeitpunkt Pestizide im Zentrum xcenter, ycenter hinzukommt
function q=QuelleZeitpunkt(py,px,ycenter,xcenter,t)
if ((sqrt((px-xcenter)^2+(py-ycenter)^2)<=1) && (mod(t,14)==0))
q=0.3675;
else q=0;
endif
endfunction
Die zweite Funktion (Ausbreitungsfunktion fDiff) simuliert über das explizite Eulerverfahren die Ausbreitung des Pflanzenschutzmittels auf dem gesamten Gebiet. Hier wird zunächst aus einer Matrix fneu, auf die später näher eingegangen wird, mit dem Befehl reshape ein Spaltenvektor ualt konstruiert. Das bedeutet, dass spaltenweise vorgegangen wird und jede Spalte der Matrix untereinander in einen Vektor überführt wird. Anschließend wird durch das explizite Eulerverfahren aus dem alten Spaltenvektor ein neuer Spaltenvektor u, indem nach Rechenvorschrift ∆t mit der Systemmatrix und dem alten Spaltenvektor multipliziert wird. Im letzten Schritt wird aus dem Spaltenvektor wieder eine Matrix konstruiert, die dieselbe Größe wie die Ausgangsmatrix hat.
%Explizites Eulerverfahren
function fDiff=Ausbreitung(fneu, Ly, Bx, deltat, A); %fneu ist eine Matrix
u_alt=reshape(fneu,Ly*Bx,1); %aus der Matrix f_neu wird ein Spaltenvektor u_alt
u=u_alt+deltat*(A*u_alt); %Spaltenvektor u
fDiff=reshape(u,Ly,Bx); %mit reshape wird fDiff wieder zu einer LyxBx Matrix
endfunction
Am Ende dieses Modellierungszyklus werden nun alle definierten Parameter und Funktionen in einer for-Schleife zusammengetragen. Die Ausgangssituation ist ein Feld ohne Pestizide, also ist falt ist die Null-Matrix (zeros). Es wird nun eine große for-Schleife programmiert in der alle Zeitschritte von 0 bis 120 in Zweierschritte abgelaufen werden. Danach werden alle Gitterpunkte abgelaufen und zu dem betrachteten Punkt in der Matrix die Pestizide der zuvor programmierten Quellfunktion hinzu addiert, falls die Bedingung in der Funktion, wie zuvor beschrieben, erfüllt ist. So entsteht die neue Matrix fneu. Anschließend wird auf diese Matrix das explizite Eulerverfahren angewandt und es entsteht die neue Matrix fDiff. In einer if-Schleife wird die so entstandene Matrix geplottet, falls der betrachtete Zeitpunkt ein Vielfaches von 14 ist. Danach werden wieder alle Gitterpunkte durchlaufen, um die berechnete Matrix fDiff auf die Matrix falt zu übertragen und der gesamte Prozess beginnt von vorne bis der letzte Zeitpunkt T erreicht ist.
falt=zeros(Ly,Bx); %am Anfang ist noch keine Pestizide
for t=[0:deltat:T] %Zeitschritte
for ix=1:Bx
for iy=1:Ly
fneu(iy,ix)=falt(iy,ix)+QuelleZeitpunkt(iy,ix,ycenter,xcenter,t);
endfor
endfor
fDiff=Ausbreitung(fneu, Ly, Bx, t, A); %Explizites Eulerverfahren
%alle 14 Tage wird ein Bild geplottet
if (mod(t,14)==0)
figure(t+2);
surfc(x,y,fDiff, 'EdgeColor', 'none')
colormap ('jet')
colorbar
axis([1 Breite 1 Laenge 0 4])
title (['Loesung in t=', num2str(t)]);
xlabel('x')
ylabel('y')
zlabel('Fungizidkonzentration')
endif
for ix=1:Bx %alle Gitterpunkte werden durchlaufen
for iy=1:Ly
falt(iy,ix)=fDiff(iy,ix); %Übergabe der Matrix
endfor
endfor
endfor
Die so entstandenen Plots zeigen die Auftragung der Pestizide alle 14 Tage in der Mitte des Weinbaugebiets. Sie wurden im nachfolgenden Video zusammengefasst.
An den rechten Skalen ist zu erkennen, dass die Pflanzenschutzmittel im Laufe der Zeitschritte in der Mitte hinzukommen. Zum Zeitpunkt t=0 ist der Zylinder schwer zu erkennen. Es wurde jedoch bewusst eine einheitliche Skalierung gewählt, um die Abbildungen untereinander besser vergleichen zu können. Zwischen dem ersten und letzten Bild ist jedoch kaum eine Ausbreitung zu den umliegenden Gitterpunkten zu erkennen, sondern lediglich eine Erhöhung der Fungizidkonzentration. Dies hängt mit dem kleinen Diffusionskoeffizienten von Kupfer zusammen. Um nun zu überprüfen, ob wirklich eine Ausbreitung auf die umliegenden Gitterpunkte stattgefunden hat, werden Ausschnitte der Matrizen falt und falt, vor und nach der Verbreitung, angezeigt. Zur besseren Übersicht wird wieder eine 3x5 Matrix (Breite=10, Laenge=6) verwendet.
In diesen beiden Abbildungen ist zu erkennen, dass sich die Pflanzenschutzmittel auf die be-nachbarten Zellen durch die Funktion fDiff ausbreiten. Des Weiteren ist zu sehen, dass im Zeitpunkt t=28 Pestizide durch die Quellfunktion in den fünf mittleren Punkte im Sternchenschema hinzukommen. Jedoch sind die Fungizidkonzentrationen im Vergleich zur Mitte sehr gering und deshalb nicht gut in der Verbreitung in den Bildern zuvor zu erkennen.
In diesem Modellierungszyklus ist es gelungen, das Weinbaugebiet durch entsprechende Abmessungen zu modellieren, eine sinnvolle Menge an Pflanzenschutzmittel in einem regelmäßigen Rhythmus aufzutragen und anschließend durch das explizite Eulerverfahren zu verteilen. Allerdings entsprach das Verspritzen in der Mitte des Gebietes noch sehr wenig der Realität. Daher wird im folgenden Zyklus ein räumlich gleichmäßigeres Verspritzen des Pflanzenschutzmittels angestrebt.
Zweiter Modellierungszyklus
[Bearbeiten]In diesem Modellierungszyklus wird auf dem gleichen Weinbaugebiet das Pflanzenschutzmittel durch Simulation eines Traktors gleichmäßig auf die Gitterpunkte aufgetragen. Somit bleiben die Anfangsbedingungen des ersten Modellierungszyklus, bis auf xcenter und ycenter, die hier nicht benötigt werden, bestehen. Auch die Systemmatrix wird wieder verwendet. Zunächst muss die Pestizide in jedem Gitterpunkt ausgerechnet werden.
%im ganzen Gebiet werden pro Spritzung 0,3675 kg ausgetragen
P=(0.3675/(Ly*Bx)) %pro Gitterpunkt wird P Pestizide verspritzt
Anschließend wird eine Spraymatrix definiert. Diese 3x2 Matrix gibt den Anteil des Pflanzenschutzmittels in Prozent an, der vom Traktor aus auf die umliegenden Weinreben (Gitterpunkte) gespritzt wird. Da das Gitter mit Schrittweite zwei gewählt wurde, liegen die Weinreben auf den ungeraden Gitterpunkten. Somit fährt der Traktor quasi entlang der geraden Gitterpunkte und spritzt dabei seitlich auf die ungeraden Gitterpunkte, die Weinreben. Bei der prozentualen Verteilung wurde angenommen, dass der Traktor in der Mitte links und rechts jeweils 40 % verspritzt und an den Eckpunkten noch jeweils 5 % der Pestizide die Weinreben erreicht. In GNU Octave wurde diese Spraymatrix so konstruiert, dass ein Eintrag der Nullmatrix durch die jeweilig gewünschte Prozentzahl ersetzt wird. Hier muss wieder beachtet werden, dass die x-Koordinate der Punkte des Weinbaugebietes der Spaltennummer und die y-Koordinate der Zeilennummer der Spraymatrix entspricht.
%Spraymatrix Verteilung der Pestizide in einer 2x3 Matrix
Spraymatrix=zeros(3,2);
Spraymatrix(2,1)=0.4; %Angaben in Prozent
Spraymatrix(2,2)=0.4;
Spraymatrix(1,1)=0.05;
Spraymatrix(3,1)=0.05;
Spraymatrix(1,2)=0.05;
Spraymatrix(3,2)=0.05;
Spraymatrix;
Im nächsten Schritt wird mit Hilfe einer for-Schleife die zuvor konstruierte 3x2 Matrix durch das gesamte Weinbaugebiet durchlaufen. Dazu wird eine Nullmatrix der Größe des betrachteten Gebietes mit dem Befehl zeros für die Pestizid-Matrix und für GP (Gitterpunkt) erzeugt. In der for-Schleife wird nun jeder Gitterpunkt in y- Richtung (vy) und in x-Richtung (vx) durchlaufen. Dabei ist es wichtig, dass Start-und Endpunkt innerhalb des Gebietes liegen. Um festzulegen, welche Punkte des Gitters pro Schritt des Traktors einen prozentualen Anteildurch die Spraymatrix erhalten, wird jeweils ein Gitterpunkt (vy, vx) als Bezugspunkt mit dem Eintrag (2,1) der Spraymatrix identifiziert. Anschaulich liegt die Spraymatrix beim jeweiligen Schritt innerhalb des Gitters, wodurch sechs Gitterpunkte durch die Spraymatrix einen Anteil des Pflanzenschutzmittels erhalten. Die übrigen Gitterpunkte, die in diesem Schritt nicht erfasst wurden, den Anteil 0 %. Der Bezugspunkt (vy, vx) wandert nun mit der Spraymatrix das gesamte Gitter entlang und bildet so für jeden Schritt eine nächste GP-Matrix. Dabei ist darauf zu achten, dass nur solche Gitterpunkte als Bezugspunkte fungieren, dass sich die Spraymatrix nie außerhalb des Gitters befindet, da nur Gitterpunkte gespritzt werden sollen. Da die GP-Matrix nach ihrer Entstehung nur Prozentwerte als Einträge besitzen, werden sie mit dem berechneten Pestizidwert P (Abbildung 26) multipliziert, so entstehen Matrizen GPneu, die nun tatsächliche Pestizidwerte enthalten. Diese werden nach jedem Schritt zur Pestizid-Matrix addiert, wodurch diese am Ende der Schleife die gesamte Pestizidverteilung für einen Traktordurchlauf repräsentiert.
%for-Schleife zum Erstellen einer Matrix für die Verspritzung von Pestizide an einem Tag
PestizidMatrix=zeros(Ly,Bx); %Gitterpunkte Ausgangslage
for vy=[1+1:1:Ly-1]; %Ränder wurden ausgenommen
for vx=[1:1:Bx-1];
%der betrachtete Gitterpunkt (GP) wird durch den Punkt der Spraymatrix ersetzt
GP=zeros(Ly,Bx); %zu Beginn sind alle Gitterpunkte Null
GP(vy,vx)=Spraymatrix(2,1); %links Mitte ist mein Bezugspunkt
GP(vy+1,vx)=Spraymatrix(3,1);
GP(vy-1,vx)=Spraymatrix(1,1);
GP(vy,vx+1)=Spraymatrix(2,2);
GP(vy-1,vx+1)=Spraymatrix(1,2);
GP(vy+1,vx+1)=Spraymatrix(3,2);
GPneu=P.*GP; %berechnete Pestizide in jedem Punkt verteilen
PestizidMatrix=PestizidMatrix+GPneu;
endfor
endfor
PestizidMatrix;
Um den erläuterten Entstehungsprozess der Pestizid-Matrix nochmals zu visualisieren, sind in der nachfolgenden Abbildung die Matrizen der einzelnen Entstehungsschritte zu sehen. Dafür wurde wieder ein kleineres Gebiet gewählt, um die entstehende Struktur übersichtlich darstellen zu können. An den Werten lässt sich außerdem ablesen, dass die einzelnen Gitterpunkte eine unterschiedliche Menge an Pflanzenschutzmittel durch den Sprayprozess erhalten. Die Stellen am Rand erhalten mehr Pflanzenschutzmittel als die im Zentrum des Gebietes und die Eckpunkte erhalten den geringsten Anteil. Dies liegt an der Konstruktion der Spraymatrix, durch die der simulierte Traktor immer auf zwei Seiten hin versprüht und somit die Stellen in der Mitte häufiger besprüht werden. Der dargstellte Wert für die Konzentration des Pflanzenschutzmittels lässt sich nach der obigen Formel berechnen.
Wie schon im Modellierungszyklus zuvor wird über den Befehl function das explizite Eulerverfahren definiert. Es wird nun durch eine weitere Funktion Traktor-Matrix-Zeitpunkt miteinbezogen, dass das Verspritzen durch die Pestizid-Matrix nur alle 14 Tage geschieht. Die Variable ist gleich der Pestizid-Matrix, falls der betrachtete Zeitpunkt ein Vielfaches von 14 (Zeile 81) ist und ansonsten ist sie null.
%in dieser Funktion soll der Traktor alle 14 Tage die zuvor erstellte PestizidMatrix auftragen
function u=TraktorMatrixZeitpunkt(t,PestizidMatrix);
if(mod(t,14)==0); %Pestizide wird nur alle 14 Tage gespritzt
u=PestizidMatrix;
else u=0;
endif
endfunction
Analog zum letzten Modellierungszyklus werden am Ende des Modellierungszyklus wesentliche Bestandteile des Codes in for-Schleife zusammengetragen. Dazu wird zunächst eine Nullermatrix falt konstruiert. Die Schleife durchläuft nun alle Zeitpunkte von 0 bis T im Abstand von ∆t . Da nun eine komplette Pestizid-Matrix erstellt wurde, die alle 14 Tage das Verspritzen simuliert, also zur vorherigen Matrix addiert wird, müssen in dieser for-Schleife, im Gegensatz zum vorherigen Modellierungszyklus, nicht alle Gitterpunkte durchlaufen werden. In Zeile 98 wird deshalb pro Zeitpunkt die Traktor-Matrix-Zeitpunkt dazu addiert, also die Funktion zu dem betrachteten Zeitpunkt ausgewertet. Dadurch wird garantiert, dass nur in den Zeitpunkten, die ein Vielfaches von 14 sind, neue Pestizide dazukommt. In Zeile 100 wird dann das explizite Eulerverfahren auf die neu konstruierte Matrix fneu angewandt, um wieder die Pestizide zu verteilen. Im letzten Schritt wird eine if-Schleife konstruiert, die alle 14 Tage die Matrix plottet. In Zeile 113 wird anschließend die neue Matrix der alten übergeben.
falt=zeros(Ly,Bx); %am Anfang ist noch keine Pestizide
for t=[0:deltat:T] %Zeitschritte
fneu=falt.+TraktorMatrixZeitpunkt(t,PestizidMatrix);
%nun nicht alle Punkte wie zuvor durchlaufen, da es zwei Matrizen sind, die addiert werden
fDiff=Ausbreitung(fneu, Ly, Bx, deltat, A); %Explizites Eulerverfahren
%alle 14 Tage wird ein Bild geplottet
if (mod(t,14)==0)
figure(t+2);
surfc(x,y,fDiff, 'EdgeColor', 'none')
colormap ('jet')
colorbar
axis([1 Breite 1 Laenge 0 0.0014])
title (['Loesung in t=', num2str(t)]);
xlabel('x')
ylabel('y')
zlabel('Fungizidkonzentration')
endif
falt=fDiff;%neue Matrix wird der alten übergeben
endfor
In der for-Schleife wird in jedem Zeitschritt eine fneu und eine fDiff Matrix vor, bzw. nach der Verteilung der Pestizide, berechnet. Die Ausgabe dieser Matrizen ist in den nachfolgenden Abbildungen für die Zeitpunkte t=14 und t=28 für die übersichtlichere 3x5 Matrix ausgeführt.
In diesen Abbildungen ist zu erkennen, dass die Pestizid-Matrix nach 14 Tagen hinzu addiert wird. Außerdem sind, wie zuvor beschrieben, deutlich die Unterschiede zwischen den inneren Punkten, Rand- und Eckpunkten zu erkennen. Zudem ist weder zum Zeitpunkt t=14 noch bei t=28 ein Unterschied zwischen den Matrizen vor und nach der Verteilung der Pestizide zuerkennen. Dies liegt an dem geringen Diffusionskoeffizienten und an den groß gewählten Gitterpunktabständen. Diese Erkenntnisse sind in den 3D Plots für die übersichtlichere Matrix in dem nachfolgenden Video leichter zu erkennen. Zudem sind die Ergebnisse für das betrachtete Weinbaugebiet auch in Video zusammengefasst.
An der Skalierung auf der rechten Seite der Plots ist gut zu erkennen, dass die Pestizidkonzentration mit der Zeit zunimmt. Da die inneren Punkte eine höhere Pestizidkonzentration besitzen, ist dieser Plot konsistent zu den Matrizen davor. Die größere Matrix besitzt mehr innere Punkte und somit hat der Plot keine spitze Form sondern ist nahezu ebenerdig erhöht. Auch hier ist zu erkennen, dass die Pestizidkonzentration zwar steigt, jedoch sehr wenig verteilt wird. Des Weiteren ist an der Skalierung zu erkennen, dass sie viel geringer ist als die Skalierung bei den kleineren Matrizen. Dies liegt daran, dass nun mehr Gitterpunkte vorhanden sind und die Pestizidkonzentration pro Gitterpunkt mit P= 0.00015 deutlich kleiner ist als zuvor.
Insgesamt entspricht der zweite Modellierungszyklus eher der Realität als der erste, da nun ein Traktor modelliert wurde, der durch die einzelnen Spalten fährt und rechts und links jeweils die Pestizide auf die Weinreben verspritzt. Ein Nachteil dieser Modellierung ist, dass die Pestizide zwar alle 14 Tage neu aufgetragen wird, sich die Pestizide aber nicht wie gewünscht verteilt. Jedoch hat dieses Modell auch noch Mängel, die optimiert werden könnten. Zum Beispiel wird im nächsten Modellierungszyklus eine Randbedingung betrachtet, da auch Pestizide nach dem Aufspritzen über den Rand hinweg diffundieren können. Zudem sind von außen wirkende Faktor, die dazu beitragen, dass die Pestizide in dem betrachteten Gebiet abnimmt, noch nicht näher betrachtet.
Dritter Modellierungszyklus
[Bearbeiten]In diesem Modellierungszyklus wird zunächst die Randbedingung näher betrachtet, um anschließend auf den prozentualen Anteil einzugehen, den vom Boden aufgenommen wird. Dazu wird die resultierende Matrix des zweiten Modellierungszyklus verwendet, weswegen zunächst sämtliche Befehle des vorherigen Zyklus übernommen werden.
Im ersten Teil dieser Modellierung soll überlegt werden, wie viel Pestizide über den Rand hinaus fließt und ob die an den Rändern gleichmäßig geschieht. Dies wird wie im zweiten Modellierungszyklus über einen prozentualen Anteil in einer Matrix berechnet. Dazu wird die zuvor übernommene kleinere Matrix nach 30 Tagen in eine größere Matrix überführt, sodass zu Beginn der Rand aus Nullen besteht.
%für die Randbedingungen wird eine Matrix erstellt,
%die um zwei Reihen und zwei Spalten größer ist als die letzte Matrix
Rand=zeros(Ly+2,Bx+2);
for n=[1:1:Bx]
Rand(:,(n+1))=[0; fDiff(:,n); 0];
endfor
Rand;
Um nun zu untersuchen, welcher Anteil an Pestizide über welchen Rand hinweg transportiert wird, wird das Höhenprofil, das mit RLP 3D (Landesamt für Vermessung und Geobasisin-formation Rheinland Pfalz) zu dem betrachteten Weinbaugebiet verwendet wurde, näher betrachtet. Zu erkennen ist, dass der Startpunkt mit 174,7 Höhenmetern der höchste Punkt ist und der niedrigste Eckpunkt diagonal gegenüber liegt und auf einer Höhe von 168,3 Höhenmetern liegt. Der markierte dritte Punkt liegt bei 172,1 Höhenmetern. Da das betrachtete Gebiet in einer Ebene liegen soll, wird angenommen, dass der vierte Eck-punkt bei 170,9 Höhenmetern liegt. Die lange Kante des Rechtecks beträgt 0,14 km und die kurz 0,07 km.
Durch die Lage der Punkte ist zu erkennen, dass ein Höhengefälle von links oben nach rechts unten diagonal verläuft. Da es sich um eine ebene Fläche handelt, wird der Verlauf des Gefälle von links nach rechts und von oben nach unten berechnet. Dazu wird eine verschobene Sinus-Kurve verwendet. Es wird der Winkel zwischen einer oberen Kante und der Senkrechten zum Boden untersucht. Der gesuchte obere Winkel liegt zwischen 0° und 360°. Es werden zwei Steigungsdreiecke betrachtet.
Über den Arkuskosinus ergibt sich ein Winkel von 88,44° des linken Steigungsdreiecks (Gefälle von links nach rechts) und ein Winkel von 87,87° des rechten Steigungdreiecks (Gefälle von oben nach unten).
Im nächsten Schritt wird mit GNU Octave eine verschobene Sinus-Kurve erstellt, mit der jeder möglichen Steigung ein zugehöriger prozentualer Anteil an auszutre-tenden Pestiziden über den Rand hinaus zugeordnet werden kann. Mit dieser Zuordnung wird erreicht, dass bei einem exakt waagerechten Gebiet genau 50 % der Pestizide über den Rand austreten. Existiert dagegen ein Gefälle mit einem gewissen Steigungswinkel, tritt über den höher liegenden Rand ein kleiner Anteil aus als über den niedriger gegenüberliegendem Rand. Die x-Koordinate gibt dabei den Winkel zwischen 0° und 360° an, der zuvor berechnet wurde. Die y-Koordinate verläuft zwischen 0 und 1. Die Schrittweite der y-Koordinate (Zeile 5) wurde so gewählt, dass beide Vektoren gleich groß sind, damit die Funktion später geplottet werden kann.
Anschließend wird eine Funktion s definiert, die die verschobene Sinuskurve beschreibt. Die-se Funktion wird benötigt, um später in GNU Octave in der Befehlsleiste die Winkel in die Funktion einzusetzen. In Zeile 14 wurde zudem dieselbe Funktion nun als f definiert, damit der anschließende Plot problemlos verläuft.
%Variablen definieren
x=[0:1:360]; %in x-Richtung von 0 bis 360 mit Schrittweite 1
y=[0:(1/360):1]; %in y-Richtung von 0 bis 1
%Schrittweite so gewählt, dass x und y gleich lange Vektoren sind
%Funktion s definieren, um später das Gefälle x pro Rand zu berechnen
function s=Sinus(x);
s=(sin((x-90)*2*(pi/360))+1)*0.5; %Sinuskurve zur Bestimmung des Gefälles
endfunction
%Funktion f wie s definieren, damit man sie später plotten kann
f=(sin((x-90)*2*(pi/360))+1)*0.5;
%Funktion plotten
plot(f)
xlabel('x')
ylabel('y')
title('Ausschnitt einer Sinuskurve')
Die zuvor berechneten Winkel der Steigungsdreiecke können nun in die Sinusfunktion einge-setzt werden. Dadurch entsteht ein berechneter prozentualer Anteil von 48,26 % von links nach rechts und ein prozentualer Anteil von 47,38 % von oben nach unten.
Da die Sinuskurve so erstellt wurde, dass immer der obere Winkel betrachtet wurde, entsprechen die prozentual berechneten Werte dem linken Rand (48,26 %), bzw. dem oberen Rand (47,38 %). Daraus ergibt sich nun, dass am rechten Rand 51,74 % und am unteren Rand 52,62 % der vorhandenen Pestizide heraustreten.
Da die prozentualen Gefälle in beiden Richtungen nicht gleich sind, müssen in der folgenden Programmierung ein Unterschied in der betrachtete Reihenfolge der Zeilen und Spalten un-terschieden werden. Zunächst werden die Spalten und anschließend die Zeilen betrachtet, so wird zuerst die Pestizide links und rechts und anschließend oben und unten über den Rand hinweg verteilt. In der for-Schleife werden nun alle Zeilen abgelaufen. In der ersten Spalte der Rand-Matrix (entspricht dem linken Rand) werden die Nullen durch den prozentualen Anteil von 48,26 % des Eintrags in derselben Zeile und in der zweiten Spalte ersetzt. In Zeile 115 wird dann die zweite Spalte von der zuvor konstruierten ersten Spalte abgezogen. In Zeile 117 und 118 wird dasselbe für den rechten Rand, also für die letzte und vorletzte Zeile berechnet.
Analog werden in einer zweiten for-Schleife alle Spalten durchlaufen und zunächst der oberen Rand mit 47,38 % in der ersten Spalte und anschließend der untere Rand mit 52,62 % in der letzten Zeile berechnet.
%Matrizen erstellen, die die Verteilung über den Rand modellieren
%berechneter Anteil der Pestizide: linken Rand 48,26% und rechter Rand 51,74%
for i=[1:1:Ly]
%erste und zweite Spalte
Rand((i+1),1)=0.4826*Rand((i+1),2); %linker Rand 48,26%
Rand((i+1),2)=Rand((i+1),2)-Rand((i+1),1);
%vorletzte und letzte Spalte
Rand((i+1),(Bx+2))=0.5174*Rand((i+1),(Bx+1)); %rechter Rand 51,74%
Rand((i+1),(Bx+1))=Rand((i+1),(Bx+1))-Rand((i+1),(Bx+2));
endfor
Rand;
%berechneter Anteil der Pestizide: oberer Rand 47,38% und unterer Rand 52,62%
for i=[1:1:Bx]
%erste und zweite Zeile
Rand(1,(i+1))=0.4738*Rand(2,(i+1)); %oberer Rand 47,38%
Rand(2,(i+1))=Rand(2,(i+1))-Rand(1,(i+1));
%vorletzte und letzte Zeile
Rand((Ly+2),(i+1))=0.5262*Rand((Ly+1),(i+1)); %unterer Rand 52,62%
Rand((Ly+1),(i+1))=Rand((Ly+1),(i+1))-Rand((Ly+2),(i+1));
endfor
Rand
Um die erstellte Rand-Matrix zu plotten, wird zunächst ein größeres Gitter analog zu den Anfangsbedingungen in jedem der Modellierungszyklen konstruiert.
BreiteR=Breite+2*h1; %Breite des Gitters in Meter
LaengeR=Laenge+2*h2; %Länge des Gitters in Meter
BR=[1:h1:BreiteR]; %Gitterindizes in x-Richtung sind es 70 Punkte, die einen Abstand von 2m haben
LR=[1:h2:LaengeR]; %Gitterindizes in y-Richtung sind es 35 Punkte, die einen Abstand von 2m haben
[xR,yR]=meshgrid(BR,LR); %Aufbau eine Gitters
figure(1);
surfc(xR,yR,Rand, 'EdgeColor', 'none')
colormap ('jet')
colorbar
axis([1 BreiteR 1 LaengeR 0 0.0014])
title ('Matrix mit Rand');
xlabel('x')
ylabel('y')
zlabel('Fungizidkonzentration')
Für die Vertauschung der Spalten-Zeilen-Verteilung wird analog zu dem Code zuvor die Rand1-Matrix konstruiert. Sie wird anders bezeichnet, damit es später nicht zu einer Überschreibung der Matrix kommt, aber sie entspricht zu Beginn der Rand-Matrix. Anschließend müssen die beiden for-Schleifen vertauscht und mit Hilfe der Befehle analog zu oben geplottet werden. So entstehen am Ende die beiden 5x7 Rand- und Rand1-Matrizen nach 30 Tagen. Werden jeweils die beiden „Eckpunkte“ links oben, bzw. rechts unten addiert, ist das gewünschte, leichte Gefälle in der Diagonalen erkennbar. Zudem unterschieden sich die beiden Matrizen leicht in den berechneten Werten.
In den beiden nachfolgenden Abbildung sind die geplotteten Pestizidkonzentrationen mit der Rand-Matrix, bzw. der Rand1-Matrix zu sehen. Ein Vergleich der beiden Abbildungen wird ein minimaler Unterschied an den Eckpunkten deutlich.
Auch wenn die kleine 5x7 Matrix nach 30 Tagen geplottet wird, ist nur ein kleiner Unterschied zwischen der Rand- Matrix und der Rand1-Matrix zu sehen.
Im zweiten Teil des Modellierungszyklus wird der prozentuale Anteil, der in den Boden versickert subtrahiert. Dazu wird angenommen, dass der Boden 10 % der Pestizide aufnehmen kann. In der nachfolgenden Abbildung wird die Endmatrix, bzw. Endmatrix1, durch Subtraktion des prozentualen Anteils der Rand-Matrix, bzw. Rand1-Matrix bestimmt und anschließend geplottet.
%Annahme: der Boden nimmt 10% der Pestizide auf dem Weinbaugebiet auf
a=0.1;
Endmatrix=Rand-a.*Rand
Endmatrix1=Rand1-a.*Rand1
figure(3);
surfc(xR,yR,Endmatrix, 'EdgeColor', 'none')
colormap ('jet')
colorbar
axis([1 BreiteR 1 LaengeR 0 0.0014])
title ('Endmatrix mit Rand');
xlabel('x')
ylabel('y')
zlabel('Fungizidkonzentration')
figure(4);
surfc(xR,yR,Endmatrix1, 'EdgeColor', 'none')
colormap ('jet')
colorbar
axis([1 BreiteR 1 LaengeR 0 0.0014])
title ('Endmatrix mit Rand1');
xlabel('x')
ylabel('y')
zlabel('Fungizidkonzentration')
In der nachfolgenden Abbildung sind die beiden geplotteten 5x7 Endmatrizen zu sehen.
Die endgültigen Matrizen des betrachteten Weinbaugebiets sind in den beiden nachfolgenden zu sehen.
Insgesamt wurden in diesem letztem Modellierungszyklus die Ränder mit Hilfe von Matrizen und einer modifizierten Sinuskurve bestimmt, um am Ende den versickerten Anteil subtra-hiert. Es zeigte sich, dass es ein leicht unterschiedliches Gefälle des betrachteten Gebietes gibt, wodurch zwischen der Reihenfolge der Berechnung der Zeilen und Spalten unterschieden werden musste. In den Plots ist kaum ein Unterschied in den Konzentrationen zu sehen. Zudem sollte noch erwähnt werden, dass Ungenauigkeit der Gefälle durch Rundungsfehler entstanden sind. Falls das Gefälle in beiden Richtungen gleich ist, genügt es, nur einen Fall zu berechnen. Außerdem wurde hier ein Gefälle von links oben nach rechts unten betrachtet. Falls das Gebiet jedoch eine andere Richtung des Gefälles aufweist, müssen die Rechnungen der Steigungsdreiecke angepasst werden.
Ausblick
[Bearbeiten]Die vorliegenden Modellierungszyklen können weiterhin über folgende Betrachtungen optimiert werden:
- Unterschiede von ökologischen und herkömmlichen Weinbau durch Verwendung anderer Pestiziden
- andersartige Gebiete, wie solche an Steilhängen
- andere Sorte von Traktor oder Hubschrauber, bei denen sich eine andere Verteilung durch die Pestizid-Matrix ergibt
- Wetterlage berücksichtigen, z.B. Regenschauer und somit einen anderen Pestizidverspritzungsrhythmus
- Recyclingspritzen, die die Pestizidausbreitung auf umliegende unerwünschte Gebiete minimieren
- verschiedene Bodentypen
- Übergang der Pestizide vom Weinbau in anliegende Gewässer und somit der Einfluss auf das Grund- und Trinkwassersystem
Quellen
[Bearbeiten]Landesamt für Vermessung und Geobasisinformation Rheinland Pfalz. RLP in 3D. Zugriff am 24.04.2021. Verfügbar unter http://www.rheinland-pfalz-in-3d.rlp.de/
Kalka, H. (2020, 05. Dezember). Tabelle der Diffusionskoeffizienten. Zugriff am 16.06.2021. Verfügbar unter https://www.aqion.de/site/diffusionskoeffizienten