diff options
| -rw-r--r-- | buch/papers/erdbeben/teil0.tex | 126 | ||||
| -rw-r--r-- | buch/papers/erdbeben/teil1.tex | 190 |
2 files changed, 183 insertions, 133 deletions
diff --git a/buch/papers/erdbeben/teil0.tex b/buch/papers/erdbeben/teil0.tex index f012225..9b33e4b 100644 --- a/buch/papers/erdbeben/teil0.tex +++ b/buch/papers/erdbeben/teil0.tex @@ -3,11 +3,12 @@ % % (c) 2020 Prof Dr Andreas Müller, Hochschule Rapperswil %% -\section{Was ist ein Erdbeben? \label{erdbeben:section:teil0}} -\rhead{Erdbeben} -Für das Verständnis möchten wir zuerst erklären, was ein Erdbeben genau ist. -Das soll uns helfen, eine Verknüpfung zwischen dem Naturphänomen und der mathematischen Problemstellung herzustellen. +%\section{Was ist ein Erdbeben? \label{erdbeben:section:teil0}} +\rhead{Erdbeben} +%Für das Verständnis möchten wir zuerst erklären, was ein Erdbeben genau ist. +%Das soll uns helfen, eine Verknüpfung zwischen dem Naturphänomen und der mathematischen Problemstellung herzustellen. +\noindent Unter einem Erdbeben verstehen wir eine Erschütterung des Erdkörpers. Dabei reiben zwei tektonische Platten aneinander, welche sich durch die Gesteinsverzahnung gegenseitig blockieren. Diese Haftreibung durch die Steine wird so lange aufgebaut, bis sie nicht mehr gehalten werden kann. @@ -17,23 +18,31 @@ Vergleichbar ist, wenn man einen Stein in einen Teich wirft und die Wellen beoba \subsection{Funktion eines Seismograph} Um ein Erdbeben kenntlich zu machen, werden in der Regel Seismographen mit vielen Sensoren verwendet. -Ein Seismograph besteht im Grunde aus einer federgelagerten Masse. Wirkt eine Bodenerregung auf das Gerät ein, schwing das Gehäuse und dadurch auch die gekoppelte Masse. -Stoppt das Erdbeben, schwingt das Gehäuse nicht mehr. -Die Masse schwing jedoch in seiner Eigendynamik weiter. -Eine Relativbewegung des Bodens kann damit als Auslenkung im Zeitverlauf gemessen werden. -In modernen Seismographen wird die Bodenbewegung in alle Richtungen gemessen, sowohl Horizontal als auch Vertikal. - - -Wir konstruieren uns eine einfachere Version eines Seismographen mit einem Gehäuse, an dem zwei Federn und eine Masse befestigt sind. -Der Seismograph ist in Abbildung ~\ref{erdbeben:Seismograph} ersichtlich. -Ein Sensor unter der Masse misst die Position, bzw. die Auslenkung der Feder und der Masse. -Dies bedeutet, unser Seismograph kann nur in eine Dimension Messwerte aufnehmen. - -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. -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. - +Ein Seismograph besteht im Grunde aus einer federgelagerten Masse. +Bei einem Erdbeben folgt das Gehäuse direkt der Bewegung des Erdbebens. +Die federgelagerte Masse wird jedoch erst durch die Feder bewegt und folgt verzögert. +Zudem schwingt die Masse auch ohne Erdbeben weiter -- das System besitzt eine Eigendynamik. +Eine Relativbewegung des Bodens kann folglich als Auslenkung im Zeitverlauf gemessen werden. +Allerdings misst man so nicht direkt das Erbeben, sondern eine Überlagerung der Effekte aus Erdbeben- und Federkraft. + + +In modernen Seismographen wird die Bodenbewegung in alle Richtungen gemessen, +sowohl Horizontal als auch Vertikal. +Wir konstruieren hier eine einfachere Version eines Seismographen mit einem Gehäuse, +an dem zwei Federn und eine Masse befestigt sind. +Abbildung~\ref{erdbeben:Seismograph} zeigt eine schematische Darstellung unseres Systems. +Ein Sensor unter der Masse misst die Position der Masse relativ zum Gehäuse. +Unser Seismograph misst also nur eindimensional. + +Für mehrere Dimensionen würde der Satz von Pythagoras für die Auslenkung der Federn benötigt. +Die benötigten Quadrate und Wurzeln brechen jedoch die Linearität des Systems. +Die Systembeschreibung wird dann deutlich komplexer. +Der Einfachheit halber beschränken wir uns deshalb auf den linearen Fall, +welcher bereits alle wesentlichen Punkte aufgezeigen kann. + +Wir werden sehen, dass diese Art der Problemstellung effektiv mittels Kalman-Filter gelöst werden kann. +Für ein nicht-lineares System werden Extended Kalman-Filter benötigt, +bei denen die System-Matrix $A$ durch die Jacobi-Matrix ersetzt wird. \begin{figure} \begin{center} @@ -45,58 +54,67 @@ Für ein nicht-lineares System werden Extended Kalman-Filter benötigt, bei dene \subsection{Ziel} Unser Seismograph misst nur die Position der Masse über die Zeit. -Wir wollen jedoch die Beschleunigung $a(t)$ des Boden, bzw. die Kraft $f(t)$, welche auf das Gehäuse wirkt, bestimmten. -Anhand dieser Beschleunigung, bzw. der Krafteinwirkung durch die Bodenbewegung, wird später das Bauwerk bemessen. +Wir wollen jedoch die Beschleunigung $a(t)$ des Boden, +respektive die Kraft $f(t)$, +welche auf das Gehäuse wirkt, bestimmten. +Anhand dieser Beschleunigung, +beziehungsweise der Krafteinwirkung durch die Bodenbewegung, +wird später das Bauwerk bemessen. Dies bedeutet, die für uns interessante Grösse $f(t)$ wird nicht durch einen Sensor erfasst. Jedoch können wir durch zweifaches ableiten der Positionsmessung $s(t)$ die Beschleunigung der Masse berechnen. -Das heisst: Die Messung ist zweifach Integriert die Kraft $f(t)$ inklusive der Eigendynamik der Masse. -Um die Krafteinwirkung der Masse zu berechnen, müssen wir Gleichungen für unser System finden. +Die Messung entspricht also dem zweiten Integral der Kraft $f(t)$, +wobei diese einerseits durch das Erdbeben, und andererseits durch die Federn zustande kommt. +Im Folgenden möchten wir die Erdbeben- und Federkräfte trennen. +Dafür benötigen wir zuerst eine mathematische Beschreibung unseres Systems. \subsection{Systemgleichung} Im Paper~\cite{erdbeben:mendezmueller} wurde das System gleich definiert und vorgegangen. Im Fall unseres Seismographen, handelt es sich um ein Feder-Masse-Pendel. -Dieser kann durch die Differentialgleichung zweiter Ordnung einer gedämpften Schwingung am harmonischen Oszillator beschrieben werden. -Die Gleichung lautet: +Dieses kann als gedämpfter harmonischer Oszillator beschrieben werden. +Die zugehörige Differentialgleichung lautet: \begin{equation} \label{erdbeben:Systemgleichung} m\ddot s + 2k \dot s + Ds = f. \end{equation} wobei $m$ die Masse, $k$ die Dämpfungskonstante und $D$ die Federkonstante bezeichnet. - -Da die Differentialgleichung linear ist möchten wir diese Gleichung in die Darstellung $\dot x = Ax$ überführen, wobei $x$ der Zustandsvektor und $A$ die Systemmatrix bezeichnet. Dazu verwenden wir die Subsitution: -\[ -s_1 = s -\qquad \text{und} \qquad -s_2 = \dot s. -\] -Somit entstehen die Gleichungen für die Geschwindigkeit $ \dot s_1(t)$ der Masse : -\[ \dot {s_1} = {s_2}\] -und -\[ \dot s_2 = -\frac{D}{m} {s_1} -\frac{2k}{m} {s_2} + \frac{f} {m} \] -für die Beschleunigung $\dot s_2(t)$ der Masse. - -Die für uns relevanten Zustände sind die Position der Masse, die Geschwindigkeit der Masse und die äussere Beschleunigung des ganzen Systems. -Unüblich ist nun, dass der Stör-Term $f$ in Gleichung~\eqref{erdbeben:Systemgleichung} gerade das ist, was wir eigentlich bestimmen möchten. -In unserem Fall wird die äusseren Beschleunigung gesucht, da diese der Erdbebenanregung gleich kommt. +Für lineare Systeme ist eine Matrix-Darstellung handlicher. +Wir möchten diese Gleichung folglich in die Darstellung $\dot x = Ax$ überführen, +wobei $x$ der Zustandsvektor und $A$ die Systemmatrix bezeichnet. +Wir subsituieren $\dot s = v$ für die Geschwindigkeit und erhalten das Gleichungssystem +\begin{align} + \dot s &= v \\ + \dot v &= -\frac{D}{m} {s} -\frac{2k}{m} {v} + \frac{f} {m}. + \label{erdbenen:systemgleichungen} +\end{align} + +Die relevanten Zustände sind also die Position $s$ und die Geschwindigkeit $v$. +Die für uns eigentlich interessante Grösse ist jedoch der Stör-Term $f$. +Dieser entspricht der Kraft durch das Erdbeben. Deshalb nehmen wir $f$ als dritte Grösse in den Zustandsvektor auf und definieren: - \[ -x= \left( \begin{array}{c} {s_1}\\ {s_2}\\{f}\end{array}\right)^T + x = \begin{pmatrix} {s} \\ {v} \\ {f} \end{pmatrix} \] -Für die Standard-Form $\dot x = Ax$ brauchen wir als nächstes die Ableitungen aller Elemente von $x$. Für $\dot s_1$ und $\dot s_2$ folgen diese direkt aus Gleichung ~\eqref{erdbeben:Systemgleichung}, aber über $\dot f$ wissen wir nichts. -Wir müssen also eine Annahme treffen: $\dot f = 0$, somit verändert sich die Kraft über die Betrachtungszeit nicht. Diese Annahme ist im Allgemeinen falsch, aber etwas Besseres haben wir zurzeit nicht zur Verfügung. +Für die Standard-Form $\dot x = Ax$ brauchen wir als nächstes die Ableitungen aller Elemente von $x$. +Für $s$ und $v$ haben wir diese in Gleichung~\eqref{erdbenen:systemgleichungen} bereits gefunden. +Über die Kraft $f$ wissen wir jedoch nichts. +Wir müssen also eine Annahme treffen: Die Kraft ändert sich nicht, $\dot f = 0$. +Diese Annahme ist im Allgemeinen natürlich falsch, aber etwas Besseres haben wir nicht zur Verfügung. Wir werden dies in einem späteren Schritt kompensieren müssen. -Durch Rücksubstituion ergibt sich uns folgende Systemgleichung in Matrix schreibweise, wobei $\dot {s_1}= v$ ist. Damit haben wir nun alles, was wir für die Matrix-Darstellung von Gleichung (20.1) benötigen. Diese lautet: +Wir haben nun alles für die Matrix-Form von Gleichung~\eqref{erdbeben:Systemgleichung} zusammen. +Sie lautet: \begin{equation} -\frac{d}{dt} \left(\begin{array}{c} s(t) \\ v(t) \\ f(t) \end{array}\right) = \left( - \begin{array}{ccc} -0 & 1& 0 \\ -- \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m}\\ -0 & 0 & 0\\ -\end{array}\right) \left(\begin{array}{c} s(t)\\ v(t)\\ f(t) \end{array}\right). + \frac{d}{dt} \begin{pmatrix} s(t) \\ v(t) \\ f(t) \end{pmatrix} + = + \begin{pmatrix} + \phantom- 0 & \phantom-1& 0 \\ + - \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m} \\ + \phantom-0 & \phantom-0 & 0\\ + \end{pmatrix} + \begin{pmatrix} s(t) \\ v(t) \\ f(t) \end{pmatrix}. + \label{erdbeben:systemmatrix} \end{equation} diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index 4d6786f..6c2539a 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -12,52 +12,70 @@ \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. -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. -Durch Kenntnis über den aktuellen Zustand und der Eigendynamik des Systems berechnen wir eine Vorhersage des nächsten Zustandes. -Die für das Filter relevante Grösse ist dann nicht mehr die eigentliche Messung, sondern die Differenz aus Messung und Vorhersage, da diese Differenz, die Innovation, eine Aussage über die nicht-deterministischen, externen Einflüsse auf das System ermöglicht. - -Die genau Herleitung des Kalman-Filter befindet sich für Interessierte im Wahrscheinlichkeit und Statistik Skript von A. Müller. -Im folgenden Abschnitt werden die Resultate zitiert. - -\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. +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, im statistischen Sinne 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 den wahrscheinlichsten Wert zwischen Normalverteilungen, in unserem Fall sind dies die Messung und die Vorhersage. -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:Gauss3} 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:Gauss3} sehen kann, ist sowohl der geschätzte Zustand als auch der gemessene Zustand normalverteilt und haben dementsprechend unterschiedliche Standardabweichungen $\sigma$ und Erwartungswerte $\mu$. Das folgende wird in~\cite{erdbeben:aragher_understanding_2012}beschrieben. +\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 wissen wir die Wahrscheinlichkeiten der beiden Aussagen. +Jedoch kennen wir auch die Wahrscheinlichkeiten der beiden Aussagen. + +\begin{figure} + \begin{center} + \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} +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}) &= @@ -69,27 +87,23 @@ Diesen werden nun multipliziert und durch deren Fläche geteilt um sie wieder zu &= \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*} -Durch geschicktes umformen resultiert aus der Kombination der beiden Verteilungen wiederum in einer Normalverteilung -mit Erwartungswert +Die genaue Berechnung ist eine reine Fingerübung. +Nach einigem Aufwand 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$ das gewichtete Mittel der beiden $\mu_{1,2}$ und die Varianzen sind die Gewichte. -Das Interessante an $\mu_{f}$ ist, dass ${\mu_2}$ ${\sigma_1}$ beeinflusst. -Somit beeinflusst die Messung die Schätzung und umgekehrt. -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 Mess- und der orangen Schätz-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. @@ -99,27 +113,35 @@ Da wir nun ein Werkzeug besitzen, dass die Beschleunigung, welche auf das Gehäu Um den 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. -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 reechts eines Vertikalstrichs. +Dies bedeutet, dass die Notation $x_{n|m}$ die Schätzung von $x$ zum Zeitpunkt $n$ aufgrund des Wissens bis zum und und mit dem Zeitpunkt $m$ repräsentiert. \subsubsection*{Vorhersage} Im Filterschritt Vorhersage wird der nächste Zustand anhand des Anfangszustand und der Systemmatrix berechnet. -Die Systemmatrix $A$ beschreibt ein kontinuierliches System $\dot x = Ax$. +Die Systemmatrix $A$ aus Gleichung~\eqref{erdbeben:systemmatrix} beschreibt jedoch 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 mittels Exponentialfunktion: +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}}. \] +Damit haben wir die Systemdynamik nun in der für unser Kalman-Filter notwendigen Form und können Vorhersagen berechnen. -Im Abschnitt ~\ref{erdbeben:Wahrscheindlichkeit} benötigten wir die Varianzen der Normalverteilungen. +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. -Dazu wird noch die Prozessunsicherheit $Q$ addiert, somit entsteht die Gleichung +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}}. \] @@ -127,14 +149,14 @@ 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. +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}} \] -beschreibt, wie genau die Vorhersage den aktuellen Messwert $z$ mittels der Systemmatrix $\Phi$ beschreiben kann. -Die Innovation ist der Teil der Messung, der nicht im Modell erfasst ist. +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 werden dann gross bzw. nur gering ausfallen. @@ -148,8 +170,9 @@ In der Literatur findet man für eine optimales Korrektur die Gleichungen: \\ {K_{k}} &= {P_{k|k-1}} {H^T}{S_{k}^{-1}} \end{align*} -Dabei ist $K$ das Kalman-Gain. -Mit dessen Hilfe erhalten wir die optimale Schätzung des nächsten Zustandes +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}}. \] @@ -159,27 +182,37 @@ Dazu kommt eine neue Kovarianz $P$ für den nächste Vorhersageschritt: \] 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{Parameter und Anfangsbedingungen} -\subsubsection*{Anfangszustand $x$} +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. -Das Filter benötigt eine Anfangsbedingung. +\subsubsection*{Anfangszustand $x$} +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$} -Für unseren Seismographen haben wir die entsprechende Matrixdarstellung in Gleichung ~\eqref{erdbeben:Systemgleichung} bereits gefunden. +\subsubsection*{Systemmatrix $A$ und $\Phi$} +Für unseren Seismographen haben wir die entsprechende Matrixdarstellung in Gleichung ~\eqref{erdbeben:Systemgleichung} bereits gefunden. +Zudem haben wir weiter oben bereits entdeckt, 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. +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( @@ -190,14 +223,13 @@ Als Initialwert für die Kovarianzmatrix ergibt sich \end{array} \right). \] -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*{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 die Ungenauigkeiten im Modell, wie die Annahme, dass sich die Kraft nicht ändert. Für uns wäre dies: \[ Q = \left( @@ -211,8 +243,8 @@ Die Standabweichungen müssten statistisch ermittelt werden, da der Fehler nicht \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. +$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) \] @@ -246,7 +278,7 @@ Das Filter beginnt mit dem Anfangszustand für $k=0$ {w_{k}}={z_{k}}-{H}{x_{k|k-1}}. \] -\item Das Kalman Filter anwenden +\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}} @@ -262,6 +294,6 @@ Das Filter beginnt mit dem Anfangszustand für $k=0$ {P_{k|k}}=(I-{K_{k}}{H}){P_{k|k-1}} \] -\item Die Outputs von $k$ werden die Inputs für ${k-1}$ und werden wieder im Schritt 1 verwendet +\item Die Outputs von $k$ werden die Inputs für ${k-1}$ und werden wieder in Schritt 1 verwendet \end{itemize} |