Ausgleichende Kugel im Gauß-Helmert-Modell

Von
wkiss
– 2011-10-16Eingestellt unter: Gauß-Helmert-Modell

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

Kreis, beschrieben durch seinen Mittelpunkt und Radius

Kreis, beschrieben durch seinen Mittelpunkt und Radius

 

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

(1)   \begin{equation*} \sqrt{\left(x_i-x_m\right)^2+\left(y_i-y_m\right)^2 +\left(z_i-z_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 & z_m & r \end{bmatrix}^T \end{equation*}

und den Beobachtungen \mathbf{l}

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

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

(4)   \begin{equation*} \mathbf{v}= \left[ \begin{array}{cccccccccccc} v_{x_1} & v_{x_2} & ... & v_{x_n} & v_{y_1} &v_{y_2} & ... & v_{y_n} & v_{z_1} & v_{z_2} & ... & v_{z_n} \end{array} \right]^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 +\left(z_i +v_{z_i} -z_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 z_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 z_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 z_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} & \frac{z_1+v_{z_1}-z_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} & \frac{z_2+v_{z_2}-z_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} & \frac{z_n+v_{z_n}-z_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+\left(z_i+v_{z_i}-z_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 drei 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{\partial \psi _1(\mathbf{x,v})}{\partial v_{z_1}} & 0 & \cdots & 0 \\ 0 & \frac{\partial \psi _2(\mathbf{x,v})}{\partial v_{z_2}} & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & \frac{\partial \psi _n(\mathbf{v,x})}{\partial v_{z_n}} \end{matrix} \right] = \end{equation*}

    \begin{equation*} \left[ \begin{matrix} \frac{x_1+v_{x_1}-y_m^0}{r_1^0} & 0 & \cdots & 0 \\ 0 & \frac{x_2+v_{x_2}-y_m^0}{r_2^0} & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & \frac{x_n+v_{x_n}-y_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_{y_2}-y_m^0}{r_2^0} & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & \frac{y_n+v_{zy_n}-y_m^0}{r_n^0} \end{matrix} \right. \end{equation*}

    \begin{equation*} \left. \begin{matrix} \frac{z_1+v_{z_1}-y_m^0}{r_1^0} & 0 & \cdots & 0 \\ 0 & \frac{z_2+v_{z_2}-y_m^0}{r_2^0} & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & \frac{z_n+v_{z_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) = \end{equation*}

    \begin{equation*} \begin{bmatrix} \frac{(x_1+v_{x_1})(x_1-x_m) + (y_1+v_{y_1})(y_1-y_m) + (z_1+v_{z_1})(z_1-z_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) + (z_2+v_{z_2})(z_2-z_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) + (z_n+v_{z_n})(z_n-z_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 + v_{x_i} - x_m \right)^2 + \left( y_i + v_{y_i} - y_m\right)^2+\left( z_i -z_m\right)^2-r^2 \end{equation*}

Diese ist zwar im streng mathematischen Sinne äquivalent zu (5), 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 & z_1+v_{z_1}^0-z_m^0 & r_0 \\ x_2+v_{x_2}^0-x_m^0 & y_2+v_{y_2}^0-y_m^0 & z_2+v_{z_2}^0-z_m^0 & r_0 \\ \vdots & \vdots & \vdots & \vdots \\ x_n+v_{x_n}^0-x_m^0 & y_n+v_{y_n}^0-y_m^0 & z_n+v_{z_n}^0-z_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_{x1}^0-x_m^0 & 0 & \cdots & 0 \\ 0 & x_2+v_{x2}^0-x_m^0 & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & x_n+v_{xn}^0-x_m^0 \end{matrix} \right. \end{equation*}

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

    \begin{equation*} \left. \begin{matrix} z_1+v_{z_1}^0-z_m^0 & 0 & \cdots & 0 \\ 0 & z_1+v_{z_2}^0-z_m^0 & \text{} & \vdots \\ \vdots & \text{} & \ddots & 0 \\ 0 & \cdots & 0 & z_n+v_{z_n}^0-z_m^0 \end{matrix} \right] \end{equation*}

und der Widerspruchspruchvektor

(13)   \begin{equation*} \mathbf{w_\phi} = \end{equation*}

    \begin{equation*} \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_{y1}^0\right)^2 +\left(y_z-z_m^0\right)^2-\left(v_{z_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 +\left(z_2-y_m^0\right)^2-\left(v_{z_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 +\left(z_n-z_m^0\right)^2-\left(v_{z_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 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)   \begin{equation*} -2x_i x_m - 2y_i y_m - 2z_i z_m + x_i^2 + y_i^2 + z_i^2 + x_m^2 + y_m^2 + z_m^2=r^2 \end{equation*}

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

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

Betrachtet man im Folgenden l_i als Beobachtungen, x_i, y_i, z_i als Konstanten und x_m,y_m,z_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 + 2z_i z_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 + z_i z_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 & z_1 & \frac{1}{2} \\ x_2 & y_2 & z_2 & \frac{1}{2} \\ \vdots & \vdots & \vdots & \vdots \\ x_n & y_n & z_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+z_1^2 & x_2^2+y_2^2 +z_3^2 & \hdots & x_n^2+y_n^2 +z_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 + z_0^2 } \end{equation*}

 

Quellcode

Download hier: mbest_fit_sphere - v1.01 - m 

%% 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

  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 [↩]
Schlagwörter: , , , , , ,