diff options
Diffstat (limited to 'buch/papers/erdbeben')
-rw-r--r-- | buch/papers/erdbeben/Teil_Fabio.tex | 90 | ||||
-rw-r--r-- | buch/papers/erdbeben/images/Makefile | 15 | ||||
-rw-r--r-- | buch/papers/erdbeben/images/messrauschen_geaendert.pdf | bin | 0 -> 481138 bytes | |||
-rw-r--r-- | buch/papers/erdbeben/images/messrauschen_geaendert.tex | 61 | ||||
-rw-r--r-- | buch/papers/erdbeben/images/prozessrauschen_geaendert.pdf | bin | 0 -> 319822 bytes | |||
-rw-r--r-- | buch/papers/erdbeben/images/prozessrauschen_geaendert.tex | 61 | ||||
-rw-r--r-- | buch/papers/erdbeben/images/standard.pdf | bin | 0 -> 362347 bytes | |||
-rw-r--r-- | buch/papers/erdbeben/images/standard.tex | 64 | ||||
-rw-r--r-- | buch/papers/erdbeben/teil0.tex | 23 | ||||
-rw-r--r-- | buch/papers/erdbeben/teil1.tex | 183 |
10 files changed, 324 insertions, 173 deletions
diff --git a/buch/papers/erdbeben/Teil_Fabio.tex b/buch/papers/erdbeben/Teil_Fabio.tex index afb2110..fa7d836 100644 --- a/buch/papers/erdbeben/Teil_Fabio.tex +++ b/buch/papers/erdbeben/Teil_Fabio.tex @@ -10,64 +10,40 @@ und weil es digital simuliert wird, haben wir keine Bauschäden zu beklagen. \subsection{Wahl der Schwingung} Wir müssen uns überlegen, mit welcher Schwingung wir ein realitätsnahes Beben erzeugen können. - Mit einer ungedämpften harmonischen Schwingung können wir zwar die meisten Vorgänge in der Physik erklären. Da aber unser Erdbeben irgendwann abklingen muss, wählen wir die gedämpfte harmonische Schwingung. Die dazugehörige Schwingungsgleichung lautet - -\begin{equation} - y = A e^{-\lambda t} \sin(\omega t) -\end{equation} - -Für die Variablen der harmonisch gedämpften Schwingung setzen wir die Werte - \begin{equation} -A = 5 + y = A e^{-\lambda t} \sin(\omega t). \end{equation} - -ein. - -$A$ ist die Amplitude der Schwingung, die uns die Heftigkeit des Erdebebens beschreibt. +Dabei ist $A=5$ die anfängliche Amplitude der Schwingung, +die uns die Heftigkeit des Erdebebens beschreibt. Sie ist vergleichbar mit der Magnitude. - -$\omega$ definiert sich durch - +$\lambda$ bezeichnet die Bodendämpfung, für die wir $0.2$ wählen. +Sie ist dafür verantwortlich, dass unser Erdbeben abklingt +und kreiert die bei gedämpften Schwingungen typische Hüllkurve. +Wir nehmen an, dass $\lambda$ ein Materialparameter von geologischen Böden ist. +Die Kreisfrequenz $\omega$ ist durch \begin{equation} \omega = 2 \pi f \end{equation} - -wobei die Frequenz $f$ mit - +gegeben, +wobei die Momentanfrequenz $f = \mathcal N(\mu_f, \sigma_f) $ einer Normalverteilung mit \begin{equation} - f = E(\mathrm{Frequenz}) + \sigma^2(\mathrm{Frequenz}) + \mu_f = \SI{15}{\hertz} + \qquad \text{und} \qquad + \sigma_f = \SI{10}{\hertz} \end{equation} +folgt. -erzeugt wird. - -Zusätzlich haben wir $f$ mit dem Savitzky-Golay-Filter gefiltert. +Zusätzlich haben wir $f$ mit einem Savitzky-Golay-Filter gefiltert. Das Savitzky-Golay-Filter schaut sich immer eine definierte Anzahl von Datenpunkte an -und bildet ein Polynom $n$-ter Ordnung. -In unserer Anwendung schaut sich das Filter, im Sinne eines verschieblichen Fensters, -jeweils zehn aufeinanderfolgende Datenpunkte an und bildet ein Polynom $0$-ter Ordnung. -Da wir den Grad $0$ gewählt haben, erhalten wir pro zehn Punkte eine Gerade mit der Steigung $0$. -Diese Art von der Filterung nennt sich gleitender Mittelwert. - -Für den Erwartungswert und die Standardabweichung setzen wir die Zahlen - -\begin{equation} -E(f) = \SI{15}{\hertz} -\end{equation} - -und -\begin{equation} -\sigma^2 = \SI{10}{\hertz} -\end{equation} - -ein. - -$\lambda$ ist die Bodendämpfung, für die wir $0.2$ wählen. -Sie ist dafür verantwortlich, dass unser Erdbeben abklingen wird und kreiert bei der gedämpften Schwingung die typische Hüllkurve der Amplitude. -Wir nehmen an, dass $\lambda$ ein Materialparameter von geologischen Böden ist. +und bildet darüber ein Polynom $n$-ter Ordnung. +In unserer Anwendung schaut sich das Filter, im Sinne eines verschiebbaren Fensters, +jeweils elf aufeinanderfolgende Datenpunkte an +und approximiert diese mit ein Polynom $0$-ter Ordnung, +also einer Konstanten. +Somit erhalten wir mit Matlab-Standardfunktionen einen gleitenden Mittelwert. \subsection{Versuch im Standardfall} Im nächsten Schritt müssen wir sinnvolle Systemparameter für unseren Seismographen definieren. @@ -90,14 +66,14 @@ Wir nehmen an, dass {0.00001}^2& 0& 0 \\ 0 & {0.00001}^2& 0\\ 0 & 0& {1 }^2\\ - \end{array}\right) + \end{array}\right). \end{equation} Auch für die Messung setzen wir ein Rauschen voraus und definieren \begin{equation} R= ({\sigma_x}^2)= -({0.00001}^2) +({0.00001}^2). \end{equation} Sind nun die benötigten Systemparameter und das Rauschen definiert, erzeugen wir das Erdbeben und schauen, wie gut das Kalman-Filter die äussere Beschleunigung schätzen kann. @@ -131,21 +107,25 @@ Das Filter war imstande die Eigenfrequenz zu eliminieren und die tatsächliche K \end{figure} \subsection{Veränderung der Systemparameter} -Was wir nun austesten möchten, sind die Auswirkungen wenn z.B. der Seismograph andere Systemparameter aufweist. -Wir nehmen an, dass sich im Vergleich zum Standardfall die Masse erhöht, die Federkonstante schwächer und die Bodendämpfung doppelt so stark wirkt. +Was wir nun testen möchten, sind die Auswirkungen wenn zum Beispiel der Seismograph andere Systemparameter aufweist. +Wir nehmen an, dass sich im Vergleich zum Standardfall die Masse erhöht, die Federkonstante schwächer und die Federdämpfung doppelt so stark wirkt. Somit gilt neu \[ -m = 0.05 -\qquad \qquad +m = 0.05, +\qquad D = 0.5 \qquad \text{und} \qquad k = 0.02. \] -Da wir mit dieser Anpassung die Trägheit des Seismogrammes erhöht haben, erwarten wir sicher eine langsamere Bewegung der Masse, das heisst die Frequenz wird sich reduzieren. +Da wir mit dieser Anpassung die Trägheit des Seismogrammes erhöht haben, +erwarten wir eine langsamere Bewegung der Masse, +das heisst die Frequenz wird kleiner. -Betrachten wir die Abbildung~\ref{erdbeben:fig:systemparameter-geaendert} können wir diese Erwartung bestätigen. -Nebst dem bemerken wir eine grössere Auslenkung der Position, die wir auf die höhere Energie der Masse und geringeren Rücklenkkraft der Feder begründen können. +Betrachten wir Abbildung~\ref{erdbeben:fig:systemparameter-geaendert}, +können wir diese Erwartung bestätigen. +Zudem bemerken wir eine grössere Auslenkung der Position, +die wir mir durch die höhere Energie der Masse und geringeren Rücklenkkraft der Feder erklären können. \begin{figure} \begin{center} @@ -181,7 +161,7 @@ Die Auswertung in Abbildung~\ref{erdbeben:fig:prozessrauschen-geaendert} zeigt a \subsection{Verstärkung des Messrauschens} Als letztes verstärken wir das Messrauschen um den Faktor $100$ und belassen wieder den Rest wie im Standardfall. Wie man eigentlich schon erwarten kann, zeigt uns die Abbildung~\ref{erdbeben:fig:messrauschen-geaendert}, dass das Signal des Messsensors vom Messrauschen gestört wird. -Weil die Messung somit ungenau wird, kann das Kalman-Filter nicht mehr genau arbeiten und produziert einen ungenauen Output. +Weil die Messung somit ungenau wird, kann das Kalman-Filter nicht mehr genau arbeiten und produziert eine ungenaues Resultat. \begin{figure} \begin{center} diff --git a/buch/papers/erdbeben/images/Makefile b/buch/papers/erdbeben/images/Makefile new file mode 100644 index 0000000..3c82cb7 --- /dev/null +++ b/buch/papers/erdbeben/images/Makefile @@ -0,0 +1,15 @@ +# +# Makefile +# +# (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule +# +all: standard.pdf messrauschen_geaendert.pdf prozessrauschen_geaendert.pdf + +standard.pdf: standard.tex + pdflatex standard.tex + +messrauschen_geaendert.pdf: messrauschen_geaendert.tex + pdflatex messrauschen_geaendert.tex + +prozessrauschen_geaendert.pdf: prozessrauschen_geaendert.tex + pdflatex prozessrauschen_geaendert.tex diff --git a/buch/papers/erdbeben/images/messrauschen_geaendert.pdf b/buch/papers/erdbeben/images/messrauschen_geaendert.pdf Binary files differnew file mode 100644 index 0000000..5056dd0 --- /dev/null +++ b/buch/papers/erdbeben/images/messrauschen_geaendert.pdf diff --git a/buch/papers/erdbeben/images/messrauschen_geaendert.tex b/buch/papers/erdbeben/images/messrauschen_geaendert.tex new file mode 100644 index 0000000..a5a7509 --- /dev/null +++ b/buch/papers/erdbeben/images/messrauschen_geaendert.tex @@ -0,0 +1,61 @@ +% +% messrauschen_geaendert.tex -- Kombination der Messrauschen-Bilder +% +% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule +% +\documentclass[tikz]{standalone} +\usepackage{amsmath} +\usepackage{times} +\usepackage{txfonts} +\usepackage{pgfplots} +\usepackage{csvsimple} +\usetikzlibrary{arrows,intersections,math} +\begin{document} +\newboolean{showgrid} +\setboolean{showgrid}{false} +\def\skala{1} +\begin{tikzpicture}[>=latex,thick,scale=\skala] + +\node at (0,0) {\includegraphics[width=14cm]{../Messrauschen_geaendert.PNG}}; + +\def\yten{-2.06} +\def\yfifteen{0.16} + +\def\links{12.2} +\def\rechts{14.2} + +\pgfmathparse{(\yfifteen-(\yten))/5} +\xdef\m{\pgfmathresult} +\xdef\b{\yten} + +\pgfmathparse{\m*(\links-10)+(\b)} +\xdef\Links{\pgfmathresult} + +\pgfmathparse{\m*(\rechts-10)+(\b)} +\xdef\Rechts{\pgfmathresult} + +\begin{scope}[yshift=-9cm] +\node at (0,0) {\includegraphics[width=14cm]{../Messrauschen_geaendert_zoom.PNG}}; +\end{scope} + +\def\breite{7} +\def\hoehe{2} + +% Gitter +\ifthenelse{\boolean{showgrid}}{ +\draw[step=0.1,line width=0.1pt] (-\breite,-\hoehe) grid (\breite, \hoehe); +\draw[step=0.5,line width=0.4pt] (-\breite,-\hoehe) grid (\breite, \hoehe); +\draw (-\breite,-\hoehe) grid (\breite, \hoehe); +\fill (0,0) circle[radius=0.05]; +}{} + +\def\unten{-4.15} + +\draw[line width=0.7pt] (\Links,\unten) rectangle (\Rechts,4.15); +\draw[line width=0.7pt] (\Links,\unten) -- (\Links,{\unten-0.05}) + -- (-6.45,-4.6) -- (-6.45,-4.9); +\draw[line width=0.7pt] (\Rechts,\unten) -- (\Rechts,{\unten-0.05}) + -- (6.86,-4.6) -- (6.86,-4.9); +\end{tikzpicture} +\end{document} + diff --git a/buch/papers/erdbeben/images/prozessrauschen_geaendert.pdf b/buch/papers/erdbeben/images/prozessrauschen_geaendert.pdf Binary files differnew file mode 100644 index 0000000..e0bf605 --- /dev/null +++ b/buch/papers/erdbeben/images/prozessrauschen_geaendert.pdf diff --git a/buch/papers/erdbeben/images/prozessrauschen_geaendert.tex b/buch/papers/erdbeben/images/prozessrauschen_geaendert.tex new file mode 100644 index 0000000..ad4dcf9 --- /dev/null +++ b/buch/papers/erdbeben/images/prozessrauschen_geaendert.tex @@ -0,0 +1,61 @@ +% +% prozessrauschen_geaendert.tex -- Kombination der Prozessrauschen-Bilder +% +% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule +% +\documentclass[tikz]{standalone} +\usepackage{amsmath} +\usepackage{times} +\usepackage{txfonts} +\usepackage{pgfplots} +\usepackage{csvsimple} +\usetikzlibrary{arrows,intersections,math} +\begin{document} +\newboolean{showgrid} +\setboolean{showgrid}{false} +\def\skala{1} +\begin{tikzpicture}[>=latex,thick,scale=\skala] + +\node at (0,0) {\includegraphics[width=14cm]{../Prozessrauschen_geaendert.PNG}}; + +\def\yten{-2.1} +\def\yfifteen{0.12} + +\def\links{13.27} +\def\rechts{14.2} + +\pgfmathparse{(\yfifteen-(\yten))/5} +\xdef\m{\pgfmathresult} +\xdef\b{\yten} + +\pgfmathparse{\m*(\links-10)+(\b)} +\xdef\Links{\pgfmathresult} + +\pgfmathparse{\m*(\rechts-10)+(\b)} +\xdef\Rechts{\pgfmathresult} + +\begin{scope}[yshift=-9cm] +\node at (0,0) {\includegraphics[width=14cm]{../Prozessrauschen_geaendert_zoom.PNG}}; +\end{scope} + +% Gitter +\def\breite{7} +\def\hoehe{2} +\ifthenelse{\boolean{showgrid}}{ +\draw[step=0.1,line width=0.1pt] (-\breite,-\hoehe) grid (\breite, \hoehe); +\draw[step=0.5,line width=0.4pt] (-\breite,-\hoehe) grid (\breite, \hoehe); +\draw (-\breite,-\hoehe) grid (\breite, \hoehe); +\fill (0,0) circle[radius=0.05]; +}{} + +\def\unten{-4.15} + +\draw[line width=0.7pt] (\Links,\unten) rectangle (\Rechts,4.18); +\draw[line width=0.7pt] (\Links,\unten) -- (\Links,{\unten-0.05}) + -- (-6.62,-4.6) -- (-6.62,-4.9); +\draw[line width=0.7pt] (\Rechts,\unten) -- (\Rechts,{\unten-0.05}) + -- (6.80,-4.6) -- (6.80,-4.9); + +\end{tikzpicture} +\end{document} + diff --git a/buch/papers/erdbeben/images/standard.pdf b/buch/papers/erdbeben/images/standard.pdf Binary files differnew file mode 100644 index 0000000..f2ca85e --- /dev/null +++ b/buch/papers/erdbeben/images/standard.pdf diff --git a/buch/papers/erdbeben/images/standard.tex b/buch/papers/erdbeben/images/standard.tex new file mode 100644 index 0000000..74ac7a1 --- /dev/null +++ b/buch/papers/erdbeben/images/standard.tex @@ -0,0 +1,64 @@ +% +% standard.tex -- Kombination der Standard-Bilder +% +% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule +% +\documentclass[tikz]{standalone} +\usepackage{amsmath} +\usepackage{times} +\usepackage{txfonts} +\usepackage{pgfplots} +\usepackage{csvsimple} +\usetikzlibrary{arrows,intersections,math} +\begin{document} +\newboolean{showgrid} +\setboolean{showgrid}{false} +\def\skala{1} +\begin{tikzpicture}[>=latex,thick,scale=\skala] + +\node at (0,0) {\includegraphics[width=14cm]{../Standard_alles.PNG}}; + +\def\yten{-2.12} +\def\yfifteen{0.1} + +\def\links{12.5} +\def\rechts{13.4} + +\pgfmathparse{(\yfifteen-(\yten))/5} +\xdef\m{\pgfmathresult} +\xdef\b{\yten} + +%\node at (0,7) {$m=\m$}; \node at (0,8) {$b=\b$}; + +\pgfmathparse{\m*(\links-10)+(\b)} +\xdef\Links{\pgfmathresult} + +\pgfmathparse{\m*(\rechts-10)+(\b)} +\xdef\Rechts{\pgfmathresult} + +\begin{scope}[yshift=-9cm] +\node at (0,0) {\includegraphics[width=14cm]{../Standard_Zoom.PNG}}; +\end{scope} + +\def\breite{7} +\def\hoehe{2} + +% Gitter +\ifthenelse{\boolean{showgrid}}{ +\draw[step=0.1,line width=0.1pt] (-\breite,-\hoehe) grid (\breite, \hoehe); +\draw[step=0.5,line width=0.4pt] (-\breite,-\hoehe) grid (\breite, \hoehe); +\draw (-\breite,-\hoehe) grid (\breite, \hoehe); +\fill (0,0) circle[radius=0.05]; +}{} + +\def\unten{-4.15} + +\draw[line width=0.7pt] (\Links,\unten) rectangle (\Rechts,4.18); +\draw[line width=0.7pt] (\Links,\unten) -- (\Links,{\unten-0.05}) + -- (-6.36,-4.6) -- (-6.36,-4.9); +\draw[line width=0.7pt] (\Rechts,\unten) -- (\Rechts,{\unten-0.05}) + -- (6.87,-4.6) -- (6.87,-4.9); + +\end{tikzpicture} +\end{document} + diff --git a/buch/papers/erdbeben/teil0.tex b/buch/papers/erdbeben/teil0.tex index d32b316..f012225 100644 --- a/buch/papers/erdbeben/teil0.tex +++ b/buch/papers/erdbeben/teil0.tex @@ -22,11 +22,19 @@ 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. + + \begin{figure} \begin{center} \includegraphics[width=5cm]{papers/erdbeben/Apperatur} @@ -50,6 +58,7 @@ 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: \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. @@ -66,25 +75,19 @@ Somit entstehen die Gleichungen für die Geschwindigkeit $ \dot s_1(t)$ der Mass 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. -Diese können wir nun in der Form -\[ \ddot f =-\frac{D}{m} {s_1} -\frac{2k}{m} {s_2} + \frac{f} {m} \] -als skalare Gleichung darstellen. 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 (20.1) gerade das ist, was wir eigentlich bestimmen möchten. +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. Deshalb nehmen wir $f$ als dritte Grösse in den Zustandsvektor auf und definieren: \[ - x = (s_1, s_2, f)^T. +x= \left( \begin{array}{c} {s_1}\\ {s_2}\\{f}\end{array}\right)^T \] -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 (20.1), aber über $\dot f$ wissen wir nichts. -Wir müssen also eine Annahme treffen: $\dot f = 0$. Diese Annahme ist im Allgemeinen falsch, aber etwas Besseres haben wir zurzeit nicht zur Verfügung. -Zudem treffen wir die Annahme, das sich die Kraft über die Beobachtungszeit nicht verändert. +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. Wir werden dies in einem späteren Schritt kompensieren müssen. -Da die Kraft unbekannt ist, wird die letzte Zeile mit Nullen gefüllt, denn genau diese Werte wollen wir. - 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: \begin{equation} diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index 014b53e..4d6786f 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -20,38 +20,30 @@ Da wir die äussere Kraft nicht direkt messen können, benötigen wir ein Werkze 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. +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. -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. +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. +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. +\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: Zwei Normalverteilungen} sind zwei Funktionen dargestellt. +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: 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. -\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} - \end{center} -\end{figure} +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. + 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. @@ -67,14 +59,17 @@ und der Messung: \] 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 +Durch geschicktes umformen resultiert aus der Kombination der beiden Verteilungen wiederum in einer Normalverteilung mit Erwartungswert \[ \mu_f = \frac{\mu_1\sigma_2^2 + \mu_2 \sigma_1^2}{\sigma_1^2 + \sigma_2^2} \] und Varianz @@ -83,30 +78,27 @@ und Varianz \] 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! +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 und der orangen Verteilung entsteht die die rote, optimale Funktion} + \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} + 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} +\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 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. -\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. @@ -114,79 +106,77 @@ Dies bedeutet, dass die Notation $x_{n|m}$ die Schätzung von $x$ zum Zeitpunkt \subsubsection*{Vorhersage} Im Filterschritt Vorhersage wird der nächste Zustand anhand des Anfangszustand und der Systemmatrix berechnet. -Dies funktioniert mit dem Rechenschritt: +Die Systemmatrix $A$ 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 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}$ 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 + +Im Abschnitt ~\ref{erdbeben:Wahrscheindlichkeit} benötigten wir die Varianzen der Normalverteilungen. +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 \[ {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. +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 den aktuellen Messwert $z$ mittels der Systemmatrix $\Phi$ beschreiben kann. +Die Innovation ist der 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. +Mit dessen Hilfe erhalten wir die optimale Schätzung des nächsten Zustandes \[ -{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}} \] -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}}. \] -\subsection{Anfangsbedingungen} +\subsection{Parameter und Anfangsbedingungen} \subsubsection*{Anfangszustand $x$} + Das Filter benötigt eine Anfangsbedingung. 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*{Anfangsfehler / Kovarianzmatrix $P$} Da auch der Anfangszustand fehlerhaft sein kann, wird für das Filter ein Anfangsfehler 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) - 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 @@ -200,30 +190,9 @@ 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. @@ -258,43 +227,41 @@ 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. \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$ + +\begin{itemize} +\item Nächster Zustand vorhersagen \[ {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}}. \] -3. Zustand wird gemessen +\item Innovation (= Messung - Vorhersage) \[ {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 Das Kalman Filter anwenden +\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 +\item Die Outputs von $k$ werden die Inputs für ${k-1}$ und werden wieder im Schritt 1 verwendet +\end{itemize} |