From b2cb974151f897ba92e0f7fbc153b6aa96a30176 Mon Sep 17 00:00:00 2001 From: Roy Seitz Date: Sat, 11 Sep 2021 22:11:07 +0200 Subject: =?UTF-8?q?=C3=9Cberarbeitung=20Teil=203.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- buch/papers/erdbeben/Gausskurve3.pdf | Bin 15413 -> 15762 bytes buch/papers/erdbeben/Gausskurve3.tex | 6 ++-- buch/papers/erdbeben/Teil_Fabio.tex | 22 ++++++++----- buch/papers/erdbeben/references.bib | 4 +-- buch/papers/erdbeben/teil0.tex | 26 +++++++-------- buch/papers/erdbeben/teil1.tex | 62 +++++++++++++++++------------------ 6 files changed, 61 insertions(+), 59 deletions(-) (limited to 'buch/papers') diff --git a/buch/papers/erdbeben/Gausskurve3.pdf b/buch/papers/erdbeben/Gausskurve3.pdf index b86023f..4f81c42 100644 Binary files a/buch/papers/erdbeben/Gausskurve3.pdf and b/buch/papers/erdbeben/Gausskurve3.pdf differ diff --git a/buch/papers/erdbeben/Gausskurve3.tex b/buch/papers/erdbeben/Gausskurve3.tex index 032d6de..6e6207b 100644 --- a/buch/papers/erdbeben/Gausskurve3.tex +++ b/buch/papers/erdbeben/Gausskurve3.tex @@ -10,8 +10,8 @@ \begin{axis}[ - xmin = -1, xmax = 4, - ymin = -0.5, ymax = 2.50, + xmin = -1.5, xmax = 3.75, + ymin = -0.25, ymax = 1.75, axis lines = center, xlabel = $\sigma$, ylabel = {$\mu$}, @@ -43,4 +43,4 @@ \end{tikzpicture} -\end{document} \ No newline at end of file +\end{document} diff --git a/buch/papers/erdbeben/Teil_Fabio.tex b/buch/papers/erdbeben/Teil_Fabio.tex index b51d1f0..76dfaa7 100644 --- a/buch/papers/erdbeben/Teil_Fabio.tex +++ b/buch/papers/erdbeben/Teil_Fabio.tex @@ -22,7 +22,7 @@ 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 +Für die Frequenz $f$ wählen wir eine Zufallssquenz mit Erwartungswert und Standardabweichung \begin{equation} \mu = \SI{15}{\hertz} \qquad\text{und}\qquad @@ -64,12 +64,16 @@ Für das Prozessrauschen werden die Bedingungen \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$. +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} Damit sind nun die benötigten Systemparameter und das Rauschen definiert. Als nächstes erzeugen wir ein Erdbeben und schauen, @@ -82,7 +86,7 @@ Leiten wir die Position einmal ab, erhalten wir die Geschwindigkeit. Die zweite Ableitung ergibt uns die Kraft, die in unserer Aufgabenstellung gesucht ist. Zoomen wir näher ran, erkennen wir wieder im Positions-Diagramm eine Überlagerung der Massen-Eigenschwingung mit der Erdbebenschwingung. -Die Masse schwingt mit einer tiefer Frequenz und hoher Amplitude, hingegen das Erdbeben mit einer hohen Frequenz und tiefer Amplitude. +Die Masse schwingt mit einer tiefen 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 wiederzugeben. @@ -96,7 +100,7 @@ Das Filter war imstande die Eigenfrequenz zu eliminieren und die tatsächliche K 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. + In der Kraft ist die Eigendynamik nicht mehr ersichtlich. Unser Filter funktioniert. } \label{erdbeben:fig:standard-alles} \end{center} @@ -107,7 +111,7 @@ Wir möchten nun testen, was die Auswirkungen sind, wenn zum Beispiel der Seismo 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 \[ -m = 0.05 +m = 0.05, \qquad \qquad D = 0.5 \qquad \text{und} \qquad @@ -139,7 +143,7 @@ Dieses Rauschen beeinflusst die Varianzen der Position und Geschwindigkeit in de 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. +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} diff --git a/buch/papers/erdbeben/references.bib b/buch/papers/erdbeben/references.bib index 9bcf57d..c0a211e 100644 --- a/buch/papers/erdbeben/references.bib +++ b/buch/papers/erdbeben/references.bib @@ -6,8 +6,6 @@ %% Saved with string encoding Unicode (UTF-8) - - @article{erdbeben:aragher_understanding_2012, author = {Faragher, Ramsey}, date-added = {2021-07-17 16:44:00 +0200}, @@ -28,7 +26,7 @@ title = {Kalmanfilter}, urldate = {2021-07-0}} -@article{erdbeben:skript:wrstat, +@online{erdbeben:skript:wrstat, author = {Andreas Müller}, url = {https://github.com/AndreasFMueller/WrStat}, date = {2021-09-11}, diff --git a/buch/papers/erdbeben/teil0.tex b/buch/papers/erdbeben/teil0.tex index a3fa6a5..0fd6173 100644 --- a/buch/papers/erdbeben/teil0.tex +++ b/buch/papers/erdbeben/teil0.tex @@ -41,7 +41,7 @@ bei denen die System-Matrix $A$ durch die Jacobi-Matrix ersetzt wird. \begin{figure} \begin{center} - \includegraphics[width=5cm]{papers/erdbeben/Apperatur} + \includegraphics[width=\linewidth,keepaspectratio]{papers/erdbeben/Apperatur} \caption{Aufbau des Seismographen mit Gehäuse, Masse, Federn und Sensor} \label{erdbeben:Seismograph} \end{center} @@ -54,9 +54,9 @@ welche auf das Gehäuse wirkt, bestimmen. 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. +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. -Die Messung entspricht also dem zweiten Integral der Kraft $f(t)$, +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. @@ -68,48 +68,48 @@ 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. +m\ddot s + 2k \dot s + Ds = F. \end{equation} wobei $m$ die Masse, $k$ die Dämpfungskonstante und $D$ die Federkonstante bezeichnet. 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 +Wir substituieren $\dot s = v$ für die Geschwindigkeit und erhalten das Gleichungssystem \begin{align} \begin{split} \dot s &= v \\ - \dot v &= -\frac{D}{m} {s} -\frac{2k}{m} {v} + \frac{f} {m}. + \dot v &= -\frac{D}{m} {s} -\frac{2k}{m} {v} + \frac{F} {m}. \end{split} \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$. +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: +Deshalb nehmen wir $F$ als dritte Grösse in den Zustandsvektor auf und definieren: \[ - x = \begin{pmatrix} {s} \\ {v} \\ {f} \end{pmatrix} + 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 $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$. +Ü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. Wir haben nun alles für die Matrix-Form von Gleichung~\eqref{erdbeben:Systemgleichung} zusammen. Sie lautet: \begin{equation} - \frac{d}{dt} \begin{pmatrix} s(t) \\ v(t) \\ f(t) \end{pmatrix} + \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}. + \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 1a893dd..7b58028 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -14,12 +14,13 @@ 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. +Das ist genau das, was Kalman-Filter tun: +Anhand 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. +Kalman-Filter wurde 1960 von Rudolf Emil Kalman erfunden und direkt von der NASA für die Apollo 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. +Dies kommt heutzutage in jedem Satellit, Navigationssystem, Smartphone und Videospiel vor. Kalman-Filter funktionieren nach folgendem Zwei-Schritt-Verfahren: Zuerst wird, @@ -48,7 +49,7 @@ Jedoch kennen wir auch die Wahrscheinlichkeiten der beiden Aussagen. \begin{figure} \begin{center} - \includegraphics[width=5cm]{papers/erdbeben/Gausskurve3.pdf} + \includegraphics[width=.5\linewidth,keepaspectratio]{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. @@ -73,7 +74,7 @@ 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: +$\odot$~beschreibt dabei die Multiplikation und die Normierung auf den Flächeninhalt eins: \begin{align*} {y_f}(x; {\mu_f}, {\sigma_f}) &= @@ -94,26 +95,26 @@ für den neuen Mittelwert und \] für die Varianz. -Interessant daran ist, dass sich die fusionierte Kurve der genauere Normal-Verteilung anpasst. +Interessant daran ist, dass sich die fusionierte Kurve der genaueren Normalverteilung 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. +Diese neue Funktion ist die bestmö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. +Dieses Prinzip mach sich das Kalman Filter zu nutze, und wird von uns für die Erdbebenberechnung 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 das Kalman Filter zu starten, müssen gewisse Bedingungen definiert werden. +wird dies 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. +Wir verwenden die Standardnotation, wie sie auch im Artikel~\cite{erdbeben:wikipedia} zu finden ist. Sie ist an die Notation der bedingten Wahrscheinlichkeiten angelehnt. Hierbei steht der betrachtete Zeitschritt links und der gegenwärtige rechts eines Vertikalstrichs. Dies bedeutet, dass die Notation $x_{n|m}$ die Schätzung von $x$ zum Zeitpunkt $n$ @@ -135,14 +136,14 @@ Damit haben wir die Systemdynamik nun in der für unser Kalman-Filter notwendige 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$. +Im mehrdimensionalen Fall übernimmt dies die Kovarianzmatrix $P$. Sie wird in jedem Schritt aktualisiert. Hinzu kommt die Prozessunsicherheit $Q$, welche als Parameter in unser Modell einfliesst. $Q$ beschreibt Unsicherheiten im Modell, wie etwa unsere Annahme, dass die Kraft sich nicht ändert, aber auch nicht-modellierbare Einflüsse wie Vibrationen. $P$ wird dabei laufend aktuallisiert. -Die optimale Gleichung lautet +Die Gleichung lautet \[ {P_{k|k-1}}=\Phi {P_{k-1|k-1}} {\Phi _{k}}^T + {Q_{k-1}}. \] @@ -150,7 +151,7 @@ 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 die Messwerte $z_k$ für unser Filter. +Der Sensor wurde noch nicht benutzt, 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 \[ @@ -190,13 +191,14 @@ Der ganze Algorithmus ist nun vollständig und beginnt wieder mit der Vorhersage \subsection{Parameter und Anfangsbedingungen} Die Grössen $P$, $Q$, $R$ und $\Phi$ können grundsätzlich in jedem Zeitschritt ändern. Für die meisten Anwendungen sind sie jedoch konstant und fliessen als Parameter ins Modell ein. +So auch in unserem Falle. Aufgrund der iterativen Arbeitsweise von Kalman-Filtern benötigen wir zudem ein paar Anfangswerte. \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) \] +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 @@ -217,14 +219,13 @@ 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( -\begin{array}{ccc} +P_0 = +\begin{pmatrix} 0 & 0 &0 \\ 0 &0 & 0 \\ 0 & 0 &0 \\ -\end{array} -\right). +\end{pmatrix} +. \] @@ -235,12 +236,12 @@ Bei unserem Modell könnte das beispielsweise ein Windstoss an die Masse sein oder auch Ungenauigkeiten im Modell, wie die Annahme, dass sich die Kraft nicht ändert. Für uns wäre dies: \[ -Q = \left( - \begin{array}{ccc} -{\sigma_s }^2& 0& 0 \\ -0 & {\sigma_v }^2& 0\\ -0 & 0& {\sigma_f }^2\\ -\end{array}\right) +Q = + \begin{pmatrix} +\sigma_s^2& 0& 0 \\ +0 & \sigma_v ^2& 0\\ +0 & 0& \sigma_f^2\\ +\end{pmatrix} . \] Die Standabweichungen müssten statistisch ermittelt werden, da der Fehler nicht vom Sensor kommt und somit nicht vom Hersteller gegeben ist. @@ -249,25 +250,24 @@ Die Messmatrix gibt an, welche Zustände gemessen werden. $H$ ist die Matrix, welche aus der Vorhersage des Zustand eine Vorhersage der Messung erzeugt. In unserem Falle messen wir nur die Position der Massen und verwenden deshalb \[ -H = (1, 0, 0) +H = \begin{pmatrix} 1 & 0 & 0 \end{pmatrix} . \] \subsubsection*{Messrauschkovarianz $R$} 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. -\clearpage \subsection{Zusammenfassung } Das Filter beginnt mit dem Anfangszustand für $k=0$. Anschliessend werden folgende Schritte iterativ ausgeführt: \begin{enumerate} -\item Nächster Zustand vorhersagen +\item Nächsten Zustand vorhersagen \[ {x_{k|k-1}}=\Phi{x_{k-1|k-1}}= \exp(A\Delta t){x_{k-1|k-1}} \] -- cgit v1.2.1