Contents
Vorwort
Im Rahmen dieses Artikels soll die Anwendung des Gauß-Helmert-Modells anhand der Bestimmung ausgleichender Geraden vorgeführt werden. Grundlegendes zur Methode der kleinsten Quadrate finden Sie (hier). Die allgemeine Herleitung des Gauß-Helmert-Modells kann (hier) nachgelesen werden.
Modellierung der Gerade im Gauß-Helmert Modell
Es gibt mehrere Möglichkeiten eine Gerade in einem ebenen Koordinatensystem zu beschreiben. Mit dem Blick auf die anschließende Parameterschätzung wird im weiteren Verlauf zunächst eine Gerade in der Punkt-Steigung-Form und danach in der Normalenform besprochen.
Punkt-Steigungs-Form
Die Gerade wird in der euklidischen (ebenen) Geometrie als die kürzeste Verbindung zwischen zwei Punkten definiert. Aus ihrem mathematischen Zusammenhang
(1)
kann über algebraische Umformung zunächst der Wertebereich von in Abhängigkeit der übrigen Parameter angegeben werden.
(2)
Betrachtet man nun den Winkel zwischen der Abszisse und der Gerade so kann
auch mit Hilfe von trigonometrischen Funktionen ausgedrückt werden
(3)
Nach dieser Substitution entsteht die Punkt-Steigung-Form
(4)
Da eine Gerade in jedem Punkt eine konstante Steigung besitzt ist auch der Winkel und damit auch der Ausdruck
konstant und kann zusammenfassend als Parameter
geführt werden so dass nun die folgende Form entsteht
(5)
Betrachtet man diesen Term für einen Punkt im Ursprung der Abszisse so dass , dann kann
auch als der Schnitt mit der Ordinate betrachtet und als Parameter
geführt werden. In (5) wurde bereits gezeigt dass die Steigung der Gerade überall konstant ist, auch an der Stelle
womit
auch vernachlässigt werden kann. So kann nun eine lineare Funktion in Abhängigkeit der beiden Parameter
und
, der Steigung der Gerade und dem Schnitt mit der Ordinate angegeben werden
(6)
Überführt man diese Funktion in das Gauß-Helmert-Modell so entsteht unter Berücksichtung der Verbesserungen die endgültige Bedingungsfunktion. Vergleiche hierzu NEITZEL, F., PETROVIC, S. (2008)1
(7)
Diese Form genießt den Vorteil dass ihre Parameter ohne jegliche Umrechnung intuitiv nachvollziehbar sind. Eine ausführliche Herleitung sowie ein nummerisches Beispiel zu dieser Bedingungsfunktion kann der o.g. Publikation von NEITZEL, F., u. PETROVIC, S. entnommen werden. Dieser Artikel ist besonders empfehlenswert da unter anderem auch auf die Tücken und Folgen des falsch angewandten Gauß-Helmert-Modells eingegangen wird.
Hessesche Normalenform
Die Gerade kann in der Ebene auch durch zwei nicht identische Punkte und
eindeutig beschrieben werden. Betrachtet man nun
als denjenigen Vektor, der senkrecht auf dem Richtungsvektor
–
der Gerade steht, so kann durch Projektion eines beliebigen auf der Gerade liegenden Punktes auf
, der kürzeste Abstand
der Gerade zum Ursprung berechnet werden.
(8)
In der Hesseschen Normalenform wird die Lage der Gerade eindeutig mit dem normierten Normalenvektor und dem Abstand der Gerade zum Ursprung
beschrieben.
(9)
Überführt man diese Funktion in das Gauß-Helmert-Modell so entsteht unter Berücksichtigung der Verbesserungen die endgültige Bedingungsfunktion
(10)
Beschrieben in dieser Form sind die Parameter nicht so intuitiv wie in der Punkt-Steigung-Form. Allerdings hat diese Form keinerlei nummerische Probleme wenn die Steigung sehr groß ist, da hier die Steigungen in der Abszisse und Ordinate getrennt betrachtet werden.
In vielen Ausgleichungsfällen wie z. B. der Zentrumsbestimmung von Zielmarken können beliebig gelagerte (auch senkrechte) Geraden im Datenmatrial vorliegen, die meistens als eine stützende Geometrie in die Ausgleichung eingeführt werden. In diesen Fällen eignet sich die Normalenform, weil sie durch die getrennte Betrachtung der Steigungen in der Abszisse und Ordinate keine nummerischen Probleme verursacht. Auch können Zusatzbedingungen wie z. B. senkrecht aufeinander stehende Geraden oder Gerade durch einen erzwungenen Punkt (Mittelpunkt der Zielmarke) mittels einfacher Vektoralgebra in die Ausgleichung eingeführt werden.
Umrechnung aus der Normalenform
Im ersten Schritt wird Winkel , zwischen dem normierten Normalenvektor
und der Y-Achse bestimmt
(11)
Aus der Winkelsumme des ebenen Dreiecks und dem nun bekannten Winkel kann wiederum der Winkel
berechnen
(12)
Da nun alle Winkel des Dreiecks bestimmt sind, lässt sich aus dem Sinussatz auch der Schnitt mit der Ordinate im Nullpunkt der Abszisse angeben
(13)
wobei die folgenden Substitutionen verwendet werden
(14)
(15)
Die Steigung der Gerade (3) lässt sich aus der Winkelsumme berechnen. Aus
(16)
folgt durch Substitution gemäß (12)
(17)
und (11) die gesuchte Steigung der Gerade
(18)
welche algebraisch vereinfacht auch wie folgt angegeben werden kann
(19)
Strenge Gauß-Helmert-Lösung der Normalenform
Die Jacobi- oder Modellmatrizen und
werden an der Stelle
differenziert und erhalten die folgende Form:
(20)
Wir erhalten bei drei Unbekannten für eine Matrix mit drei Spalten. Die Anzahl der Zeilen hängt von der Anzahl gegebener Punkte ab. Da Verbesserungen berücksichtigt wurden, erscheinen diese in der Modellmatrix
(21)
Die Modellmatrix wird als Bandmatrix mit zwei horizontal zusammengesetzten und nicht überlappenden Diagonalmatrizen definiert. Jedes Band beinhaltet nur einen Wert. Beide Bänder haben n > 2 Elemente und beginnen jeweils in der ersten Zeile.
(22)
Durch die Bildung der Ableitung entfallen die Verbesserungen. Dies ist ein Sonderfall und darf nicht verallgemeinert ohne eine explizite Überprüfung auf andere Modelle übertragen werden.
(23)
Explizit ausgeschrieben, lautet der Widerspruch des Gleichungssystems,
(24)
Dieser Vektor kann nach der Auflösung der Klammern algebraisch so zusammengefasst werden, dass eine reduzierte Form entsteht, die keinerlei Verbesserungen beinhaltet.
(25)
An dieser Stelle muss noch einmal darauf hingewiesen werden, dass es im Allgemeinen nicht üblich ist, dass die Verbesserungen in der strengen Lösung im Widerspruchsvektor entfallen. Auf diesen Sonderfall weisen bereits NEITZEL, PETROVIC (2008) am Beispiel der ausgleichenden Geraden hin.
Einführung der Zwangsbedingung
Da man den Normalenvektor der Bedingungsgleichung (10) mit einer beliebigen Zahl ungleich Null multiplizieren kann, gibt es unendlich viele solcher Gleichungen, die jeweils die gleiche Gerade beschreiben. Das Normalgleichungssystem bleibt in diesem Fall singulär. Geometrisch gesehen entspricht dies, einer Verlängerung des Normalenvektors mit gleichbleibender Orientierung.
Um diese Mehrdeutigkeit zu lösen fordern wir dass der Betrag des Normalenvektors der Gerade die Länge 1 annimmt und führen hierzu eine Zwangsbedingung (engl. constraint, siehe Lösung mit Nebenbedingungen ) ein
(26)
Diese wird nur an der Stelle differenziert, weil hiermit eine Zwangsbedingung zwischen den Parametern beschreiben. Somit lautet unser Constraint-Vektor
(27)
Der dazugehörige skalare Widerspruch beschreibt, wie groß die Abweichung von der Länge 1 ist
(28)
Quellcode
Nachfolgend sei der Octave Quellcode für die strenge Lösung der vier Parameter Ebene gelistet. Für die Bildung der Matrizen wurde auf Octave-Matrizenfunktionen zurückgegriffen, da die Bildung über Summen, die in Octave z. B. als for()-Schleifen realisiert werden können, wesentlich ineffizienter ist. Generell sollte man in Octave, soweit es möglich ist, auf Matrizenoperationen zurückgreifen, da diese intern über hocheffiziente, nummerisch stabile und teilweise mehrprozessoroptimierte BLAS2 und LAPACK Bibliotheken realisiert und nicht über die wesentlich langsamere M-Scriptsprache berechnet werden. Diese Behauptung kann natürlich nicht allgemein auf andere Programmiersprachen wie C/C++ übertragen werden, die auch mit for()-Schleifen sehr effizient eingesetzt werden können.
Quellcode als Octave Datei herunterladen:
Download hier: best_fit_line
%% Berechnung ausgleichender Geradenparameter im Gauss-Helmert-Modell
%
% v1.03
%
% Kisser Waldemar (2011)
% http://kisser.online/ausgleichung/ghm/gerade
clc;
clear all;
%% Datensatz waehlen
% Z.B. aus Datei laden: "Punkte = dlmread( 'C:/daten/linefit/daten.xy' );"
% oder
Punkte = ...
[ ...
0.0 0.0
1.0 1.0
2.0 4.0
3.0 9.0
];
%% Hilfsvariablen
% Anzahl der Beobachtungen , Dimension
[ n , d ] = size( Punkte );
% Initialisierung der Naeherungswerte
% nx, ny, d
X_dach = rand( 3 , 1 );
X_dach( 1 : 2 ) = X_dach( 1 : 2 ) ./ norm( X_dach( 1 : 2 ) );
X_dach( 3 ) = norm( mean( Punkte ) );
%% Schaetzung der exakten Minimierer X0,V0
% Vektor der Verbesserungen
v = zeros( d * n , 1 );
% Anzahl der Unbekannten
u = size( X_dach , 1 );
% Bedingungsvektor
c = 1;
C = zeros( c , u );
% Anzahl der Bedingungen
b = n * d;
% Legt den gewuenschten Konvergenzbetrag fest ab dem die
% Iterationen abgebrochen werden koennen.
EPSILON = 1e-12;
% Standardabweichung der Gewichtseinheit a priori
s0_prio = 1.0;
% Kofaktormatrix (hier gleichgenaue unkorrelierte Beobachtungen)
Cl = eye(2*n);
% Kovarianzmatrix
Ql = 1 / s0_prio^2 * Cl;
% Maximale Iterationsanzahl
maxit = 20;
% Beginn des Iterationsvorgangs
for iteration = 0 : maxit
% Modellmatrix A
A = [ Punkte + reshape( v , n , d ) , ones( n , 1 ) ];
% Bildung des Bedingungsvektors
C( 1 : u - 1 ) = X_dach( 1 : u - 1 ) ./ norm( X_dach( 1 : u - 1 ) );
% Bildung der B-Modellmatrix
% Da diese 2(n^2-n) Nullelemente hat, eignet sich der Sparse
% Datentyp fuer duennbesetzte Matrizen.
B = spdiags( repmat( X_dach( 1 : 2 ).' , n , 1 ) , [ 0 , n ] , n , 2 * n );
% Bildung der Inversen zur B*Ql*B^T Matrix
BQBT = B * Ql * B.';
% Bildung der Normalgleichungsmatrix
N = [ -A.' *( BQBT \ A ) , C.' ; C , zeros( c , c ) ];
Ncond = rcond( N );
Ndet = abs( det( N ) );
if ( ( Ndet < EPSILON ) || ( Ncond < EPSILON ) )
fprintf( 'Fehler: Normalgleichungsmatrix singul?r oder\n' );
fprintf( 'schlecht konditioniert\n' );
fprintf( 'Determinante : %.3e\n' , Ndet );
fprintf( 'Konditionierung : %.3e\n\n' , Ncond );
return
end
% Bilde den Widerspruchsvektor
wm = [ Punkte , ones( n , 1 ) ] * X_dach;
% Bilde den Widerspruchsskalar
wc = 1 - C * X_dach;
% Berechne die Differenzloesung
dx = N \ [ A.' * ( BQBT \ wm ) ; wc ];
% Berechne die Verbesserungen dieser Iteration
v = Ql * B.' *( BQBT \ ( wm - A * dx( 1 : u ) ) );
% Addiere die Differenzloesung
X_dach = X_dach + dx( 1 : u );
% Zweite Abbruchbedingung
if max( abs( dx( 1 : u ) ) ) < EPSILON
break
end
end
%% Normalenvektor gemaess Hesseschen Normalenform
vz = 1.0;
if X_dach( end ) < 0.0
X_dach = - X_dach;
vz = -1.0;
end
%% Ausgabe
fprintf( '--- AUSGLEICHUNGSPROTOKOLL ---' );
fprintf( '\n\n-- Ausgleichungsvorgang --' );
fprintf( '\nModell\t\t\t\t: Gauss-Helmert-Modell' );
fprintf( '\nVersion\t\t\t\t: 1.03' );
fprintf( '\nForm\t\t\t\t: 3 Parameter Gerade' );
if iteration < maxit
fprintf( '\nKonvergenz\t\t\t: Erfolgt' );
fprintf( '\nKonvergenzgrenze\t: %.1e' , EPSILON );
fprintf( '\nAnzahl Iterationen\t: %d' , iteration );
fprintf( '\n\n-- Statistik --' );
redundanz = n - u + c;
vTPv = v.' * ( Ql \ v );
s0_post = sqrt( vTPv / redundanz );
Qxdach = inv( - N );
fprintf( '\nAnzahl Beobachtungen\t: %d' , n );
fprintf( '\nAnzahl Parameter\t\t: %d' , size( X_dach , 1 ) );
fprintf( '\nAnzahl Bedingungen\t\t: %d' , c );
fprintf( '\nGesamtredundanz\t\t\t: %d' , redundanz );
fprintf( '\n\nvTPv\t: %12.8f' , vTPv );
fprintf( '\ns0_prio\t: %12.8f' , s0_prio );
fprintf( '\ns0_post\t: %12.8f' , s0_post );
fprintf( '\n\n-- Parameter --' );
fprintf( [ '\nnx\t\t: %12.8f\t' , char(177) , '%12.8f' ] , X_dach( 1 ) , s0_post * sqrt( Qxdach( 1 , 1 ) ) );
fprintf( [ '\nny\t\t: %12.8f\t' , char(177) , '%12.8f' ] , X_dach( 2 ) , s0_post * sqrt( Qxdach( 2 , 2 ) ) );
fprintf( [ '\nd\t\t: %12.8f\t' , char(177) , '%12.8f' ] , X_dach( 3 ) , s0_post * sqrt( Qxdach( 3 , 3 ) ) );
%% Berechnung der Punkt-Steigung-Form
if abs( X_dach( 2 ) ) > EPSILON
a0 = tan( acos( X_dach( 2 ) ) );
b0 = vz * X_dach( 3 ) / X_dach( 2 );
%% Varianzfortpflanzung
A2 = ...
[ ...
0 , - vz / X_dach( 2 )^2 / sqrt( 1 - X_dach( 2 )^2 ) , 0 ;
0 , - vz * X_dach( 3 ) / X_dach( 2 )^2 , vz / X_dach( 2 )
];
sab = s0_post * sqrt( diag( A2 * Qxdach( 1 : u , 1 : u ) * A2.' ) );
fprintf( '\n\n-- 2D-Geradenparameter --' );
fprintf( [ '\na\t\t: %12.8f\t' , char(177) , '%12.8f' ] , a0 , sab( 1 ) );
fprintf( [ '\nb\t\t: %12.8f\t' , char(177) , '%12.8f' ] , b0 , sab( 2 ) );
end
fprintf( '\n' );
else
fprintf( '\nKonvergenz: Nicht erfolgt\n' );
end
Quellen und Fußnoten
- NEITZEL, F., PETROVIC, S. (2008): Total Least Squares (TLS) im Kontext der Ausgleichung nach kleinsten Quadraten am Beispiel der ausgleichenden Geraden. Zeitschrift für Geodäsie, Geoinformation und Landmanagement 3/2008 S. 141-148, Wißner-Verlag, Augsburg [↩]
- Basic Linear Algebra Subprograms ist eine Softwarebibliothek für elementare Matrizenoperationen. Dank der starken Unterstützung von Hardwareentwicklern wie Intel oder AMD ist BLAS derzeit ein Quasi-Standard und bildet damit eine feste Teilmenge von erweiterten nummerischen Bibliotheken wie LAPACK/Intel MKL/GSL/CUDA SDK und vielen mehr. [↩]