Partikelfilter zur Parameterschätzung in der Hydrogeophysik

Über meine eigene Arbeit habe ich bisher noch wenig geschrieben. Aber ich denke, jetzt wo die EGU in Wien vor der Tür steht, die größte Geowissenschaften-Konferenz in Europa, ist eine gute Gelegenheit mein letztes Paper einmal vorzustellen, wie ich es auch auf der Konferenz tun werde.

Gehen wir doch zuerst mal ins Feld, Daten sammeln. Aber eigentlich ist es mehr ein Labor als das Feld, seht mal das Foto, es zeigt ein Deichmodell, das sich in Karlsruhe an der Bundesanstalt für Wasserbau befindet. Modell heißt aber nicht, dass es klein wäre, es ist ein vollmaßstäblicher Deich von etwa 3.5 m Höhe und 22 m Länge. Man kann auf dem Bild erkennen, dass er in einer Vertiefung aufgebaut ist und auf einer schwarzen, wasserundurchlässigen Folie steht. Auf der linken Seite befinden sich Rohre, mit denen man Wasser in das Becken einfüllen kann, um so einen Einstauversuch durchzuführen. Damit bildet man ein Hochwasser nach, man sieht im Bild wie das Wasser am Deich ansteht. Es wird auch in den Deich eindringen und auf der Landseite, am Fuß des Deiches, in einen Kiesdrän abfließen.

Ich habe das einen Laborversuch genannt. Der Grund ist einerseits, dass man durch den kontrollierten Einstau die Möglichkeit hat, in wenigen Tagen einen Boden von trocken nach nass zu bringen. Das wird noch besonders wichtig, denn das liefert uns später die Information die wir suchen. Außerdem ist der Deich homogen aus einem Sand aufgebaut. Böden sind eigentlich so ziemlich das komplexeste System, das sich finden lässt. Da erleichtert ein homogenes Material erst einmal sehr die Untersuchungen mit neuen Methoden.

Gemessen habe ich an diesem Deich während meiner Doktorarbeit. Ich habe mit “Electrical Resistivity Tomography” (ERT) gearbeitet, was ich hier schon erklärt habe. Nur kurz, mit ERT bestimmt man von der Oberfläche aus die elektrische Widerstandsverteilung im Untergrund. Und diese hängt vom Wassergehalt ab, also lässt sich der Wassergehalt (vielmehr seine Änderung) mit ERT messen und überwachen! Die Ergebnisse vom Deich hatte ich auch schon einmal vorgestellt.

Es gibt aber noch eine weitergehende Fragestellung, nämlich ob man das Verhalten eines Deiches, oder allgemein eines Bodens, auch simulieren kann. Dazu muss man die hydrologischen oder hydraulischen Parameter des Bodens kennen, beispielsweise die (hydraulische) Leitfähigkeit für Wasser, wenn alle Poren des Bodens mit Wasser gefüllt sind. Hier in Jülich befasse ich mich mit gekoppelter hydrogeophysikalischer Inversion, der Schätzung dieser Parameter “aus der Ferne” mit Geophysik, also ohne den Boden aufgraben und ins Labor tragen zu müssen. Das macht den Deich zu so einem interessanten Versuchsobjekt: Da er kontrolliert von trocken bis gesättigt beobachtbar ist, kann man gerade aus der Beobachtung der Wasserfront die hydraulischen Parameter ableiten.

Gekoppelte Inversion

Unser Ansatz dazu versucht, den Schritt über die Bestimmung des Wassergehaltes wegzulassen – da er mit einigen Fehlern und Ungenauigkeiten behaftet sein kann. Das Bild hilft ein wenig, den Weg zu verdeutlichen:

Wir haben ein bißchen Welt, die wir untersuchen. Darauf machen wir ein hydrologisches Experiment, den Einstau, und messen mit einer geophysikalischen Methode. Diese misst “scheinbare Widerstände” und nicht die zweidimensionale Widerstandsverteilung. Letzere müsste man mit einem numerischen Verfahren ermitteln, das aber nicht eindeutig und ggf. quantitativ mit sehr großen Fehlern behaftet ist. Daher kommen wir von der anderen Seite: Wir bauen uns ein Modell des Deiches und simulieren die Reaktion auf den Einstauversuch. Zu verschiedenen Zeitpunkten haben wir eine geophysikalische Messung gemacht, und zu diesen Zeiten ermitteln wir, welches Messergebnis des Computermodell denn geliefert hätte.
Und das geht so: Wir nehmen die Wassergehaltsverteilung, die das hydrologische Modell simuliert. Darauf wenden wir eine petrophysikalische Gleichung an – diese Gleichungen verbinden Wassergehalt und eine geophysikalisch interessante Größe, bei uns der elektrische Widerstand. Aus dem Bild der Verteilung des Wassergehaltes entsteht also ein Bild der Verteilung des el. Widerstandes. Diese Verteilung wiederum nimmt man für ein elektrisches Vorwärtsmodell, das dann berechnet was man bei solch einer Verteilung gemessen hätte.
Dieser ganzer Vorgang ist ein Durchlauf innerhalb einer Optimierung. Nach jedem solchen Schritt vergleicht man die modellierte mit der tatsächlichen Messung. Dann dreht man an den Schrauben der Modelle. Das sind die Parameter, und zwar einerseits die hydraulischen, die uns interessieren, aber auch die der petrophysikalischen Gleichung, die wir brauchen.

Die Anzahl an Parametern wird schnell hoch, wir haben in einer früheren Veröffentlichung 12 bestimmt. In dieser Arbeit haben wir zur Parameterschätzung einen “state-of-the-art” Monte Carlo Markoc Chain (MCMC)-Algorithmus eingesetzt. Dieser hat den schönen Namen DREAM (DiffeRential Evolutionary Adaptive Metropolis). Der Einsatz von MCMC-Methoden trägt vor allem einem Rechnung: Dass wir die Parameter gar nicht auf den Punkt bestimmen können, sondern mit Unsicherheiten arbeiten und stattdessen lieber Wahrscheinlichkeitsverteilungen für die Parameter angeben. Das hat gut geklappt, auch wenn man damit einen PC-Cluster einige Wochen beschäftigen kann.

Sequentielle Optimierung

Zu jedem Optimierungsschritt in dieser Methode gehört es, alle Messungen gleichzeitig zu betrachten, um die Parameter zu optimieren. Dann bekommt man nach vielen Zehntausend Schritten (bei denen also jeweils das hydrologische Modell einmal und das geoelektrische Modell einmal je Messung durchlaufen muss) eine Konvergenz und kann aus einigen Tausend Schritten im konvergierten Zustand die Wahrscheinlichskeitsverteilung ableiten. Vorher gibt man Grenzen für jeden Parameter an, aus denen der Algorithmus dann zufällig fischt – die sogenannte a-priori-Verteilung. Durch den Vergleich mit den Messungen ermittelt man die a-posteriori-Verteilung.
Aber hier ist ein Gedanke – was wäre denn wenn man mit einer solchen ERT-Methode einen Deich permanent überwachen möchte? Dann hätte man gar kein Ende – wann startet man dann den MCMC-Algorithmus? Und wichtiger noch – man wäre doch sehr daran interessiert, nach jeder Messung Informationen über die Unsicherheit in den Parametern, und vor allem auch im Zustand (der Wassergehaltsverteilung) zu bekommen.
Auftritt sequentielles Bayes-Verfahren – Partikelfilter!

Particle Filtering

Filtern stammt eigentlich aus der Signalverarbeitung. Es geht hier darum, den Zustand eines Systems zu beobachten. Der Zustand ist aber nicht direkt beobachtbar (wie bei uns der Wassergehalt IM Deich). Man kann aber eine Messung des Zustandes beobachten und hat ein Modell, das die zeitliche Entwicklung des Zustandes simuliert. Dann passt man bei jeder Messung den Zustand immer so gut wie möglich an – wobei die Hauptaufgabe in der Signalverarbeitung des Herausfiltern des Signals aus dem Rauschen ist. Wir verwenden Filtern hier aber hauptsächlich zur Parameterschätzung, was man so “nebenbei” mitlaufen lassen kann. Der wichtigeste und verbreitetste Filter ist der Kalman-Filter, der in erweiterten Varianten auch in der Hydrologie und Meteorologie fleißig zum Einsatz kommt. Der Kalman-Filter lebt aber in einer linearen Gauss-Welt. Auch die erweiterten Varianten, die nicht-lineare Modelle berücksichtigen können, setzen Gauss-Verteilungen voraus. Das ist nicht so toll für uns – alleine schon der Gedanke, dass unsere Parameterschätzung am Anfang nur eine Ober- und Untergrenze vorgibt. Wir wollen annehmen können, dass der Parameter irgendwo dort drin liegen kann, mit gleicher Wahrscheinlichkeit. Die Lösung bietet sequentielles Bayes-Filtern, oder einfach als Partikelfilter bekannt. Hier sieht man die Grundidee:

In blau ist die Wahrscheinlichkeitsverteilung für einen Parameter angedeutet. Am Anfang wissen wir, dass er irgendwo zwischen zwei Grenzen liegt. Der Filter soll diese Verteilung jetzt durch ständiges Vergleichen mit Messungen zu der “echten” Verteilung umbiegen, die unser Wissen über den Parameter eingebaut hat.
Da wir so eine komplizierte Verteilung schlecht als Funktion angeben können, versuchen wir sie einfach an vielen diskreten Stellen abzugucken, als gute Näherung an die wirkliche Form. Ein Partikel trägt eine Möglichkeit des Zustandes und ein für jeden Parameter. Die Form der Wahrscheinlichkeit drücken wir durch das Gewicht eines Partikels aus – dicke Partikel zeigen an, dass dieses Partikel den Zustand und die Parameter nahe an der Wirklichkeit trägt.

Rennen der Partikel

Nachdem wir die Partikel zufällig aus der Anfangsverteilung mit Wassergehalt und Parametern versorgt haben, tragen die Partikel ein Rennen aus!

Eine “Runde” im Rennen ist quasi der Zeitraum zwischen zwei Messungen – und dann zeigt sich welche Partikel schneller sind. Die “schnelleren”, deren Parameter die Wirklichkeit besser wiedergeben, kommen früher rund und bekommen ein höheres Gewicht.

Das ganze lässt sich auch noch “seriöser” darstellen:

Wir beginnen mit einer Verteilung der Parameter mit gleichem Gewicht – mitten im Bild vor dem dritten Kästchen “Propagation”. Im Schritt der Zeitpropagation stecken wir den Zustand (Wassergehalt) als Anfangszustand in ein hydrologisches Modell. Die Parameter die das Partikel mitbringt, bestimmen dann den Verlauf der Simulation, und angetrieben wird diese durch die veränderten Randbedingungen. Das sind in unserem Fall das Ansteigen und Absinken des Wasserstands im Becken (ihr erinnert euch noch?). Wir rechnen vorwärts bis ein Messzeitpunkt erreicht ist. Wie oben beschrieben, bestimmen wir dann, was wir bei einer solchen Wassergehaltsverteilung gemessen hätten. Ist diese simulierte Messung nahe bei unserer echten Messung, erhält das Partikel ein hohes Gewicht (Jetzt sind wir oben im Bild, bei “Observation”). Wir machen diese ganze Simuliererei für viele Partikel mit verschieden guten Parametern. Wenn wir ein paar Runden dieses Rennens durchführen, werden ein paar das Fett ansammeln, weil sie gute Parameter rundtragen. Damit nicht irgendwann ein paar Partikel alles Gewicht angesaugt haben wie ein Schwarzes Loch, gibt es dann einen “Resampling”-Schritt. Nach diesem haben alle Partikel wieder gleiches Gewicht, sitzen aber an neuen Positionen. Wo vorher ein fettes Partikel saß, hocken jetzt mehrere Partikel, und manch Leichtgewicht fliegt vielleicht völlig raus.

Durchführung

Durchgeführt haben wir dieses Particle Filtering für die Messungen mit ERT; aber zum Vergleich haben wir auch die elektrische Modellierung weggelassen und einen zweiten Filter laufen lassen, der direkt mit (unabhängig) gemessenem Wassergehalt verglichen hat.
Implementiert habe ich den Filter, in dem ich um die bestehenden Modelle (HYDRUS für die Hydrologie und CRMOD für die elektrische Modellierung) drumrum mit Python eine Software gebaut habe, die die Partikel verwaltet und die Modelle parallel auf mehrere Prozessoren verteilt startet. Ich habe ihr den dämlichen Namen PFIFF (Particle filtering inversion-friendly framework) gegeben.
In dieser einfachen Variante eines Partikelfilters, die man SIR (Sampling Importance Resampling) nennt, sitzen nach dem Resampling mehrere Partikel auf der gleichen Stelle. Da das natürlich unnötig ist und man außerdem den Parameterraum durchsuchen möchte, gibt man etwas Gauss-Rauschen drauf (bewegt sie also zufällig ein klein wenig). Diese Methode erlaubt, die besten Parameter zu suchen, ist aber nicht systematisch wie ein Metropolis-Algorithmus und braucht daher sehr viele Partikel. In Zukunft werden wir das verbessern, aber jetzt konnten wir erstmal nur 3 Parameter schätzen und haben auch dafür 4096 Partikel gebraucht und 50 Stunden parallel auf 26 CPUs gerechnet. Der Vergleich mit dem Wassergehalt ging etwas schneller, aber darauf will ich nicht weiter eingehen.
Aber um nochmal den großen Vorteil des Partikelfilters zu betonen – zu jedem Zeitpunkt, an dem eine Messung gemacht wurde bekommen wir eine verbesserte Wahrscheinlichkeitsverteilung in Parameter und Zustand – und nicht erst nachdem alles durch ist wie beim MCMC.

Ergebnisse

Ich beschränke mich auf das Ergebnis der gekoppelten Inversion, und auf den Parameter k^S, das ist die hydraulische Leitfähigkeit bei Sättigung – salopp wie schnell das Wasser durch den Deich fließt.

Im Bild ist auf der y-Achse der Wertbereich (logarithmisch) des Parameters dargestellt. Es kommt also so etwas wie 0,001 m/s heraus. Gezeigt wird im Bild die zeitliche Entwicklung der Wahrscheinlichkeitsverteilung des Parameters. Die blaue Linie zeigt den Median der Verteilung, der blaue und graue Bereich sind Konfidenzintervalle. Das heißt, dass im blauen Bereich 50% aller Partikel um den Median herum liegen, und im grauen 90%. Die grünen Punkte geben einerseits die Messzeitpunkte an und andererseits das gewichtete Mittel der Parameterschätzung.
Klar ist zu erkennen, wie zu der Zeit, bei der das ERT beginnt die steigende Wasserfront im Deich zu sehen (erinnert euch an die Bilder hier), die Breite der Wahrscheinlichkeitsverteilung abnimmt und der Median zu einem Wert hinstrebt. Dieser passt übrigens gut mit Labormesswerten zusammen. Heureka!

Die beiden anderen Parameter waren die petrophysikalischen Parameter, und wurden ebenfalls gut bestimmt. Die gesamten Ergebnisse lassen sich im Paper nachlesen, das übrigens als Open Access Paper von jedem kostenlos angesehen werden kann.

Partikelfilter bieten einige Möglichkeiten, vor allem da man damit den Zustand eines Systems überwachen und immer verbesserte Informationen über die Parameter bekommen kann. Für die Zukunft wollen wir die Partikelfilter verbessern, indem eine systematische Durchsuchung des Parameterraumes eingebaut wird. Außerdem untersuchen wir den Umgang mit zeitlich variablen Parametern.