... und mein Script funktioniert mit Deinen Werten oder nicht?

X=t0+mRx
und ja dieser Ansatz beruht auf eine Affintransformation

In dieser Schreibweise aber ggf. missverständlich, da es eigentlich Rmx lauten müsste. Bei der 7-Parameter-Transformation ist dies aber egal. Schätzt Du hingegen drei unterschiedliche Maßstäbe, ist die Reihenfolge zu beachten, da i.A. MR != RM - siehe die Literaturstellen im letzten Posting. Ferner kannst Du über Deinen Ansatz nicht so einfach einzelne Scherparameter fixieren.

Und daher liegen kleine Rotationen und ein kleiner Maßstab vor.

Wozu dann den Aufwand, wenn Du Dir keinen Mehrwert versprichst?

Viele Grüße
Micha

also es handelt sich hierbei um Transformationen von GNSS-Koordinaten ins Landessystem. Und daher liegen kleine Rotationen und ein kleiner Maßstab vor.
Mit dem linearen Bursa-Wolf-Modell habe ich schon gute Ergebnisse erzielt und die lineare Affintransformation erzielte auch gute Ergebnisse.

Bei diesem Ansatz:

X=t0+mRx

mit...
t0=(X0,Y0,Z0)^T
R=[r11,r12,r13;r21,r22,r23;r31,r32,r33]

besitze ich 13 Unbekannte (3 Translationen, 9 Kooefizienten der Rotationsmatrix und 1 Maßstab)

und ja dieser Ansatz beruht auf eine Affintransformation und mit hinzufügen der 6 Bedingungsgleichungen (3 Orthogonalitätsbedingungen, 3 Normierungsbedingungen) für die Rotationsmatrix, vermindert sich die Anzahl der Unbekannten auf 7-Parameter. Denn Wir wissen ja das die Rotationsmatrix orthonormal ist.

Deine Umsetzung ähnelt sehr der von Lehmann (2014): Transformation model selection by multiple hypotheses testing., welcher auf Andrei (2006): 3D affine coordinate transformations. verweist. Das Paper von Lehmann habe ich bzgl. der Wahl der Restriktionen kurz überflogen und bei meiner Implementierung herangezogen. Dennoch konnte ich auf die schnell keine zufriedenstellende Lösung finden, die auch bei "großen" Koordinaten numerisch stabil bleibt. Im Anhang findest Du meinen Versuch aus Octave/Matlab. Vielleicht hilft es aber weiter in dieser Diskussion.

Dennoch zwei Hinweise, die Dir ggf. bei Deinem Problem weiterhelfen können, wenn Du auch für andere Lösungswege offen bist. Zum Bestimmen der 7 Parameter der Helmert-Transformation existiert eine geschlossene Lösung, die keine Iterationen benötigt. Siehe hierzu auch den Artikel: Bestimmung der 7 Parameter einer Helmerttransformation mit Quaternionen.

Solltest Du vorhaben, Dein Modell auf mehr als 7 Parameter zu erweitern bspw. Affintransformation, beachte bitte, dass in der Geodäsie üblicherweise folgendes Modell verwendet wird:

X = T + RSMx

worin T die Translation, R die Drehung, S die Scherung und M die Matrix der Maßstäbe ist. Um dieses Modell zu reduzieren, könnte bspw. S = I gesetzt werden, sodass keine Scherungen auftreten. Die Maßstäbe in diesem Modell sind dem Startsystem zugeordnet und stehen somit nicht vor der Drehung, vgl. auch Koch & Fröhlich (2000): Transformation räumlicher variabler Koordinaten. Für dieses allg. Transformationsmodell und die nötigen Bedingungen zur Reduktion der Parameter findest Du entsprechende Gleichungen in Lösler & Eschelbach (2014): Zur Bestimmung der Parameter einer räumlichen Affintransformation.

Ich hoffe, dass hilft Dir weiter.

Viele Grüße
Micha

clear all
close all
format long g
clc
 
 
%% Punkte im Startsystem
PunkteS = [
    304.507 2262.732 8226.374
    3267.980 5992.254 4529.311
    568.348 4081.212 6567.109
];
 
%% Punkte im Zielsystem
PunkteZ = [
    275.525 1668.545 7810.831
    3239.034 5398.100 4113.819
    539.378 3487.020 6151.622
];
 
schwerpunktS = mean(PunkteS);
schwerpunktZ = mean(PunkteZ);
PunkteS = PunkteS - (schwerpunktS'*ones(1,size(PunkteS,1)))';
PunkteZ = PunkteZ - (schwerpunktZ'*ones(1,size(PunkteZ,1)))';
 
nop = size(PunkteS,1);
iter = 0;
%x=(f11 f12 f13 f21 f22 f23 f31 f32 f33)
X0=[1 0 0   0 1 0   0 0 1 ]';
dx = inf;
while(iter < 200 && max(abs(dx)) > eps^(0.75))
    iter = iter + 1;
    A = zeros(size(PunkteS,1)*size(PunkteS,2), length(X0));
    w = zeros(size(A,1),1);
    R = zeros(5, length(X0));
    r = zeros(size(R,1), 1);
 
    f11 = X0(1);
    f12 = X0(2);
    f13 = X0(3);
 
    f21 = X0(4);
    f22 = X0(5);
    f23 = X0(6);
 
    f31 = X0(7);
    f32 = X0(8);
    f33 = X0(9);
 
    for i=1:nop
        xP = PunkteS(i,1);
        yP = PunkteS(i,2);
        zP = PunkteS(i,3);
 
        XP = PunkteZ(i,1);
        YP = PunkteZ(i,2);
        ZP = PunkteZ(i,3);
 
        A(i,:)       = [xP yP zP   0  0  0    0  0  0 ];
        A(i+nop, :)  = [0  0  0   xP yP zP    0  0  0 ];
        A(i+2*nop,:) = [0  0  0    0  0  0    xP yP zP];
 
        w(i, 1)      = (f11*xP+f12*yP+f13*zP) - XP;
        w(i+nop, 1)  = (f21*xP+f22*yP+f23*zP) - YP;
        w(i+2*nop,1) = (f31*xP+f32*yP+f33*zP) - ZP;
    end
    R(1,:) = [f21 f22 f23   f11 f12 f13   0 0 0 ];
    R(2,:) = [0 0 0   f31 f32 f33   f21 f22 f23 ];
    R(3,:) = [f31 f32 f33   0 0 0   f11 f12 f13 ];
 
    r(1,1) = f11*f21 + f12*f22 + f13*f23;
    r(2,1) = f21*f31 + f22*f32 + f23*f33;
    r(3,1) = f11*f31 + f12*f32 + f13*f33;
 
    R(4,:) = [2*f11 2*f12 2*f13   -2*f21 -2*f22 -2*f23   0 0 0 ];
    R(5,:) = [-f11 -f12 -f13   -f21 -f22 -f23   2*f31 2*f32 2*f33 ];
 
    r(4,1) = f11*f11 + f12*f12 + f13*f13 - (f21*f21 + f22*f22 + f23*f23);
    r(5,1) = f31*f31 + f32*f32 + f33*f33 - 0.5*(f11*f11 + f12*f12 + f13*f13  + f21*f21 + f22*f22 + f23*f23);
 
    N = [A'*A  R'
        R     zeros(length(r))];
 
    n = -[A'*w; r];
 
    dxk = N\n;
    dx = dxk(1:length(X0));
    X0 = X0 + dx;
end
 
F     = [f11 f12 f13; f21 f22 f23; f31 f32 f33];
scale = mean([norm(F(1,:)) norm(F(2,:)) norm(F(3,:))]);
R     = 1/scale.*F;
 
% Pruefung auf orth. Matrix R'*R == I
disp(R'*R);
 
disp('Lösung')
disp(R);
disp(scale);
 

ich möchte gern eine räumliche Helmerttransformation mit Matlab programmieren, aber ständig kommt bei mir die Fehlermeldung das die N-Matrix singulär ist, ich habe sie dann mit einer Restriktionsmatrix G gerändert aber trotzdem ist diese geränderte Matrix singulär. Vlt habe ich ein Fehler in funktionalen Modell.

Helmert-Transformation mit analytischer Rotationsmatrix.

mein Ansatz lautet: X_z = t0 + mRX_a (oder einfacher: X = t0 + mRx)

mit t0…Translationsvektor
m…Maßstab
R…Rotationsmatrix (analytisch mit 9-Parameter: r11,r12,r13,…,r33)
X_z…Zielsystem
X_a…Ausgangssystem

ich führe dann Schwerpunktreduzierte Koordinaten ein, damit ist mein neuer Ansatz.

X_red,z = mRX_red,a (oder einfacher: X = mRx)

Meine Fehlergleichung der Form (l+v=f(x)…Gauß-Markov-Modell) lautet:

X + v = m(r11*x+r12*y+r13*z)
Y + v = m(r21*x+r22*y+r23*z)
Z + v = m(r31*x+r32*y+r33*z)

für den linearisierten Fall gilt dann die Teildesignmatrix A_1 (erster identischer Punkt) mit

A_1 =

|m*x m*y m*z 0 0 0 0 0 0 r11*x+r12*y+r13*z |
| 0 0 0 m*x m*y m*z 0 0 0 r21*x+r22*y+r23*z |
| 0 0 0 0 0 0 m*x m*y m*z r31*x+r32*y+r33*z |

für jeden identischen Punkt gilt eine Teildesignmatrix und diese werden in der Designmatrix A festgehalten (mind 3 Punkte)

A = [A_1 A_2 A_3 ... A_n]^T

Vektor der Unbekannten: x=(r11 r12 r13 r21 r22 r23 r31 r32 r33 m)^T

Vektor der Näherungswerte x_0= (1 0 0 0 1 0 0 0 1 1)^T

Da die Rotationsmatrix Orthonormal ist, gelten 3 Normierungsbedingungen und 3 Orthogonalitätsbedingungen

g1: r11^2 + r12^2 + r13^2 - 1 = 0
g2: r21^2 + r22^2 + r23^2 - 1 = 0
g3: r31^2 + r32^2 + r33^2 - 1 = 0
g4: r11*r21 + r12*r22 + r13*r23 = 0
g5: r21*r31 + r22*r32 + r23*r33 = 0
g6: r11*r31 + r12*r32 + r13*r33 = 0

für den linearisierten Fall stehen diese Bedingungen in der G-Matrix

G =

| 2*r11 2*r12 2*r13 0 0 0 0 0 0 0 |
| 0 0 0 2*r21 2*r22 2*r23 0 0 0 0 |
| 0 0 0 0 0 0 2*31 2*r32 2*r33 0 |
| r21 r22 r23 r11 r12 r13 0 0 0 0 |
| 0 0 0 r31 r32 r33 r21 r22 r23 0 |
| r31 r32 r33 0 0 0 r11 r12 r13 0 |

Wenn ich nun die Designmatrix mit der G-Matrix ränder…und nach GMMB (Gauß-Markow-Modell mit Bedingungen) löse, ist die inverse der Rändern nicht möglich. Habe ich noch irgendwo einen kleinen Denkfehler?

Danke
LG Anni