Ausgleichender Kreis im Gauß-Helmert-Modell

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

Kreis, beschrieben durch seinen Mittelpunkt und Radius

Kreis, beschrieben durch seinen Mittelpunkt und Radius

 

Ein Kreis kann unter anderem auch durch seinen Mittelpunkt und Radius eindeutig beschrieben werden. Jeder Punkt P_i(X,Y) der auf dem Kreis liegt, erfüllt die Gleichung

(1)   \begin{equation*} \sqrt{\left(x_i-x_m\right)^2+\left(y_i-y_m\right)^2}=r \end{equation*}

womit dieser Term den funktionalen Zusammenhang zwischen den Unbekannten \mathbf{x}

(2)   \begin{equation*} \mathbf{x}= \begin{bmatrix} x_m & y_m & r \end{bmatrix}^T \end{equation*}

und den Beobachtungen \mathbf{l}

(3)   \begin{equation*} \mathbf{l}= \begin{bmatrix} x_1 & y_1 & x_2 & y_2 & \ldots & x_n & y_n \end{bmatrix}^T \end{equation*}

beschreibt. Überführt man die Kreisgleichung (1) unter Berücksichtigung der Verbesserungen \mathbf{v}

(4)   \begin{equation*} \mathbf{v}= \begin{bmatrix} v_{x_1} & v_{x_2} & \ldots & v_{x_n} & v_{y_1} & v_{y_2} & \ldots & v_{y_n} \end{bmatrix}^T \end{equation*}

in das GHM so entsteht die folgende Bedingungsfunktion

(5)   \begin{equation*} \psi _i (\mathbf{x,v}) = \sqrt{\left(x_i + v_{x_i} - x_m \right)^2+\left(y_i + v_{y_i} - y_m\right)^2}-r \end{equation*}

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 \mathbf{v}_0 und Parameter \mathbf{x}_0 (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)   \begin{equation*} \mathbf{A_\psi(x,v)} = \left.\frac{\partial \mathbf{\phi}}{\partial \mathbf{x}}\right|_{\mathbf{x_0,v_0}} = \end{equation*}

    \begin{equation*} \begin{bmatrix} \frac{\partial \psi _1}{\partial x_m^0} & \frac{\partial \psi _1}{\partial y_m^0} & \frac{\partial \psi _1}{\partial r^0} \\ \frac{\partial \psi _2}{\partial x_m^0} & \frac{\partial \psi _2}{\partial y_m^0} & \frac{\partial \psi _2}{\partial r^0} \\ \vdots & \vdots & \vdots \\ \frac{\partial \psi _n}{\partial x_m^0} & \frac{\partial \psi _n}{\partial y_m^0} & \frac{\partial \psi _n}{\partial r^0} \end{bmatrix} = - \begin{bmatrix} \frac{x_1+v_{x_1}-x_m^0}{r^0_1} & \frac{y_1+v_{y_1}-y_m^0}{r^0_1} & 1 \\ \frac{x_2+v_{x_2}-x_m^0}{r^0_2} & \frac{y_2+v_{y_2}-y_m^0}{r^0_2} & 1 \\ \vdots & \vdots & \vdots \\ \frac{x_n+v_{x_n}-x_m^0}{r^0_n} & \frac{y_n+v_{y_n}-x_m^0}{r^0_n} & 1 \end{bmatrix} \end{equation*}

mit der Hilfsvariable

(7)   \begin{equation*} r_i^0 = \sqrt{\left(x_i+v_{x_i}-x_m^0\right)^2+\left(y_i+v_{y_i}-y_m^0\right)^2} \end{equation*}

Weil die Lösung iterativ bestimmt wird, können Taylor Glieder höherer Ordnung vernachlässigt werden. Die Modelllmatrix \mathbf{B} wird als eine Bandmatrix mit zwei an einander gehängten nicht überlappenden Diagonalen gebildet.

(8)   \begin{equation*} \mathbf{B_\psi(x,v)} = \left.\frac{\partial \mathbf{\phi}}{\partial \mathbf{v}}\right|_{\mathbf{x_0,v_0}} = \end{equation*}

    \begin{equation*} \left[ \begin{matrix} \frac{\partial \psi _1(\mathbf{x,v})}{\partial v_{x_1}} & 0 & \cdots & 0 & \frac{\partial \psi _1(\mathbf{x,v})}{\partial v_{y_1}} & 0 & \cdots & 0 \\ 0 & \frac{\partial \psi _2(\mathbf{x,v})}{\partial v_{x_2}} & \text{} & \vdots & 0 & \frac{\partial \psi _2(\mathbf{x,v})}{\partial v_{y_2}} & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 & \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & \frac{\partial \psi _n(\mathbf{v,x})}{\partial v_{x_n}} & 0 & \cdots & 0 & \frac{\partial \psi _n(\mathbf{v_x})}{\partial v_{y_n}} \end{matrix} \right] \end{equation*}

    \begin{equation*} = \left[ \begin{matrix} \frac{x_1+v_{x_1}-x_m^0}{r_1^0} & 0 & \cdots & 0 \\ 0 & \frac{x_2+v_{x_2}-x_m^0}{r_2^0} & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & \frac{x_n+v_{x_n}-x_m^0}{r_n^0} \end{matrix} \right. \end{equation*}

    \begin{equation*} \left. \begin{matrix} \frac{y_1+v_{y_1}-y_m^0}{r_1^0} & 0 & \cdots & 0 \\ 0 & \frac{y_2+v_{y2}-y_m^0}{r_2^0} & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & \frac{y_n+v_{y_n}-y_m^0}{r_n^0} \end{matrix} \right] \end{equation*}

Da man eine Linearisierung vornimmt muss auch ein Widerspruchsvektor eingeführt werden.

(9)   \begin{equation*} \mathbf{w_\psi} = -\mathbf{Bv_0} + \psi \left( \mathbf{x_0},\mathbf{v_0} \right) = \begin{bmatrix} \frac{(x_1+v_{x_1})(x_1-x_m) +(y_1+v_{y_1})(y_1-y_m) -r^0\cdot r^0_1}{r^0_1} \\ \frac{(x_2+v_{x_2})(x_2-x_m) +(y_2+v_{y_2})(y_2-y_m) -r^0\cdot r^0_2}{r^0_2} \\ \vdots \\ \frac{(x_n+v_{x_n})(x_n-x_m) +(y_n+v_{y_n})(y_n-y_m) -r^0\cdot r^0_n}{r^0_n} \\ \end{bmatrix} \end{equation*}

Alternative Bedingungsfunktion im GHM

Alternativ kann auch eine eine andere Form der Kreisgleichung verwendet werden

(10)   \begin{equation*} \phi _i = \left(x_i-x_m\right)^2+\left(y_i-y_m\right)^2-r^2 \end{equation*}

Diese ist zwar im streng mathematischen Sinne äquivalent zu (1), nummerisch betrachtet haben die Matrizen der partiellen Ableitungen

(11)   \begin{equation*} \mathbf{A_\phi(x,v)} = \\ -2\begin{bmatrix} x_1+v_{x_1}^0-x_m^0 & y_1+v_{y_1}^0-y_m^0 & r_0 \\ x_2+v_{x_2}^0-x_m^0 & y_2+v_{y_2}^0-y_m^0 & r_0 \\ \vdots & \vdots & \vdots \\ x_n+v_{x_n}^0-x_m^0 & y_n+v_{y_n}^0-y_m^0 & r_0 \end{bmatrix} \end{equation*}

(12)   \begin{equation*} \mathbf{B_\phi(x,v)} = \\ \end{equation*}

    \begin{equation*} 2\left[ \begin{matrix} x_1+v_{x_1}^0-x_m^0 & 0 & \cdots & 0 \\ 0 & x_2+v_{x_2}^0-x_m^0 & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & x_n+v_{x_n}^0-x_m^0 \end{matrix} \right. \end{equation*}

    \begin{equation*} \left. \begin{matrix} y_1+v_{y_1}^0-y_m^0 & 0 & \cdots & 0 \\ 0 & y_1+v_{y_2}^0-y_m^0 & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & y_n+v_{y_n}^0-y_m^0 \end{matrix} \right] \end{equation*}

und der Widerspruchspruchvektor

(13)   \begin{equation*} \mathbf{w_\phi} = \begin{bmatrix} \left(x_1-x_m^0\right)^2-\left(v_{x_1}^0\right)^2+\left(y_1-y_m^0\right)^2-\left(v_{y_1}^0\right)^2-r_0^2 \\ \left(x_2-x_m^0\right)^2-\left(v_{x_2}^0\right)^2+\left(y_2-y_m^0\right)^2-\left(v_{y_2}^0\right)^2-r_0^2 \\ \vdots \\ \left(x_n-x_m^0\right)^2-\left(v_{x_n}^0\right)^2+\left(y_n-y_m^0\right)^2-\left(v_{y_n}^0\right)^2-r_0^2 \end{bmatrix} \end{equation*}

deutlich größere Beträge, welche je nach verwendetem Datentyp zu einer geringeren Rechenschärfe führen können, da nicht alle Reele Zahlen \mathbb{R} 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

Tabelle 1. Kreisdatensatz
X_i Y_i R_i
P_1 1,49 2,29 0,416
P_2 1,69 2,12 0,154
P_2 1,99 1,99 0,180
P_2 2,09 1,72 0,414
S_x S_y \bar{r}
1,815 2,030 0,291

 

welcher in der nachfolgenden Abbildung visualisiert wird

Visualisierung Tabelle 1

Visualisierung Tabelle 1

 

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)   \begin{equation*} -2x_i x_m - 2y_i y_m + y_i^2 + x_i^2 + x_m^2 + y_m^2=r^2 \end{equation*}

(15)   \begin{equation*} \underbrace{y_i^2 + x_i^2}_{\mathbf{l}} = 2x_i x_m + 2y_i y_m + \underbrace{r^2 - x_m^2 - y_m^2}_{a} \end{equation*}

(16)   \begin{equation*} l_i = 2x_i x_m + 2y_i y_m + a \end{equation*}

Betrachtet man im Folgenden l_i als Beobachtungen, x_i, y_i als Konstanten und x_m,y_m,a als Unbekannte der Ausgleichungsaufgabe so können die neuen Parameter geschlossen berechnet werden.

Überführt in das GHM entsteht die nun folgende Bedingungsgleichung

(17)   \begin{equation*} \tau _i (\mathbf{x,v}) = 2x_i x_m + 2y_i y_m + a - l_i - v_i \left. \right| : 2 \end{equation*}

    \begin{equation*} \tau _i (\mathbf{x,v}) = x_i x_m + y_i y_m + \frac{1}{2}\left(a - l_i - v_i\right) \end{equation*}

Die Jacobi-Matrizen lauten

(18)   \begin{equation*} \mathbf{A_{\tau}(x,v)}= \begin{bmatrix} x_1 & y_1 & \frac{1}{2} \\ x_2 & y_2 & \frac{1}{2} \\ \vdots & \vdots & \vdots \\ x_n & y_n & \frac{1}{2} \end{bmatrix} \end{equation*}

(19)   \begin{equation*} \mathbf{B_{\tau}(x,v)}= \frac{1}{2} \left[ \begin{array}{rrrr} -1 & 0 & \cdots & 0 \\ 0 & -1 & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & -1 \end{array} \right] \end{equation*}

Der neue Beobachtungsvektor lautet

(20)   \begin{equation*} \mathbf{L} = \frac{1}{2} \begin{bmatrix} x_1^2+y_1^2 & x_2^2+y_2^2 & \hdots & x_n^2+y_n^2 \end{bmatrix}^T \end{equation*}

Aus \mathbf{A_{\tau}, B_{\tau}} und \mathbf{L} entsteht die nun folgende Blockmatrizenform.

(21)   \begin{equation*} \begin{bmatrix} \mathbf{BQB}^T & \mathbf{A} \\ \mathbf{A}^T & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{k} \\ \mathbf{\tilde X } \end{bmatrix} -\begin{bmatrix} \mathbf{L} \\ \mathbf{0} \end{bmatrix} = \mathbf{0} \end{equation*}

sowie dessen Lösung

(22)   \begin{equation*} \begin{bmatrix} \mathbf{BQB}^T & \mathbf{A} \\ \mathbf{A}^T & \mathbf{0} \end{bmatrix}^{-1} \begin{bmatrix} \mathbf{L} \\ \mathbf{0} \end{bmatrix} = \begin{bmatrix} \mathbf{k} \\ \mathbf{\tilde X } \end{bmatrix} \end{equation*}

Durch die Substitution in (15) gelingt es deutlich stabilere Schätzwerte \tilde X für die exakten Parameter \hat x zu bekommen. Aus der geschätzten Hilfvariable a erhält man mittels Rücksubstitution den geschätzen Radius

(23)   \begin{equation*} r^0 = \sqrt{a + x_0^2+y_0^2} \end{equation*}

 

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

  1. 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 [↩]
  2. Tim Mattson, Ken Strandberg, Parallelization and Floating Point Numbers, 2008-12-17, PDF, http://software.intel.com/file/1645 [↩]