diff options
Diffstat (limited to 'buch/papers/erdbeben/teil1.tex')
| -rw-r--r-- | buch/papers/erdbeben/teil1.tex | 341 |
1 files changed, 172 insertions, 169 deletions
diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index 014b53e..1a893dd 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -9,187 +9,213 @@ % (c) 2020 Prof Dr Andreas Müller, Hochschule Rapperswil %% - - \rhead{Kalman-Filter} - \section{Kalman-Filter} -Die interessante Grösse ist also das Integral der Überlagerung zweier Kräfte. -Wir brauchen also die zweite Ableitung der Messung, ohne deren Eigendynamik. -Da wir die äussere Kraft nicht direkt messen können, benötigen wir ein Werkzeug, welches aus der gemessenen Position, die Krafteinwirkung auf unsere System schätzt. -Dies ist eine typische Anwendung für das Kalman-Filter. - -Das Filter schätzt den Zustand eines Systems anhand von Messungen und kann den nächsten Zustand errechnen und aus dieser Schätzung auch eine erwartete Messung herleiten. -Die für das Filter relevante Grösse ist dann nicht mehr die eigentliche Messung, sondern die Differenz aus Messung und Erwartung, da diese Differenz, die Innovation, eine Aussage über die nicht-deterministischen, externen Einflüsse auf das System ermöglicht. -Das Filter berücksichtigt dazu nicht nur die Messung und den Zustand, sondern auch die Unsicherheiten dieser beiden Grössen, welche als Parameter in das Modell des Systems einfliessen. - -Unser Ziel ist es, anhand der Messung die eigentlich interessante Grösse $f$ zu bestimmen. -Dabei wird durch eine deterministische Vorhersage, in dem der Zustand mit der Eigendynamik des Systems multipliziert wird. -Die Idee dahinter ist, dass das Kalman-Filter die nicht-deterministische Grösse $f$ anhand der Messung und der Vorhersage zu bestimmen. - -Für mehrere Dimensionen (x,y,z) würde der Satz von Pythagoras für das System benötigt. -Da sich der Pythagoras bekanntlich nicht linear verhält, kann kein lineares Kalman-Filter implementiert werden. -Da das Kalman-Filter besonders effektiv und einfach für lineare Abläufe geeignet ist, würde eine zweidimensionale Betrachtung den Rahmen dieser Arbeit sprengen. -Einfachheitshalber beschränken wir uns auf den linearen Fall, da dadurch die wesentlichen Punkte bereits aufgezeigt werden. -Für ein nicht-lineares System werden Extended Kalman-Filter benötigt, bei denen die System-Matrix (A) durch die Jacobi-Matrix des System ersetzt wird. - -\subsection{Geschichte} -Das Kalman-Filter wurde 1960 von Rudolf Emil Kalman entdeckt und direkt von der NASA für die Appollo Mission benutzt. -Das Filter kommt mit wenig Rechenleistung aus und war somit dafür geeignet die Rakete bei der Navigation zu unterstützen.
Eine typische Anwendungen des Kalman-Filters ist Glättung von verrauschten Daten und die Schätzung von Parametern. Dies kommt heutzutage in jedem Satellit, Navigationssystem, Smartphones und Videospielen vor. - -\subsection{Wahrscheinlichkeit} -Das Kalman-Filter schätzt den wahrscheinlichsten Wert zwischen Normalverteilungen. -Dies bedeutet, das Filter schätzt nicht nur den Mittelwert, sondern auch die Standartabweichung. -Da Normalverteilungen dadurch vollständig definiert sind, schätzt ein Kalman-Filter die gesamte Verteilungsfunktion des Zustandes. -In der Abbildung~\ref{erdbeben: Zwei Normalverteilungen} sind zwei Funktionen dargestellt. -Die eine Funktion zeigt die errechnete Vorhersage des Zustands, bzw. deren Normalverteilung. -Die andere Funktion zeigt die verrauschte Messung des nächsten Zustand, bzw. deren Normalverteilung. -Wie man am Beispiel der Gauss-Verteilungen in Abblidung~\ref{erdbeben: Zwei Normalverteilungen} sehen kann, ist sowohl der geschätzte Zustand als auch der gemessene Zustand normalverteilt und haben dementsprechend unterschiedliche Standardabweichungen $\sigma$ und Erwartungswerte $\mu$. Dies wird in~\cite{erdbeben:aragher_understanding_2012}beschrieben. +Im letzten Abschnitt haben wir Gleichungen für unser System gefunden. +Als nächstes brauchen wir also ein Werkzeug, +um aus der Messung der Position $s(t)$ den gesammten Zustand $x(t)$ zu schätzen. +Das ist genau das, was Kalman-Filter tun: Ahand von Messungen den Zustand eines Systems schätzen. + +Kalman-Filter wurde 1960 von Rudolf Emil Kalman erfunden und direkt von der NASA für die Appollo Mission benutzt. +Diese Filter kommen mit wenig Rechenleistung aus und waren somit geeignet, die Rakete bei der Navigation zu unterstützen. +Heutige, typische Anwendungen von Kalman-Filtern sind die Glättung verrauschter Daten und die Schätzung von Parametern. +Dies kommt heutzutage in jedem Satellit, Navigationssystem, Smartphones und Videospielen vor. + +Kalman-Filter funktionieren nach folgendem Zwei-Schritt-Verfahren: +Zuerst wird, +ausgehend von der aktuellen Schätzung des Zustands und der Eigendynamik des Systzems, +eine Vorhersage berechnet. +Daraus lässt sich eine erwartete Messung ableiten. +Anschliessend wird diese Vorhersage korrigiert, +wobei die Korrektur abhänging von der Differenz zwischen erwarteter und effektiver Messung ist. + +Dabei sind sowohl die Vorhersage als auch die Messung nur Schätzungen und unweigerlich fehlerbehaftet. +Unter der Annahme, dass die Fehler normalverteilt sind, +lassen sich beide Schätzungen zu einer neuen, optimalen Schätzung kombinieren. +Die genaue Herleitung des Kalman-Filters ist relativ aufwendig +und kann unter Anderem in \cite{erdbeben:skript:wrstat} nachgelesen werden. + +\subsection{Exkurs Wahrscheinlichkeit} +\label{erdbeben:Wahrscheindlichkeit} +Das Kalman-Filter schätzt also den wahrscheinlichsten Wert zwischen zwei Normalverteilungen, +genauer gesagt zwischen einer Messung und einer Vorhersage. +In diesem Abschnitt wollen wir auffrischen, wie dies genau passiert. + +Das Folgende wird in \cite{erdbeben:aragher_understanding_2012} beschrieben. +Wir haben eine Vorhersage aus der Systemdynamik und eine Messung des Zustandes. +Diese widersprechen sich im Allgemeinen. +Jedoch kennen wir auch die Wahrscheinlichkeiten der beiden Aussagen. + \begin{figure} \begin{center} - \includegraphics[width=5cm]{papers/erdbeben/Gausskurve2.pdf} - \caption{Zwei Normalerteilungen; Die eine Funktion zeigt die Vorhersage, die andere die Messung} - \label{erdbeben: Zwei Normalverteilungen} + \includegraphics[width=5cm]{papers/erdbeben/Gausskurve3.pdf} + \caption{ + Seien blau und orange zwei normalverteilte Schätzungen eines Zustandes, etwa eine Vorhersage und eine Messung. + Dann ist die rote Kurve die optimale Schätzung. + Sie entspricht bis auf Normierung dem Produkt von blau und orange.} + \label{erdbeben:Gauss3} \end{center} \end{figure} -Wir haben eine Vorhersage aus der Systemdynamik und eine Messung des Zustandes. -Diese widersprechen sich im Allgemeinen. -Jedoch wissen wir die Wahrscheinlichkeiten der beiden Aussagen. +Abbildung~\ref{erdbeben:Gauss3} zeigt in blau und rot zwei Normalverteilungen, +je eine für die Vorhersage und eine für die Messung. +Diese unterscheiden sich sowohl in ihren Mittelwerten $\mu_{1,2}$, als auch in ihren Standardabweichungen $\sigma_{1,2}$. Um eine genauere Schätzung des Zustandes zu machen, wird nun ein Wert zwischen den beiden Verteilungen berechnet. -Nun wird eine Eigenschaft der Normalverteilung ausgenutzt. Durch das Multiplizieren zweier Normalverteilungen entsteht eine neue Normalverteilung. -Wir haben eine Normalverteilung der Vorhersage: +Nun wird eine Eigenschaft der Normalverteilung ausgenutzt: +Durch das Multiplizieren zweier Normalverteilungen entsteht eine neue Normalverteilung. + +Wir haben also eine Normalverteilung der Vorhersage \[ {y_1}(x;{\mu_1},{\sigma_1})=\frac{1}{\sqrt{2\pi\sigma_1^2}}\quad e^{-\frac{(x-{\mu_1})^2}{2{\sigma_1}^2}} \] -und der Messung: +und der Messung \[ {y_2}(x;{\mu_2},{\sigma_2})=\frac{1}{\sqrt{2\pi\sigma_2^2}}\quad e^{-\frac{(x-{\mu_2})^2}{2{\sigma_2}^2}}. \] -Diesen werden nun multipliziert und durch deren Fläche geteilt um sie wieder zu normieren, $\odot$ beschreibt dabei die Multiplikation und die Normierung auf den Flächeninhalt eins : +Diesen werden nun multipliziert und durch deren Fläche geteilt, +um sie wieder zu normieren. +$\odot$ beschreibt dabei die Multiplikation und die Normierung auf den Flächeninhalt eins: \begin{align*} - {y_f}(x; {\mu_f}, {\sigma_f}) = {y_1}(x;{ \mu_1},{ \sigma_1}) \odot {y_2}(x; {\mu_2}, {\sigma_2}) + {y_f}(x; {\mu_f}, {\sigma_f}) + &= + {y_1}(x;{ \mu_1},{ \sigma_1}) \odot {y_2}(x; {\mu_2}, {\sigma_2}) + \\ &= \frac{1}{\sqrt{2\pi\sigma_1^2}}\quad e^{-\frac{(x-{\mu_1})^2}{2{\sigma_1}^2}} \odot \frac{1}{\sqrt{2\pi\sigma_2^2}}\quad e^{-\frac{(x-{\mu_2})^2}{2{\sigma_2}^2}} \\ &= \frac{ \frac{1}{\sqrt{2\pi\sigma_1^2}}e^{-\frac{(x-{\mu_1})^2}{2{\sigma_1}^2}} \cdot \frac{1}{\sqrt{2\pi\sigma_2^2}}e^{-\frac{(x-{\mu_2})^2}{2{\sigma_2}^2}}}{\int {y_1} {y_2} dx}. \end{align*} -Diese Kombination der beiden Verteilungen resultiert wiederum in einer Normalverteilung -mit Erwartungswert +Die genaue Berechnung ist nicht schwierig aber aufwendig und wird hier deshalb ausgelassen. +Nach einigem Rechnen findet man die Ausdrücke \[ \mu_f = \frac{\mu_1\sigma_2^2 + \mu_2 \sigma_1^2}{\sigma_1^2 + \sigma_2^2} \] -und Varianz +für den neuen Mittelwert und \[ \sigma_f^2 = \frac{\sigma_1^2 \sigma_2^2}{\sigma_1^2 + \sigma_2^2}. \] -Dadurch gleicht sich die neue Kurve den anderen an. Interessant daran ist, dass die fusionierte Kurve sich der genauere Normal-Verteilung anpasst. -Ist ${\sigma_2}$ klein und ${\sigma_1}$ gross, so wird sich die fusionierte Kurve näher an ${y_2}(x;{\mu_2},{\sigma_2})$ begeben. -Somit ist $\mu_f$ ist das gewichtete Mittel der beiden $\mu_{1,2}$, und die Varianzen sind die Gewichte! -Die neue Funktion ist die best mögliche Schätzung für zwei Verteilungen, welche den selben Zustand beschreiben. -Dies ist in der Abbildung~\ref{erdbeben:Gauss3} anhand der rote Funktion ersichtlich. -\begin{figure} - \begin{center} - \includegraphics[width=5cm]{papers/erdbeben/Gausskurve3.pdf} - \caption{Durch das Multiplizieren der blauen und der orangen Verteilung entsteht die die rote, optimale Funktion} - \label{erdbeben:Gauss3} - \end{center} -\end{figure} +für die Varianz. + +Interessant daran ist, dass sich die fusionierte Kurve der genauere Normal-Verteilung anpasst. +Ist ${\sigma_2}$ klein und ${\sigma_1}$ gross, +so wird sich die fusionierte Kurve näher an ${y_2}(x;{\mu_2},{\sigma_2})$ begeben. +$\mu_f$ ist das gewichtete Mittel der beiden $\mu_{1,2}$, und die Varianzen $\sigma_{1,2}$ sind die Gewichte. +Das Interessante an $\mu_{f}$ ist, dass ${\mu_2}$ das Gewicht für ${\sigma_1}$ ist. +Somit ist die Unsicherheit der Messung das Gewicht der Vorhersage und umgekehrt. +Diese neue Funktion ist die best mögliche Schätzung für zwei Verteilungen, welche den selben Zustand beschreiben. +Dies ist in der Abbildung~\ref{erdbeben:Gauss3} anhand der roten Funktion ersichtlich. + Was in zwei Dimensionen erklärt wurde, funktioniert auch in mehreren Dimensionen. Dieses Prinzip mach sich das Kalman Filter zu nutze, und wird von uns für die Erdbeben Berechnung genutzt. -\section{Filter-Matrizen} -Da wir nun ein Werkzeug besitzen, dass die Beschleunigung, welche auf das Gehäuse wirkt, ermitteln kann, wird dieses nun Schritt für Schritt erklärt. -Um den Kalman Filter zu starten, müssen gewisse Bedingungen definiert werden. +\subsection{Filter-Matrizen} +Da wir nun ein Werkzeug besitzen, dass die Beschleunigung, welche auf das Gehäuse wirkt, ermitteln kann, +wird dieses nun Schritt für Schritt erklärt. +Um das Kalman Filter zu starten, müssen gewisse Bedingungen definiert werden. In diesem Abschnitt werden die einzelnen Parameter und Matrizen erklärt und erläutert, wofür sie nützlich sind. -\subsection{Fiter-Agorithmus} -Nachdem alle Parameter aufgestellt sind, wird das Filter initialisiert. -Zuerst wird der nächste Zustand der Masse vorhergesagt, danach wird die Messung präzisiert und laufend aktualisiert. -Das Filter berechnet aufgrund der aktuellen Schätzung eine Vorhersage. -Diese wird, sobald verfügbar, mit der Messung verglichen. -Aus dieser Differenz und den Unsicherheiten des Prozesses ($Q$) und der Messung ($R$) wird der wahrscheinlichste, neue Zustand geschätzt. -Dabei muss genau auf den Index geachtet werden. Nach dem Artikel~\cite{erdbeben:wikipedia} ist die Indexierung so genormt: -Der Zeitschritt wird mit $k$ definiert, $k-1$ ist somit ein Zeitschritt vor $k$. -Auf der linken Seite von | wird der aktuelle Zustand verlangt, bzw. ausgegeben, auf der rechten Seiten den bisherigen Zustand. -Dies bedeutet, dass die Notation $x_{n|m}$ die Schätzung von $x$ zum Zeitpunkt $n$ bis und mit zur Zeitpunkt $m \leq \ n$ präsentiert. +Dabei muss genau auf den Index geachtet werden. +Wir verwenden die Standard-Notation, wie sie auch im Artikel~\cite{erdbeben:wikipedia} zu finden ist. +Sie ist an die Notation der bedingten Wahrscheinlichkeiten angelehnt. +Hierbei steht der betrachtete Zeitschritt links und der gegenwärtige rechts eines Vertikalstrichs. +Dies bedeutet, dass die Notation $x_{n|m}$ die Schätzung von $x$ zum Zeitpunkt $n$ +aufgrund des Wissens bis zum und mit dem Zeitpunkt $m$ repräsentiert. \subsubsection*{Vorhersage} -Im Filterschritt Vorhersage wird der nächste Zustand anhand des Anfangszustand und der Systemmatrix berechnet. -Dies funktioniert mit dem Rechenschritt: +Im Filterschritt Vorhersage wird anhand des aktuellen Zustands und der Systemmatrix eine Schätzung für den nächsten Zustand berechnet. +Die Systemmatrix $A$ aus Gleichung~\eqref{erdbeben:systemmatrix} beschreibt ein kontinuierliches System $\dot x = Ax$. +Wir benötigen jedoch ein Zeit-diskretes System $x_{k+1} = \Phi x_k$. + +Die Exponentialfunktion $\exp(At)$ beschreibt die Entwicklung eine Zustandes im Laufe der Zeit. +Die Übergangs-Matrix $\Phi$ erhalten wir folglich aus der Systemdynamikmatrix $A$ durch die Exponentialfunktion +\[\Phi = \exp(A\Delta t). \] +Die Matrix $\Phi$ beschreibt die Übergänge zwischen zeitlich aufeinanderfolgenden Zuständen $x_{k-1}$ und $x_{k}$ anhand folgender Gleichung: \[ {x_{k|k-1}}=\Phi{x_{k-1|k-1}}= \exp(A\Delta t){x_{k-1|k-1}}. \] -Die Kovarianz $P_{k|k-1}$ wird ebenfalls neu berechnet. Zudem kommt noch die Prozessunsicherheit $Q$ dazu, so dass die Unsicherheit des Anfangsfehlers $P$ laufend verändert. -Dies funktioniert durch multiplizieren der Systemmatrix mit dem aktualisierten Anfangsfehler. -Dazu wird noch die Prozessunsicherheit addiert, somit entsteht die Gleichung +Damit haben wir die Systemdynamik nun in der für unser Kalman-Filter notwendigen Form und können Vorhersagen berechnen. + +Als nächstes benötigen wir die Unsicherheit der Vorhersage. +Im Abschnitt ~\ref{erdbeben:Wahrscheindlichkeit} haben wir dafür die Varianzen der Normalverteilungen verwendet. +Im mehrdimensionalen Fall übernimmt dies die Kovarinanzmatrix $P$. +Sie wird in jedem Schritt aktualisiert. +Hinzu kommt die Prozessunsicherheit $Q$, welche als Parameter in unser Modell einfliesst. +$Q$ beschreibt Unsicherheiten im Modell, +wie etwa unsere Annahme, dass die Kraft sich nicht ändert, +aber auch nicht-modellierbare Einflüsse wie Vibrationen. +$P$ wird dabei laufend aktuallisiert. +Die optimale Gleichung lautet \[ {P_{k|k-1}}=\Phi {P_{k-1|k-1}} {\Phi _{k}}^T + {Q_{k-1}}. \] -Es vergeht genau $\Delta t$ Zeit, und dieser Vorgang wird wiederholt. -Das hochgestellte T bezeichnet die transponierte Matrix. -Dabei wird in den späteren Schritten überprüft, wie genau die letzte Anpassung von $P$ zur Messung stimmt. -Ist der Unterschied klein, wird die Kovarianz $P$ kleiner, ist der Unterschied gross, wird auch die Kovarianz grösser. +Es vergeht genau $\Delta t$ Zeit, und dieser Vorgang wird wiederholt. Das Filter passt sich selber an und korrigiert sich bei grosser Abweichung. \subsubsection*{Messen} -Der Sensor wurde noch nicht benutz, doch genau der liefert Werte für das Filter. -Die aktuellen Messwerte $z$ werden die Innovation $w$ mit dem Zustandsvektor $x$ und der Messmatrix $H$ zusammengerechnet. -Hier bei wird lediglich die Messung mit dem Fehler behaftet, und die Messmatrix $H$ mit der Vorhersage multipliziert. +Der Sensor wurde noch nicht benutz, doch genau der liefert die Messwerte $z_k$ für unser Filter. +Aus der Vorhersage des Zustandes $x_{k|k-1}$ und der Messmatrix $H$ erhalten wird eine Vorhersage der Messung. +Die Innovation \[ -{w_{k}}={z_{k}}-{H}{x_{k|k-1}}. +{w_{k}}={z_{k}}-{H}{x_{k|k-1}} \] -Die Innovation ist der Teil der Messung, die nicht durch die Systemdynamik erklärt werden kann. -Die Hilfsgröße Innovation beschreibt, wie genau die Vorhersage den aktuellen Messwert mittels der Systemmatrix $\Phi$ beschreiben kann. +beschreibt, wie genau die Vorhersage $x_{k|k-1}$ zur aktuellen Messung $z_k$ passt. +Die Innovation ist also derjenige Teil der Messung, der nicht im Modell erfasst ist. +Dies leuchtet ein, eine Innovation von $0$ bedeutet, dass die Messung nichts Neues hervorbrachte. Für eine schlechte Vorhersage wird die dazugehörige Innovation gross, für eine genaue Vorhersage dagegen klein sein. -Entsprechende Korrekturen müssen dann gross bzw. nur gering ausfallen. -Innovation = Messung - Vorhersage. Dies leuchtet ein, eine Innovation von 0 bedeutet, dass die Messung nichts Neues hervorbrachte. - -Im nächsten Schritt wir analysiert, mit welcher Kovarianz weiter gerechnet wird. -Hierbei wird die Unsicherheit $P$, die Messmatrix $H$ und die Messunsicherheit $R$ miteinander verrechnet. -\[ -{S_{k}}={H}{P_{k|k-1}}{H}^T+{R_{k}} -\] +Entsprechende Korrekturen werden dann gross bzw. nur gering ausfallen. \subsubsection*{Aktualisieren} -Im nächsten Schritt kommt nun die Wahrscheinlichkeit dazu. -\[{K_{k}}= {P_{k|k-1}} {H^T}{S_{k}^{-1}}\] -Die Grösse $K$ wird Kalman-Gain genannt. -Das Kalman-Gain gibt dem Zustand die Gewichtung, bzw. wie die Vorhersage auf den Zustand passt. -Vereinfacht gesagt: Es wird das das Verhältnis zwischen der Unsicherheit der Vorhersage $P_k$ zu der zugehörigen Messunsicherheit $R_k$ gebildet. -In unserem Fall wird werden die Elemente der Kalman-Matrix vorweg berechnet, da das Kalman-Gain ohne Messungen auskommt. - -Anhand der Informationen aus der Innovation wird das Kalman-Gain $K$ gebildet. Dabei beschreibt das Kalman-Gain die Wirkung der Innovation auf den geschätzten Zustand. So wird das System aktualisiert. + +Für eine optimale Schätzung des Zustandes muss die Vorhersage entsprechend der Innovation korrigiert werden. +In der Literatur findet man für eine optimales Korrektur die Gleichungen +\begin{align*} +{S_{k}} &={H}{P_{k|k-1}}{H}^T+{R_{k}} +\\ +{K_{k}} &= {P_{k|k-1}} {H^T}{S_{k}^{-1}} +\end{align*} +Dabei ist $K$ das Kalman-Gain. +$K$ beschreibt, wie die Vorhersage korrigiert werden muss. +Die optimale Schätzung des neuen Zustandes wird dann zu \[ -{x_{k|k}}={x_{k|k-1}}+{K_{k}}{w_{k}} +{x_{k|k}}={x_{k|k-1}}+{K_{k}}{w_{k}}. \] -Dabei wird der Unterschied zwischen dem erwarteten, errechneten, Zustand und dem gemessenen Zustand berechnet. - -Dazu kommt eine neue Kovarianz für den nächste Vorhersageschritt: +Dazu kommt eine neue Kovarianz $P$ für den nächste Vorhersageschritt: \[ -{P_{k|k}}=(I-{K_{k}}{H}){P_{k|k-1}} +{P_{k|k}}=(I-{K_{k}}{H}){P_{k|k-1}}. \] -Der ganze Algorithmus und beginnt wieder mit der Vorhersage +Der ganze Algorithmus ist nun vollständig und beginnt wieder mit der Vorhersage \[ -{x_{k|k-1}}=\Phi{x_{k-1|k-1}}= \exp(A\Delta t){x_{k|k-1}}. +{x_{k+1|k}}=\Phi{x_{k|k}}= \exp(A\Delta t){x_{k|k}}. \] -\subsection{Anfangsbedingungen} +\subsection{Parameter und Anfangsbedingungen} +Die Grössen $P$, $Q$, $R$ und $\Phi$ können grundsätzlich in jedem Zeitschritt ändern. +Für die meisten Anwendungen sind sie jedoch konstant und fliessen als Parameter ins Modell ein. +Aufgrund der iterativen Arbeitsweise von Kalman-Filtern benötigen wir zudem ein paar Anfangswerte. + \subsubsection*{Anfangszustand $x$} -Das Filter benötigt eine Anfangsbedingung. +Für die erste Vorhersage benötigt das Filter einen Anfangszustand. In unserem Fall ist es die Ruhelage, die Masse bewegt sich nicht. Zudem erfährt die Apparatur keine äussere Kraft. \[ {x_0 }= \left( \begin{array}{c} {s_0}\\ {v_0}\\{f_0}\end{array}\right) = \left( \begin{array}{c} 0\\ 0\\ 0\end{array}\right) \] +\subsubsection*{Systemmatrix $A$ und $\Phi$} +Für unseren Seismographen haben wir die entsprechende Matrixdarstellung +in Gleichung~\eqref{erdbeben:systemmatrix} bereits gefunden. +Zudem haben wir weiter oben bereits beschrieben, +wie wir mittels Exponentialfunktion zu einer zeitdiskreten Beschreibung für das Kalman-Filter kommen. +Es gilt +\[ \Phi = \exp(A \Delta t) .\] + \subsubsection*{Anfangsfehler / Kovarianzmatrix $P$} -Da auch der Anfangszustand fehlerhaft sein kann, wird für das Filter ein Anfangsfehler verwendet. +Da auch der Anfangszustand fehlerhaft sein kann, wird für das Filter eine Anfangsunsicherheit verwendet. Auf der Diagonalen werden die Varianzen eingesetzt, in den restlichen Felder stehen die Kovarianzen. -Zur Erinnerung: Die Varianz ist ein Mass für die Streuung eines Wertes, die Kovarianz hingegen beschreibt die Abhängigkeit der Streuungen zweier Werte. - -Kovarianz: Cov(x, y) und Varianz: Var(x) = Cov(x, x) +Für einen gut bekannten Zustandsvektor können kleine Werte eingesetzt werden, für ungenaue Anfangsbedingungen sollten grosse Werte verwendet werden. +Grosse Werte ermöglichen dem Filter sich schnell einzupendeln. In unserem Fall ist der Anfangszustand gut bekannt. -Wir gehen davon aus, dass das System in Ruhe und in Abwesenheit eines Erdbeben startet, somit kann die Matrix mit Nullen bestückt werden. -Als Initialwert für die Kovarianzmatrix ergibt sich +Wir gehen davon aus, +dass das System in Ruhe und in Abwesenheit eines Erdbeben startet. +Somit kann die Matrix mit Nullen bestückt werden und wir starten mit \[ {P_0 }= \left( @@ -200,35 +226,13 @@ Als Initialwert für die Kovarianzmatrix ergibt sich \end{array} \right). \] -Diese Matrix beschreibt die Unsicherheit des geschätzten Zustandes und wird sowohl für die Vorhersage als auch die Korrektur benötigt. -Sie wird nach jeder Schätzung aktualisiert. -Für einen gut bekannten Zustandsvektor können kleine Werte eingesetzt werden, für ungenaue Anfangsbedingungen sollten grosse Werte verwendet werden. -Grosse Werte ermöglichen dem Filter sich schnell einzupendeln. -\subsubsection*{Dynamikmatrix $A$} -Das Kalman-Filter benötigt für die Vorhersage des nächsten Zustandes eine Beschreibung der Systemdynamik. -Die Dynamikmatrix bildet den Kern des Filters. Diese wurde weiter oben bereits beschrieben. -Dabei wollen wird die äussere Kraft des Systems ermitteln. -Da nichts über die äussere Kraft bekannt ist, müssen wir annehmen das deren Ableitung 0 ist. -Die System-Matrix lautet daher: -\[ -A = \left( - \begin{array}{ccc} -0 & 1& 0 \\ -- \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m}\\ -0 & 0& 0\\ -\end{array}\right) - \] -Dabei soll der Kalman-Filter in diskreten Zeitschritten $\Delta t$ arbeiten. -$A$ beschreibt ein kontinuierliches System ($\dot x = Ax$), wir benötigen jedoch ein Zeit-diskretes System $x_{k+1} = \Phi x_k$. -Die Übergangs-Matrix erhalten wir aus der Systemdynamikmatrix mittels Exponentialfunktion: -\[\Phi = \exp(A\Delta t). \] -Die Matrix $\Phi$ beschreibt die Übergänge zwischen zeitlich aufeinanderfolgenden Zuständen $x_{k-1}$ und $x_{k}$ \subsubsection*{Prozessrauschkovarianzmatrix $Q$} Die Prozessrauschmatrix teilt dem Filter mit, wie sich der Prozess verändert. Die Matrix $Q$ beschreibt die Unsicherheit, die der Prozess mit sich bringt. -Bei unserem Modell könnte das beispielsweise ein Windstoss an die Masse sein oder auch die Ungenauigkeiten im Modell, wie die Annahme das dich die Kraft nicht ändert. +Bei unserem Modell könnte das beispielsweise ein Windstoss an die Masse sein +oder auch Ungenauigkeiten im Modell, wie die Annahme, dass sich die Kraft nicht ändert. Für uns wäre dies: \[ Q = \left( @@ -241,9 +245,9 @@ Q = \left( Die Standabweichungen müssten statistisch ermittelt werden, da der Fehler nicht vom Sensor kommt und somit nicht vom Hersteller gegeben ist. \subsubsection*{Messmatrix $H$} -Die Messmatrix gibt an, welche Parameter gemessen werden. -$H$ ist die Matrix für die Vorhersage der Messung. -In unserem Falle ist es die Position der Massen. +Die Messmatrix gibt an, welche Zustände gemessen werden. +$H$ ist die Matrix, welche aus der Vorhersage des Zustand eine Vorhersage der Messung erzeugt. +In unserem Falle messen wir nur die Position der Massen und verwenden deshalb \[ H = (1, 0, 0) \] @@ -252,49 +256,48 @@ H = (1, 0, 0) Die Messrauschkovarianzmatrix beinhaltet, wie der Name schon sagt, das Rauschen der Messung. In unserem Fall wird nur die Position der Masse gemessen. Da wir keine anderen Sensoren haben ist $R$ lediglich: \[ -R= ({\sigma_\mathrm{sensor}}^2). +R= (\sigma_\mathrm{sensor}^2). \] Diese Messrauchen wird meistens vom Sensorhersteller angegeben. -Für unsere theoretische Apparatur wird hier ein kleiner Fehler eingesetzt da heutige Sensoren sehr genau messen können. +Für unsere theoretische Apparatur wird hier ein kleiner Fehler eingesetzt, +da heutige Sensoren sehr genau messen können. +\clearpage \subsection{Zusammenfassung } -Zusammenfassend kann das Kalman-Filter in offizieller Typus dargestellt werden. -Dabei beginnt das Filter mit dem Anfangszustand für $k=0$ - -1. Nächster Zustand vorhersagen +Das Filter beginnt mit dem Anfangszustand für $k=0$. +Anschliessend werden folgende Schritte iterativ ausgeführt: +\begin{enumerate} +\item Nächster Zustand vorhersagen \[ -{x_{k|k-1}}=\Phi{x_{k-1|k-1}}= \exp(A\Delta t){x_{k-1|k-1}}. +{x_{k|k-1}}=\Phi{x_{k-1|k-1}}= \exp(A\Delta t){x_{k-1|k-1}} \] -2. Nächste Fehlerkovarianz vorhersagen + \item Nächste Fehlerkovarianz vorhersagen \[ -{P_{k|k-1}}=\Phi {P_{k-1|k-1}} {\Phi _{k}}^T + {Q_{k-1}}. +{P_{k|k-1}}=\Phi {P_{k-1|k-1}} {\Phi _{k}}^T + {Q_{k-1}} \] -3. Zustand wird gemessen +\item Innovation (= Messung - Vorhersage) \[ -{w_{k}}={z_{k}}-{H}{x_{k|k-1}}. +{w_{k}}={z_{k}}-{H}{x_{k|k-1}} \] -4. Innovation (= Messung - Vorhersage) -\[ -{S_{k}}={H}{P_{k|k-1}}{H}^T+{R_{k}} -\] - -5. Das Kalman Filter anwenden -\[ -{K_{k}}= {P_{k|k-1}} {H^T}{S_{k}^{-1}} -\] +\item Optimales Kalman-Gain berechnen +\begin{align*} +{S_{k}} &={H}{P_{k|k-1}}{H}^T+{R_{k}}\\ +{K_{k}} &= {P_{k|k-1}} {H^T}{S_{k}^{-1}} +\end{align*} -6. Schätzung aktualisieren +\item Schätzung aktualisieren \[ {x_{k|k}}={x_{k|k-1}}+{K_{k}}{w_{k}} \] -7. Fehlerkovarianz aktualisieren +\item Fehlerkovarianz aktualisieren \[ {P_{k|k}}=(I-{K_{k}}{H}){P_{k|k-1}} \] -8. Die Outputs von $k$ werden die Inputs für ${k-1}$ und werden wieder im Schritt 1 verwendet +\end{enumerate} +Die Outputs von $k$ werden die Inputs für ${k+1}$ und werden wieder in Schritt 1 verwendet. |