From 30eb23bfa08b5119e29efd1e4a3cbfc6ff9d81c7 Mon Sep 17 00:00:00 2001 From: Lukaszogg <82384106+Lukaszogg@users.noreply.github.com> Date: Thu, 8 Jul 2021 19:56:31 +0200 Subject: Update teil1.tex --- buch/papers/erdbeben/teil1.tex | 228 ++++++++++++++++++++++++++++++----------- 1 file changed, 168 insertions(+), 60 deletions(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index 0d21f84..d6f5638 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -11,12 +11,15 @@ +\begin{document} + + \section{Kalman Filter} \subsection{Geschichte} Das Kalman Filter wurde 1960 von Rudolf Emil Kalman entdeckt und direkt von der NASA für die Appollo Mission benutzt. Der Filter kommt mit wenig Rechenleistung aus und war somit dafür geeignet die Rakete bei der Navigation zu unterstützen. Das Filter schätzt den Zustand eines Systems anhand von Messungen und kann den nächsten Zustand errechnen. Typische Anwendungen des Kalman-Filters sind die Glättung von verrauschten Daten und die Schätzung von Parametern und kommt heutzutage in jedem Satellit, Navigationssystem, Smartphones und Videospielen vor. \subsection{Wahrscheinlichkeit} -Das Kalman Filter versucht nichts anderes, als ein geeigneter Wert zwischen zwei Normalverteilungen zu schätzen. Die eine Kurve zeigt die errechnete Vorhersage des Zustands, bzw. deren Normal- Gauss-Verteilung. Die andere Kurve zeigt die verrauschte Messung des nächsten Zustand, bzw. deren Normal-Verteilung. Wie man in am Beispiel dieser zwei Gauss-Verteilungen sehen kann, ist sowohl der geschätzte Zustand als auch der gemessene Zustand nicht am selben Punkt. +Das Kalman-Filter schätzt den wahrscheinlichsten Wert zwischen zwei Normalverteilungen oder auch Gauss-Verteilung. Die eine Kurve zeigt die errechnete Vorhersage des Zustands, bzw. deren Normalverteilung. Die andere Kurve zeigt die verrauschte Messung des nächsten Zustand, bzw. deren Normalverteilung. Wie man am Beispiel dieser zwei Gauss-Verteilungen sehen kann, ist sowohl der geschätzte Zustand als auch der gemessene Zustand verteilt und haben unterschiedliche Standardabweichungen $\sigma$ und Erwartungswerte $\mu$. @@ -43,10 +46,10 @@ und für die Messung: Diesen werden nun Multipliziert und durch deren Fläche geteilt um sie wieder zu Normieren: \begin{equation} -{y_f}(x;{\mu_f},{\sigma_f})=\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}\,} +{y_f}(x;{\mu_f},{\sigma_f})=\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{equation} -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. Sie ist also Gewichtet und die best mögliche Schätzung. +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. Sie ist also gewichtet und die best mögliche Schätzung. \begin{figure} @@ -59,14 +62,13 @@ Dadurch gleicht sich die neue Kurve den anderen an. Interessant daran ist, dass Was in 2 Dimensionen erklärt wurde, funktioniert auch in mehreren Dimensionen. Dieses Prinzip mach sich der Kalman Filter zu nutze, und wird von uns für die Erdbeben Berechnung genutzt. -\subsection{Anwendungsgrenzen} -Nicht lineare Systeme %Noch nicht Fertig + \section{Aufbau} -Um ein Erdbeben kenntlich zumachen 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, bleibt die gekoppelte Masse in der regel stehen und das Gehäuse schwingt mit.Relativbewegung des Bodens kann damit als Längenänderung 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, welcher rein mechanisch funktioniert. Zudem kann er nur in eine Dimension Messwerte aufnehmen. Würde das System ausgebaut werden, um alle Horizontalbewegungen aufzunehmen, würde der Verwendung des Kalman-Filters zu kompliziert werden. Für zwei Dimensionen (x,y) würde der Pythagoras für das System benötigt werden. 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. 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. +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, bleibt die gekoppelte Masse stehen und das Gehäuse schwingt mit.Relativbewegung des Bodens kann damit als Längenänderung 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, welcher rein mechanisch funktioniert. Zudem kann er nur in eine Dimension Messwerte aufnehmen. Würde das System ausgebaut werden, um alle Horizontalbewegungen aufzunehmen, würde der Verwendung des Kalman-Filters zu kompliziert werden. Für zwei Dimensionen (x,y) würde der Pythagoras für das System benötigt werden. 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. 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. Einfachheitshalber beschränken wir uns aber auf den linearen Fall, da dadurch die wesentlichen punkte bereits aufgezeigt werden. \begin{figure} \begin{center} @@ -76,36 +78,20 @@ Wir konstruieren uns eine einfachere Version eines Seismographen, welcher rein m \end{figure} -\subsection{Optionen} -Wollte man einen 2D Seismographen aufbauen, ohne den Pythagroas zu verwenden, kann dies mit der Annahme, das die Feder sehr lang sind erfolgen. Da sich bei langen Federn die Auslenkungen verkleiner...!!Noch nicht fertig! \section{Systemgleichung} -Da das Kalman-Filter zum Schätzen des nächsten Zustand verwendet wird, wird eine Gleichung, welche das System beschreibt. Das Kalman-Filter benötigt eine Beschreibung der Systemdynamik. Im Fall unseres Seismographen, kann die Differentialgleichung zweiter Ordnung einer gedämpften Schwingung am harmonischen Oszillator verwendet werden. Diese lautet: +Da das Kalman-Filter zum Schätzen des nächsten Zustand verwendet wird, benötigt das Kalman-Filter eine Beschreibung der Systemdynamik. Im Fall unseres Seismographen, kann die Differentialgleichung zweiter Ordnung einer gedämpften Schwingung am harmonischen Oszillator verwendet werden. Diese lautet: \begin{equation} m\ddot x + 2k \dot x + Dx = f \end{equation} mit den Konstanten $m$ = Masse, $k$ = Dämpfungskonstante und $D$ = Federkonstante. -Um diese nun in die Systemmatrix umzuwandeln, wird aus der Differentialgleichung zweiter Ordnung durch eine Substitution eine DGL erster Ordnung: - -\begin{equation} -{x_1}=x, \qquad +Um diese nun in die Systemmatrix umzuwandeln, wird aus der Differentialgleichung zweiter Ordnung durch die Substitution \[ {x_1}=x, \qquad {x_2}=\dot x, \qquad -{x_3}=\ddot x\qquad \mid \quad \text {Substitution} -\end{equation} +{x_3}=\ddot x\qquad\] erhalten wir die Differentialgleichung \[ m{x_3}+ 2k{x_2} + D{x_1} = f.\] Diese können wir nun in der Form \[ {x_3}=-\frac{D}{m} {x_1} -\frac{2k}{m} {x_2} + \frac{f} {m} \] auch als Matrix-Vektor-Gleichung darstellen. -\begin{equation} -m{x_3}+ 2k{x_2} + D{x_1} = f\qquad \mid \quad \text {DGL 1. Ordnung} -\end{equation} - -\begin{equation} -{x_3}=-\frac{D}{m} {x_1} -\frac{2k}{m} {x_2} + \frac{f} {m} \qquad \mid \quad \text {nach} \quad{x_3} -\end{equation} -auch als Matrix-Vektor-Gleichung schreiben. -Hierbei beschreibt die Matrix $A$ die gesamte Systemdynamik in der Form, wie sie ein Kalman-Filter benötigt. - -Um die lineare Differentialgleichung in das Kalman-Filter zu implementieren, muss dieses als Vektor-Gleichung umgewandelt werden. Dafür wird die Gleichung in die Zustände aufgeteilt. Die für uns relevanten Zustände sind die Position der Masse, die Geschwindigkeit der Masse und äussere Beschleunigung des ganzen System. Dabei muss unterschieden werden. um welche Beschleunigung es sich handelt. Das System beinhaltet sowohl eine Beschleunigung der Masse bzw. Feder (innere Beschleunigung), als auch eine Beschleunigung der ganzen Apparatur (äussere Beschleunigung). In unserem Fall wird die äusseren Beschleunigung gesucht, da diese der Erdbeben Anregung gleich kommt. +Dafür wird die Gleichung in die Zustände aufgeteilt. Die für uns relevanten Zustände sind die Position der Masse, die Geschwindigkeit der Masse und die äussere Beschleunigung des ganzen System. Dabei muss unterschieden werden, um welche Beschleunigung es sich handelt. Das System beinhaltet sowohl eine Beschleunigung der Masse bzw. Feder (innere Beschleunigung), als auch eine Beschleunigung der ganzen Apparatur (äussere Beschleunigung). In unserem Fall wird die äusseren Beschleunigung gesucht, da diese der Erdbeben Anregung gleich kommt. \begin{equation} @@ -134,7 +120,7 @@ Um den Kalman Filter zu starten, müssen gewisse Bedingungen definiert werden. I \subsection{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. +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. \begin{equation} {x_0 }= \left( \begin{array}{c} 0\\ 0\\ 0\end{array}\right) @@ -142,6 +128,8 @@ Das Filter benötigt eine Anfangsbedingung. In unserem Fall ist es die Ruhelage, \subsubsection*{Anfangsfehler / Kovarianzmatrix $P$} Da auch der Anfangszustand fehlerhaft sein kann, wird für den Filter einen Anfangsfehler eingeführt. 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) undVarianz: 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. Somit ergibt sich für die Kovarianzmatrix \begin{equation} @@ -183,16 +171,16 @@ Q = \left( \end{array}\right) \end{equation} -Die Standabweichungen müssten Statistisch ermittelt werden, da der Fehler nicht vom Sensor kommt und somit nicht vom Hersteller gegeben ist. Das Bedeutet wiederum dass $Q$ die Unsicherheit des Prozesses beschreibt, und die Messung. +Die Standabweichungen müssten Statistisch ermittelt werden, da der Fehler nicht vom Sensor kommt und somit nicht vom Hersteller gegeben ist. Das Bedeutet wiederum dass $Q$ die Unsicherheit des Prozesses beschreibt, und die der Messung. \subsubsection*{Messmatrix $H$} -Die Messmatrix gibt an, welcher Parameter gemessen werden soll. In unserem Fall ist es nur die Position der Massen. +Die Messmatrix gibt an, welche Parameter gemessen werden soll. In unserem Falle ist es nur die Position der Massen. \[ H = (1, 0, 0) \] \subsubsection*{Messrauschkovarianz $R$} -Die Messrauschkovarianzmatrix beinhaltet, wie der Name es schon sagt, das Rauschen der Positionssensoren. In unserem Fall wird nur die Position der Masse gemessen. Da wir keine anderen Sensoren haben ist dies lediglich: +Die Messrauschkovarianzmatrix beinhaltet, wie der Name es schon sagt, das Rauschen der Positionssensoren. In unserem Fall wird nur die Position der Masse gemessen. Da wir keine anderen Sensoren haben ist $R$ lediglich: \begin{equation} R= ({\sigma_x}^2). \end{equation} @@ -203,17 +191,14 @@ Nachdem alle Parameter aufgestellt sind, wird der Filter initialisiert und wird \subsubsection*{Vorhersage} -Im Filterschritt Vorhersage wird der nächste Zustand anhand des Anfangszustand und der Systemmatrix berechnet. Dies funktioniert ganz Trivial mit dem Rechenschritt: +Im Filterschritt Vorhersage wird der nächste Zustand anhand des Anfangszustand und der Systemmatrix berechnet. Dies funktioniert mit dem Rechenschritt: \begin{equation} {x_{t+1}}=A\cdot{x_t}. \end{equation} -Die Kovarianz $P_{pred}$ wird ebenfalls neu berechnet, da die Unsicherheit im Vorhersage grösser wird als im Aktuellen. Da wir ein mehrdimensionales System haben, kommt noch die Messunsicherheit $Q$ dazu, so dass die Unsicherheit des Anfangsfehlers $P$ immer grösser wird. Dies funktioniert durch multiplizieren der Systemmatrix, deren Ableitung und mit dem aktualisierten Anfangsfehler. Dazu wird noch die Messunsicherheit addiert, somit entsteht die Gleichung - -\begin{equation} -{P_{pred}}=APA^T+Q. -\end{equation} +Die Kovarianz $P_{pred}$ wird ebenfalls neu berechnet. Da wir ein mehrdimensionales System haben, 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 +\[ P_\mathrm{pred} = A P A^T + Q . \] wird dieser Vorgang wiederholt, schaut der Filter 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. Das Filter passt sich selber an und korrigiert sich bei grosser Abweichung. @@ -223,41 +208,164 @@ Hier bei wird lediglich die Messung mit dem Fehler behaftet, und die Messmatrix \begin{equation} w=Z-(H\cdot x) \end{equation} -Die Innovation ist der Teil der Messung, die nicht durch die Systemdynamik erklärt werden kann. -Innovation = Messung - Vorhersage. Dies ist Intuitiv logisch, eine Innovation von 0 bedeutet, dass die Messung nichts Neues hervorbrachte. +Die Innovation ist der Teil der Messung, die nicht durch die Systemdynamik erklärt werden kann. Die Hilfsgröße Innovation beschreibt, wie genau der vorhergesagte Mittelwert den aktuellen Messwert mittels der Beobachtungsgleichung beschreiben kann. 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 ist intuitiv logisch, 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. +\begin{equation} +S=Z-(H\cdot P\cdot H`+R) +\end{equation} + + +\subsubsection*{Aktualisieren} +Im nächsten Schritt kommt nun die Wahrscheinlichkeit nach Gauss dazu. + +\begin{equation} +K= \frac{P \cdot H`}S +\end{equation} +Dieser Vorgang wird Kalman-Gain genannt. Er sagt aus, welcher Kurve mehr Vertraut werden soll, dem Messwert oder der Systemdynamik. +Das Kalman-Gain wird geringer wen der Messwert dem vorhergesagten Systemzustand entspricht. Sind die Messwerte komplett anders als die Vorhersage, wo werden die Elemente in der Matrix $K$ grösser. -\subsubsection*{Korrigieren} -Udpdate -\section{Anfügen der Schwingung} +Anhand der Informationen aus dem Kalman-Gain $K$ wird das System geupdated. + +\begin{equation} +x=x+(K \cdot w) +\end{equation} -Ein Erdbeben breitet sich im Boden wellenartig aus und bringt Objekte, wie zum Beispiel ein Gebäude, in Schwingung. -Diese Schwingungen pflanzen sich im Gebäude mit gleicher Amplitude, Geschwindigkeit und Beschleunigung in horizontaler und vertikaler Bewegung fort. -Wir möchten herauszufinden, wie gross die Massenbeschleunigung infolge eines Erdbeben ist. -Mit Hilfe von fiktiven Sensoren, die eine Ortsveränderung des Gebäude messen, können wir mit Anwendung von Matrizen und dem Kalman-Filter die Beschleunigung berechnen. +Dazu kommt eine neue Kovarianz für den nächste Vorhersageschritt: \begin{equation} -\int_a^b x^2\, dx -= -\left[ \frac13 x^3 \right]_a^b -= -\frac{b^3-a^3}3. -\label{erdbeben:equation1} +P=(I-(K \cdot H)) \cdot P \end{equation} -\section{Erreger-Schwingung} -Wir möchten mit einer gedämpften harmonischen Schwingung ein einfaches Erdbeben simulieren, die im Kalman Filter eingespeist wird. -Die Gleichung lautet +Der ganze Ablauf wird nun zum Algorithmus und beginnt wieder mit der Vorhersage + +\begin{equation} +{x_{t+1}}=e^{A\Delta t}{ x_t}. +\end{equation} + + +\subsection{Zusammenfassung des Filters} + + +1. Nächster Zustand vorhersagen +\begin{equation} +{x_{k|k-1}}={A_{k-1}}{x_{k-1}}+{B_{k-1}}{u_{k-1}} +\end{equation} +2. Nächste Fehlerkovarianz vorhersagen \begin{equation} -x(t)=Ae^{t/2}sin(t). +{P_{k|k-1}}={A_{k-1}}{P_{k-1}}{A_{k-1}^T}+{Q_{k-1}} +\end{equation} + + +3. Das Kalman Filter anwenden +\begin{equation} +{K_k}={P_{k-1}}{H_{k}^T({H_k}{P_{k|k-1}}{H_k}^T}+{R_k})^{-1} \end{equation} -Mit dieser Schwingung können wir ein einachsiger Seismograph simulieren, der eine Ortsverschiebung auf der x-Achse durchführt. -Die Dämpfung der Schwingung ist relevant, da das System beim Schwingungsvorgang durch die Federkonstante und der Reibung, Energie verliert. +4. Schätzung aktualisieren +\begin{equation} +{x_k}={x_{k|k-1}}+{K_k}({z_k}-{H_k}{x_{k|k-1}}) +\end{equation} + +5. Fehlerkovarianz aktualisieren +\begin{equation} +{P_k}=(I-{K_k}{H_k}){P_{k|k-1}} +\end{equation} + +6. Die Outputs von $k$ werden die Inputs für ${k-1}$ und werden wieder im Schritt 1 verwendet + + +\section{Matlab-Code} +Um das simulierte Erdbeben auf die theoretische Apparatur zu bringen und mit dem Kalman-Filter Resultate zu generieren, wurde in Matlab der Algorithmus programmiert. +\begin{lstlisting} +%% Initialisierte Werte +t0 = 0.00; % Anfangszeit +deltat = 0.01; % Zeitschritt +tend = 50.00; % Endzeit + +% Standard-Abweichungen Prozess +sigmax = 0.05e-3; % Position +sigmav = 0.01e-3; % Geschwindigkeit +sigmaf = 1; % (Äussere) Kraft + +% Standard-Abweichung Messung +sigmam = 0.01e-3; + +% Systemparameter +m = 1.00; % Masse +D = 0.30; % Federkonstante +k = 0.10; % Dämpfung + + +%% Kalmanfilter +% Initialisierung + + +% Anfangszustand (Position, Geschwindigkeit, Kraft) +x0 = [0; 0; 0]; + +% Unsicherheit des Anfangszustand +P0 = [0, 0, 0; ... + 0, 0, 0; ... + 0, 0, 0]; + +% Systemmatrizen +A = [0, 1, 0;... % Dynamikmatrix + -D/m, -2*k/m, 1;... + 0, 0, 0]; % Ableitungen von f(t) unbekant. Annahme: 0 +A = expm(A * deltat); + +Q = [sigmax^2, 0, 0;... + 0, sigmav^2, 0;... + 0, 0, sigmaf^2]; % Prozessrauschen (Covarianz) + +% Messprozess +H = [1, 0, 0]; % Messmatrix +R = sigmam^2; % Messrauschen (Könnte durch Versuche bestimmt werden) + +I = eye(3); % Identity matrix (Einheitsmatrix) + +% Filterprozess + +% Initialisieren der Variablen +N = length(t); % Anzahl Punkte im Einheitsvektor (= Anzahl Messwerte) +xhat = zeros(3, N); % Matrix mit geschätzten Zuständen + +% Index ':' bedeutet: 'alles' +% Index '(1, :)' bedeutet: 'alles aus der 1. Zeile' + +% Anfangszustand setzen +xhat(:, 1) = x0; +P = P0; + +% Kalman-Matrizen konvergiert. Vorab-Berechnung in 'genügenden' Iterationen +for idx = 1:100 + Ppred = A * P * A' + Q; % Prädizieren der Kovarianz + S = (H * Ppred * H' + R); % Innovationskovarianz + K = Ppred * H' / S; % Filter-Matrix (Kalman-Gain) + P = (I - K * H) * Ppred; % Aktualisieren der Kovarianz +end + +% Anfangszustand gegeben +% Erster zu berechnender Wert ist der zweite +for idx = 2:N + % Vorhersage + xpred = A * xhat(:, idx-1); % Prädizierter Zustand aus Bisherigem und System + % Ppred = A * P * A' + Q; % Prädizieren der Kovarianz + + % Korrektur + y = xt(idx) - H * xpred; % Messungen/ Kraft aus System - Vohersage + % S = (H * Ppred * H' + R); % Innovationskovarianz + % K = Ppred * H' / S; + + xhat(:, idx) = xpred + K * y; % Aktualisieren des Systemzustands + % P = (I - K * H) * Ppred; % Aktualisieren der Kovarianz +end +\end{lstlisting} + -Die Ergebnisse dieser Schwingung setzen wir in die Messmatrix ein und können den Kalman-Filter starten. -- cgit v1.2.1 From 7ad1f062b8b34fb5ba4b99fc0a34aba0567090d6 Mon Sep 17 00:00:00 2001 From: Lukaszogg <82384106+Lukaszogg@users.noreply.github.com> Date: Thu, 8 Jul 2021 20:04:43 +0200 Subject: Update teil1.tex --- buch/papers/erdbeben/teil1.tex | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index d6f5638..98bbd9e 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -245,8 +245,8 @@ Der ganze Ablauf wird nun zum Algorithmus und beginnt wieder mit der Vorhersage \end{equation} -\subsection{Zusammenfassung des Filters} - +\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 \begin{equation} -- cgit v1.2.1 From 6fdb95a0d5ec6c8f8bd016e5343a397a40a18c9e Mon Sep 17 00:00:00 2001 From: Roy Seitz Date: Sun, 11 Jul 2021 12:27:32 +0200 Subject: Remove '\begin{document}' from chapter part file. --- buch/papers/erdbeben/teil1.tex | 4 ---- 1 file changed, 4 deletions(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index 98bbd9e..bb3bdd4 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -10,10 +10,6 @@ % - -\begin{document} - - \section{Kalman Filter} \subsection{Geschichte} Das Kalman Filter wurde 1960 von Rudolf Emil Kalman entdeckt und direkt von der NASA für die Appollo Mission benutzt. Der Filter kommt mit wenig Rechenleistung aus und war somit dafür geeignet die Rakete bei der Navigation zu unterstützen. Das Filter schätzt den Zustand eines Systems anhand von Messungen und kann den nächsten Zustand errechnen. Typische Anwendungen des Kalman-Filters sind die Glättung von verrauschten Daten und die Schätzung von Parametern und kommt heutzutage in jedem Satellit, Navigationssystem, Smartphones und Videospielen vor. -- cgit v1.2.1 From 84b45ebb3d4788db314796095ddcca2c901d406d Mon Sep 17 00:00:00 2001 From: fabioviecelli <80270098+fabioviecelli@users.noreply.github.com> Date: Fri, 16 Jul 2021 16:13:27 +0200 Subject: Einleitung und Schwingungsgleichung --- buch/papers/erdbeben/Teil_Fabio.tex | 110 ++++++++++++++++++++++++++++++++++++ buch/papers/erdbeben/main.tex | 5 +- 2 files changed, 113 insertions(+), 2 deletions(-) create mode 100644 buch/papers/erdbeben/Teil_Fabio.tex (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/Teil_Fabio.tex b/buch/papers/erdbeben/Teil_Fabio.tex new file mode 100644 index 0000000..63b9648 --- /dev/null +++ b/buch/papers/erdbeben/Teil_Fabio.tex @@ -0,0 +1,110 @@ +\section{Kalman Filter} +\subsection{Was ist ein Erdbeben?} +Für das Verständnis möchten wir zuerst klären, was ein Erdbeben genau ist. +Das soll uns helfen, eine Verknüpfung zwischen dem Naturphänomen und der mathematischen Lösungsfindung herzustellen. + +Unter einem Erdbeben verstehen wir eine Erschütterung des Erdkörpers. +Dabei reiben zwei tektonische Platten aneinander, welche aber sich durch die Gesteinsverzahnung gegenseitig blockieren. +Aufgrund dieser Haftreibung entstehen Spannungen, die sich immer mehr bis zum Tipping Point aufbauen. +Irgendwann ist der Punkt erreicht, in dem die Scherfestigkeit der Gesteine überwunden wird. +Wenn dies passiert, entladet sich die aufgebaute Spannung und setzt enorme Energien frei, die wir als Erdbeben wahrnehmen. + +Ein Erdbeben pflanzt 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. + +Wir möchten nun mittels Kalman-Filter die Erdbebenbeschleunigung herausfinden. +Die Erdbebenbeschleunigung ist in der Praxis zur Entwicklung von Erdbebengefährdungskarten, sowie der Ausarbeitung von Baunormen für erdbebengerechte Bauweise von Bedeutung. + + +\subsection{Künstliche Erdbebendaten} +Nun möchten wir anhand eines eigenen Beispiels das Kalman-Filter anwenden. +Wir müssen Erdbebendaten künstlich erzeugen, um sie in das Filter zu geben und somit den Prozess zu starten. +Dafür nehmen wir die Formel für harmonische gedämpfte Schwingungen, die + +\begin{equation} + y = A \sin(\omega t e^{-lambda t}) +\end{equation} + +lautet. + + + + +A ist die Amplitude der Schwingung und beschreibt die Heftigkeit eines Erdbebens, die Magnitude. +Omega repräsentiert die Erdbebenfrequenz, die in der Realität zwischen 1 Hz und 30 Hz betragen kann. +Wir wählen als Erwartungswert 15 Herz und für die Standardabweichung 1 Hz. +Lambda ist die Bodendämpfung, für die wir 0.2 wählen. +Wir haben diese Zahl aus der Literatur entnommen und ist für das Bauwesen bedeutend. +Je grösser Lambda gewählt wird, desto stärker wirkt die Dämpfung der Massenschwingung. +Die Funktion ist zeitabhängig und wir lassen pro Sekunde zehn Messwerte generieren. + +Die Frequenz soll im Matlab als Zufallszahl generiert werden. +Mit dem Golay-Filter glätten wir unsere Werte, um unser Output näher an die Realität zu bringen. +Zusätzlich werden Ausreisser nicht vernachlässigt und wirken geglättet in unsere Datenmenge. + +Grafik einfügen + +In der Grafik erkennen wir in den Sekunden 0 bis 10, dass die Sinuskurve gezackt ist. +Das deutet darauf hin, dass die Frequenz des Erdbebens einen hohen Einfluss auf die Masse des Seismographen hat. +Ab der 10. Sekunde bis zu tend, pendelt sich die Masse in ihre Eigenfrequenz ein und verhält sich unabhängiger vom Erdbeben. + +\subsection{Versuch} +Um den Kalman-Filter auszuprobieren, setzen wir nun Werte ein. +Für die Systemparameter wählen wir m=1.0, D = 0.3 und k = 0.1 und fügen es in die Differentialgleichung + +\begin{equation} + m\ddot x + 2k \dot x + Dx = f +\end{equation} + +ein und erhalten + +\begin{equation} + 1\ddot x + 0.1 \dot x + 0.3x = f +\end{equation} + + +\subsubsection*{Prozessrauschkovarianzmatrix $Q$} + + + + + +\begin{equation} + Q = \left( + \begin{array}{ccc} + (5 \cdot 10^{-5})^2 & 0 & 0 \\ + 0 & (1 \cdot 10^{-5})^2 & 0\\ + 0 & 0& ( 1 )^2\\ + \end{array}\right) +\end{equation} + + + + + +\subsection{Resultate} + +Vergleichen wir die künstlichen Messdaten mit der geschätzten Schwingung des Kalman-Filters, stellen wir fest, dass wir eine gute Methode gefunden haben, die Erdbebenbeschleunigung zu schätzen. +Obwohl die künstlichen Daten mit einer random-Funktion erzeugt werden, kann das Kalman-Filter präzise Vorhersagungen bilden. + +Für die Differentialgleichung zweiter Ordnung brauchen wir im Matlab die Funktion ode45. +Mit dieser Funktion können wir Differentialgleichungen auflösen. + + + + + + + + + + + + + +In Matlab fügen wir die Formel und unsere definierten Werte ein. +Die Frequenz generieren wir mit einem Zufallscode, +Mit einem Zufallscode und einen Zeitraum + +Matlabcode einfügen + diff --git a/buch/papers/erdbeben/main.tex b/buch/papers/erdbeben/main.tex index 83ef295..8f9c8d5 100644 --- a/buch/papers/erdbeben/main.tex +++ b/buch/papers/erdbeben/main.tex @@ -29,8 +29,9 @@ Bilden Sie auch für Formeln kurze Zeilen, einerseits der besseren \input{papers/erdbeben/teil0.tex} \input{papers/erdbeben/teil1.tex} -\input{papers/erdbeben/teil2.tex} -\input{papers/erdbeben/teil3.tex} +%\input{papers/erdbeben/teil2.tex} +%\input{papers/erdbeben/teil3.tex} +\input{papers/erdbeben/Teil_Fabio.tex} \printbibliography[heading=subbibliography] \end{refsection} -- cgit v1.2.1 From 04b1cb248d2e559814dd6551cb331d95b9df9fdf Mon Sep 17 00:00:00 2001 From: fabioviecelli <80270098+fabioviecelli@users.noreply.github.com> Date: Sat, 17 Jul 2021 15:19:58 +0200 Subject: Einleitung, Schwingungsgleichung, Matlab-code --- buch/papers/erdbeben/Teil_Fabio.tex | 216 +++++++++++++++++++++++++----------- 1 file changed, 154 insertions(+), 62 deletions(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/Teil_Fabio.tex b/buch/papers/erdbeben/Teil_Fabio.tex index 63b9648..9f5d092 100644 --- a/buch/papers/erdbeben/Teil_Fabio.tex +++ b/buch/papers/erdbeben/Teil_Fabio.tex @@ -1,15 +1,16 @@ -\section{Kalman Filter} +\section{Kalman-Filter} \subsection{Was ist ein Erdbeben?} -Für das Verständnis möchten wir zuerst klären, was ein Erdbeben genau ist. -Das soll uns helfen, eine Verknüpfung zwischen dem Naturphänomen und der mathematischen Lösungsfindung herzustellen. +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. + Unter einem Erdbeben verstehen wir eine Erschütterung des Erdkörpers. -Dabei reiben zwei tektonische Platten aneinander, welche aber sich durch die Gesteinsverzahnung gegenseitig blockieren. +Dabei reiben zwei tektonische Platten aneinander, welche sich durch die Gesteinsverzahnung gegenseitig blockieren. Aufgrund dieser Haftreibung entstehen Spannungen, die sich immer mehr bis zum Tipping Point aufbauen. Irgendwann ist der Punkt erreicht, in dem die Scherfestigkeit der Gesteine überwunden wird. -Wenn dies passiert, entladet sich die aufgebaute Spannung und setzt enorme Energien frei, die wir als Erdbeben wahrnehmen. +Wenn dies passiert, entlädt sich die aufgebaute Spannung und setzt enorme Energien frei, die wir als Erdbeben wahrnehmen. -Ein Erdbeben pflanzt sich vom Erdbebenherd in allen Richtungen gleich aus. +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. Wir möchten nun mittels Kalman-Filter die Erdbebenbeschleunigung herausfinden. @@ -18,29 +19,47 @@ Die Erdbebenbeschleunigung ist in der Praxis zur Entwicklung von Erdbebengefähr \subsection{Künstliche Erdbebendaten} Nun möchten wir anhand eines eigenen Beispiels das Kalman-Filter anwenden. -Wir müssen Erdbebendaten künstlich erzeugen, um sie in das Filter zu geben und somit den Prozess zu starten. -Dafür nehmen wir die Formel für harmonische gedämpfte Schwingungen, die +Da wir keine Rohdaten über vergangene Erdbeben zur Hand haben, müssen wir künstliche Daten erzeugen, um sie in das Filter einzugeben und somit den Prozess starten. +Dafür nehmen wir die Formel für harmonisch gedämpfte Schwingungen, die \begin{equation} y = A \sin(\omega t e^{-lambda t}) \end{equation} -lautet. - - - +lautet. A ist die Amplitude der Schwingung und beschreibt die Heftigkeit eines Erdbebens, die Magnitude. -Omega repräsentiert die Erdbebenfrequenz, die in der Realität zwischen 1 Hz und 30 Hz betragen kann. +Omega repräsentiert die Erdbebenfrequenz, die in der Realität zwischen 1 Hz und 30 Hz beträgt. Wir wählen als Erwartungswert 15 Herz und für die Standardabweichung 1 Hz. Lambda ist die Bodendämpfung, für die wir 0.2 wählen. -Wir haben diese Zahl aus der Literatur entnommen und ist für das Bauwesen bedeutend. +Wir haben diese Zahl aus der Literatur entnommen, denn sie ist für das Bauwesen bedeutend. +Lambda ist ein Materialparameter von Böden. + Je grösser Lambda gewählt wird, desto stärker wirkt die Dämpfung der Massenschwingung. Die Funktion ist zeitabhängig und wir lassen pro Sekunde zehn Messwerte generieren. -Die Frequenz soll im Matlab als Zufallszahl generiert werden. -Mit dem Golay-Filter glätten wir unsere Werte, um unser Output näher an die Realität zu bringen. -Zusätzlich werden Ausreisser nicht vernachlässigt und wirken geglättet in unsere Datenmenge. +Die Frequenz basiert auf einer random-Funktion, da wir das Erdbeben unberechenbar gestalten möchten. +Mit dem Golay-Filter können wir hohe Frequenz-Anteile in die Berechnung mit einfliessen lassen, anstatt sie abzuschneiden. +Die Bildung eines üblichen Mittelwerts wäre hier weniger geeignet. + +\begin{lstlisting} +freq = sgolayfilt(randn(size(Time)),0,11)*freqstd... ++freqmean; +\end{lstlisting} + +Mit der Frequenz erhalten wir die Winkelbeschleunigung und damit können wir die Amplitude berechnen. + + +\begin{lstlisting} +w = 2 * pi * freq; +a = Amplitude*sin(cumsum(w.*[0;diff(Time)])).*exp(-lambda*Time); +\end{lstlisting} + +Mit der Matlab-Funktion ode45 haben wir eine Funktion gefunden, um die Differentialgleichung aufzulösen. ode45 basiert auf dem Runge-Kutta-Verfahren, einem Einschrittverfahren, bei dem die Lösung ausgehend von einem gegebenen Anfangswert, in einer Näherung gesucht wird. + +\begin{lstlisting} +[T,Y] = ode45(@(t,x)ErzeugteSchwingung(t,x,m,k,d,a,Time),[0 tend], IC, SolverOptions); +\end{lstlisting} Grafik einfügen @@ -48,7 +67,8 @@ In der Grafik erkennen wir in den Sekunden 0 bis 10, dass die Sinuskurve gezackt Das deutet darauf hin, dass die Frequenz des Erdbebens einen hohen Einfluss auf die Masse des Seismographen hat. Ab der 10. Sekunde bis zu tend, pendelt sich die Masse in ihre Eigenfrequenz ein und verhält sich unabhängiger vom Erdbeben. -\subsection{Versuch} +\subsection{Versuch (bin noch dran)} + Um den Kalman-Filter auszuprobieren, setzen wir nun Werte ein. Für die Systemparameter wählen wir m=1.0, D = 0.3 und k = 0.1 und fügen es in die Differentialgleichung @@ -62,49 +82,121 @@ ein und erhalten 1\ddot x + 0.1 \dot x + 0.3x = f \end{equation} - -\subsubsection*{Prozessrauschkovarianzmatrix $Q$} - - - - - -\begin{equation} - Q = \left( - \begin{array}{ccc} - (5 \cdot 10^{-5})^2 & 0 & 0 \\ - 0 & (1 \cdot 10^{-5})^2 & 0\\ - 0 & 0& ( 1 )^2\\ - \end{array}\right) -\end{equation} - - - - +\subsection{Matlab Code} + + +\begin{lstlisting} + %% Initialisierte Werte + t0 = 0.00; % Anfangszeit + deltat = 0.01; % Zeitschritt + tend = 50.00; % Endzeit +\end{lstlisting} +Ein natürliches Erdbeben dauert zwischen wenigen Sekunden bis etwa eine Minute an. +50 Sekunden genügen für unsere Daten. +Pro Sekunde erhalten wir 100 Messpunkte, die für den Prozess des Filters eine präzise Anwendung ermöglichen. + +\begin{lstlisting} + % Standard-Abweichungen Prozess + sigmax = 0.05e-3; % Position + sigmav = 0.01e-3; % Geschwindigkeit + sigmaf = 1; % (Äussere) Kraft + + % Standard-Abweichung Messung + sigmam = 0.01e-3; +\end{lstlisting} + +Wir vertrauen dem System und geben kleine Standardabweichungen für die Position, Geschwindigkeit und Kraft ein. +Bei der Messung erwarten wir auch, dass die Sensoren genau funktionieren. +Jedoch hängt das vom Hersteller ab oder muss statistisch ermittelt werden. + + +\begin{lstlisting} + % Systemparameter +m = 1.00; % Masse +D = 0.30; % Federkonstante +k = 0.10; % Dämpfung +\end{lstlisting} +Hier werden die Spezifikationen des Seismographen definiert. + +\begin{lstlisting} +%% Kalmanfilter +% Initialisierung + +% Anfangszustand (Position, Geschwindigkeit, Kraft) +x0 = [0; 0; 0]; + +% Unsicherheit des Anfangszustand +P0 = [0, 0, 0; ... +0, 0, 0; ... +0, 0, 0]; + +% Systemmatrizen +A = [0, 1, 0;... % Dynamikmatrix +-D/m, -2*k/m, 1;... +0, 0, 0]; % Ableitungen von f(t) unbekant. Annahme: 0 +A = expm(A * deltat); + +Q = [sigmax^2, 0, 0;... +0, sigmav^2, 0;... +0, 0, sigmaf^2]; % Prozessrauschen (Covarianz) + + +\begin{lstlisting} +% Messprozess +H = [1, 0, 0]; % Messmatrix +R = sigmam^2; % Messrauschen (Könnte durch Versuche bestimmt werden) +\end{lstlisting} +Tritt ein Erdbeben ein, wird die Position der Masse in die Messmatrix eingetragen. + + +I = eye(3); % Identity matrix (Einheitsmatrix) + +\begin{lstlisting} +% Filterprozess + +% Initialisieren der Variablen +N = length(t); % Anzahl Punkte im Einheitsvektor (= Anzahl Messwerte) +xhat = zeros(3, N); % Matrix mit geschätzten Zuständen + +% Index ':' bedeutet: 'alles' +% Index '(1, :)' bedeutet: 'alles aus der 1. Zeile' + +% Anfangszustand setzen +xhat(:, 1) = x0; +P = P0; +\end{lstlisting} + +\begin{lstlisting} + +% Kalman-Matrizen konvergiert. Vorab-Berechnung in 'genügenden' Iterationen +for idx = 1:100 +Ppred = A * P * A' + Q; % Prädizieren der Kovarianz +S = (H * Ppred * H' + R); % Innovationskovarianz +K = Ppred * H' / S; % Filter-Matrix (Kalman-Gain) +P = (I - K * H) * Ppred; % Aktualisieren der Kovarianz +end +\end{lstlisting} + +In diesem Schritt wird die Kovarianz vorhergesagt, mit der Messung verglichen und nach jeder Berechnung aktualisiert. + +\begin{lstlisting} +% Anfangszustand gegeben +% Erster zu berechnender Wert ist der zweite +for idx = 2:N +% Vorhersage +xpred = A * xhat(:, idx-1); % Prädizierter Zustand aus Bisherigem und System +% Ppred = A * P * A' + Q; % Prädizieren der Kovarianz + +% Korrektur +y = xt(idx) - H * xpred; % Messungen/ Kraft aus System - Vohersage +% S = (H * Ppred * H' + R); % Innovationskovarianz +% K = Ppred * H' / S; + +xhat(:, idx) = xpred + K * y; % Aktualisieren des Systemzustands +% P = (I - K * H) * Ppred; % Aktualisieren der Kovarianz +end +\end{lstlisting} \subsection{Resultate} - -Vergleichen wir die künstlichen Messdaten mit der geschätzten Schwingung des Kalman-Filters, stellen wir fest, dass wir eine gute Methode gefunden haben, die Erdbebenbeschleunigung zu schätzen. -Obwohl die künstlichen Daten mit einer random-Funktion erzeugt werden, kann das Kalman-Filter präzise Vorhersagungen bilden. - -Für die Differentialgleichung zweiter Ordnung brauchen wir im Matlab die Funktion ode45. -Mit dieser Funktion können wir Differentialgleichungen auflösen. - - - - - - - - - - - - - -In Matlab fügen wir die Formel und unsere definierten Werte ein. -Die Frequenz generieren wir mit einem Zufallscode, -Mit einem Zufallscode und einen Zeitraum - -Matlabcode einfügen - +Grafik einfügen +Wir erkennen, dass wir mit dem Kalman-Filter eine gute Methode gefunden haben, die äussere Beschleunigung zu schätzen. Die Schätzung der nächsten Position der Federmasse liegt immer ziemlich nahe der tatsächlichen Messung. Man muss aber auch berücksichtigen, dass die Federschwingung ziemlich kontrolliert verläuft und das Kalman-Filter somit präzise Vorhersagen treffen kann. -- cgit v1.2.1 From 20dda3d0a116d5b4126f9b1210c3f7b2c1bab097 Mon Sep 17 00:00:00 2001 From: Lukaszogg <82384106+Lukaszogg@users.noreply.github.com> Date: Sat, 17 Jul 2021 16:02:04 +0200 Subject: Teil1+Teil0 --- buch/papers/erdbeben/teil0.tex | 97 +++++++-- buch/papers/erdbeben/teil1.tex | 439 ++++++++++++++++++----------------------- 2 files changed, 269 insertions(+), 267 deletions(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/teil0.tex b/buch/papers/erdbeben/teil0.tex index 6e89821..8ac5d6d 100644 --- a/buch/papers/erdbeben/teil0.tex +++ b/buch/papers/erdbeben/teil0.tex @@ -2,21 +2,88 @@ % einleitung.tex -- Beispiel-File für die Einleitung % % (c) 2020 Prof Dr Andreas Müller, Hochschule Rapperswil -% +%% \section{Teil 0\label{erdbeben:section:teil0}} -\rhead{Teil 0} -Lorem ipsum dolor sit amet, consetetur sadipscing elitr, sed diam -nonumy eirmod tempor invidunt ut labore et dolore magna aliquyam -erat, sed diam voluptua \cite{erdbeben:bibtex}. -At vero eos et accusam et justo duo dolores et ea rebum. -Stet clita kasd gubergren, no sea takimata sanctus est Lorem ipsum -dolor sit amet. - -Lorem ipsum dolor sit amet, consetetur sadipscing elitr, sed diam -nonumy eirmod tempor invidunt ut labore et dolore magna aliquyam -erat, sed diam voluptua. -At vero eos et accusam et justo duo dolores et ea rebum. Stet clita -kasd gubergren, no sea takimata sanctus est Lorem ipsum dolor sit -amet. +\rhead{Erdbeben} +\section{Erdbebenmessung} +\subsection{Was ist ein Erdbeben} +Fabio +\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, bleibt die gekoppelte Masse stehen aber das Gehäuse schwingt mit. +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 eine Gehäuse, an dem zwei Federn und eine Masse befestigt ist. +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. + +\begin{figure} + \begin{center} + \includegraphics[width=5cm]{papers/erdbeben/Apperatur} + \caption{Aufbau des Seismographen mit Gehäuse, Masse, Federn und Sensor} + \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 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. +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)$ + der Eigendynamik der Masse. +Um die Bewegung der Masse zu berechnen, müssen wir Gleichungen für unser System finden. + +\subsection{Systemgleichung} +Im Fall unseres Seismographen, kann die Differentialgleichung zweiter Ordnung einer gedämpften Schwingung am harmonischen Oszillator verwendet werden. +Diese lautet: +\begin{equation} +m\ddot s + 2k \dot s + Ds = f +\end{equation} +mit den Konstanten $m$ = Masse, $k$ = Dämpfungskonstante und $D$ = Federkonstante. +Um diese nun in die Systemmatrix umzuwandeln, wird die Differentialgleichung zweiter Ordnung substituiert: +\[ {x_1}=s \qquad +{x_2}=\dot s, \qquad\] +Somit entstehen die Gleichungenür die Position $s(t)$ der Masse : +\[ \dot {x_1} = {x_2}\] +und +\[ \dot x_2 = -\frac{D}{m} {x_1} -\frac{2k}{m} {x_2} + \frac{f} {m} \] für die Geschwindigkeit $v(t)$ der Masse. + +Diese können wir nun in der Form +\[ {x_3}=-\frac{D}{m} {s_1} -\frac{2k}{m} {s_2} + \frac{f} {m} \] +auch als Matrix-Vektor-Gleichung darstellen. +Dafür wird die Gleichung in die Zustände aufgeteilt. +Die für uns relevanten Zustände sind die Position der Masse, die Geschwindigkeit der Masse und die äussere Beschleunigung des ganzen System. +Dabei muss unterschieden werden, um welche Beschleunigung es sich handelt. +Das System beinhaltet sowohl eine Beschleunigung der Masse (innere Beschleunigung), als auch eine Beschleunigung der ganzen Apparatur (äussere Beschleunigung). +In unserem Fall wird die äusseren Beschleunigung gesucht, da diese der Erdbebenanregung gleich kommt. +\begin{equation} +\frac{d}{dt} \left(\begin{array}{c} {s_1} \\ {s_2} \end{array}\right) = \left( + \begin{array}{ccc} +0 & 1& 0 \\ +- \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m}\\ +\end{array}\right) \left(\begin{array}{c} {s_1} \\ {s_2} \\ {s_3} \end{array}\right). +\end{equation} + +Durch Rücksubstituion ergibt sich: +\begin{equation} +\frac{d}{dt} \left(\begin{array}{c} s(t) \\ v(t) \end{array}\right) = \left( + \begin{array}{ccc} +0 & 1& 0 \\ +- \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m}\\ +\end{array}\right) \left(\begin{array}{c} s(t)\\ v(t)\\ f(t) \end{array}\right). +\end{equation} +Wir wissen nicht wie sich die Kraft verhält. +Deshalb treffen wir die Annahme, das sich die Kraft über die Beobachtungszeit nicht verändert. +Diese unzutreffende Annahme wird später durch einen grossen Systemfehler kompensiert. +Da die Kraft unbekannt ist, wird die letzte Zeile mit Nullen gefüllt, denn genau diese Werte wollen wir. + + + + + + + + + diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index 98bbd9e..a4e2220 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -7,132 +7,107 @@ % teil2.tex -- Beispiel-File für teil2 % % (c) 2020 Prof Dr Andreas Müller, Hochschule Rapperswil -% +%% + +\rhead{Kalman-Filter} -\begin{document} +\section{Kalman-Filter} +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 den linearen Kalman-Filter. +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 \cdot Eigendynamik des Systems gerechnet. +Die Idee dahinter ist, dass das Kalman-Filter die nicht-deterministische Grösse $f$ anhand der Messung und der Vorhersage zu bestimmen. +Für mehrere Dimensionen (x,y,z) würde der Pythagoras für das System benötigt werden. +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. +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. +Einfachheitshalber beschränken wir uns auf den linearen Fall, da dadurch die wesentlichen Punkte bereits aufgezeigt werden. -\section{Kalman Filter} \subsection{Geschichte} -Das Kalman Filter wurde 1960 von Rudolf Emil Kalman entdeckt und direkt von der NASA für die Appollo Mission benutzt. Der Filter kommt mit wenig Rechenleistung aus und war somit dafür geeignet die Rakete bei der Navigation zu unterstützen. Das Filter schätzt den Zustand eines Systems anhand von Messungen und kann den nächsten Zustand errechnen. Typische Anwendungen des Kalman-Filters sind die Glättung von verrauschten Daten und die Schätzung von Parametern und kommt heutzutage in jedem Satellit, Navigationssystem, Smartphones und Videospielen vor. +Das Kalman-Filter wurde 1960 von Rudolf Emil Kalman entdeckt und direkt von der NASA für die Appollo Mission benutzt. Der Filter kommt mit wenig Rechenleistung aus und war somit dafür geeignet die Rakete bei der Navigation zu unterstützen. Das Filter schätzt den Zustand eines Systems anhand von Messungen und kann den nächsten Zustand errechnen. 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 zwei Normalverteilungen oder auch Gauss-Verteilung. Die eine Kurve zeigt die errechnete Vorhersage des Zustands, bzw. deren Normalverteilung. Die andere Kurve zeigt die verrauschte Messung des nächsten Zustand, bzw. deren Normalverteilung. Wie man am Beispiel dieser zwei Gauss-Verteilungen sehen kann, ist sowohl der geschätzte Zustand als auch der gemessene Zustand verteilt und haben unterschiedliche Standardabweichungen $\sigma$ und Erwartungswerte $\mu$. - - +Das Kalman-Filter schätzt den wahrscheinlichsten Wert zwischen Normalverteilungen. +Dies bedeutet, das Filter schätzt nicht nur den Mittelwert, sondern auch die Standartabweichung. +Da Normalverteilungen dadurch vollständig definiert sind, schätzt ein Kalman-Filter die gesamte Verteilungsfunktion des Zustandes. +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 unten sehen kann, ist sowohl der geschätzte Zustand als auch der gemessene Zustand normalverteilt und haben dementsprechend unterschiedliche Standardabweichungen $\sigma$ und Erwartungswerte $\mu$. \begin{figure} \begin{center} \includegraphics[width=5cm]{papers/erdbeben/Gausskurve2.pdf} - \caption{System} + \caption{Zwei Normalerteilungen; Die eine Funktion zeigt die Vorhersage, die andere die Messung} \end{center} \end{figure} - -Um eine genauere Schätzung des Zustandes zu machen, wird nun ein Wert zwischen den beiden Verteilungen gesucht. An diesem Punkt wird nun eine Eigenschaft ausgenutzt. Durch das Multiplizieren zweier Normalverteilungen entsteht eine neue Normalverteilung. - +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: -\begin{equation} -{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}} -\end{equation} -und für die Messung: - -\begin{equation} -{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}}. -\end{equation} - -Diesen werden nun Multipliziert und durch deren Fläche geteilt um sie wieder zu Normieren: -\begin{equation} -{y_f}(x;{\mu_f},{\sigma_f})=\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{equation} -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. Sie ist also gewichtet und die best mögliche Schätzung. +\[ {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: +\[ {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}}. \] -\begin{figure} - \begin{center} - \includegraphics[width=5cm]{papers/erdbeben/Gausskurve3.pdf} - \caption{System} - \end{center} -\end{figure} - - -Was in 2 Dimensionen erklärt wurde, funktioniert auch in mehreren Dimensionen. Dieses Prinzip mach sich der Kalman Filter zu nutze, und wird von uns für die Erdbeben Berechnung genutzt. +Diesen werden nun Multipliziert und durch deren Fläche geteilt um sie wieder zu Normieren: +\[ +{y_f}(x;{\mu_f},{\sigma_f})=\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}\cdot{y_2} dx\,} + \] +Diese Kombination der beiden Verteilungen resultiert wiederum in einer Normalverteilung +\[ {y_f}(x; {\mu_f}, {\sigma_f}) = {y_1}(x;{ \mu_1},{ \sigma_1}) {\cdot y_2}(x; {\mu_2}, {\sigma_2}), \] +mit Erwartungswert +\[ \mu_f = \frac{\mu_1\sigma_2^2 + \mu_2 \sigma_1^2}{\sigma_1^2 + \sigma_2^2} \] +und Varianz +\[ \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. +Sie ist also gewichtet und die best mögliche Schätzung. -\section{Aufbau} -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, bleibt die gekoppelte Masse stehen und das Gehäuse schwingt mit.Relativbewegung des Bodens kann damit als Längenänderung 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, welcher rein mechanisch funktioniert. Zudem kann er nur in eine Dimension Messwerte aufnehmen. Würde das System ausgebaut werden, um alle Horizontalbewegungen aufzunehmen, würde der Verwendung des Kalman-Filters zu kompliziert werden. Für zwei Dimensionen (x,y) würde der Pythagoras für das System benötigt werden. 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. 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. Einfachheitshalber beschränken wir uns aber auf den linearen Fall, da dadurch die wesentlichen punkte bereits aufgezeigt werden. \begin{figure} \begin{center} - \includegraphics[width=5cm]{papers/erdbeben/Apperatur} - \caption{System} + \includegraphics[width=5cm]{papers/erdbeben/Gausskurve3.pdf} + \caption{Durch das Multiplizieren der blauen und der orangen Verteilung entsteht die die rote, optimale Funktion} \end{center} \end{figure} +Was in 2 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{Systemgleichung} -Da das Kalman-Filter zum Schätzen des nächsten Zustand verwendet wird, benötigt das Kalman-Filter eine Beschreibung der Systemdynamik. Im Fall unseres Seismographen, kann die Differentialgleichung zweiter Ordnung einer gedämpften Schwingung am harmonischen Oszillator verwendet werden. Diese lautet: -\begin{equation} -m\ddot x + 2k \dot x + Dx = f -\end{equation} -mit den Konstanten $m$ = Masse, $k$ = Dämpfungskonstante und $D$ = Federkonstante. - -Um diese nun in die Systemmatrix umzuwandeln, wird aus der Differentialgleichung zweiter Ordnung durch die Substitution \[ {x_1}=x, \qquad -{x_2}=\dot x, \qquad -{x_3}=\ddot x\qquad\] erhalten wir die Differentialgleichung \[ m{x_3}+ 2k{x_2} + D{x_1} = f.\] Diese können wir nun in der Form \[ {x_3}=-\frac{D}{m} {x_1} -\frac{2k}{m} {x_2} + \frac{f} {m} \] auch als Matrix-Vektor-Gleichung darstellen. - - -Dafür wird die Gleichung in die Zustände aufgeteilt. Die für uns relevanten Zustände sind die Position der Masse, die Geschwindigkeit der Masse und die äussere Beschleunigung des ganzen System. Dabei muss unterschieden werden, um welche Beschleunigung es sich handelt. Das System beinhaltet sowohl eine Beschleunigung der Masse bzw. Feder (innere Beschleunigung), als auch eine Beschleunigung der ganzen Apparatur (äussere Beschleunigung). In unserem Fall wird die äusseren Beschleunigung gesucht, da diese der Erdbeben Anregung gleich kommt. - - -\begin{equation} -\frac{d}{dt} \left(\begin{array}{c} {x_1} \\ {x_2} \end{array}\right) = \left( - \begin{array}{ccc} -0 & 1& 0 \\ -- \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m}\\ -\end{array}\right) \left(\begin{array}{c} {x_1} \\ {x_2} \\ {x_3} \end{array}\right). -\end{equation} - -Durch die Rücksubstituion ergibt sich: -\begin{equation} -\frac{d}{dt} \left(\begin{array}{c} x(t) \\ v(t) \end{array}\right) = \left( - \begin{array}{ccc} -0 & 1& 0 \\ -- \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m}\\ -\end{array}\right) \left(\begin{array}{c} x(t)\\ v(t)\\ f(t) \end{array}\right). -\end{equation} - - -Da die Kraft unbekannt ist, wird die letzte Zeile später mit Nullen bestückt, denn genau diese Werte wollen wir. - -\section{Kalman Filter} -Um den Kalman Filter zu starten, müssen gewisse Bedingungen definiert werden. In diesem Abschnitt werden die einzelnen Parameter/Matrizen erläutert und Erklärt, wofür sie nützlich sind. - +\section{Filter-Matrizen} +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{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. +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. -\begin{equation} -{x_0 }= \left( \begin{array}{c} 0\\ 0\\ 0\end{array}\right) -\end{equation} + +\[ {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*{Anfangsfehler / Kovarianzmatrix $P$} -Da auch der Anfangszustand fehlerhaft sein kann, wird für den Filter einen Anfangsfehler eingeführt. 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) undVarianz: Var(x) = Cov(x, x) +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. Somit ergibt sich für die Kovarianzmatrix +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 für die Kovarianzmatrix ergibt sich -\begin{equation} +\[ {P_0 }= \left( \begin{array}{ccc} @@ -141,229 +116,189 @@ In unserem Fall ist der Anfangszustand gut bekannt. Wir gehen davon aus, dass da 0 & 0 &0 \\ \end{array} \right). -\end{equation} -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 (1 Million) verwendet werden. Grosse Werte ermöglichen dem Filter sich schnell einzupendeln. - + \] +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$} -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. +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 Vektor-Gleichung lautet daher: - - -\begin{equation} +\[ A = \left( \begin{array}{ccc} 0 & 1& 0 \\ - \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m}\\ 0 & 0& 0\\ \end{array}\right) -\end{equation} + \] +Dabei soll der Kalman-Filter in diskreten Zeitschritten $\Delta t$ arbeiten. +Die Übergangs-Matrix erhalten wir aus der Systemdynamikmatrix mittels Exponentialfunktion: +\[\Phi = \exp(A\Delta t). \] \subsubsection*{Prozessrauschkovarianzmatrix $Q$} -Die Prozessrauschmatrix teilt dem Filter mit, wie sich der Systemzustand verändert. Kalman-Filter berücksichtigen Unsicherheiten wie Messfehler und -rauschen. Bei unserem Modell könnte das beispielsweise ein Windstoss an die Masse sein. Für uns wäre dies: -\begin{equation} +Die Prozessrauschmatrix teilt dem Filter mit, wie sich der Prozess verändert. +Kalman-Filter berücksichtigen Unsicherheiten wie Messfehler und -rauschen. +Bei unserem Modell könnte das beispielsweise ein Windstoss an die Masse sein. +Für uns wäre dies: +\[ Q = \left( \begin{array}{ccc} -{\sigma_x }^2& 0& 0 \\ +{\sigma_s }^2& 0& 0 \\ 0 & {\sigma_v }^2& 0\\ 0 & 0& {\sigma_f }^2\\ \end{array}\right) -\end{equation} + \] -Die Standabweichungen müssten Statistisch ermittelt werden, da der Fehler nicht vom Sensor kommt und somit nicht vom Hersteller gegeben ist. Das Bedeutet wiederum dass $Q$ die Unsicherheit des Prozesses beschreibt, und die der Messung. +Die Standabweichungen müssten statistisch ermittelt werden, da der Fehler nicht vom Sensor kommt und somit nicht vom Hersteller gegeben ist. +Das Bedeutet wiederum dass $Q$ die Unsicherheit des Prozesses beschreibt und nicht die der Messung. \subsubsection*{Messmatrix $H$} -Die Messmatrix gibt an, welche Parameter gemessen werden soll. In unserem Falle ist es nur die Position der Massen. +Die Messmatrix gibt an, welche Parameter gemessen werden +In unserem Falle ist es die Position der Massen. \[ H = (1, 0, 0) \] - \subsubsection*{Messrauschkovarianz $R$} -Die Messrauschkovarianzmatrix beinhaltet, wie der Name es schon sagt, das Rauschen der Positionssensoren. In unserem Fall wird nur die Position der Masse gemessen. Da wir keine anderen Sensoren haben ist $R$ lediglich: -\begin{equation} -R= ({\sigma_x}^2). -\end{equation} -Diese Messrauchen wird meistens vom Sensorhersteller angegeben. Für unsere Theoretische Apparatur wird hier ein kleiner Fehler eingesetzt. - -\subsection{Fiter Algorithmus} -Nachdem alle Parameter aufgestellt sind, wird der Filter initialisiert und wird den Zustand der Feder vorherzusagen, die Messung zu präzisieren und laufend zu aktualisieren. 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. - +Die Messrauschkovarianzmatrix beinhaltet, wie der Name es schon sagt, das Rauschen der Positionsmessung. +In unserem Fall wird nur die Position der Masse gemessen. Da wir keine anderen Sensoren haben ist $R$ lediglich: +\[ R= ({\sigma_{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. + +\subsection{Fiter-Agorithmus} +Nachdem alle Parameter aufgestellt sind, wird das Filter initialisiert. +Zuerst wird der nächste Zustand der Feder vorhergesagt, danach wird die Messung präzisiert und laufend zu aktualisieren. +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. \subsubsection*{Vorhersage} -Im Filterschritt Vorhersage wird der nächste Zustand anhand des Anfangszustand und der Systemmatrix berechnet. Dies funktioniert mit dem Rechenschritt: -\begin{equation} -{x_{t+1}}=A\cdot{x_t}. -\end{equation} - +Im Filterschritt Vorhersage wird der nächste Zustand anhand des Anfangszustand und der Systemmatrix berechnet. +Dies funktioniert mit dem Rechenschritt: +\[ +{x_{k|k-1}}=\Phi \cdot {x_{k-1|k-1}}= \exp(A\Delta t)\cdot{x_{k|k-1}}. + \] + +Die Kovarianz $P_{pred}$ wird ebenfalls neu berechnet. Da wir ein mehrdimensionales System haben, 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 +\[ {P_{k|k-1}} = {\Phi_k} {P_{k-1|k-1}} {\Phi_k} ^T + {Q_{k-1}} .\] +Es vergeht genau $dt$ Zeit, und dieser Vorgang wird wiederholt. +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. +Das Filter passt sich selber an und korrigiert sich bei grosser Abweichung. -Die Kovarianz $P_{pred}$ wird ebenfalls neu berechnet. Da wir ein mehrdimensionales System haben, 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 -\[ P_\mathrm{pred} = A P A^T + Q . \] +\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 -wird dieser Vorgang wiederholt, schaut der Filter 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. Das Filter passt sich selber an und korrigiert sich bei grosser Abweichung. +\[{w_{k}}={z_{k}}-{H_{k}}\cdot{x_{k|k-1}}.\] -\subsubsection*{Messen} -Der Sensor wurde noch nicht benutz, doch genau der liefert Werte für den 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$ -\begin{equation} -w=Z-(H\cdot x) -\end{equation} -Die Innovation ist der Teil der Messung, die nicht durch die Systemdynamik erklärt werden kann. Die Hilfsgröße Innovation beschreibt, wie genau der vorhergesagte Mittelwert den aktuellen Messwert mittels der Beobachtungsgleichung beschreiben kann. 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 ist intuitiv logisch, eine Innovation von 0 bedeutet, dass die Messung nichts Neues hervorbrachte. +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. +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 ist intuitiv logisch, 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. -\begin{equation} -S=Z-(H\cdot P\cdot H`+R) -\end{equation} - +\[ +{S_{k}}={H_{k}}{P_{k|k-1}}{H_{k}}^T+{R_{k}} + \] \subsubsection*{Aktualisieren} Im nächsten Schritt kommt nun die Wahrscheinlichkeit nach Gauss dazu. - -\begin{equation} -K= \frac{P \cdot H`}S -\end{equation} -Dieser Vorgang wird Kalman-Gain genannt. Er sagt aus, welcher Kurve mehr Vertraut werden soll, dem Messwert oder der Systemdynamik. -Das Kalman-Gain wird geringer wen der Messwert dem vorhergesagten Systemzustand entspricht. Sind die Messwerte komplett anders als die Vorhersage, wo werden die Elemente in der Matrix $K$ grösser. - +\[ +{K_{k}}= {{P_{k|k-1}} \cdot {H_{k}^T}}\cdot {S_{k}}^{-1} + \] +Dieser Vorgang wird Kalman-Gain genannt. +Er sagt aus, welcher Kurve mehr Vertraut werden soll, dem Messwert oder der Systemdynamik. +Das Kalman-Gain wird geringer wen der Messwert dem vorhergesagten Systemzustand entspricht. +Sind die Messwerte komplett anders als die Vorhersage, wo werden die Elemente in der Matrix $K$ grösser. Anhand der Informationen aus dem Kalman-Gain $K$ wird das System geupdated. -\begin{equation} -x=x+(K \cdot w) -\end{equation} +\[ +{x_{k|k}}={x_{k|k-1}}+({K_{k}}\cdot {w_{k}}) + \] Dazu kommt eine neue Kovarianz für den nächste Vorhersageschritt: -\begin{equation} -P=(I-(K \cdot H)) \cdot P -\end{equation} +\[ +{P_{k|k}}=(I-({K_{k}} \cdot {H_{k}})) \cdot {P_{k|k-1}} + \] Der ganze Ablauf wird nun zum Algorithmus und beginnt wieder mit der Vorhersage -\begin{equation} -{x_{t+1}}=e^{A\Delta t}{ x_t}. -\end{equation} +\[ +{x_{k|k-1}}=\Phi \cdot {x_{k-1|k-1}}= \exp(A\Delta t)\cdot{x_{k|k-1}}. + \] \subsection{Zusammenfassung } -Zusammenfassend kann das Kalman-Filter in offizieller Typus dargestellt werden. Dabei beginnt das Filter mit dem Anfangszustand für $k=0$ +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 -\begin{equation} -{x_{k|k-1}}={A_{k-1}}{x_{k-1}}+{B_{k-1}}{u_{k-1}} -\end{equation} +\[{x_{k|k-1}}=\Phi \cdot {x_{k-1|k-1}}= \exp(A\Delta t)\cdot{x_{k|k-1}}.\] 2. Nächste Fehlerkovarianz vorhersagen -\begin{equation} -{P_{k|k-1}}={A_{k-1}}{P_{k-1}}{A_{k-1}^T}+{Q_{k-1}} -\end{equation} - +\[{P_{k|k-1}}={\Phi _{k}} {P_{k-1|k-1}} {\Phi _{k}}^T + {Q_{k-1}}.\] 3. Das Kalman Filter anwenden -\begin{equation} -{K_k}={P_{k-1}}{H_{k}^T({H_k}{P_{k|k-1}}{H_k}^T}+{R_k})^{-1} -\end{equation} +\[{K_{k}}= {P_{k|k-1}} \cdot {H_{k}^T}\cdot {S_{k}^{-1}}\] 4. Schätzung aktualisieren -\begin{equation} -{x_k}={x_{k|k-1}}+{K_k}({z_k}-{H_k}{x_{k|k-1}}) -\end{equation} +\[{x_{k|k}}={x_{k|k-1}}+({K_{k}}\cdot {w_{k}}) \] 5. Fehlerkovarianz aktualisieren -\begin{equation} -{P_k}=(I-{K_k}{H_k}){P_{k|k-1}} -\end{equation} - -6. Die Outputs von $k$ werden die Inputs für ${k-1}$ und werden wieder im Schritt 1 verwendet - - -\section{Matlab-Code} -Um das simulierte Erdbeben auf die theoretische Apparatur zu bringen und mit dem Kalman-Filter Resultate zu generieren, wurde in Matlab der Algorithmus programmiert. -\begin{lstlisting} -%% Initialisierte Werte -t0 = 0.00; % Anfangszeit -deltat = 0.01; % Zeitschritt -tend = 50.00; % Endzeit - -% Standard-Abweichungen Prozess -sigmax = 0.05e-3; % Position -sigmav = 0.01e-3; % Geschwindigkeit -sigmaf = 1; % (Äussere) Kraft - -% Standard-Abweichung Messung -sigmam = 0.01e-3; +\[{P_{k|k}}=(I-({K_{k}}\cdot {H_{k}})) \cdot {P_{k|k-1}} \] -% Systemparameter -m = 1.00; % Masse -D = 0.30; % Federkonstante -k = 0.10; % Dämpfung - -%% Kalmanfilter -% Initialisierung - - -% Anfangszustand (Position, Geschwindigkeit, Kraft) -x0 = [0; 0; 0]; - -% Unsicherheit des Anfangszustand -P0 = [0, 0, 0; ... - 0, 0, 0; ... - 0, 0, 0]; - -% Systemmatrizen -A = [0, 1, 0;... % Dynamikmatrix - -D/m, -2*k/m, 1;... - 0, 0, 0]; % Ableitungen von f(t) unbekant. Annahme: 0 -A = expm(A * deltat); - -Q = [sigmax^2, 0, 0;... - 0, sigmav^2, 0;... - 0, 0, sigmaf^2]; % Prozessrauschen (Covarianz) - -% Messprozess -H = [1, 0, 0]; % Messmatrix -R = sigmam^2; % Messrauschen (Könnte durch Versuche bestimmt werden) - -I = eye(3); % Identity matrix (Einheitsmatrix) - -% Filterprozess - -% Initialisieren der Variablen -N = length(t); % Anzahl Punkte im Einheitsvektor (= Anzahl Messwerte) -xhat = zeros(3, N); % Matrix mit geschätzten Zuständen - -% Index ':' bedeutet: 'alles' -% Index '(1, :)' bedeutet: 'alles aus der 1. Zeile' - -% Anfangszustand setzen -xhat(:, 1) = x0; -P = P0; - -% Kalman-Matrizen konvergiert. Vorab-Berechnung in 'genügenden' Iterationen -for idx = 1:100 - Ppred = A * P * A' + Q; % Prädizieren der Kovarianz - S = (H * Ppred * H' + R); % Innovationskovarianz - K = Ppred * H' / S; % Filter-Matrix (Kalman-Gain) - P = (I - K * H) * Ppred; % Aktualisieren der Kovarianz -end - -% Anfangszustand gegeben -% Erster zu berechnender Wert ist der zweite -for idx = 2:N - % Vorhersage - xpred = A * xhat(:, idx-1); % Prädizierter Zustand aus Bisherigem und System - % Ppred = A * P * A' + Q; % Prädizieren der Kovarianz - - % Korrektur - y = xt(idx) - H * xpred; % Messungen/ Kraft aus System - Vohersage - % S = (H * Ppred * H' + R); % Innovationskovarianz - % K = Ppred * H' / S; - - xhat(:, idx) = xpred + K * y; % Aktualisieren des Systemzustands - % P = (I - K * H) * Ppred; % Aktualisieren der Kovarianz -end -\end{lstlisting} +6. Die Outputs von $k$ werden die Inputs für ${k-1}$ und werden wieder im Schritt 1 verwendet +\end{document} + + +@article{faragher_understanding_2012, + title = {Understanding the Basis of the Kalman Filter Via a Simple and Intuitive Derivation [Lecture Notes]}, + volume = {29}, + issn = {1053-5888}, + url = {http://ieeexplore.ieee.org/document/6279585/}, + doi = {10.1109/MSP.2012.2203621}, + pages = {128--132}, + number = {5}, + journaltitle = {{IEEE} Signal Processing Magazine}, + shortjournal = {{IEEE} Signal Process. Mag.}, + author = {Faragher, Ramsey}, + urldate = {2021-07-09}, + date = {2012-09} +} +@article{Wikipedia, + title = Kalmanfilter}, + url = {https://de.wikipedia.org/wiki/Kalman-Filter}, + pages = {128--132}, + number = {5}, + urldate = {2021-07-09}, +} + +@article{mueller_deconvolving_2008, + title = {Deconvolving oscillatory transients with a {Kalman} filter}, + url = {http://arxiv.org/abs/0809.4676}, + abstract = {This paper describes a method to filter oscillatory transients from measurements of a time series which were at least an order of magnitude larger than the signal to be measured. Based on a Kalman filter, it has an optimality property and a natural scaling parameter that allows to tune it to high resolution or low noise.}, + urldate = {2021-07-09}, + journal = {arXiv:0809.4676 [math]}, + author = {Mueller, Andreas}, + month = sep, + year = {2008}, + note = {arXiv: 0809.4676}, + keywords = {93E11, Mathematics - Optimization and Control}, + annote = {Comment: 12 pages, 9 figures}, -- cgit v1.2.1 From 5a77dba064b9f1a841930dcc823cb2e83a759839 Mon Sep 17 00:00:00 2001 From: Lukaszogg <82384106+Lukaszogg@users.noreply.github.com> Date: Sat, 17 Jul 2021 16:39:12 +0200 Subject: =?UTF-8?q?L=C3=B6schen=20von=20\end?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- buch/papers/erdbeben/teil1.tex | 44 +----------------------------------------- 1 file changed, 1 insertion(+), 43 deletions(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index a4e2220..52872f6 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -17,7 +17,7 @@ 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 den linearen Kalman-Filter. 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 \cdot Eigendynamik des Systems gerechnet. +Dabei wird durch eine deterministische Vorhersage, in dem der Zustand * Eigendynamik des Systems gerechnet. Die Idee dahinter ist, dass das Kalman-Filter die nicht-deterministische Grösse $f$ anhand der Messung und der Vorhersage zu bestimmen. Für mehrere Dimensionen (x,y,z) würde der Pythagoras für das System benötigt werden. @@ -262,46 +262,4 @@ Dabei beginnt das Filter mit dem Anfangszustand für $k=0$ 6. Die Outputs von $k$ werden die Inputs für ${k-1}$ und werden wieder im Schritt 1 verwendet -\end{document} - - -@article{faragher_understanding_2012, - title = {Understanding the Basis of the Kalman Filter Via a Simple and Intuitive Derivation [Lecture Notes]}, - volume = {29}, - issn = {1053-5888}, - url = {http://ieeexplore.ieee.org/document/6279585/}, - doi = {10.1109/MSP.2012.2203621}, - pages = {128--132}, - number = {5}, - journaltitle = {{IEEE} Signal Processing Magazine}, - shortjournal = {{IEEE} Signal Process. Mag.}, - author = {Faragher, Ramsey}, - urldate = {2021-07-09}, - date = {2012-09} -} -@article{Wikipedia, - title = Kalmanfilter}, - url = {https://de.wikipedia.org/wiki/Kalman-Filter}, - pages = {128--132}, - number = {5}, - urldate = {2021-07-09}, -} - -@article{mueller_deconvolving_2008, - title = {Deconvolving oscillatory transients with a {Kalman} filter}, - url = {http://arxiv.org/abs/0809.4676}, - abstract = {This paper describes a method to filter oscillatory transients from measurements of a time series which were at least an order of magnitude larger than the signal to be measured. Based on a Kalman filter, it has an optimality property and a natural scaling parameter that allows to tune it to high resolution or low noise.}, - urldate = {2021-07-09}, - journal = {arXiv:0809.4676 [math]}, - author = {Mueller, Andreas}, - month = sep, - year = {2008}, - note = {arXiv: 0809.4676}, - keywords = {93E11, Mathematics - Optimization and Control}, - annote = {Comment: 12 pages, 9 figures}, - - - - - -- cgit v1.2.1 From 51d891e986fa62c66cb18c6d558458cec41dd540 Mon Sep 17 00:00:00 2001 From: Lukaszogg <82384106+Lukaszogg@users.noreply.github.com> Date: Sat, 17 Jul 2021 16:49:21 +0200 Subject: Update references.bib --- buch/papers/erdbeben/references.bib | 70 ++++++++++++++++++++++++------------- 1 file changed, 46 insertions(+), 24 deletions(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/references.bib b/buch/papers/erdbeben/references.bib index aef5de9..56ca24b 100644 --- a/buch/papers/erdbeben/references.bib +++ b/buch/papers/erdbeben/references.bib @@ -1,35 +1,57 @@ -% -% references.bib -- Bibliography file for the paper erdbeben -% -% (c) 2020 Autor, Hochschule Rapperswil -% +%% This BibTeX bibliography file was created using BibDesk. +%% https://bibdesk.sourceforge.io/ + +%% Created for lukas zogg at 2021-07-17 16:48:19 +0200 + + +%% Saved with string encoding Unicode (UTF-8) + + + +@article{aragher_understanding_2012, + author = {Faragher, Ramsey}, + date-added = {2021-07-17 16:44:00 +0200}, + date-modified = {2021-07-17 16:45:54 +0200}, + journal = { Signal Processing Magazine}, + month = {09}, + number = {5}, + pages = {128--132}, + title = {Understanding the Basis of the Kalman Filter Via a Simple and Intuitive Derivation }, + volume = {29}, + year = {2012}, + Bdsk-File-1 = {YnBsaXN0MDDSAQIDBFxyZWxhdGl2ZVBhdGhZYWxpYXNEYXRhXxByLi4vLi4vLi4vLi4vLi4vLi4vRG93bmxvYWRzL1VuZGVyc3RhbmRpbmcgdGhlIEJhc2lzIG9mIHRoZSBLYWxtYW4gRmlsdGVyIFZpYSBhIFNpbXBsZSBhbmQgSW50dWl0aXZlIERlcml2YXRpb24ucGRmTxECbgAAAAACbgACAAAMTWFjaW50b3NoIEhEAAAAAAAAAAAAAAAAAAAAAAAAAEJEAAH/////H1VuZGVyc3RhbmRpbmcgdGhlICNGRkZGRkZGRi5wZGYAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAP////8AAAAAAAAAAAAAAAAABgACAAAKIGN1AAAAAAAAAAAAAAAAAAlEb3dubG9hZHMAAAIAci86VXNlcnM6bHVrYXN6b2dnOkRvd25sb2FkczpVbmRlcnN0YW5kaW5nIHRoZSBCYXNpcyBvZiB0aGUgS2FsbWFuIEZpbHRlciBWaWEgYSBTaW1wbGUgYW5kIEludHVpdGl2ZSBEZXJpdmF0aW9uLnBkZgAOAK4AVgBVAG4AZABlAHIAcwB0AGEAbgBkAGkAbgBnACAAdABoAGUAIABCAGEAcwBpAHMAIABvAGYAIAB0AGgAZQAgAEsAYQBsAG0AYQBuACAARgBpAGwAdABlAHIAIABWAGkAYQAgAGEAIABTAGkAbQBwAGwAZQAgAGEAbgBkACAASQBuAHQAdQBpAHQAaQB2AGUAIABEAGUAcgBpAHYAYQB0AGkAbwBuAC4AcABkAGYADwAaAAwATQBhAGMAaQBuAHQAbwBzAGgAIABIAEQAEgBwVXNlcnMvbHVrYXN6b2dnL0Rvd25sb2Fkcy9VbmRlcnN0YW5kaW5nIHRoZSBCYXNpcyBvZiB0aGUgS2FsbWFuIEZpbHRlciBWaWEgYSBTaW1wbGUgYW5kIEludHVpdGl2ZSBEZXJpdmF0aW9uLnBkZgATAAEvAAAVAAIAEP//AAAACAANABoAJACZAAAAAAAAAgEAAAAAAAAABQAAAAAAAAAAAAAAAAAAAws=}} + +@url{erdbeben:wikipedia, + author = {https://de.wikipedia.org/wiki/Kalman-Filter}, + date-added = {2021-07-17 16:42:22 +0200}, + date-modified = {2021-07-17 16:43:53 +0200}, + title = {Kalmanfilter}, + urldate = {2021-07-0}} @online{erdbeben:bibtex, + date = {2020-02-06}, + day = {6}, + month = {2}, title = {BibTeX}, url = {https://de.wikipedia.org/wiki/BibTeX}, - date = {2020-02-06}, year = {2020}, - month = {2}, - day = {6} -} + Bdsk-Url-1 = {https://de.wikipedia.org/wiki/BibTeX}} @book{erdbeben:numerical-analysis, - title = {Numerical Analysis}, author = {David Kincaid and Ward Cheney}, - publisher = {American Mathematical Society}, - year = {2002}, - isbn = {978-8-8218-4788-6}, inseries = {Pure and applied undegraduate texts}, - volume = {2} -} + isbn = {978-8-8218-4788-6}, + publisher = {American Mathematical Society}, + title = {Numerical Analysis}, + volume = {2}, + year = {2002}} @article{erdbeben:mendezmueller, - author = { Tabea Méndez and Andreas Müller }, - title = { Noncommutative harmonic analysis and image registration }, - journal = { Appl. Comput. Harmon. Anal.}, - year = 2019, - volume = 47, - pages = {607--627}, - url = {https://doi.org/10.1016/j.acha.2017.11.004} -} - + author = {Tabea M{\'e}ndez and Andreas M{\"u}ller}, + journal = {Appl. Comput. Harmon. Anal.}, + pages = {607--627}, + title = {Noncommutative harmonic analysis and image registration}, + url = {https://doi.org/10.1016/j.acha.2017.11.004}, + volume = 47, + year = 2019, + Bdsk-Url-1 = {https://doi.org/10.1016/j.acha.2017.11.004}} -- cgit v1.2.1 From 07e868a4275c0cecd4d743e9bd0c1e4f2b7c7be1 Mon Sep 17 00:00:00 2001 From: Lukaszogg <82384106+Lukaszogg@users.noreply.github.com> Date: Tue, 20 Jul 2021 17:19:47 +0200 Subject: Korrektur und Anpassungen --- buch/papers/erdbeben/teil0.tex | 47 ++++++++++++++++--------- buch/papers/erdbeben/teil1.tex | 78 ++++++++++++++++++++++++------------------ 2 files changed, 75 insertions(+), 50 deletions(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/teil0.tex b/buch/papers/erdbeben/teil0.tex index 8ac5d6d..844245c 100644 --- a/buch/papers/erdbeben/teil0.tex +++ b/buch/papers/erdbeben/teil0.tex @@ -6,16 +6,29 @@ \section{Teil 0\label{erdbeben:section:teil0}} \rhead{Erdbeben} \section{Erdbebenmessung} -\subsection{Was ist ein Erdbeben} -Fabio +subsection{Was ist ein 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. + +Unter einem Erdbeben verstehen wir eine Erschütterung des Erdkörpers. +Dabei reiben zwei tektonische Platten aneinander, welche sich durch die Gesteinsverzahnung gegenseitig blockieren. +Aufgrund dieser Haftreibung entstehen Spannungen, die sich immer mehr bis zum Tipping Point aufbauen. +Irgendwann ist der Punkt erreicht, in dem die Scherfestigkeit der Gesteine überwunden wird. +Wenn dies passiert, entlädt sich die aufgebaute Spannung und setzt enorme Energien frei, die wir als Erdbeben wahrnehmen. + +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} 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, bleibt die gekoppelte Masse stehen aber das Gehäuse schwingt mit. +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. 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 eine Gehäuse, an dem zwei Federn und eine Masse befestigt ist. +Wir konstruieren uns eine einfachere Version eines Seismographen mit eine Gehäuse, an dem zwei Federn und eine Masse befestigt sind. 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. +Dies bedeutet, unser Seismograph kann nur in eine Dimension Messwerte aufnehmen. \begin{figure} \begin{center} @@ -30,7 +43,7 @@ Wir wollen jedoch die Beschleunigung $a(t)$ des Boden bzw. die Kraft $f(t)$ welc Anhand dieser Beschleunigung bzw. 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)$ + der Eigendynamik der Masse. +Das heisst: Die Messung ist zweifach Integriert die Kraft $f(t)$ inklusive der Eigendynamik der Masse. Um die Bewegung der Masse zu berechnen, müssen wir Gleichungen für unser System finden. \subsection{Systemgleichung} @@ -40,21 +53,21 @@ Diese lautet: m\ddot s + 2k \dot s + Ds = f \end{equation} mit den Konstanten $m$ = Masse, $k$ = Dämpfungskonstante und $D$ = Federkonstante. -Um diese nun in die Systemmatrix umzuwandeln, wird die Differentialgleichung zweiter Ordnung substituiert: -\[ {x_1}=s \qquad -{x_2}=\dot s, \qquad\] -Somit entstehen die Gleichungenür die Position $s(t)$ der Masse : -\[ \dot {x_1} = {x_2}\] +Da die DGL linear ist, kann sie in die kompaktere und einfachere Matrix-Form umgewandelt werden. Dazu wird die Differentialgleichung zweiter Ordnung substituiert: +\[ {s_1}=s \qquad +{s_2}=\dot s, \qquad\] +Somit entstehen die Gleichungen für die Position $s(t)$ der Masse : +\[ \dot {s_1} = {s_2}\] und -\[ \dot x_2 = -\frac{D}{m} {x_1} -\frac{2k}{m} {x_2} + \frac{f} {m} \] für die Geschwindigkeit $v(t)$ der Masse. +\[ \dot s_2 = -\frac{D}{m} {s_1} -\frac{2k}{m} {s_2} + \frac{f} {m} \] für die Beschleunigung $a(t)$ der Masse. Diese können wir nun in der Form -\[ {x_3}=-\frac{D}{m} {s_1} -\frac{2k}{m} {s_2} + \frac{f} {m} \] +\[ {s_3}=-\frac{D}{m} {s_1} -\frac{2k}{m} {s_2} + \frac{f} {m} \] auch als Matrix-Vektor-Gleichung darstellen. Dafür wird die Gleichung in die Zustände aufgeteilt. Die für uns relevanten Zustände sind die Position der Masse, die Geschwindigkeit der Masse und die äussere Beschleunigung des ganzen System. Dabei muss unterschieden werden, um welche Beschleunigung es sich handelt. -Das System beinhaltet sowohl eine Beschleunigung der Masse (innere Beschleunigung), als auch eine Beschleunigung der ganzen Apparatur (äussere Beschleunigung). +Das System beinhaltet sowohl eine Beschleunigung der Masse, innere Beschleunigung, als auch eine Beschleunigung der ganzen Apparatur, äussere Beschleunigung. In unserem Fall wird die äusseren Beschleunigung gesucht, da diese der Erdbebenanregung gleich kommt. \begin{equation} \frac{d}{dt} \left(\begin{array}{c} {s_1} \\ {s_2} \end{array}\right) = \left( @@ -70,11 +83,13 @@ Durch Rücksubstituion ergibt sich: \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). \end{equation} Wir wissen nicht wie sich die Kraft verhält. -Deshalb treffen wir die Annahme, das sich die Kraft über die Beobachtungszeit nicht verändert. -Diese unzutreffende Annahme wird später durch einen grossen Systemfehler kompensiert. +Deshalb treffen wir die Annahme, das sich die Kraft über die Beobachtungszeit nicht verändert. +Diese Annahme ist nicht zulässig, jedoch ist dies das beste, was wir Annehmen können. +Diese unzutreffende Annahme wird späteren Berechnungen berücksichtigen werden Da die Kraft unbekannt ist, wird die letzte Zeile mit Nullen gefüllt, denn genau diese Werte wollen wir. diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index 52872f6..e07800f 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -15,7 +15,7 @@ \section{Kalman-Filter} 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 den linearen Kalman-Filter. +Dies ist eine typische Anwendung für das Kalman-Filter. 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 * Eigendynamik des Systems gerechnet. Die Idee dahinter ist, dass das Kalman-Filter die nicht-deterministische Grösse $f$ anhand der Messung und der Vorhersage zu bestimmen. @@ -27,7 +27,9 @@ Für ein nicht-lineares System werden Extended Kalman-Filter benötigt, bei dene Einfachheitshalber beschränken wir uns auf den linearen Fall, da dadurch die wesentlichen Punkte bereits aufgezeigt werden. \subsection{Geschichte} -Das Kalman-Filter wurde 1960 von Rudolf Emil Kalman entdeckt und direkt von der NASA für die Appollo Mission benutzt. Der Filter kommt mit wenig Rechenleistung aus und war somit dafür geeignet die Rakete bei der Navigation zu unterstützen. Das Filter schätzt den Zustand eines Systems anhand von Messungen und kann den nächsten Zustand errechnen. 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 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. +Das Filter schätzt den Zustand eines Systems anhand von Messungen und kann den nächsten Zustand errechnen. 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. @@ -80,7 +82,7 @@ Sie ist also gewichtet und die best mögliche Schätzung. \end{figure} -Was in 2 Dimensionen erklärt wurde, funktioniert auch in mehreren Dimensionen. +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} @@ -105,7 +107,7 @@ 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 für die Kovarianzmatrix ergibt sich +Als Initialwert für die Kovarianzmatrix ergibt sich \[ {P_0 }= @@ -127,7 +129,7 @@ Das Kalman-Filter benötigt für die Vorhersage des nächsten Zustandes eine Bes 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 Vektor-Gleichung lautet daher: +Die System-Matrix lautet daher: \[ A = \left( \begin{array}{ccc} @@ -139,10 +141,12 @@ A = \left( Dabei soll der Kalman-Filter in diskreten Zeitschritten $\Delta t$ arbeiten. 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. -Kalman-Filter berücksichtigen Unsicherheiten wie Messfehler und -rauschen. +Kalman-Filter berücksichtigen sowohl Unsicherheiten wie Messfehler und -rauschen. +In der Matrix $Q$ geht es jedoch im die Unsicherheit die der Prozess mit sich bringt. Bei unserem Modell könnte das beispielsweise ein Windstoss an die Masse sein. Für uns wäre dies: \[ @@ -158,22 +162,23 @@ Die Standabweichungen müssten statistisch ermittelt werden, da der Fehler nicht Das Bedeutet wiederum dass $Q$ die Unsicherheit des Prozesses beschreibt und nicht die der Messung. \subsubsection*{Messmatrix $H$} -Die Messmatrix gibt an, welche Parameter gemessen werden +Die Messmatrix gibt an, welche Parameter gemessen werden. +$H$ ist die Gleichung die für die Vorhersage der Messung. In unserem Falle ist es die Position der Massen. \[ H = (1, 0, 0) \] \subsubsection*{Messrauschkovarianz $R$} -Die Messrauschkovarianzmatrix beinhaltet, wie der Name es schon sagt, das Rauschen der Positionsmessung. +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_{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. \subsection{Fiter-Agorithmus} Nachdem alle Parameter aufgestellt sind, wird das Filter initialisiert. -Zuerst wird der nächste Zustand der Feder vorhergesagt, danach wird die Messung präzisiert und laufend zu aktualisieren. +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. @@ -182,14 +187,14 @@ Aus dieser Differenz und den Unsicherheiten des Prozesses ($Q$) und der Messung Im Filterschritt Vorhersage wird der nächste Zustand anhand des Anfangszustand und der Systemmatrix berechnet. Dies funktioniert mit dem Rechenschritt: \[ -{x_{k|k-1}}=\Phi \cdot {x_{k-1|k-1}}= \exp(A\Delta t)\cdot{x_{k|k-1}}. +{x_{k-1}}=\Phi \cdot {x_{k-1}}= \exp(A\Delta t)\cdot{x_{k-1}}. \] Die Kovarianz $P_{pred}$ wird ebenfalls neu berechnet. Da wir ein mehrdimensionales System haben, 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 -\[ {P_{k|k-1}} = {\Phi_k} {P_{k-1|k-1}} {\Phi_k} ^T + {Q_{k-1}} .\] -Es vergeht genau $dt$ Zeit, und dieser Vorgang wird wiederholt. +\[ {P_{k-1}} = {\Phi_k} {P_{k-1}} {\Phi_k} ^T + {Q_{k-1}} .\] +Es vergeht genau $t$ Zeit, und dieser Vorgang wird wiederholt. 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. Das Filter passt sich selber an und korrigiert sich bei grosser Abweichung. @@ -199,10 +204,10 @@ 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 -\[{w_{k}}={z_{k}}-{H_{k}}\cdot{x_{k|k-1}}.\] +\[{w_{k}}={z_{k}}-{H}\cdot{x_{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. +Die Hilfsgröße Innovation beschreibt, wie genau die Vorhersage den aktuellen Messwert mittels der Systemmatrix $\Phi$ beschreiben kann. 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 ist intuitiv logisch, eine Innovation von 0 bedeutet, dass die Messung nichts Neues hervorbrachte. @@ -210,34 +215,34 @@ Innovation = Messung - Vorhersage. Dies ist intuitiv logisch, eine Innovation vo 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_{k}}{P_{k|k-1}}{H_{k}}^T+{R_{k}} +{S_{k}}={H}{P_{k-1}}{H}^T+{R_{k}} \] \subsubsection*{Aktualisieren} -Im nächsten Schritt kommt nun die Wahrscheinlichkeit nach Gauss dazu. +Im nächsten Schritt kommt nun die Wahrscheinlichkeit dazu. \[ -{K_{k}}= {{P_{k|k-1}} \cdot {H_{k}^T}}\cdot {S_{k}}^{-1} +{K_{k}}= {{P_{k-1}} \cdot {H_{k}^T}}\cdot {S_{k}}^{-1} \] Dieser Vorgang wird Kalman-Gain genannt. Er sagt aus, welcher Kurve mehr Vertraut werden soll, dem Messwert oder der Systemdynamik. -Das Kalman-Gain wird geringer wen der Messwert dem vorhergesagten Systemzustand entspricht. -Sind die Messwerte komplett anders als die Vorhersage, wo werden die Elemente in der Matrix $K$ grösser. -Anhand der Informationen aus dem Kalman-Gain $K$ wird das System geupdated. +Das Kalman-Gain wird geringer, wenn der Messwert dem vorhergesagten Systemzustand entspricht. +Sind die Messwerte komplett anders als die Vorhersage, werden die Elemente in der Matrix $K$ grösser. +Anhand der Informationen aus dem Kalman-Gain $K$ wird das System aktualisiert. \[ -{x_{k|k}}={x_{k|k-1}}+({K_{k}}\cdot {w_{k}}) +{x_{k|k}}={x_{k-1}}+({K_{k}}\cdot {w_{k}}) \] Dazu kommt eine neue Kovarianz für den nächste Vorhersageschritt: \[ -{P_{k|k}}=(I-({K_{k}} \cdot {H_{k}})) \cdot {P_{k|k-1}} +{P_{k}}=(I-({K_{k}} \cdot {H})) \cdot {P_{k-1}} \] -Der ganze Ablauf wird nun zum Algorithmus und beginnt wieder mit der Vorhersage +Der ganze Algorithmus und beginnt wieder mit der Vorhersage \[ -{x_{k|k-1}}=\Phi \cdot {x_{k-1|k-1}}= \exp(A\Delta t)\cdot{x_{k|k-1}}. +{x_{k-1}}=\Phi \cdot {x_{k-1}}= \exp(A\Delta t)\cdot{x_{k-1}}. \] @@ -246,20 +251,25 @@ 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 -\[{x_{k|k-1}}=\Phi \cdot {x_{k-1|k-1}}= \exp(A\Delta t)\cdot{x_{k|k-1}}.\] +\[{x_{k-1}}={\Phi} \cdot {x_{k-1}}= \exp(A\Delta t)\cdot{x_{k-1}}.\] 2. Nächste Fehlerkovarianz vorhersagen -\[{P_{k|k-1}}={\Phi _{k}} {P_{k-1|k-1}} {\Phi _{k}}^T + {Q_{k-1}}.\] +\[{P_{k-1}}={\Phi} {P_{k-1}} {\Phi _{k}}^T + {Q_{k-1}}.\] -3. Das Kalman Filter anwenden -\[{K_{k}}= {P_{k|k-1}} \cdot {H_{k}^T}\cdot {S_{k}^{-1}}\] +3. Zustand wird gemessen +\[{w_{k}}={z_{k}}-{H}\cdot{x_{k-1}}.\] -4. Schätzung aktualisieren -\[{x_{k|k}}={x_{k|k-1}}+({K_{k}}\cdot {w_{k}}) \] +4. Innovation (= Messung - Vorhersage) +\[ {S_{k}}={H}{P_{k-1}}{H}^T+{R_{k}}\] -5. Fehlerkovarianz aktualisieren -\[{P_{k|k}}=(I-({K_{k}}\cdot {H_{k}})) \cdot {P_{k|k-1}} \] +5. Das Kalman Filter anwenden +\[{K_{k}}= {P_{k-1}} \cdot {H^T}\cdot {S_{k}^{-1}}\] +6. Schätzung aktualisieren +\[{x_{k}}={x_{k-1}}+({K_{k}}\cdot {w_{k}}) \] -6. Die Outputs von $k$ werden die Inputs für ${k-1}$ und werden wieder im Schritt 1 verwendet +7. Fehlerkovarianz aktualisieren +\[{P_{k}}=(I-({K_{k}}\cdot {H})) \cdot {P_{k-1}} \] + +8. Die Outputs von $k$ werden die Inputs für ${k-1}$ und werden wieder im Schritt 1 verwendet -- cgit v1.2.1 From 98b291860c3de86df4d35a7ef3d58315722c2c18 Mon Sep 17 00:00:00 2001 From: Lukaszogg <82384106+Lukaszogg@users.noreply.github.com> Date: Thu, 22 Jul 2021 18:47:29 +0200 Subject: Update teil0.tex --- buch/papers/erdbeben/teil0.tex | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/teil0.tex b/buch/papers/erdbeben/teil0.tex index 844245c..ba6552b 100644 --- a/buch/papers/erdbeben/teil0.tex +++ b/buch/papers/erdbeben/teil0.tex @@ -12,10 +12,8 @@ Das soll uns helfen, eine Verknüpfung zwischen dem Naturphänomen und der mathe Unter einem Erdbeben verstehen wir eine Erschütterung des Erdkörpers. Dabei reiben zwei tektonische Platten aneinander, welche sich durch die Gesteinsverzahnung gegenseitig blockieren. -Aufgrund dieser Haftreibung entstehen Spannungen, die sich immer mehr bis zum Tipping Point aufbauen. -Irgendwann ist der Punkt erreicht, in dem die Scherfestigkeit der Gesteine überwunden wird. +Diese Haftreibung durch die Steine wird so lange aufgebaut, bis sie nicht mehr gehalten werden kann. Wenn dies passiert, entlädt sich die aufgebaute Spannung und setzt enorme Energien frei, die wir als Erdbeben wahrnehmen. - 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. -- cgit v1.2.1 From ca9b453a796c6fb80a563d0bd979b5393acec373 Mon Sep 17 00:00:00 2001 From: Lukaszogg <82384106+Lukaszogg@users.noreply.github.com> Date: Fri, 23 Jul 2021 11:00:22 +0200 Subject: Anpassungen teil0 und main --- buch/papers/erdbeben/main.tex | 25 +++---------------------- buch/papers/erdbeben/teil0.tex | 4 +--- 2 files changed, 4 insertions(+), 25 deletions(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/main.tex b/buch/papers/erdbeben/main.tex index 8f9c8d5..95f1f4b 100644 --- a/buch/papers/erdbeben/main.tex +++ b/buch/papers/erdbeben/main.tex @@ -3,30 +3,11 @@ % % (c) 2020 Hochschule Rapperswil % -\chapter{Thema\label{chapter:erdbeben}} +\chapter{Erdbebenmessung\label{chapter:erdbeben}} \lhead{Thema} \begin{refsection} -\chapterauthor{Hans Muster} - -Ein paar Hinweise für die korrekte Formatierung des Textes -\begin{itemize} -\item -Absätze werden gebildet, indem man eine Leerzeile einfügt. -Die Verwendung von \verb+\\+ ist nur in Tabellen und Arrays gestattet. -\item -Die explizite Platzierung von Bildern ist nicht erlaubt, entsprechende -Optionen werden gelöscht. -Verwenden Sie Labels und Verweise, um auf Bilder hinzuweisen. -\item -Beginnen Sie jeden Satz auf einer neuen Zeile. -Damit ermöglichen Sie dem Versionsverwaltungssysteme, Änderungen -in verschiedenen Sätzen von verschiedenen Autoren ohne Konflikt -anzuwenden. -\item -Bilden Sie auch für Formeln kurze Zeilen, einerseits der besseren -Übersicht wegen, aber auch um GIT die Arbeit zu erleichtern. -\end{itemize} - +\chapterauthor{Lukas Zogg und +Fabio Veicelli} \input{papers/erdbeben/teil0.tex} \input{papers/erdbeben/teil1.tex} %\input{papers/erdbeben/teil2.tex} diff --git a/buch/papers/erdbeben/teil0.tex b/buch/papers/erdbeben/teil0.tex index ba6552b..8ce8ff2 100644 --- a/buch/papers/erdbeben/teil0.tex +++ b/buch/papers/erdbeben/teil0.tex @@ -3,10 +3,8 @@ % % (c) 2020 Prof Dr Andreas Müller, Hochschule Rapperswil %% -\section{Teil 0\label{erdbeben:section:teil0}} +\section{Was ist ein Erdbeben? \label{erdbeben:section:teil0}} \rhead{Erdbeben} -\section{Erdbebenmessung} -subsection{Was ist ein 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. -- cgit v1.2.1 From d11655b383e154e8ad5bb7006e33383f99e8c62c Mon Sep 17 00:00:00 2001 From: Lukaszogg <82384106+Lukaszogg@users.noreply.github.com> Date: Wed, 28 Jul 2021 15:02:47 +0200 Subject: Anpassungen nach Besprechung --- buch/papers/erdbeben/Gausskurve2.pdf | Bin 26978 -> 14941 bytes buch/papers/erdbeben/Gausskurve2.tex | 5 +- buch/papers/erdbeben/Gausskurve3.pdf | Bin 27445 -> 15413 bytes buch/papers/erdbeben/Gausskurve3.tex | 5 +- buch/papers/erdbeben/main.tex | 2 +- buch/papers/erdbeben/references.bib | 8 +- buch/papers/erdbeben/teil0.tex | 57 ++++++------ buch/papers/erdbeben/teil1.tex | 168 +++++++++++++++++++---------------- 8 files changed, 131 insertions(+), 114 deletions(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/Gausskurve2.pdf b/buch/papers/erdbeben/Gausskurve2.pdf index bee3bc0..5e4afdf 100644 Binary files a/buch/papers/erdbeben/Gausskurve2.pdf and b/buch/papers/erdbeben/Gausskurve2.pdf differ diff --git a/buch/papers/erdbeben/Gausskurve2.tex b/buch/papers/erdbeben/Gausskurve2.tex index 44319c3..2441766 100644 --- a/buch/papers/erdbeben/Gausskurve2.tex +++ b/buch/papers/erdbeben/Gausskurve2.tex @@ -1,13 +1,12 @@ \documentclass{standalone} \usepackage{pgfplots} - +\usepackage{txfonts} \pgfplotsset{compat = newest} \begin{document} - -\begin{tikzpicture} +\begin{tikzpicture}[>=latex,thick] \begin{axis}[ diff --git a/buch/papers/erdbeben/Gausskurve3.pdf b/buch/papers/erdbeben/Gausskurve3.pdf index e86a403..b86023f 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 85455ef..032d6de 100644 --- a/buch/papers/erdbeben/Gausskurve3.tex +++ b/buch/papers/erdbeben/Gausskurve3.tex @@ -1,13 +1,12 @@ \documentclass{standalone} \usepackage{pgfplots} - +\usepackage{txfonts} \pgfplotsset{compat = newest} \begin{document} - -\begin{tikzpicture} +\begin{tikzpicture}[>=latex,thick] \begin{axis}[ diff --git a/buch/papers/erdbeben/main.tex b/buch/papers/erdbeben/main.tex index 95f1f4b..4167475 100644 --- a/buch/papers/erdbeben/main.tex +++ b/buch/papers/erdbeben/main.tex @@ -4,7 +4,7 @@ % (c) 2020 Hochschule Rapperswil % \chapter{Erdbebenmessung\label{chapter:erdbeben}} -\lhead{Thema} +\lhead{Erdbeben} \begin{refsection} \chapterauthor{Lukas Zogg und Fabio Veicelli} diff --git a/buch/papers/erdbeben/references.bib b/buch/papers/erdbeben/references.bib index 56ca24b..444c82d 100644 --- a/buch/papers/erdbeben/references.bib +++ b/buch/papers/erdbeben/references.bib @@ -1,22 +1,22 @@ %% This BibTeX bibliography file was created using BibDesk. %% https://bibdesk.sourceforge.io/ -%% Created for lukas zogg at 2021-07-17 16:48:19 +0200 +%% Created for lukas zogg at 2021-07-27 17:56:45 +0200 %% Saved with string encoding Unicode (UTF-8) -@article{aragher_understanding_2012, +@article{erdbeben:aragher_understanding_2012, author = {Faragher, Ramsey}, date-added = {2021-07-17 16:44:00 +0200}, date-modified = {2021-07-17 16:45:54 +0200}, - journal = { Signal Processing Magazine}, + journal = {Signal Processing Magazine}, month = {09}, number = {5}, pages = {128--132}, - title = {Understanding the Basis of the Kalman Filter Via a Simple and Intuitive Derivation }, + title = {Understanding the Basis of the Kalman Filter Via a Simple and Intuitive Derivation}, volume = {29}, year = {2012}, Bdsk-File-1 = {YnBsaXN0MDDSAQIDBFxyZWxhdGl2ZVBhdGhZYWxpYXNEYXRhXxByLi4vLi4vLi4vLi4vLi4vLi4vRG93bmxvYWRzL1VuZGVyc3RhbmRpbmcgdGhlIEJhc2lzIG9mIHRoZSBLYWxtYW4gRmlsdGVyIFZpYSBhIFNpbXBsZSBhbmQgSW50dWl0aXZlIERlcml2YXRpb24ucGRmTxECbgAAAAACbgACAAAMTWFjaW50b3NoIEhEAAAAAAAAAAAAAAAAAAAAAAAAAEJEAAH/////H1VuZGVyc3RhbmRpbmcgdGhlICNGRkZGRkZGRi5wZGYAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAP////8AAAAAAAAAAAAAAAAABgACAAAKIGN1AAAAAAAAAAAAAAAAAAlEb3dubG9hZHMAAAIAci86VXNlcnM6bHVrYXN6b2dnOkRvd25sb2FkczpVbmRlcnN0YW5kaW5nIHRoZSBCYXNpcyBvZiB0aGUgS2FsbWFuIEZpbHRlciBWaWEgYSBTaW1wbGUgYW5kIEludHVpdGl2ZSBEZXJpdmF0aW9uLnBkZgAOAK4AVgBVAG4AZABlAHIAcwB0AGEAbgBkAGkAbgBnACAAdABoAGUAIABCAGEAcwBpAHMAIABvAGYAIAB0AGgAZQAgAEsAYQBsAG0AYQBuACAARgBpAGwAdABlAHIAIABWAGkAYQAgAGEAIABTAGkAbQBwAGwAZQAgAGEAbgBkACAASQBuAHQAdQBpAHQAaQB2AGUAIABEAGUAcgBpAHYAYQB0AGkAbwBuAC4AcABkAGYADwAaAAwATQBhAGMAaQBuAHQAbwBzAGgAIABIAEQAEgBwVXNlcnMvbHVrYXN6b2dnL0Rvd25sb2Fkcy9VbmRlcnN0YW5kaW5nIHRoZSBCYXNpcyBvZiB0aGUgS2FsbWFuIEZpbHRlciBWaWEgYSBTaW1wbGUgYW5kIEludHVpdGl2ZSBEZXJpdmF0aW9uLnBkZgATAAEvAAAVAAIAEP//AAAACAANABoAJACZAAAAAAAAAgEAAAAAAAAABQAAAAAAAAAAAAAAAAAAAws=}} diff --git a/buch/papers/erdbeben/teil0.tex b/buch/papers/erdbeben/teil0.tex index 8ce8ff2..c099340 100644 --- a/buch/papers/erdbeben/teil0.tex +++ b/buch/papers/erdbeben/teil0.tex @@ -23,6 +23,7 @@ Die Masse schwing jedoch in seiner Eigendynamik weiter. 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 eine 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. @@ -30,52 +31,52 @@ Dies bedeutet, unser Seismograph kann nur in eine Dimension Messwerte aufnehmen. \begin{center} \includegraphics[width=5cm]{papers/erdbeben/Apperatur} \caption{Aufbau des Seismographen mit Gehäuse, Masse, Federn und Sensor} + \label{erdbeben:Seismograph} \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 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, 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. 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 Bewegung der Masse zu berechnen, müssen wir Gleichungen für unser System finden. +Um die Krafteinwirkung der Masse zu berechnen, müssen wir Gleichungen für unser System finden. \subsection{Systemgleichung} -Im Fall unseres Seismographen, kann die Differentialgleichung zweiter Ordnung einer gedämpften Schwingung am harmonischen Oszillator verwendet werden. -Diese lautet: +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: \begin{equation} -m\ddot s + 2k \dot s + Ds = f +m\ddot s + 2k \dot s + Ds = f. \end{equation} -mit den Konstanten $m$ = Masse, $k$ = Dämpfungskonstante und $D$ = Federkonstante. -Da die DGL linear ist, kann sie in die kompaktere und einfachere Matrix-Form umgewandelt werden. Dazu wird die Differentialgleichung zweiter Ordnung substituiert: -\[ {s_1}=s \qquad -{s_2}=\dot s, \qquad\] -Somit entstehen die Gleichungen für die Position $s(t)$ der Masse : +wobei $m$ die Masse, $k$ die Dämpfungskonstante und $D$ die Federkonstante bezeichnet. +Da die Differentialgleichung linear ist, kann sie in die kompaktere und einfachere Matrix-Form umgewandelt werden. +Dazu verwenden wir die Subsitution: +\[ s_1 = s \qquad \text{und} \qquad s_2 = \dot s . \] +Somit entstehen die Gleichungen für die Position $ \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 $a(t)$ der Masse. - +\[ \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 -\[ {s_3}=-\frac{D}{m} {s_1} -\frac{2k}{m} {s_2} + \frac{f} {m} \] +\[ f =-\frac{D}{m} {s_1} -\frac{2k}{m} {s_2} + \frac{f} {m} \] auch als Matrix-Vektor-Gleichung darstellen. Dafür wird die Gleichung in die Zustände aufgeteilt. -Die für uns relevanten Zustände sind die Position der Masse, die Geschwindigkeit der Masse und die äussere Beschleunigung des ganzen System. -Dabei muss unterschieden werden, um welche Beschleunigung es sich handelt. -Das System beinhaltet sowohl eine Beschleunigung der Masse, innere Beschleunigung, als auch eine Beschleunigung der ganzen Apparatur, äussere Beschleunigung. -In unserem Fall wird die äusseren Beschleunigung gesucht, da diese der Erdbebenanregung gleich kommt. -\begin{equation} -\frac{d}{dt} \left(\begin{array}{c} {s_1} \\ {s_2} \end{array}\right) = \left( - \begin{array}{ccc} -0 & 1& 0 \\ -- \frac{D}{m} &-\frac{2k}{m} & \frac{1} {m}\\ -\end{array}\right) \left(\begin{array}{c} {s_1} \\ {s_2} \\ {s_3} \end{array}\right). -\end{equation} - -Durch Rücksubstituion ergibt sich: +Die für uns relevanten Zustände sind die Position der Masse, die Geschwindigkeit der Masse und die äussere Beschleunigung des ganzen Systems. + +Dabei muss unterschieden werden, um welche Beschleunigung es sich handelt. +Das System beinhaltet sowohl eine Beschleunigung der Masse (innere Beschleunigung) als auch eine Beschleunigung der ganzen Apparatur (äussere Beschleunigung). +In unserem Fall wird die äusseren Beschleunigung gesucht, da diese der Erdbebenanregung gleich kommt. +Dazu wird ein Zustandsvektor definiert: +\[ + \left(\begin{array}{c} {s_1} \\ {s_2} \\ {f} \end{array}\right). + \] +Durch Rücksubstituion ergibt sich uns folgende Systemgleichung in Matrix schreibweise, , wobei $\sot {s_1}= v$ ist: \begin{equation} -\frac{d}{dt} \left(\begin{array}{c} s(t) \\ v(t) \end{array}\right) = \left( +\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}\\ diff --git a/buch/papers/erdbeben/teil1.tex b/buch/papers/erdbeben/teil1.tex index e07800f..6c334bf 100644 --- a/buch/papers/erdbeben/teil1.tex +++ b/buch/papers/erdbeben/teil1.tex @@ -14,6 +14,8 @@ \rhead{Kalman-Filter} \section{Kalman-Filter} +Interessante Grösse ist also Integral von Überlagerung zweier Kräfte. +Wir brauchen also dir zweite Ableitung von 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. Unser Ziel ist es, anhand der Messung die eigentlich interessante Grösse $f$ zu bestimmen. @@ -23,8 +25,8 @@ Die Idee dahinter ist, dass das Kalman-Filter die nicht-deterministische Grösse Für mehrere Dimensionen (x,y,z) würde der Pythagoras für das System benötigt werden. 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. -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. Einfachheitshalber beschränken wir uns auf den linearen Fall, da dadurch die wesentlichen Punkte bereits aufgezeigt werden. +Für ein nicht-lineares System werden Extended Kalman-Filter benötigt, bei denen die System-Matrix (A) durch die Jacobi-Matrix des System ersetzt wird. \subsection{Geschichte} Das Kalman-Filter wurde 1960 von Rudolf Emil Kalman entdeckt und direkt von der NASA für die Appollo Mission benutzt. @@ -35,57 +37,60 @@ Das Filter schätzt den Zustand eines Systems anhand von Messungen und kann den Das Kalman-Filter schätzt den wahrscheinlichsten Wert zwischen Normalverteilungen. Dies bedeutet, das Filter schätzt nicht nur den Mittelwert, sondern auch die Standartabweichung. Da Normalverteilungen dadurch vollständig definiert sind, schätzt ein Kalman-Filter die gesamte Verteilungsfunktion des Zustandes. +In der Abbildung~\ref{erdbeben: Zwei Normalverteilungen} sind zwei Funktionen dargestellt. Die eine Funktion zeigt die errechnete Vorhersage des Zustands, bzw. deren Normalverteilung. Die andere Funktion zeigt die verrauschte Messung des nächsten Zustand, bzw. deren Normalverteilung. -Wie man am Beispiel der Gauss-Verteilungen unten sehen kann, ist sowohl der geschätzte Zustand als auch der gemessene Zustand normalverteilt und haben dementsprechend unterschiedliche Standardabweichungen $\sigma$ und Erwartungswerte $\mu$. - +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} - - +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. 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: - -\[ {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}} \] +\[ +{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: -\[ {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: -\[ -{y_f}(x;{\mu_f},{\sigma_f})=\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}\cdot{y_2} dx\,} - \] - +\[ +{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 : +\begin{align*} {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 -\[ {y_f}(x; {\mu_f}, {\sigma_f}) = {y_1}(x;{ \mu_1},{ \sigma_1}) {\cdot y_2}(x; {\mu_2}, {\sigma_2}), \] mit Erwartungswert \[ \mu_f = \frac{\mu_1\sigma_2^2 + \mu_2 \sigma_1^2}{\sigma_1^2 + \sigma_2^2} \] und Varianz -\[ \sigma_f^2 = \frac{\sigma_1^2 \sigma_2^2}{\sigma_1^2 + \sigma_2^2}. \] - +\[ +\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. -Sie ist also gewichtet und die best mögliche Schätzung. - - +Somit ist $\mu_f$ ist das gewichtete Mittel der beiden $\mu_{1,2}$, und die Varianzen sind die Gewichte! +Die neue Funktion ist die best mögliche Schätzung für zwei Verteilungen, welche den selben Zustand beschreiben. +Dies ist in der Abbildung~\ref{erdbeben:Gauss3} anhand der rote Funktion ersichtlich. \begin{figure} \begin{center} \includegraphics[width=5cm]{papers/erdbeben/Gausskurve3.pdf} \caption{Durch das Multiplizieren der blauen und der orangen Verteilung entsteht die die rote, optimale Funktion} + \label{erdbeben:Gauss3} \end{center} \end{figure} - - Was in zwei Dimensionen erklärt wurde, funktioniert auch in mehreren Dimensionen. Dieses Prinzip mach sich das Kalman Filter zu nutze, und wird von uns für die Erdbeben Berechnung genutzt. \section{Filter-Matrizen} +Da wir nun ein Werkzeug besitzen, dass die Beschleunigung, welche auf das Gehäuse wirkt, ermitteln kann, wird dieses nun Schritt für Schritt erklärt. Um den Kalman Filter zu starten, müssen gewisse Bedingungen definiert werden. In diesem Abschnitt werden die einzelnen Parameter und Matrizen erklärt und erläutert, wofür sie nützlich sind. @@ -94,8 +99,6 @@ In diesem Abschnitt werden die einzelnen Parameter und Matrizen erklärt und erl 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*{Anfangsfehler / Kovarianzmatrix $P$} @@ -108,7 +111,6 @@ 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 - \[ {P_0 }= \left( @@ -145,9 +147,9 @@ Die Matrix $\Phi$ beschreibt die Übergänge zwischen zeitlich aufeinanderfolgen \subsubsection*{Prozessrauschkovarianzmatrix $Q$} Die Prozessrauschmatrix teilt dem Filter mit, wie sich der Prozess verändert. -Kalman-Filter berücksichtigen sowohl Unsicherheiten wie Messfehler und -rauschen. -In der Matrix $Q$ geht es jedoch im die Unsicherheit die der Prozess mit sich bringt. -Bei unserem Modell könnte das beispielsweise ein Windstoss an die Masse sein. +Kalman-Filter berücksichtigen Unsicherheiten wie Messfehler und -rauschen. +In der Matrix $Q$ geht es jedoch um 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. Für uns wäre dies: \[ Q = \left( @@ -157,7 +159,6 @@ Q = \left( 0 & 0& {\sigma_f }^2\\ \end{array}\right) \] - Die Standabweichungen müssten statistisch ermittelt werden, da der Fehler nicht vom Sensor kommt und somit nicht vom Hersteller gegeben ist. Das Bedeutet wiederum dass $Q$ die Unsicherheit des Prozesses beschreibt und nicht die der Messung. @@ -165,13 +166,15 @@ Das Bedeutet wiederum dass $Q$ die Unsicherheit des Prozesses beschreibt und nic Die Messmatrix gibt an, welche Parameter gemessen werden. $H$ ist die Gleichung die für die Vorhersage der Messung. In unserem Falle ist es die Position der Massen. - -\[ H = (1, 0, 0) \] +\[ +H = (1, 0, 0) +\] \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_{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. @@ -182,19 +185,25 @@ Zuerst wird der nächste Zustand der Masse vorhergesagt, danach wird die Messung Das Filter berechnet aufgrund der aktuellen Schätzung eine Vorhersage. Diese wird, sobald verfügbar, mit der Messung verglichen. Aus dieser Differenz und den Unsicherheiten des Prozesses ($Q$) und der Messung ($R$) wird der wahrscheinlichste, neue Zustand geschätzt. +Dabei muss genau auf den Index geachtet werden. Nach dem Artikel~\cite{erdbeben:wikipedia} ist die Indexierung so genormt: +Der Zeitschritt wird mit $k$ definiert, $k-1$ ist somit ein Zeitschritt vor $k$. +Auf der linken Seite von | wird der aktuelle Zustand verlangt, bzw. ausgegeben, auf der rechten Seiten den bisherigen Zustand. +Dies bedeutet, dass die Notation $x_{n|m}$ die Schätzung von $x$ zum Zeitpunkt $n$ bis und mit zur Zeitpunkt $m \leq \ n$ präsentiert. \subsubsection*{Vorhersage} Im Filterschritt Vorhersage wird der nächste Zustand anhand des Anfangszustand und der Systemmatrix berechnet. Dies funktioniert mit dem Rechenschritt: -\[ -{x_{k-1}}=\Phi \cdot {x_{k-1}}= \exp(A\Delta t)\cdot{x_{k-1}}. - \] - -Die Kovarianz $P_{pred}$ wird ebenfalls neu berechnet. Da wir ein mehrdimensionales System haben, kommt noch die Prozessunsicherheit $Q$ dazu, so dass die Unsicherheit des Anfangsfehlers $P$ laufend verändert. +\[ +{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 -\[ {P_{k-1}} = {\Phi_k} {P_{k-1}} {\Phi_k} ^T + {Q_{k-1}} .\] -Es vergeht genau $t$ Zeit, und dieser Vorgang wird wiederholt. +\[ +{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. Das Filter passt sich selber an und korrigiert sich bei grosser Abweichung. @@ -202,74 +211,83 @@ 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 - -\[{w_{k}}={z_{k}}-{H}\cdot{x_{k-1}}.\] - +Hier bei wird lediglich die Messung mit dem Fehler behaftet, und die Messmatrix $H$ mit der Vorhersage multipliziert. +\[ +{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. 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 ist intuitiv logisch, eine Innovation von 0 bedeutet, dass die Messung nichts Neues hervorbrachte. +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-1}}{H}^T+{R_{k}} - \] +{S_{k}}={H}{P_{k|k-1}}{H}^T+{R_{k}} +\] \subsubsection*{Aktualisieren} Im nächsten Schritt kommt nun die Wahrscheinlichkeit dazu. -\[ -{K_{k}}= {{P_{k-1}} \cdot {H_{k}^T}}\cdot {S_{k}}^{-1} - \] +\[{K_{k}}= {P_{k|k-1}} {H^T}{S_{k}^{-1}}\] Dieser Vorgang wird Kalman-Gain genannt. -Er sagt aus, welcher Kurve mehr Vertraut werden soll, dem Messwert oder der Systemdynamik. -Das Kalman-Gain wird geringer, wenn der Messwert dem vorhergesagten Systemzustand entspricht. -Sind die Messwerte komplett anders als die Vorhersage, werden die Elemente in der Matrix $K$ grösser. -Anhand der Informationen aus dem Kalman-Gain $K$ wird das System aktualisiert. +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. -\[ -{x_{k|k}}={x_{k-1}}+({K_{k}}\cdot {w_{k}}) - \] +Anhand der Informationen aus dem Kalman-Gain $K$ wird das System aktualisiert. +\[ +{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: - -\[ -{P_{k}}=(I-({K_{k}} \cdot {H})) \cdot {P_{k-1}} - \] - +\[ +{P_{k|k}}=(I-{K_{k}}{H}){P_{k|k-1}} +\] Der ganze Algorithmus und beginnt wieder mit der Vorhersage - -\[ -{x_{k-1}}=\Phi \cdot {x_{k-1}}= \exp(A\Delta t)\cdot{x_{k-1}}. - \] - +\[ +{x_{k|k-1}}=\Phi{x_{k-1|k-1}}= \exp(A\Delta t){x_{k|k-1}}. +\] \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 -\[{x_{k-1}}={\Phi} \cdot {x_{k-1}}= \exp(A\Delta t)\cdot{x_{k-1}}.\] +\[ +{x_{k|k-1}}=\Phi{x_{k-1|k-1}}= \exp(A\Delta t){x_{k-1|k-1}}. +\] 2. Nächste Fehlerkovarianz vorhersagen -\[{P_{k-1}}={\Phi} {P_{k-1}} {\Phi _{k}}^T + {Q_{k-1}}.\] +\[ +{P_{k|k-1}}=\Phi {P_{k-1|k-1}} {\Phi _{k}}^T + {Q_{k-1}}. +\] 3. Zustand wird gemessen -\[{w_{k}}={z_{k}}-{H}\cdot{x_{k-1}}.\] +\[ +{w_{k}}={z_{k}}-{H}{x_{k|k-1}}. +\] 4. Innovation (= Messung - Vorhersage) -\[ {S_{k}}={H}{P_{k-1}}{H}^T+{R_{k}}\] +\[ +{S_{k}}={H}{P_{k|k-1}}{H}^T+{R_{k}} +\] 5. Das Kalman Filter anwenden -\[{K_{k}}= {P_{k-1}} \cdot {H^T}\cdot {S_{k}^{-1}}\] +\[ +{K_{k}}= {P_{k|k-1}} {H^T}{S_{k}^{-1}} +\] 6. Schätzung aktualisieren -\[{x_{k}}={x_{k-1}}+({K_{k}}\cdot {w_{k}}) \] +\[ +{x_{k|k}}={x_{k|k-1}}+{K_{k}}{w_{k}} +\] 7. Fehlerkovarianz aktualisieren -\[{P_{k}}=(I-({K_{k}}\cdot {H})) \cdot {P_{k-1}} \] +\[ +{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 -- cgit v1.2.1 From 4d8e9b6051dcd25c34b6270c1fc1938668e7df6d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20M=C3=BCller?= Date: Wed, 28 Jul 2021 18:05:37 +0200 Subject: fix files broken by JODBaer pull request --- buch/papers/erdbeben/teil0.tex | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'buch/papers/erdbeben') diff --git a/buch/papers/erdbeben/teil0.tex b/buch/papers/erdbeben/teil0.tex index c099340..c985713 100644 --- a/buch/papers/erdbeben/teil0.tex +++ b/buch/papers/erdbeben/teil0.tex @@ -74,7 +74,7 @@ Dazu wird ein Zustandsvektor definiert: \[ \left(\begin{array}{c} {s_1} \\ {s_2} \\ {f} \end{array}\right). \] -Durch Rücksubstituion ergibt sich uns folgende Systemgleichung in Matrix schreibweise, , wobei $\sot {s_1}= v$ ist: +Durch Rücksubstituion ergibt sich uns folgende Systemgleichung in Matrix schreibweise, , wobei $\dot {s_1}= v$ ist: \begin{equation} \frac{d}{dt} \left(\begin{array}{c} s(t) \\ v(t) \\ f(t) \end{array}\right) = \left( \begin{array}{ccc} -- cgit v1.2.1