Ausgleichung mit Restriktionen (Geodäsie/Vermessung)

DoreenH @, Friday, 26.01.2024, 19:22 (vor 92 Tagen)

Hallo,
kann mir jemand die Ausgleichung mit Restriktionen erklären? Den Allgemeinfall der Ausgleichung -Gauß-Helmert-Modell- sowie die vermittelnde Ausgleichung sind mir bekannt. Aber ich konnte bislang keine Literatur finden, die den Berechnungsweg der Ausgleichung mit Restriktionen (verständlich) erläutert. Das Thema wird immer schnell abgehandelt.
Die zu erstellende B-Matrix wird nicht wie im Gauß-Helmert-Modell aus den Ableitungen der Beobachtungen erstellt, sondern aus der Ableitung der Bedingung erstellt und ergibt dann eine einspaltige B-Matrix --> so mein Stand. Aber wie sieht dann die Normalengleichung aus?

Mein Beispiel: Höhennetz mit 2 Festpunkten (A+B) und 2 Neupunkten (1+2)
geg. HA=115,584 m; HB=118,460 m (fehlerfrei); ΔhA1 = 0,5619 m; ΔhA2 = 3,5699 m; ΔhB1 = -2,3152 m; ΔhB2 = 0,6970 m; Δh12 = 3,0121 m (fehlerfrei)
Die ersten vier Höhendifferenzen haben eine a priori Standardabweichung von 1 mm, die letzte
Höhendifferenz ist fehlerfrei.

Vielleicht kann mir auch jemand einen Literatur-Tipp senden, in den man auch mal Rechenbeispiele findet.

ich danke euch.
Doreen

Avatar

Ausgleichung mit Restriktionen

MichaeL ⌂, Bad Vilbel, Friday, 26.01.2024, 21:09 (vor 92 Tagen) @ DoreenH

Hallo Doreen,

wie sieht dann die Normalengleichung aus?

Die Normalgleichung $\mathbf{Nx=n}$ ändert sich nicht. Sie wird nur von den Bedingungen $\mathbf{Rx=r}$ gerändert. Hierzu nutzt man die Lagrange Funktion, die Du bereits beim Gauß-Helmert-Modell kennengelernt hast. Sei $\mathbf{k}$ der Vektor der Lagrange Multiplikatoren, dann lautet das erweiterte System nun

$\begin{bmatrix} \mathbf{N} & \mathbf{R^{\mathrm{T}}} \\ \mathbf{R} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{x}\\\mathbf{k} \end{bmatrix} = \begin{bmatrix} \mathbf{y}\\\mathbf{r} \end{bmatrix} $

und Du suchst die Lösung für $\mathbf{x}$ und $\mathbf{k}$. Ich habe mal anstelle von $\mathbf{B}$, wie Du es genannt hattest, $\mathbf{R}$ verwendet, um Verwechslungen zu vermeiden.

Mein Beispiel: Höhennetz mit 2 Festpunkten (A+B) und 2 Neupunkten (1+2)
geg. HA=115,584 m; HB=118,460 m (fehlerfrei); ΔhA1 = 0,5619 m; ΔhA2 = 3,5699 m; ΔhB1 = -2,3152 m; ΔhB2 = 0,6970 m; Δh12 = 3,0121 m (fehlerfrei)
Die ersten vier Höhendifferenzen haben eine a priori Standardabweichung von 1 mm, die letzte
Höhendifferenz ist fehlerfrei.

Dieses Beispiel kann auf verschiedene Weisen gelöst werden. Da es Dir um die Restriktionen geht, würde ich alle Deine fehlerfreien Informationen als Bedingung formulieren wollen.

Zunächst stellst Du - wie gehabt - die Beobachtungsgleichungen auf

$\Delta hA1 = H1 - HA = 0,5619\,\mathrm{m}$
$\Delta hA2 = H2 - HA = 3,5699\,\mathrm{m}$
$\Delta hB1 = H1 - HB =-2,3152\,\mathrm{m}$
$\Delta hB2 = H2 - HB = 0,6970\,\mathrm{m}$

und bildest die Designmatrix (Spalten korrespondieren mit $HA,\,HB,\,H1,\,H2$)

$\mathbf{A} = \begin{pmatrix} -1 & 0 & +1 & 0\\ -1 & 0 & 0 & +1\\ 0 & -1 & +1 & 0\\ 0 & -1 & 0 & +1 \end{pmatrix}$

und den Vektor der Beobachtungen

$\mathbf{y} = \begin{pmatrix} 0,5619 \\ 3,5699 \\ -2,3152 \\ 0,6970 \end{pmatrix}$

Anmerkung: Häufig auch als $\mathbf{l}$-Vektor bezeichnet. Da $\mathbf{l}$, $\mathbf{I}$ und $\mathbf{1}$ schwer zu unterscheiden sind, nutze ich hier lieber $\mathbf{y}$.

Das stochastische Modell können wir vernachlässigen, da alle Beobachtungen das identische Gewicht besitzen und die Gewichtsmatrix dadurch als Einheitsmatrix dargestellt werden kann. Das Normalgleichungssystem lautet damit $\mathbf{Nx=n}$, worin $\mathbf{N=A{^\mathrm{T}}A}$ und $\mathbf{n=A{^\mathrm{T}}y}$ sind. Soweit sollte es Dir bekannt sein bzw. vorkommen - mal abgesehen von einigen Bezeichnungen vielleicht.

Mit den Restriktionen verfährst Du genauso. Wenn wir alle fehlerfreien Informationen als Restriktion betrachten wollen, dann hast Du drei Restriktionen:
$HA\overset{!}{=}115,584\,\mathrm{m}$,
$HB\overset{!}{=}118,460\,\mathrm{m}$,
$\Delta h12 = H2 - H1\overset{!}{=}3,0121\,\mathrm{m}$.

Mit

$\mathbf{R} = \begin{pmatrix} +1 & 0 & 0 & 0\\ 0 & +1 & 0 & 0\\ 0 & 0 & -1 &+1 \end{pmatrix}$

und

$\mathbf{r} = \begin{pmatrix} 115,5840 \\ 118,4600 \\ 3,0121 \end{pmatrix}$

ergeben sich diese in Matrixnotation $\mathbf{Rx=r}$.

Nun hast Du alles, was Du benötigst. Einsetzen in das erweitere Modell oben und auflösen nach $\mathbf{x}$ und $\mathbf{k}$. Die Lösung lautet dann

$\mathbf{\hat{x}} = \begin{pmatrix} 115,5840\\ 118,4600\\ 116,14435\\ 119,15645 \end{pmatrix}$

und

$\mathbf{\hat{k}} = \begin{pmatrix} 0,001\\ -0,001\\ -0,002 \end{pmatrix}$

Alles klar?

Viele Grüße
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Tags:
Ausgleichung, Höhennetz, Parameterschätzung, Gauß-Markov-Modell, Restriktionen, Bedingungen

Ausgleichung mit Restriktionen

DoreenH @, Saturday, 27.01.2024, 20:52 (vor 91 Tagen) @ MichaeL

Hallo Micha,

vielen lieben Dank für deine tolle Erläuterung. Mein großes Problem war letztendlich, dass ich mir mit der Auflösung der Normalengleichung schwer getan habe. Mir wurde die Literatur von Herrn Niemeier empfohlen und hier werden x und k separat aufgelöst. Ich habe jetzt einfach x und k in einem Vektor gelassen und die Normalengleichung berechnet.
Desweiteren muss ich die a priori und a posteriori Standardabweichungen der Neupunkte bestimmen. Ich musste jetzt feststellen, dass die N-Matrix singulär ist und diese sich nicht invertieren lässt. Meine Suche nach einer Lösung dieses Problems blieb leider erfolglos. Kannst du mir hier nochmals weiterhelfen?
Noch eine Frage: Kannst du mir einen Literaturhinweis geben, indem auch mal Beispielaufgaben enthalten sind oder darf ich dich einfach fragen?

Viele Grüße
Doreen

Avatar

Ausgleichung mit Restriktionen

MichaeL ⌂, Bad Vilbel, Sunday, 28.01.2024, 12:33 (vor 90 Tagen) @ DoreenH

Hallo Doreen,

hier werden x und k separat aufgelöst.

Diese Vorgehensweise setzt jedoch voraus, dass die Normalgleichungsmatrix $\mathbf{N}$ vollen Rang hat und damit invertiert werden kann. In vielen Anwendungen werden aber Restriktionen hinzugefügt, da $\mathbf{N}$ einen Defekt aufweist. In diesem Fall kann $\mathbf{N^{\mathrm{-1}}}$ nicht gebildet werden. Eine stufenweise Auflösung nach $\mathbf{x}$ und $\mathbf{k}$ kann demnach nicht gelingen. Eine typische Anwendung in der Geodäsie ist hier die freie Netzausgleichung. Die Normalgleichung weist keinen vollen Rang auf. Zum Lösen werden üblicherweise Bedingungsgleichungen für spezifische Punkte - Datumspunkte - eingeführt. Das erweitere System ist dann wieder regulär.

Ich habe jetzt einfach x und k in einem Vektor gelassen und die Normalengleichung berechnet.

Dies ist auch der übliche Weg, der unabhängig von einem etwaigen Defekt von $\mathbf{N}$ funktioniert (wenn das um Restriktionen erweitere System regulär ist).

Desweiteren muss ich die a priori und a posteriori Standardabweichungen der Neupunkte bestimmen. Ich musste jetzt feststellen, dass die N-Matrix singulär ist und diese sich nicht invertieren lässt.

Da die Höhenunterschiede keine Information zur Lagerung des Netzes liefern, sondern nur die relativen Informationen zwischen zwei Punkten beschreiben, können keine absoluten Höhen bestimmt werden. Die Normalgleichungsmatrix $\mathbf{N}$ weist einen Datumsdefekt auf und der Rang ist $\mathrm{rg(\mathbf{N})} = 3$ und nicht 4 (= Anzahl der Spalten in $\mathbf{N}$). Geometrisch bedeutet dies, dass wir die Höhen von drei Punkten im Netz bestimmen können, wenn wir die Höhe von einem Punkt kennen. In Deinem Fall kennen wir sogar die Höhen von zwei Punkten HA, HB. Insofern ist die Voraussetzung erfüllt und die Höhen H1, H2 lassen sich in Bezug auf HA, HB bestimmen. Da HA und HB fehlerfrei sein sollen, können für diese beiden Punkte keine Unsicherheiten aus der Ausgleichung abgeschätzt werden. Restriktionen schaffen also harte Fakten. In der Kofaktormatrix werden für HA und HB also nur Nullen stehen. Um diese Kofaktormatrix zu erhalten, invertierst Du einfach das erweiterte System und erhältst dann automatisch die Kofaktormatrix von $\mathbf{x}$, $\mathbf{k}$ und deren Abhängigkeiten.


$\begin{bmatrix} \mathbf{N} & \mathbf{R^{\mathrm{T}}} \\ \mathbf{R} & \mathbf{0} \end{bmatrix}^{-1} = \begin{pmatrix} 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0,25 & 0,25 & 0,5 & 0,5 & -0,5 \\ 0 & 0 & 0,25 & 0,25 & 0,5 & 0,5 & 0,5 \\ 1 & 0 & 0,5 & 0,5 & -1 & 1 & 0 \\ 0 & 1 & 0,5 & 0,5 & 1 & -1 & 0 \\ 0 & 0 & -0,5 & 0,5 & 0 & 0 & -1 \\ \end{pmatrix}$

So, wie Du $\mathbf{x}$ und $\mathbf{k}$ gemeinsam bestimmst, bestimmst Du auch die gemeinsame Kofaktormatrix. Und so wie Du $\mathbf{\hat{x}}$ ausgeschnitten hast, musst Du auch $\mathbf{Q_{\mathbf{\hat{x}}}}$ ausschneiden. Es ist die obere linke 4x4 Blockmatrix

$\mathbf{Q_{\mathbf{\hat{x}}}} = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0,25 & 0,25 \\ 0 & 0 & 0,25 & 0,25 \\ \end{pmatrix}$

Die untere 3x3 Blockmatrix ist im übrigen die (negative) Kofaktormatrix $\mathbf{Q_{\mathbf{\hat{k}}}}$ der geschätzten Lagrange Multiplikatoren $\mathbf{\hat{k}}$. Den Nebendiagonalblöcke beschreiben wiederum die Abhängigkeiten zwischen $\mathbf{\hat{x}}$ und $\mathbf{\hat{k}}$.

Wie Du siehst, sind die Elemente für HA und HB alle Null. Es würde sich also anbieten, diese erst gar nicht im Modell $\mathbf{y = Ax+e}$ einzuführen. Hierzu müssen wir die funktionalen Zusammenhänge leicht umformulieren und auf gekürzte Beziehungen übergehen. Bisher haben wir folgende Gleichungen direkt verwendet und dabei die Kenntnis der Höhen HA und HB nicht berücksichtigt.

$\Delta hA1 = H1 - HA = 0,5619\,\mathrm{m}$
$\Delta hA2 = H2 - HA = 3,5699\,\mathrm{m}$
$\Delta hB1 = H1 - HB =-2,3152\,\mathrm{m}$
$\Delta hB2 = H2 - HB = 0,6970\,\mathrm{m}$

wenn wir diese aber bereits im funktionalen Modell berücksichtigen wollen, dann ergibt sich:

$\Delta hA1 = H1 - 115,584 = 0,5619\,\mathrm{m}$
$\Delta hA2 = H2 - 115,584 = 3,5699\,\mathrm{m}$
$\Delta hB1 = H1 - 118,460 =-2,3152\,\mathrm{m}$
$\Delta hB2 = H2 - 118,460 = 0,6970\,\mathrm{m}$

und weiter zusammengefasst indem wir die bekannten Höhen auf die andere Seite bringen

$H1 = 116,1459\,\mathrm{m}$
$H2 = 119,1539\,\mathrm{m}$
$H1 = 116,1448\,\mathrm{m}$
$H2 = 119,1570\,\mathrm{m}$

In der Designmatrix haben wir nun nur noch zwei Unbekannte, die Höhen von H1, H2.

$\mathbf{A} = \begin{pmatrix} +1 & 0\\ 0 & +1\\ +1 & 0\\ 0 & +1 \end{pmatrix}$

Diese Matrix entspricht also genau dem rechten Teil der ursprünglichen Designmatrix. Wenn Du nun aber $\mathbf{N=A{^\mathrm{T}}A}$ bildest und den Rang bestimmst, wirst Du feststellen, dass dieser $\mathrm{rg(\mathbf{N})} = 2$ ist und demnach genau der Anzahl der Spalten entspricht. $\mathbf{N^{\mathrm{-1}}}$ existiert und beide Punkte können demnach bestimmt werden, wobei der Beobachtungsvektor nun

$\mathbf{y} = \begin{pmatrix} 116,1459\\ 119,1539\\ 116,1448\\ 119,1570 \end{pmatrix}$

lautet. Die Lösung für dieses Gleichungssystem ist

$\mathbf{\hat{x}} = \begin{pmatrix} 116.14535\\ 119.15545 \end{pmatrix}$

Dies ist also exakt der jeweilige Mittelwert der beiden Höhen in $\mathbf{y}$.
Die Lösung weicht von der bereits gezeigten Lösung jedoch ab, da wir zwar die Information der Festpunkthöhen implizit berücksichtigt haben im funktionalen Modell aber die Restriktion

$\Delta h12 = H2 - H1\overset{!}{=}3,0121\,\mathrm{m}$

noch nicht. Mit

$\mathbf{R} = \begin{pmatrix} -1 &+1 \end{pmatrix}$

und

$\mathbf{r} = \begin{pmatrix} 3,0121 \end{pmatrix}$

erweitern wir wieder das Gleichungssystem zu

$\begin{bmatrix} \mathbf{N} & \mathbf{R^{\mathrm{T}}} \\ \mathbf{R} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{x}\\\mathbf{k} \end{bmatrix} = \begin{bmatrix} \mathbf{y}\\\mathbf{r} \end{bmatrix} $

und erhalten als Lösung wiederum

$\mathbf{\hat{x}} = \begin{pmatrix} 116,14435\\ 119,15645 \end{pmatrix}$

sowie

$\mathbf{\hat{k}} = \begin{pmatrix} -0,002 \end{pmatrix}$

Die Kofaktormatrizen lauten

$\mathbf{Q_{\mathbf{\hat{x}}}} = \begin{pmatrix} 0,25 & 0,25 \\ 0,25 & 0,25 \\ \end{pmatrix}$

und

$\mathbf{Q_{\mathbf{\hat{k}}}} = \begin{pmatrix} 1 \end{pmatrix}$

Da $\mathbf{N^{\mathrm{-1}}}$ existiert, kannst Du hier auch die stufenweise Auflösung nutzen, die Du bereits in der Literatur gefunden hast.

Noch eine Frage: Kannst du mir einen Literaturhinweis geben, indem auch mal Beispielaufgaben enthalten sind oder darf ich dich einfach fragen?

Ich persönlich schaue am häufigsten in das Buch Klassische und robuste Ausgleichungsverfahren von Jäger u.a., welches hoffentlich dieses Jahr in der 2. Auflage erscheint. Es enthält jedoch für die bedingte Ausgleichung auch nicht sonderlich viele Beispiele. Die meisten beziehen sich auf die freie Netzausgleichung. Ich denke aber, dass man sich leicht eigene Beispiele überlegen kann. Gerade in der Formschätzung tritt der Fall häufig auf, dass ein Vektor normiert sein soll oder dass man den Radius eines Kreis/einer Kugel a-priori bereits kennst und dies berücksichtigen will. Sei kreativ. ;-)

Viele Grüße
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Tags:
Ausgleichung, Höhennetz, Parameterschätzung, Gauß-Markov-Modell, freie Netzausgleichung, Restriktionen, Bedingungen, Stufenweise Lösung

Ausgleichung mit Restriktionen

DoreenH @, Sunday, 28.01.2024, 14:43 (vor 90 Tagen) @ MichaeL

Vielen lieben Dank. Das mit den gekürzten Parametern klappt zwar problemlos, aber für mich verständlicher ist deine erste Variante.
Total toll und verständlich sind deine Erläuterungen --> diese helfen mir unwahrscheinlich viel, mich in das Thema rein zu finden.

Viele Grüße
Doreen

Avatar

Ausgleichung mit Restriktionen

MichaeL ⌂, Bad Vilbel, Monday, 29.01.2024, 19:38 (vor 89 Tagen) @ DoreenH

Hi Doreen,

Das mit den gekürzten Parametern klappt zwar problemlos, aber für mich verständlicher ist deine erste Variante.

Du solltest Dir beide Varianten dennoch ansehen. Der Ansatz mit den gekürzten Darstellungen bietet den Vorteil, dass dieser ein kleines Gleichungssystem liefert. Du musst hier im erweiterten Modell nur eine 3x3-Matrix invertieren, während Du im allgemeinen Ansatz eine 7x7 Matrix zu invertieren hast. Man spart sich also Rechenzeit. Zum anderen wirst Du um eine gekürzte Darstellung nicht drum herumkommen, wenn Du nichtlineare Probleme bearbeiten willst.

Alles Gute weiterhin
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Ausgleichung mit Restriktionen

DoreenH @, Tuesday, 20.02.2024, 21:03 (vor 67 Tagen) @ MichaeL

Hallo Micha,
ich bin mit meinen Aufgaben ganz gut vorangekommen. Danke nochmals.
Ich habe jetzt eine Frage zur Biber-Schätzung, hier muss ich irgendein Denkfehler in meiner Berechnung haben.
Ich habe ein Lagenetz, bestehend aus 4 Punkten. (P-fehlerfrei, A-Anschluss, B+C=Neupunkte), Beobachtungen: rPA, rPB, rPC, sPA, sPB, sPC, dxBC, dyBC
gegeben sind ebenfalls Standardabweichungen für r und s, sowie die Kovarianzmatrix Basislinie BC und für A.
Dieses habe ich nach dem Allgemeinfall ausgeglichen und das passt auch.
Im Anschluss folgte die Homogenisierung der Beobachtungen. Nun habe ich l', v', A' und P ist jetzt eine Einheitsmatrix (so habe ich das verstanden) --> erst nach der Berechnung erfolgt die Rücktransformation mit den Gewichten.

Nun beginne ich ja mit der Biber-Schätzung.
Vorbereitung: 1. Berechnung der Redundanz ri = I – A' x N^-1 x A't
2. sigma0 = 1, sigmaVi = sigma0 x Wurzel ri
3. Biber-Wert (gegeben =3,5): ki = 3,5 x sigma vi
4. Berechnung Normierte Verbesserung: |vi| : sigmaVi

Problem: alle meine normierten Verbesserungen sind unter 3,5 (und vi < ki) --> aber eigentlich soll ich mindestens 3 Iterationen rechnen. Also muss ich bei der Berechnung einen Fehler gemacht haben. Ist meine Berechnung von sigma falsch?

Beste Grüße
Doreen

Avatar

Ausgleichung mit Restriktionen

MichaeL ⌂, Bad Vilbel, Wednesday, 21.02.2024, 11:49 (vor 66 Tagen) @ DoreenH

Hallo Doreen,

ich bin derzeit auf einer Messkampagne und zeitlich etwas eingespannt. Ich kann daher nicht ganz so ausführlich sein im Moment.

Ist meine Berechnung von sigma falsch?

Soweit ich das nachvollziehen konnte, sieht es stimmig aus. Die normierte Verbesserung wird mit dem Grenzwert (hier: 3,5) verglichen und dann werden die Gewichte (sofern Du diesen Weg gewählt hast) entsprechend angepasst, siehe Kapitel 12 in der Arbeit von Wicki. Wenn Du Matlab verwendest, könntest Du ggf. Deinen Quellcode zeigen, sodass man die Schritte im Detail nachvollziehen kann mit Deinen Zahlen.
Wenn Du alternativ bezüglich Deiner Implementierung sicher gehen willst, dann kannst Du eines der einfachen Beispiele (ab Seite 144ff) mal nachrechnen, welches Wicki in seiner Dissertation ausführlich abhandelt.

Viele Grüße
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Tags:
Ausgleichungsrechnung, Biber-Schätzer, Robuste Schätzung, Wicki, M-Schätzer

Ausgleichung mit Restriktionen

DoreenH @, Thursday, 22.02.2024, 22:27 (vor 65 Tagen) @ MichaeL

Hallo Micha,

lieben Dank für deine Rückmeldung. Bei meiner Studienarbeit soll ich den Weg der Modifikation der Beobachtungen gehen.
Den ersten Teil konnte ich lösen. Aber nach der 1. Iteration muss ich mein neu berechnetes li noch irgendwie anpassen. Die erste Iteration verläuft nach der Formel im Wicki 12.2.4, aber ab dem 2. Durchlauf muss "die Beobachtung des 1. Iterationsschrittes als auch die jetzige so angepasst werden, bis beide ein vi aufweisen, dass dem entsprechenden ki entspricht."
Hier bin ich noch am rätseln, was ich machen muss. Den ersten neu berechneten li-Wert einfach weiterführen oder immer nach der Formel neu berechnen führt nicht zur Lösung und die Schlussrechnung passt ebenfalls nicht.

Naja mal gucken, ob mir da noch was sinnvolles einfällt.
Beste Grüße
Doreen

Avatar

Ausgleichung mit Restriktionen

MichaeL ⌂, Bad Vilbel, Friday, 23.02.2024, 06:15 (vor 64 Tagen) @ DoreenH

Hi Doreen,

Naja mal gucken, ob mir da noch was sinnvolles einfällt.

Ich habe mal schnell das Beispiel aus der Dissertation von Wicki (Kapitel 18.2.4) in Matlab implementiert. Ich poste hier mal den Quelltext in der Hoffnung, dass es Dir hilft. Ich habe es nicht kommentiert und ggf. sind einige Schritte überflüssig. Ich hoffe dennoch, dass Du es lesen kannst.

Viele Grüße
Micha

clear variables;
close all;
format long g;
clc;
 
set(groot,'defaultAxesTickLabelInterpreter','latex');  
set(groot,'defaulttextinterpreter','latex');
set(groot,'defaultLegendInterpreter','latex');
 
% Ausgleichung ohne Fehler in den Beobachtungen
A = [
    -1  1  0  0
     0 -1  1  0
     0 -1  0  1
     0  0 -1  1
     1  0  0  0
     0  0  0  1
     0  0  1  0
     0  1  0  0
    -1  0  1  0
];
 
P = diag(1./0.001^2.*[0.1276 0.1372 0.0772 0.1041 0.0693 0.1041 0.0918 0.1372 0.0865]);
 
y = [
     32.059
     -6.556
     26.170
     32.726
    -27.809
     30.419
     -2.317
      4.246
     25.496
];
 
N = A'*P*A;
 
Qx = inv(N);
x = Qx * A'*P*y;
 
v =  A*x - y;
sigma_v = sqrt(diag(inv(P) - A*Qx*A'));
 
w = v./sigma_v;
 
% Beobachtungsfehler in 1 und 7
c = 3.5;
y([1 7]) = y([1 7]) + [+0.1; -0.1];
 
x = Qx * A'*P*y;
v = A*x - y;
 
Q_v = inv(P) - A*Qx*A';
R   = Q_v * P; 
sigma_v = sqrt(diag(Q_v));
r = diag(R);
 
w = v./sigma_v;
k = c*sigma_v;
 
while abs(max(abs(w)) - c) > 0.01
    % groesste Teststatistik
    [~, idx] = max(abs(w));
 
    d = zeros(size(v));
    d(idx) = v(idx) - sign(v(idx)).*k(idx);
    y = y - d./-r;
 
    x = Qx * A'*P*y;
    v = A*x - y;
 
    w = v./sigma_v;
end
 
disp(x);
disp(w);

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Ausgleichung mit Restriktionen

DoreenH @, Sunday, 28.01.2024, 12:27 (vor 90 Tagen) @ MichaeL

Hallo Micha,

ich habe mal die A-Matrix gekürzt und nur H1 und H2 sowie nur die Bedingung Δh12 in die Berechnung einbezogen, dann erhalte ich eine reguläre N-Matrix. Dennoch wäre es interessant zu wissen, wie die Genauigkeiten der ersten Berechnung erfolgen würde.
Lieben Dank
Gruß Doreen

RSS-Feed dieser Diskussion