Vorwort
Im Rahmen dieses Artikels soll die Anwendung des Gauß-Helmert-Modells anhand der Bestimmung ausgleichender Kreise 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 des Kreises im Gauß-Helmert-Modell
Ein Kreis kann unter anderem auch durch seinen Mittelpunkt und Radius eindeutig beschrieben werden. Jeder Punkt der auf dem Kreis 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 Kreisgleichung (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 zwei 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 (1), 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 gleichzeitig 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 zeigt dieser Datensatz
![]() |
![]() |
![]() |
|
---|---|---|---|
![]() |
1,49 | 2,29 | 0,416 |
![]() |
1,69 | 2,12 | 0,154 |
![]() |
1,99 | 1,99 | 0,180 |
![]() |
2,09 | 1,72 | 0,414 |
![]() |
![]() |
![]() |
|
1,815 | 2,030 | 0,291 |
welcher in der nachfolgenden Abbildung visualisiert wird
Es wird offensichtlich dass die Abschätzung des Kreismittelpunktes durch die Mittelwertmethode insbesondere bei Datensätzen die nur einen kleinen Kreissektor abbilden mangelhafte Ergebnisse liefert und weil der defizient approximierte Mittelpunkt auch in die Radius-Schätzung einfließt, kann dieser ebenfalls nur mit einer unzureichenden Genauigkeit ermittelt werden.
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_circle.m
%% Berechnung ausgleichender Kreisparameter im Gauss-Helmert-Modell.
%
% v1.04
%
% Kisser Waldemar (2011)
% http://kisser.online/ausgleichung/ghm/kreis
clc
clear all
format long G
%% Laden des Datensatzes
% Punkte = dlmread('c:\date\boppkrausskreis.xy');
% oder
Punkte = ...
[ ...
164.595 73.414
159.396 62.563
136.455 45.842
112.648 46.136
85.822 71.995
84.701 95.772
89.246 106.901
113.501 125.639
137.304 125.374
164.847 97.226
];
%% Hilfsvariablen
% Anzahl Punkte
[ n , d ] = size( Punkte );
%% Schaetzung der Naeherungswerte
A = [ Punkte , ones( n , 1 ) * 0.5 ];
% Normalgleichungsmatrix
N = A.' * A;
% Loesung
X_dach = N \ A.' * 0.5 * sum( Punkte.^2 , 2 );
% Ruecksubstitution von a in r0
X_dach(3) = sqrt( X_dach( 3 ) + sum( X_dach( 1 : 2 ).^2 ) );
% Anzahl Unbekannte
u = length( X_dach );
%% Schaetzung der exakten Minimierer X,V
% Vektor der Verbesserungen
v = zeros( d * n , 1 );
% Maximale Iterationsanzahl
maxit = 20;
% Betrag der Differenzloesung (Startwert)
max_dx = Inf;
% Konvergenzgrenze (kann geaendert werden)
EPSILON = 1e-12;
% Standardabweichung der Gewichtseinheit (kann geaendert werden)
s0_prio = 1.0;
% 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 );
dvx = dx + v( 1 : n );
dvy = dy + v( n + 1 : 2 * n );
% Radius an der Stelle X0,V0
r0 = hypot( dvx , dvy );
% Modellmatrix A
A(:,1:2) = -[ dvx , dvy ] ./ [ r0 , r0 ];
% Widerspruchsvektor des Modells an der Stelle X0,V0
wm = ( dx .* dvx + dy .* dvy - X_dach( 3 ) .* r0 ) ./ r0;
% Modellmatrix B an der Stelle X0,V0
B = - sparse( [ diag( A( : , 1 ) ) , diag( A( : , 2 ) ) ] );
% 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.' / ( B * Ql * B.' ) * wm;
end
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' , size( X_dach , 1 ) );
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( [ '\nr\t\t: %12.8f\t' , char(177) , '%12.8f\n' ] , X_dach( 3 ) , s0_post * sqrt( Qxdach( 3 , 3 ) ) );
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 [↩]