diff options
Diffstat (limited to 'buch/papers')
| -rw-r--r-- | buch/papers/erdbeben/Teil_Fabio.tex | 229 | ||||
| -rw-r--r-- | buch/papers/erdbeben/teil0.tex | 41 | ||||
| -rw-r--r-- | buch/papers/erdbeben/teil1.tex | 62 |
3 files changed, 170 insertions, 162 deletions
diff --git a/buch/papers/erdbeben/Teil_Fabio.tex b/buch/papers/erdbeben/Teil_Fabio.tex index 48ec054..737c1d1 100644 --- a/buch/papers/erdbeben/Teil_Fabio.tex +++ b/buch/papers/erdbeben/Teil_Fabio.tex @@ -1,106 +1,80 @@ \section{Anwendung des Kalman-Filters} -\subsection{Ziel} -Bis jetzt haben wir gelesen, was das Kalman-Filter bewirkt und wie es funktioniert. -Nun möchten wir mit einem Beispiel herausfinden, ob das Filter unsere gesuchte Grösse $f(t)$ bestimmen kann. +Bis jetzt haben wir gesehen, was das Kalman-Filter bewirkt und wie es funktioniert. +Nun möchten wir mit einem konkreten Beispiel herausfinden, +ob das Filter unsere gesuchte Grösse $f(t)$ bestimmen kann. -\subsection{Künstliche Erdbebendaten} -Da wir keine Rohdaten über vergangene Erdbeben zur Hand haben, müssen wir mittels Matlab künstliche Daten erzeugen und sie dann in das Filter eingeben. +Da wir keine Rohdaten über vergangene Erdbeben zur Hand haben, +müssen wir mittels Simulation künstliche Daten erzeugen. +Diese können wir dann mit unserem Filter verarbeiten. Diese Vorgehensweise erlaubt uns das Erdbeben beliebig zu gestalten -und weil es digital simuliert wird, haben wir keine Bauschäden zu beklagen. +und weil es digital simuliert wird, haben wir auch 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 -\end{equation} - -ein. - -$A$ ist die Amplitude der Schwingung, die uns die Heftigkeit des Erdebebens beschreibt. -Sie ist vergleichbar mit der Magnitude. - -$\omega$ definiert sich durch - +Da ein Erdbeben vorteilhafterweise irgendwann abklingen sollte, +wählen wir eine gedämpfte harmonische Schwingung \begin{equation} - \omega = 2 \pi f + y = A e^{-\lambda t} \sin(2\pi f t). \end{equation} -wobei die Frequenz $f$ mit +In unsere Simulation können wir die Parameter frei wählen. +Wir setzten $A = 5$ für die Amplitude der Schwingung. +Sie beschreibt die Heftigkeit des Erdebebens und ist vergleichbar mit der Magnitude. +Für die Frequenz $f$ wählen wir eine Zufalls-Sequenz mit Erwartungswert und Standardabweichung \begin{equation} - f = E(\mathrm{Frequenz}) + \sigma^2(\mathrm{Frequenz}) + \mu = \SI{15}{\hertz} + \qquad\text{und}\qquad + \sigma = \SI{10}{\hertz}. \end{equation} - -erzeugt wird. - -Zusätzlich haben wir $f$ mit dem Savitzky-Golay-Filter gefiltert. -Das Savitzky-Golay-Filter schaut sich immer eine definierte Anzahl von Datenpunkte an +Zusätzlich haben wir $f$ mit einem Savitzky-Golay-Filter gefiltert. +Ein 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. +In unserer Anwendung schaut sich das Filter, im Sinne eines verschiebbaren Fensters, +jeweils elf aufeinanderfolgende Datenpunkte an und bildet ein Polynom $0$-ter Ordnung, +also eine Konstante. +Somit erhalten wir mit Matlab-Standardfunktionen einen gleitenden Mittelwert, +um all zu schnelle Änderungen der Frequenz zu unterdrücken. $\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. +Sie ist dafür verantwortlich, dass unser Erdbeben abklingt +und kreiert bei der gedämpften Schwingung die typische Hüllkurve. Wir nehmen an, dass $\lambda$ ein Materialparameter von geologischen Böden ist. \subsection{Versuch im Standardfall} Im nächsten Schritt müssen wir sinnvolle Systemparameter für unseren Seismographen definieren. Eine kurze Recherche zeigt, dass die Masse ein Gewicht von ca.\ \SI{100}{\gram} hat. -Zur Federkonstante D und Dämpfung k konnten wir leider keine brauchbaren Grössen finden und treffen die Annahme, dass $D = 1$ und $k = 0.01$. +Zur Federkonstante $D$ und Dämpfung $k$ konnten wir leider keine brauchbaren Grössen finden. +Wir treffen die Annahmen $D = 1$ und $k = 0.01$. Für die Masse definieren wir $m = 0.01$. Für das Prozessrauschen werden die Bedingungen - \begin{equation} - Q = \left( - \begin{array}{ccc} - {\sigma_x }^2& 0& 0 \\ - 0 & {\sigma_v }^2& 0\\ - 0 & 0& {\sigma_f }^2\\ - \end{array}\right)= \left( - \begin{array}{ccc} - {0.00001}^2& 0& 0 \\ - 0 & {0.00001}^2& 0\\ - 0 & 0& {1 }^2\\ - \end{array}\right) + Q = + \begin{pmatrix} + \sigma_x ^2 & 0 & 0 \\ + 0 & \sigma_v^2 & 0\\ + 0 & 0 & \sigma_f^2 \\ + \end{pmatrix}= + \begin{pmatrix} + 0.00001^2& 0& 0 \\ + 0 & 0.00001^2& 0\\ + 0 & 0& 1^2 \\ + \end{pmatrix} \end{equation} - angesetzt. +Die Annahme, dass sich die Erdbebenkraft $f$ nicht ändert, +kompensieren wir hier endlich durch einen grossen Wert von $\sigma_f^2$. Auch für die Messung setzen wir ein Rauschen voraus und definieren - \begin{equation} -R= ({\sigma_x}^2)= -({0.00001}^2) +R= (\sigma_x^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. +Damit sind nun die benötigten Systemparameter und das Rauschen definiert. +Als nächstes erzeugen wir ein Erdbeben und schauen, +wie gut das Kalman-Filter die äussere Beschleunigung schätzen kann. \subsection*{Ergebnis} @@ -112,18 +86,25 @@ Zoomen wir näher ran, erkennen wir wieder im Positions-Diagramm eine Überlager Die Masse schwingt mit einer tiefer Frequenz und hoher Amplitude, hingegen das Erdbeben mit einer hohen Frequenz und tiefer Amplitude. Vergleichen wir nun die Position mit der Kraft, stellen wir fest, dass das Kalman-Filter eine Schätzung wiedergibt, die auch eine Frequenz von \SI{15}{\hertz} hat. -Das Filter war imstande die Eigenfrequenz zu eliminieren und die tatsächliche Kraft des Erdbebens zu wiedergeben. +Das Filter war imstande die Eigenfrequenz zu eliminieren und die tatsächliche Kraft des Erdbebens wiederzugeben. \begin{figure} - \begin{center} - \includegraphics[width=\linewidth,keepaspectratio]{papers/erdbeben/images/standard.PDF} - \caption{Das Position-Zeit-Diagramm zeigt uns die typische Aufzeichnung eines Seismographen während eines Erdbebens. Um die Geschwindigkeit zu erhalten müssen wir die Position einmal ableiten. Ein weiteres Ableiten erzeugt uns die Beschleunigung, respektive die Kraft. Sehr gut ersichtlich ist die Hüllkurve der Amplitude, wie wir sie bei einer gedämpften Schwingung erwarten. Erst das Vergrössern an die Datenpunkte zeigt uns auf, wie gut die Schätzung des Kalman-Filters funktioniert.} + \begin{center} + \includegraphics[width=.95\linewidth,keepaspectratio]{papers/erdbeben/images/standard.PDF} + \caption{ + Das Position-Zeit-Diagramm zeigt eine typische Aufzeichnung eines Seismographen während eines Erdbebens. + Sehr gut ersichtlich ist die Hüllkurve, wie wir sie bei einer gedämpften Schwingung erwarten. + In der Vergrösserung wird die Überlagerung aus Eigenschwingung und Erdbeben gut ersichtlich. + Die Geschwindigkeit und schliesslich die Kraft weden aus der Position durch unser Kalman-Filter geschätzt. + Erst das Vergrössern an die Datenpunkte zeigt, wie gut die Schätzung des Kalman-Filters funktioniert. + In der Kraft ist die Eigendynamim nicht mehr ersichtlich. Unser Filter funktioniert. + } \label{erdbeben:fig:standard-alles} - \end{center} + \end{center} \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 möchten nun testen, was die Auswirkungen sind, 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 Bodendämpfung doppelt so stark wirkt. Somit gilt neu \[ @@ -133,59 +114,97 @@ 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 eine langsamere Bewegung der Masse, das heisst die Eigenfrequenz wird reduziert. -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. - -Betrachten wir die Abbildung~\ref{erdbeben:fig:systemparameter-geaendert} können wir diese Erwartung bestätigen. +Betrachten wir 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. \begin{figure} - \begin{center} - \includegraphics[width=\linewidth,keepaspectratio]{papers/erdbeben/Systemparameter_geaendert_2.PNG} - \caption{Im Geschwindigkeits-Diagramm erkennen wir, dass sich im Vergleich zum Standardfall, die Auslenkung und Frequenz vergrössert hat. Dies wird mit der Erhöhung der Masse und somit der Trägheit begründet. Auch stellen wir fest, dass die Positionsmessung überwiegend die Eigenfrequenz misst.} + \begin{center} + \includegraphics[width=\linewidth,keepaspectratio]{papers/erdbeben/Systemparameter_geaendert.PNG} + \caption{ + Im Geschwindigkeits-Diagramm erkennen wir, + dass sich im Vergleich zum Standardfall die Auslenkung und Frequenz vergrössert hat. + Dies wird mit der Erhöhung der Masse und somit der Trägheit begründet. + Auch stellen wir fest, dass die Positionsmessung überwiegend die Eigenfrequenz misst. + Diese wird in der Schätzung der Kraft dennoch sehr gut kompensiert. + } \label{erdbeben:fig:systemparameter-geaendert} - \end{center} + \end{center} \end{figure} \subsection{Verstärkung des Prozessrauschens} Falls wir unseren Seismographen in der Nähe einer grösseren Stadt aufstellen, so müssen wir aufgrund der Vibrationen mit einem stärkeren Prozessrauschen rechnen. Dieses Rauschen beeinflusst die Varianzen der Position und Geschwindigkeit in der Matrix $Q$. -Aus diesem Grund erhöhen wir die Standardabweichungen in der Matrix $Q$ um den Faktor $100$. -Die Auswertung in Abbildung~\ref{erdbeben:fig:prozessrauschen-geaendert} zeigt auf, dass das Kalman-Filter die Schätzung der Kraft nur gering an den Messwerten anpasst. -Die Theorie dazu haben wir im Kapitel \nameref{Wahrscheinlichkeit} angeschaut. +Aus diesem Grund erhöhen wir die Standardabweichungen der Positions $\sigma_s$ und Geschwindigkeit $\sigma_v$ in der Matrix $Q$ um den Faktor $100$. +Die Auswertung in Abbildung~\ref{erdbeben:fig:prozessrauschen-geaendert} zeigt auf, +dass das Kalman-Filter die Schätzung der Kraft nur gering an den Messwerten anpasst, +da wir den Schätzungen für die Position nun wenig Vertrauen schenken und stärker der Modell-Annahme $\dot f = 0$ folgen. +Die Theorie dazu haben wir im Abschnitt~\ref{erdbeben:Wahrscheindlichkeit} angeschaut. \begin{figure} - \begin{center} - \includegraphics[width=\linewidth,keepaspectratio]{papers/erdbeben/images/Prozessrauschen_geaendert.PDF} - \caption{Mit dem Erhöhen des Prozessrauschens gehen wir von einer grösseren Unsicherheit der Systemmatrix aus. Aus diesem Grund folgt das Filter vor allem den Messwerten, was sichtbare Folgen für die Schätzkurve im Kraft-Zeit-Diagramm hat. Hier möchte das Filter auch den Messwerten folgen. Da wir aber für die Kraft keine Messwerte aufzeichnen, erhalten wir eine sehr schwache Kurve. Die Position kann immernoch präzise geschätzt werden und die Ableitung zur Geschwindigkeit ergibt gute Resultate. Jedoch ist die Schätzkurve der Kraft sehr weit von der idealen Kurve entfernt und nicht nutzbar.} + \begin{center} + \includegraphics[width=.95\linewidth,keepaspectratio]{papers/erdbeben/images/Prozessrauschen_geaendert.PDF} + \caption{ + Mit dem Erhöhen des Prozessrauschens gehen wir von einer grösseren Unsicherheit der Systemmatrix aus. + Aus diesem Grund folgt das Filter vor allem den Messwerten, + was sichtbare Folgen für die Schätzkurve im Kraft-Zeit-Diagramm hat. + Hier möchte das Filter auch den Messwerten folgen. + Da wir aber für die Kraft keine Messwerte aufzeichnen, + erhalten wir eine sehr schwache Kurve. + Die Position kann immernoch präzise geschätzt werden und die Ableitung zur Geschwindigkeit ergibt gute Resultate. + Jedoch ist die Schätzkurve der Kraft sehr weit von der idealen Kurve entfernt und nicht nutzbar. + } \label{erdbeben:fig:prozessrauschen-geaendert} - \end{center} + \end{center} \end{figure} \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. -Auch hier verweisen wir auf Kapitel \nameref{Wahrscheinlichkeit}. +Wie man eigentlich schon erwarten kann, +zeigt uns die Abbildung~\ref{erdbeben:fig:messrauschen-geaendert}, +dass das Signal des Messsensors vom Messrauschen start gestört wird. +Weil die Messung zu ungenau ist, +kann das Kalman-Filter nicht mehr gut arbeiten und produziert einen ungenauen Output. \begin{figure} - \begin{center} - \includegraphics[width=\linewidth,keepaspectratio]{papers/erdbeben/images/Messrauschen_geaendert.PDF} - \caption{Im Kraft-Zeit-Diagramm erhalten wir nur bis ca. $t = 10$ gute Schätzwerte. Von $t = 10$ bis $t = 30$ wirkt das Messrauschen zu stark und erhalten keine brauchbaren Werte mehr. Im Position-Zeit-Diagramm erhielten wir bis jetzt immer genaue Schätzungen. Mit einem starken Messrauschen fällt es nun dem Filter schwerer, präzise Werte zu generieren. Die Nahaufnahme im Kraft-Zeit-Diagramm bestätigt uns aber, dass die Messfehler zu gross sind, um ein klares Bild über die äussere Kraft zu erhalten.} - \label{erdbeben:fig:messrauschen-geaendert} - \end{center} + \begin{center} + \includegraphics[width=.95\linewidth,keepaspectratio]{papers/erdbeben/images/Messrauschen_geaendert.PDF} + \caption{ + Im Kraft-Zeit-Diagramm erhalten wir nur bis ca.\ $t = 10$ gute Schätzwerte. + Ab $t = 10$ wirkt das Messrauschen zu stark und wir erhalten keine brauchbaren Werte mehr. + Im Position-Zeit-Diagramm erhielten wir bis jetzt immer genaue Schätzungen. + Mit einem starken Messrauschen fällt es dem Filter nun schwerer, + präzise Schätzungen zu berechnen. + Die Nahaufnahme im Kraft-Zeit-Diagramm bestätigt uns, + dass die Messfehler zu gross sind, + um ein klares Bild über die äussere Kraft zu erhalten. + } + \label{erdbeben:fig:messrauschen-geaendert} + \end{center} \end{figure} \subsection{Zusammenfassung} -Wir haben uns zum Ziel gesetzt, die äussere Beschleunigung $a(t)$, bzw. die Kraft $f(t)$ eines Erdbebens zu ermitteln. +Wir haben uns zum Ziel gesetzt, +die äussere Beschleunigung $a(t)$, +beziehungsweise die Kraft $f(t)$ eines Erdbebens +aus den Messugnen eines Seismographen zu berechen. + +Wir haben einen Seismographen mathematisch beschrieben und +mit der Software Matlab Messresultate während eines künstlichen Erdbebens erzeugt. +Diese Messwerte haben wir mit einem Kalman-Filter bearbeitet, +um aus den Messwerten wieder das Erdbeben zu gewinnen. -Mit der Software Matlab haben wir einen virtuellen Seismographen gebaut und ein künstliches Erdbeben erzeugt. -Der Seismograph war fähig die Position der Masse während der Einwirkung des Erdbebens aufzuzeichnen. +Der Seismograph war fähig, die Position der Masse während der Einwirkung des Erdbebens aufzuzeichnen. $a(t)$ kann zwar nicht mit Sensoren gemessen werden, jedoch erhalten wir $a(t)$ durch zweifaches Ableiten. Da wir so aber die innere Beschleunigung erhalten, mussten wir das Kalman-Filter anwenden. -Das Kalman-Filter half uns die äussere Beschleunigung zu schätzen und lieferte erstaunlich genaue Werte. +Das Kalman-Filter half uns, die äussere Beschleunigung zu schätzen, und lieferte erstaunlich genaue Werte. Ausserdem hat es das Filter geschafft, die Eigenfrequenz der Masse und die Erdbebenfrequenz zu separieren. Folglich erhielten wir eine Schätzung, die nur das Erdbeben betraf. -Zuletzt haben wir aufgezeigt, das Veränderungen an den System- und Rauschparametern die Genauigkeit und Zuverlässigkeit des Kalman-Filters beeinträchtigen können. +Zuletzt haben wir aufgezeigt, +dass Veränderungen an den System- und Rauschparametern die Genauigkeit und Zuverlässigkeit +des Kalman-Filters beeinträchtigen können. +Wir haben gesehen, dass aus zu schlechten Sensordaten auch mittels Filterung keine genauen Aussagen möglich sind. diff --git a/buch/papers/erdbeben/teil0.tex b/buch/papers/erdbeben/teil0.tex index 9b33e4b..a3fa6a5 100644 --- a/buch/papers/erdbeben/teil0.tex +++ b/buch/papers/erdbeben/teil0.tex @@ -4,10 +4,7 @@ % (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. \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. @@ -16,18 +13,17 @@ Wenn dies passiert, entlädt sich die aufgebaute Spannung und setzt enorme Energ Ein Erdbeben breitet sich vom Erdbebenherd in allen Richtungen gleich aus. Vergleichbar ist, wenn man einen Stein in einen Teich wirft und die Wellen beobachten kann, die sich ausbreiten. -\subsection{Funktion eines Seismograph} +\section{Funktion eines Seismographen} Um ein Erdbeben kenntlich zu machen, werden in der Regel Seismographen mit vielen Sensoren verwendet. 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. +Eine Relativbewegung des Gehäuses 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. +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. @@ -36,9 +32,8 @@ 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. +Die Systembeschreibung wird dann deutlich komplexer, bringt aber nichts wesentlich Neues hervor. +Wir beschränken uns deshalb auf den linearen Fall. 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, @@ -52,11 +47,10 @@ bei denen die System-Matrix $A$ durch die Jacobi-Matrix ersetzt wird. \end{center} \end{figure} -\subsection{Ziel} -Unser Seismograph misst nur die Position der Masse über die Zeit. -Wir wollen jedoch die Beschleunigung $a(t)$ des Boden, -respektive die Kraft $f(t)$, -welche auf das Gehäuse wirkt, bestimmten. +Unser Seismograph misst jedoch nur die Position der Masse über die Zeit. +Wir wollen aber die Beschleunigung $a(t)$ des Boden, +respektive die Kraft, +welche auf das Gehäuse wirkt, bestimmen. Anhand dieser Beschleunigung, beziehungsweise der Krafteinwirkung durch die Bodenbewegung, wird später das Bauwerk bemessen. @@ -83,8 +77,10 @@ Wir möchten diese Gleichung folglich in die Darstellung $\dot x = Ax$ überfüh 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}. + \begin{split} + \dot s &= v \\ + \dot v &= -\frac{D}{m} {s} -\frac{2k}{m} {v} + \frac{f} {m}. + \end{split} \label{erdbenen:systemgleichungen} \end{align} @@ -117,14 +113,3 @@ Sie lautet: \label{erdbeben:systemmatrix} \end{equation} - - - - - - - - - - - diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index 6c2539a..1a893dd 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -9,8 +9,6 @@ % (c) 2020 Prof Dr Andreas Müller, Hochschule Rapperswil %% - - \rhead{Kalman-Filter} \section{Kalman-Filter} Im letzten Abschnitt haben wir Gleichungen für unser System gefunden. @@ -33,7 +31,7 @@ wobei die Korrektur abhänging von der Differenz zwischen erwarteter und effekti 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. +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. @@ -87,8 +85,8 @@ $\odot$ beschreibt dabei die Multiplikation und die Normierung auf den Flächeni &= \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*} -Die genaue Berechnung ist eine reine Fingerübung. -Nach einigem Aufwand findet man die Ausdrücke +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} \] für den neuen Mittelwert und \[ @@ -109,20 +107,23 @@ Was in zwei Dimensionen erklärt wurde, funktioniert auch in mehreren Dimensione Dieses Prinzip mach sich das Kalman Filter zu nutze, und wird von uns für die Erdbeben Berechnung genutzt. \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. +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. 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. +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. -Die Systemmatrix $A$ aus Gleichung~\eqref{erdbeben:systemmatrix} beschreibt jedoch ein kontinuierliches System $\dot x = Ax$. +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). \] @@ -164,7 +165,7 @@ Entsprechende Korrekturen werden dann gross bzw. nur gering ausfallen. \subsubsection*{Aktualisieren} 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: +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}} \\ @@ -178,7 +179,7 @@ Die optimale Schätzung des neuen Zustandes wird dann zu \] 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 ist nun vollständig und beginnt wieder mit der Vorhersage \[ @@ -198,8 +199,10 @@ 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:Systemgleichung} bereits gefunden. -Zudem haben wir weiter oben bereits entdeckt, wie wir mittels Exponentialfunktion zu einer zeitdiskreten Beschreibung für das Kalman-Filter kommen. +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) .\] @@ -229,7 +232,7 @@ Somit kann die Matrix mit Nullen bestückt werden und wir starten mit 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, dass sich die Kraft nicht ändert. +oder auch Ungenauigkeiten im Modell, wie die Annahme, dass sich die Kraft nicht ändert. Für uns wäre dies: \[ Q = \left( @@ -242,7 +245,7 @@ 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. +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 \[ @@ -253,29 +256,30 @@ 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 } - -Das Filter beginnt mit dem Anfangszustand für $k=0$ - -\begin{itemize} +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}} \] \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}} \] -\item Innovation (= Messung - Vorhersage) +\item Innovation (= Messung - Vorhersage) \[ -{w_{k}}={z_{k}}-{H}{x_{k|k-1}}. +{w_{k}}={z_{k}}-{H}{x_{k|k-1}} \] \item Optimales Kalman-Gain berechnen @@ -294,6 +298,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 in Schritt 1 verwendet -\end{itemize} +\end{enumerate} +Die Outputs von $k$ werden die Inputs für ${k+1}$ und werden wieder in Schritt 1 verwendet. |