Ausgleichende Ebene im Gauß-Helmert-Modell

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 \mathbb R^3 unter anderem durch drei nicht identische Punkte \mathbf{P}_{1}, \mathbf{P}_{2}, \mathbf{P}_{3} eindeutig definiert die nicht auf nicht auf einer Gerade liegen. Bildet man nun ein Vektorprodukt

(1)   \begin{equation*} \mathbf{n} = ( \mathbf{P}_{2} - \mathbf{P}_{1} ) \times ( \mathbf{P}_{3} - \mathbf{P}_{1} ) \end{equation*}

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 \mathbf{P}_{i} auf den normierten Normalenvektor \mathbf{n}_{0}

(2)   \begin{equation*} d=\frac{\mathbf{n} \cdot \mathbf{P}_i}{|\mathbf{n}|}=\mathbf{n}_o\mathbf{P}_i \end{equation*}

In der Hesse’schen Normalenform wird die Lage der Ebene eindeutig mit dem normierten Normalenvektor \mathbf{n}_0 und dem Abstand der Ebene zum Ursprung d beschrieben.

(3)   \begin{equation*} E: \frac{ 1 }{ \sqrt{n_x^2 + n_y^2 + n_z^2 } }\begin{bmatrix} n_x \\ n_y \\ n_z \end{bmatrix} \cdot \begin{bmatrix} x_i \\ y_i \\ z_i \end{bmatrix} +d=0 \end{equation*}

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)   \begin{equation*} \psi _i(\mathbf{v,x})=n_x\left(x_i+v_{x_i}\right)+n_y\left(y_i+v_{y_i}\right)+n_z\left(z_i+v_{z_i}\right)+d=0 \end{equation*}

Die Jacobi- oder Modellmatrizen \mathbf{A} und \mathbf{B} werden an der Stelle \mathbf{x_0,v_0} ( Näherungswerte aus der aktuellen Iteration ) differenziert und erhalten die folgende Form:

(5)   \begin{equation*} \mathbf{A}(\mathbf{x,v}) = \left. \frac{ \partial \mathbf{\phi} }{ \partial \mathbf{x} } \right|_{\mathbf{x}_0,\mathbf{v}_0 } = \begin{bmatrix} \frac{\partial \psi _1}{\partial n_x^0} & \frac{\partial \psi _1}{\partial n_y^0} & \frac{\partial \psi _1}{\partial n_z^0} & \frac{\partial \psi _1}{\partial d} \\ \frac{\partial \psi _2}{\partial n_x^0} & \frac{\partial \psi _2}{\partial n_y^0} & \frac{\partial \psi _2}{\partial n_z^0} & \frac{\partial \psi _2}{\partial d} \\ \vdots & \vdots & \vdots & \vdots \\ \frac{\partial \psi _n}{\partial n_x^0} & \frac{\partial \psi _n}{\partial n_y^0} & \frac{\partial \psi _n}{\partial n_z^0} & \frac{\partial \psi _n}{\partial d} \end{bmatrix} \end{equation*}

Wir erhalten bei vier Unbekannten für \mathbf{A} 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 \mathbf{A}

(6)   \begin{equation*} \mathbf{A}(\mathbf{x,v})= \begin{bmatrix} x_1+v^0_{x_1} & y_1+v^0_{y_1} & z_1+v^0_{z_1} & 1 \\ x_2+v^0_{x_2} & y_2+v^0_{y_2} & z_2+v^0_{z_2} & 1 \\ \vdots & \vdots & \vdots & \vdots \\ x_n+v^0_{x_n} & y_n+v^0_{y_n} & z_n+v^0_{z_n} & 1 \end{bmatrix}_{n\times 4} \end{equation*}

Die Modellmatrix \mathbf{B} 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)   \begin{equation*} \mathbf{B}(\mathbf{x,v}) = \left. \frac{ \partial \mathbf{\phi} }{ \partial \mathbf{v} }\right|_{\mathbf{x}_0,\mathbf{v}_0} = \end{equation*}

    \begin{equation*} \left[ \begin{array}{cccccccccccc} \frac{\partial \psi _1}{\partial v^0_{x_1}} & 0 & \cdots & 0 & \frac{\partial \psi _1}{\partial v^0_{y_1}} & 0 & \cdots & 0 & \frac{\partial \psi _1}{\partial v^0_{z_1}} & 0 & \cdots & 0 \\ 0 & \frac{\partial \psi _2}{\partial v^0_{x_2}} & 0 & \vdots & 0 & \frac{\partial \psi _2}{\partial v^0_{y_2}} & \text{} & \vdots & 0 & \frac{\partial \psi _2}{\partial v^0_{z_2}} & \text{} & \vdots \\ \vdots & 0 & \ddots & 0 & \vdots & \text{} & \ddots & 0 & \vdots & \text{} & \ddots & 0 \\ 0 & \hdots & 0 & \frac{\partial \psi _n}{\partial v^0_{x_n}} & 0 & \cdots & 0 & \frac{\partial \psi _n}{\partial v^0_{y_n}} & 0 & \hdots & 0 & \frac{\partial \psi _n}{\partial v^0_{z_n}} \end{array} \right] \end{equation*}

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)   \begin{equation*} \mathbf{B}(\mathbf{x,v})= \left[ \begin{array}{cccccccccccc} n_x^0 & 0 & \cdots & 0 & n_y^0 & 0 & \cdots & 0 & n_z^0 & 0 & \cdots & 0\\ 0 & n_x^0 & \text{} & \vdots & 0 & n_y^0 & \text{} & \vdots & 0 & n_z^0 & \text{} & \vdots\\ \vdots & \text{} & \ddots & 0 & \vdots & \text{} & \ddots & 0 & \vdots & \text{} & \ddots & 0\\ 0 & \cdots & 0 & n_x^0 & 0 & \cdots & 0 & n_y^0 & 0 & \cdots & 0 & n_z^0 \end{array} \right]_{n \times 3n} \end{equation*}

Explizit ausgeschrieben, lautet der Widerspruchsvektor, \mathbf{w}_m

(9)   \begin{equation*} \mathbf{w}_m= -\mathbf{Bv}_0 + \phi \left( \mathbf{x},\mathbf{v} \right) = \end{equation*}

    \begin{equation*} \left[ \begin{array}{c} -n_x^0v^0_{x_1}-n_y^0v^0_{y_1}-n_z^0v^0_{z_1}+n_x^0\left(x_1+v^0_{x_1}\right) + n_y^0\left(y_1+v^0_{y_1}\right)+n_z^0\left(z_1+v^0_{z_1}\right)+d^0 \\ -n_x^0v^0_{x_2}-n_y^0v^0_{y_2}-n_z^0v^0_{z_2}+n_x^0\left(x_2+v^0_{x_2}\right) + n_y^0\left(y_2+v^0_{y_2}\right)+n_z^0\left(z_3+v^0_{z_3}\right)+d^0 \\ \vdots \\ -n_x^0v^0_{x_n}-n_y^0v^0_{y_n}-n_z^0v^0_{z_n}+n_x^0\left(x_n+v^0_{x_n}\right) +n_y^0\left(y_n+v^0_{y_n}\right)+n_z^0\left(z_n+v^0_{z_n}\right)+d^0 \end{array} \right] \end{equation*}

Dieser Vektor kann nach der Auflösung der Klammern algebraisch so zusammengefasst werden, dass eine reduzierte Form entsteht, die keinerlei Verbesserungen beinhaltet.

(10)   \begin{equation*} \mathbf{w}_m= \left[ \begin{array}{c} n_x^0x_1+n_y^0y_1+n_z^0z_1+d^0 \\ n_x^0x_2+n_y^0y_2+n_z^0z_2+d^0 \\ \vdots \\ n_x^0x_n+n_y^0y_n+n_z^0z_n+d^0 \end{array} \right]_{n\times 1} \end{equation*}

An dieser Stelle muss noch einmal darauf hingewiesen werden, dass es im Allgemeinen nicht üblich ist, dass die Verbesserungen v_i 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)   \begin{equation*} \gamma (\mathbf{x})=\sqrt{\left(n_x^0\right){}^2+\left(n_y^0\right){}^2+\left(n_z^0\right){}^2}=1 \end{equation*}

Dieser wird an der Stelle \mathbf{x}_0 differenziert, weil er eine Zwangsbedingung zwischen den Parametern beschreibt. Somit lautet unser Constraint-Vektor \mathbf{C}

(12)   \begin{equation*} \mathbf{C}(\mathbf{x})=\left.\frac{\partial \gamma }{\partial \mathbf{x}}\right|_{\mathbf{x}_0} = \frac{1}{\gamma (\mathbf{x})}\left[ \begin{array}{ccc} n_x^0 & n_y^0 & n_z^0 \end{array} \right]_{1\times 3} \end{equation*}

Der dazugehörige skalare Widerspruch w_c beschreibt, wie groß die Abweichung von der Länge 1 ist

(13)   \begin{equation*} 1-w_c = \gamma(\mathbf{x}) \end{equation*}

 

Bildung der Blockmatrizenform

Für die Ausgleichung wird die folgende Blockmatrizenform gebildet

(14)   \begin{equation*} \begin{bmatrix} \underset{4 \times 4}{-\mathbf{A}^T\left(\mathbf{B}\mathbf{Q}\mathbf{B}^T\right)^{-1}\mathbf{A}} & \underset{3 \times 1}{\mathbf{C}^T} \\ \underset{1 \times 3}{\mathbf{C}} & \underset{1 \times 1}{0} \end{bmatrix}_{5 \times 5} \begin{bmatrix} \underset{4 \times 1}{\hat{\mathbf{x}}-\mathbf{x}^0} \\ \underset{1 \times 1}{k_c} \end{bmatrix}_{5 \times 1} = \begin{bmatrix} \underset{4 \times 1}{\mathbf{A}^T\left(\mathbf{B}\mathbf{Q}\mathbf{B}^T\right)^{-1}\mathbf{w}_m} \\ \underset{1 \times 1}{w_c} \end{bmatrix}_{5 \times 1} \end{equation*}

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)
% https://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

  1. Benannt nach dem Mathematiker Otto Hesse (1811 – 1874) [↩]
  2. 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 [↩]