Vorwort
Im Rahmen dieses Artikels soll die Anwendung des Gauß-Helmert-Modells anhand der Bestimmung ausgleichender Kugeln vorgeführt werden. Grundlegendes zur Methode der kleinsten Quadrate finden Sie (hier). Die allgemeine Herleitung des Gauß-Helmert-Modells (GHM) kann (hier) nachgelesen werden.
Modellierung der Kugel im Gauß-Helmert-Modell
Eine Kugel kann unter anderem auch durch seinen Mittelpunkt und Radius eindeutig beschrieben werden. Jeder Punkt der auf der Kugel liegt, erfüllt die Gleichung
(1)
womit dieser Term den funktionalen Zusammenhang zwischen den Unbekannten
(2)
und den Beobachtungen
(3)
beschreibt. Überführt man die Kugelgleichung (1) unter Berücksichtigung der Verbesserungen
(4)
in das GHM so entsteht die folgende Bedingungsfunktion
(5)
Dieser Term ist nichtlinear. In der Ausgleichungsrechnung ist es üblich nichtlineare Funktionen durch eine Linearisierung an der Stelle der Näherungswerte für die Verbesserungen und Parameter
(der aktuellen Iteration) zu linearisieren um anschließend in einem iterativen Prozess eine Lösung zu bestimmen.
Bildung der Blockmatrizen im GHM
Die Jacobimatrizen welche nach der Linearisierung von (1) entstehen, sehen wie folgt aus
(6)
mit der Hilfsvariable
(7)
Weil die Lösung iterativ bestimmt wird, können Taylor Glieder höherer Ordnung vernachlässigt werden. Die Modelllmatrix wird als eine Bandmatrix mit drei an einander gehängten nicht überlappenden Diagonalen gebildet.
(8)
Da man eine Linearisierung vornimmt muss auch ein Widerspruchsvektor eingeführt werden.
(9)
Alternative Bedingungsfunktion im GHM
Alternativ kann auch eine eine andere Form der Kreisgleichung verwendet werden
(10)
Diese ist zwar im streng mathematischen Sinne äquivalent zu (5), nummerisch betrachtet haben die Matrizen der partiellen Ableitungen
(11)
(12)
und der Widerspruchspruchvektor
(13)
deutlich größere Beträge, welche je nach verwendetem Datentyp zu einer geringeren Rechenschärfe führen können, da nicht alle reele Zahlen auf Fließkommazahlen gemäß IEEE 7541 abgebildet werden können2.
Näherungswertbestimmung im GHM
Das vorgestellte Gleichungssystem erfordert a priori Informationen über die Parameter welche nicht immer gegeben sein müssen. Ein sehr einfacher Ansatz wäre die Annahme dass der Kreismittelpunkt der Schwerpunkt des Datensatzes ist und der Radius der Mittelwert aller Abstände zum eben diesem Schwerpunkt sein müsste. Dass dieser Ansatz äußerst ungünstig sein kann wurde schon im Artikel zum ausgleichenden Kreis gezeigt: ( Näherungswertbestimmung für ausgleichende Kreise )
Eine stabilere Näherungslösung erhält man indem man die man Gleichung (10) wie folgt umformt
(14)
(15)
(16)
Betrachtet man im Folgenden als Beobachtungen,
,
,
als Konstanten und
,
,
,
als Unbekannte der Ausgleichungsaufgabe so können die neuen Parameter geschlossen berechnet werden.
Überführt in das GHM entsteht die nun folgende Bedingungsgleichung
(17)
Die Jacobi-Matrizen lauten
(18)
(19)
Der neue Beobachtungsvektor lautet
(20)
Aus und
entsteht die nun folgende Blockmatrizenform.
(21)
sowie dessen Lösung
(22)
Durch die Substitution in (15) gelingt es deutlich stabilere Schätzwerte für die exakten Parameter
zu bekommen. Aus der geschätzten Hilfvariable
erhält man mittels Rücksubstitution den geschätzen Radius
(23)
Quellcode
Download hier: best_fit_sphere
%% Berechnung ausgleichender Kugelparameter im Gauss-Helmert-Modell.
%
% v1.01
%
% Kisser Waldemar (2011)
% http://kisser.online/ausgleichung/ghm/kugel
clc
clear all
format short
%% Laden des Datensatzes
% Punkte = dlmread( 'c:\daten\kugel\jaeger.xyz' );
% oder
Punkte = ...
[ ...
14.551 9.495 8.433
9.071 5.002 10.892
12.865 4.306 8.775
9.170 7.898 11.930
9.993 9.575 11.745
12.152 6.017 11.054
];
%% Hilfsvariablen
% Anzahl Punkte, Dimension
[ n , d ] = size( Punkte );
%% Schaetzung der Naeherungswerte
A = [ Punkte , ones(n,1) * 0.5 ];
% Loesung
X_dach = ( A.' * A ) \ A.' * 0.5 * sum( Punkte.^2 , 2 );
% Ruecksubstitution von a in r0
X_dach(4) = sqrt( X_dach(4) + sum( X_dach( 1 : 3 ).^2 ) );
% Anzahl Unbekannte
u = length( X_dach );
%% Schaetzung der exakten Minimierer X,V
% Maximale Iterationsanzahl
maxit = 20;
% Betrag der Differenzloesung (Startwert)
max_dx = Inf;
% Aktuelle Iteration (Startwert)
iteration = 0;
% Konvergenzgrenze (kann geaendert werden)
EPSILON = 1e-12;
% Standardabweichung der Gewichtseinheit (kann geaendert werden)
s0_prio = 1.0;
% Verbesserungen
v = zeros( n * d , 1 );
% Modellmatrix A
A = [ zeros( n , u - 1 ) , - ones( n , 1 ) ] ;
% Kofaktormatrix
Cl = eye( d * n );
% Kovarianzmatrix
Ql = 1 / s0_prio^2 * Cl;
% Ausgleichungskern
for iteration = 1 : maxit
% Hilfsgroessen
dx = Punkte( : , 1 ) - X_dach( 1 );
dy = Punkte( : , 2 ) - X_dach( 2 );
dz = Punkte( : , 3 ) - X_dach( 3 );
dvx = dx + v( 1 : n );
dvy = dy + v( n + 1 : 2 * n );
dvz = dz + v( 2 * n + 1 : 3 * n );
% Radius an der Stelle X0,V0
r0 = hypot( hypot( dvx , dvy ) , dvz );
% Modellmatrix A
A(:,1:3) = -[ dvx , dvy , dvz ] ./ [r0 , r0 , r0];
% Widerspruchsvektor an der Stelle X0,V0
wm = ( dx .* dvx + dy .* dvy + dz .* dvz - X_dach( 4 ) .* r0 ) ./ r0;
% Modellmatrix B an der Stelle X0,V0
B = - sparse( [ diag( A( : , 1 ) ) , diag( A( : , 2 ) ) , diag( A( : , 3 ) ) ] );
% Zusammensetzung der Normalgleichungsmatrix
N = - A.' / ( B * Ql * B.' ) * A;
Ncond = rcond( N );
Ndet = abs( det( N ) );
if ( ( Ndet < EPSILON ) || ( Ncond < EPSILON ) )
fprintf( 'Fehler: Normalgleichungsmatrix singulär oder' );
fprintf( 'Schlecht konditioniert\n' );
fprintf( 'Determinante : %.3e\n' , Ndet );
fprintf( 'Konditionierung : %.3e\n\n' , Ncond );
return;
else
% Berechnung der Differenzloesung
x = N \ A.' * wm;
end
% Korrelaten des Modells
km = ( B * Ql * B.' ) \ ( wm - A * x );
% Berechnung der Verbesserungen
v = Ql * B.' * km;
% Addieren der Differenzloesung. Neue Naeherungswerte
X_dach = X_dach + x;
% Konvergenzkriterium
if max( abs( x ) ) < EPSILON
break
end
end
%% Ausgabe
fprintf( ' --- AUSGLEICHUNGSPROTOKOLL ---' );
fprintf( ' \n\n-- Ausgleichungsvorgang --' );
fprintf( ' \nModell\t\t\t\t: Gauss-Helmert-Modell' );
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;
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( '\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( [ '\nx\t\t: %12.8f\t' , char(177) , '%12.8f' ] , X_dach( 1 ) , s0_post * sqrt( Qxdach( 1 , 1 ) ) );
fprintf( [ '\ny\t\t: %12.8f\t' , char(177) , '%12.8f' ] , X_dach( 2 ) , s0_post * sqrt( Qxdach( 2 , 2 ) ) );
fprintf( [ '\nz\t\t: %12.8f\t' , char(177) , '%12.8f' ] , X_dach( 3 ) , s0_post * sqrt( Qxdach( 3 , 3 ) ) );
fprintf( [ '\nr\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
- The Institute of Electrical and Electronics Engineers Inc, 345 East 47th Street, New York, NY 10017, USA, IEEE Standard for Binary Floating-Point Arithmetic, 1985, PDF, http://754r.ucbtest.org/standards/754.pdf [↩]
- Tim Mattson, Ken Strandberg, Parallelization and Floating Point Numbers, 2008-12-17, PDF, http://software.intel.com/file/1645 [↩]