Contents
Vorwort
Im Rahmen dieses Artikels soll die korrekte Anwendung des Gauß-Helmert-Modells anhand der Bestimmung ausgleichender Ebenen vorgeführt werden. Grundlegendes zur Methode der kleinsten Quadrate finden Sie (hier). Eine allgemeine Herleitung des Gauß-Helmert-Modells kann (hier) nachgelesen werden.
Modellierung der Ebene im Gauß-Helmert Modell
Je nachdem ob man in den Modellmatrizen auf den Distanzparameter verzichtet, kann die Ebene durch drei oder vier Parameter modelliert werden. Die explizite Einführung des Distanzparameters ist nicht zwingend notwendig, da dieser aus dem Normalenvektor und den Punkten abgeleitet werden kann. Dennoch wird der Distanzparameter in diesem Artikel modelliert. Diese Form hat mehrere Vorteile. Ein wesentlicher Vorteil ist dass durch die Einführung einer Bedingung die zu invertierende Matrix sich im Gauß-Helmert-Modell auf die Dimension 5 reduziert und nicht mehr linear mit der Anzahl verwendeter Punkte wächst.
Hesse’sche Normalenform

Ebene in der Hesseschen Normalenform1
Die Ebene ist als zweidimensionaler Vektorraum im unter anderem durch drei nicht identische Punkte
eindeutig definiert die nicht auf nicht auf einer Gerade liegen. Bildet man nun ein Vektorprodukt
(1)
so ist n derjenige Vektor, der senkrecht auf der Ebene steht. Den kürzesten Abstand der Ebene zum Ursprung d erhält man durch Projektion eines beliebigen Punktes auf den normierten Normalenvektor
(2)
In der Hesse’schen Normalenform wird die Lage der Ebene eindeutig mit dem normierten Normalenvektor und dem Abstand der Ebene zum Ursprung d beschrieben.
(3)
Strenge Lösung der mit vier Parametern
Ausmultipliziert ergibt die Ebenengleichung (3) bei Berücksichtigung der unterstellten Verbesserungen im Gauß-Helmert-Modell die Bedingungsgleichung der Ebene
(4)
Die Jacobi- oder Modellmatrizen und
werden an der Stelle
( Näherungswerte aus der aktuellen Iteration ) differenziert und erhalten die folgende Form:
(5)
Wir erhalten bei vier Unbekannten für eine Matrix mit vier Spalten. Die Anzahl der Zeilen hängt von der Anzahl der Punkte ab. Da Verbesserungen berücksichtigt wurden, erscheinen diese in der Modellmatrix
(6)
Die Modellmatrix wird als Bandmatrix mit drei horizontal zusammengesetzten und nicht überlappenden Diagonalmatrizen definiert. Jedes Band beinhaltet nur einen Wert. Alle drei Bänder haben n > 3 Elemente und beginnen jeweils in der ersten Zeile.
(7)
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.
(8)
Explizit ausgeschrieben, lautet der Widerspruchsvektor,
(9)
Dieser Vektor kann nach der Auflösung der Klammern algebraisch so zusammengefasst werden, dass eine reduzierte Form entsteht, die keinerlei Verbesserungen beinhaltet.
(10)
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)2 an einem Beispiel der ausgleichenden Geraden im strengen Gauß-Helmert Modell hin.
Einführung der Zwangsbedingung
Da man die Bedingungsgleichung (4) mit einer beliebigen Zahl ungleich Null multiplizieren kann, gibt es unendlich viele solcher Gleichungen, die die Ebene beschreiben. Das Normalgleichungssystem bleibt in diesem Fall singulär. Um diese Mehrdeutigkeit zu lösen, wird gefordert, dass der Betrag des Normalenvektors der Ebene die Länge 1 annimmt. Hierzu wird Bedingungsvektor (engl. constraint) eingeführt
(11)
Dieser wird an der Stelle differenziert, weil er eine Zwangsbedingung zwischen den Parametern beschreibt. Somit lautet unser Constraint-Vektor
(12)
Der dazugehörige skalare Widerspruch beschreibt, wie groß die Abweichung von der Länge 1 ist
(13)
Bildung der Blockmatrizenform
Für die Ausgleichung wird die folgende Blockmatrizenform gebildet
(14)
Der nachfolgende Code verwendet den Levenberg-Marquard Algorithmus zur Lösung des Blockmatrizensystems.
Quellcode: Lösung mit vier Parametern in OCTAVE
Download hier:
best_fit_plane.m
%% Berechnung ausgleichender Ebenenparameter im Gauss-Helmert-Modell
%
% v1.05
%
% Kisser Waldemar (2011)
% http://kisser.online/ausgleichung/ghm/ebene
clc
clear all
% Liest eine Koordinatendatei im folgenden Format ein
% 0.0 0.0 0.0
% ... ... ...
% xi yi zi
% Punkte = dlmread( 'C:\daten\boppkrauss.xyz' );
% oder
Punkte = ...
[
96.993 42.615 12.749
97.320 44.299 13.899
97.451 43.588 13.487
96.497 43.323 13.098
96.704 42.596 12.605
96.859 42.915 12.838
];
%% Hilfsvariablen
% Anzahl der Beobachtungen , Dimension
[ n , d ] = size( Punkte );
% nx, ny, nz, d
X_dach = rand( 4 , 1 );
X_dach = X_dach ./ norm( X_dach( 1 : 3 ) );
X_dach(4) = 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( 3 * 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:3) = X_dach( 1 : u - 1 ) ./ norm( X_dach( 1 : u - 1 ) );
% Bildung der B-Modellmatrix
% Da diese 3(n^2-n) Nullelemente hat, eignet sich der Sparse
% Datentyp fuer duennbesetzte Matrizen.
B = spdiags( repmat( X_dach( 1 : 3 ).' , n , 1 ) , [ 0 , n , 2 * n ] , n , 3 * 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
%% Ausgabe
fprintf( '--- AUSGLEICHUNGSPROTOKOLL ---' );
fprintf( '\n\n-- Ausgleichungsvorgang --' );
fprintf( '\nModell\t\t\t\t: Gauss-Helmert-Modell' );
fprintf( '\nVersion\t\t\t\t: 1.05' );
fprintf( '\nForm\t\t\t\t: 4 Parameter Ebene' );
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' , u );
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 );
if X_dach( end ) < 0
X_dach = - X_dach;
end
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( [ '\nnz\t\t: %12.8f\t' , char(177) , '%12.8f' ] , X_dach( 3 ) , s0_post * sqrt( Qxdach( 3 , 3 ) ) );
fprintf( [ '\nd\t\t: %12.8f\t' , char(177) , '%12.8f\n' ] , X_dach( 4 ) , s0_post * sqrt( Qxdach( 4 , 4 ) ) );
else
fprintf( '\nKonvergenz: Nicht erfolgt\n' );
end
Quellen und Fußnoten
- Benannt nach dem Mathematiker Otto Hesse (1811 – 1874) [↩]
- 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 [↩]