Die Ausgleichung der erhobenen Daten in JAG3D erfolgt nach der Methode der kleinsten Verbesserungsquadratsumme im Gauß-Markov-Modell (GMM). Dieses Modell liegt immer dann vor, wenn die Beobachtungen $\mathbf{l}$ als Funktion der unbekannten Parameter $\mathbf{x}$ dargestellt werden können. Das Gauß-Markov-Modell gilt als Standardmodell in der geodätischen Netzauswertung, da sich alle Beobachtungsgleichungen nach diesem Schema formulieren lassen. Die Zielfunktion lautet
$$\Omega = \mathbf{v^TPv} = Min$$
worin $\mathbf{v}$ die Verbesserungen und $\mathbf{P}$ die Gewichtsmatrix sind. Eine detaillierte Darstellungen des Gauß-Markov-Modell findet sich u.a. bei Jäger et al. (2005).
Im Allgemeinen werden unterschiedliche Beobachtungstypen wie bspw. Strecken und Richtungen gemeinsam in die Ausgleichung als Beobachtungen eingeführt. Diese Beobachtungstypen unterscheiden sich i.d.R. in ihren Genauigkeiten. So sind bei der klassischen Tachymetrie die Winkelmessungen gegenüber der Streckenmessung häufig genauer. Aber auch Beobachtungen gleichen Typs können unterschiedlich genau ausfallen, wenn diese bspw. mit unterschiedlichen Instrumenten oder Verfahren registriert wurden. Das stochastisches Modell hat somit u.a. das Ziel, eine Balance zwischen den unterschiedlichen Beobachtungsgenauigkeiten herzustellen und somit den Einfluss von genaueren Messungen auf das Ausgleichungsergebnis zu erhöhen. Als Maßzahl für die Güte einer Beobachtung wird deren a-priori Standardunsicherheit $\sigma$ verwendet und in der sogenannten Kovarianzmatrix $\mathbf{C_{ll}}$ zusammengefasst. Wird von unabhängigen Beobachtungen ausgegangen, so ist in dieser Matrix lediglich die Hauptdiagonale besetzt.
$$ \mathbf{C_{ll}} = \begin{pmatrix} \sigma_1^2 & 0 & \dots & 0 \\ 0 & \sigma_i^2 & \dots & 0 \\ \vdots & 0 & \ddots & \vdots \\ 0 & \dots & 0 & \sigma_n^2 \end{pmatrix} $$
Durch Invertieren der Kovarianzmatrix gelangt man zur o.g. Gewichtsmatrix $\mathbf{P}$. Der vorgezogene Varianzfaktor $\sigma_0^2$ ist ein beliebig wählbarer Faktor und kann dazu verwendet werden, den Rechenaufwand zu minimieren. In JAG3D wird der Standardwert $\sigma_0 = 1$ verwendet.
$$ \mathbf{C_{ll}} = \sigma_0^2 \mathbf{Q_{ll}} = \sigma_0^2 \mathbf{P^{-1}} $$
Die Matrix $\mathbf{Q_{ll}}$ entspricht hierbei der Kofaktormatrix. Sie unterscheidet sich von der a-priori Varianz-Kovarianz-Matrix $\mathbf{C_{ll}}$ lediglich um den Vorfaktor $\sigma_0^2$ und ist mit $\mathbf{C_{ll}}$ identisch, wenn $\sigma_0 = 1$ gilt.
Neben dem stochastischen Modell existiert das funktionale oder mathematische Modell der Ausgleichung. In diesem wird der lineare bzw. linearisierte Zusammenhang zwischen den (gekürzten) Beobachtungen $\mathbf{l}$ und den zu schätzenden (Zuschlägen der) Unbekannten $\mathbf{\hat x}$ abgebildet
$$\mathbf{l} + \mathbf{v} = \mathbf{\hat{l}} = \mathbf{A\hat{x}}$$
worin $\mathbf{A}$ die Design- oder Jacobimatrix und $\mathbf{v}$ die Beobachtungsresiduen sind. Das hier verwendete funktionale Modell setzt voraus, dass jede Beobachtung $l_i$ als eigenständige Funktion der unbekannten Parameter $l_i = \mathrm{f(\mathbf{\hat x})}$ beschrieben werden kann. In der geodätischen Literatur wird dieses Modell als Gauß-Markov Modell bezeichnet.
Die Normalgleichung ergibt sich zu:
$$\mathbf{A^TPA\hat{x}} = \mathbf{A^TPl}$$
bzw.
$$\mathbf{N\hat{x}} = \mathbf{n}$$
mit der Normalgleichungsmatrix
$$\mathbf{N} = \mathbf{A^TPA}$$
und dem Absolutgliedvektor
$$\mathbf{n} = \mathbf{A^TPl}$$
Durch Auflösen des Normalgleichungssystems nach $\mathbf{\hat x}$ ergibt sich der Lösungsvektor der unbekannten Parameter
$$\mathbf{\hat{x}} = \mathbf{N^{-1}n}$$
welcher aus der lineare Funktion
$$\mathbf{\hat{x}} = \mathbf{Fl}$$
resultiert, worin $\mathbf{F} = \left(\mathbf{A^TPA}\right)^{-1} \mathbf{A^TP}$ die Funktionalmatrix der linearen Transformation ist.
Die geschätzten Beobachtungsverbesserungen lauten
$$\mathbf{v} = \mathbf{A\hat{x}-l} = \left(\mathbf{AF-I}\right)\mathbf{l}$$
und lassen sich ebenfalls als lineare Funktion von $\mathbf{l}$ darstellen.
Im linearisierten Modell enthält der geschätzte Lösungsvektor der unbekannten Parameter lediglich die Zuschläge zur Entwicklungsstelle $\mathbf{x}_0$.
Die Anwendung des Unsicherheitsfortpflanzungsgesetzes liefert neben der Kofaktormatrix der ausgeglichenen Parameter
$$\mathbf{Q_{\hat{x}\hat{x}}} = \mathbf{FQ_{ll}F^T} = \mathbf{N^{-1}}$$
auch die Kofaktormatrix der ausgeglichenen Beobachtungen
$$\mathbf{Q_{\hat{l}\hat{l}}} = \mathbf{AQ_{\hat{x}\hat{x}}A^T}$$
und die Kofaktormatrix der Verbesserungen
$$\mathbf{Q_{vv}} = \mathbf{Q_{ll} - Q_{\hat{l}\hat{l}}}$$
Die geschätzten Kovarianzmatrizen für die Unbekannten, die ausgeglichenen Beobachtungen und die Verbesserungen ergeben sich aus der Multiplikation der entsprechenden Kofaktormatrizen mit dem Varianzfaktor $\sigma^2_0$ der Gewichtseinheit. Die Schätzung für $\sigma^2_0$ wird häufig als a-posteriori Varianzfaktor bezeichnet und ergibt sich aus
$$\hat{\sigma}^2_0 = \frac{\mathbf{v^TPv}} {n-u+d} = \frac{\Omega} {\tr \mathbf{R}}$$
worin $n$ der Anzahl der Beobachtungen, $u$ der Anzahl der Unbekannten und $d$ einen möglichen Datumsdefekt beschreiben. Die Matrix $\mathbf{R}$ wird Redundanzmatrix genannt.
$$\mathbf{R} = \mathbf{Q_{vv}P}$$
Der Gesamtfreiheitsgrad $f$ der Netzausgleichung kann auch aus der Summe der Hauptdiagonalelemente - der sogenannten Spur (engl. Trace) einer Matrix - der Redundanzmatrix $\mathbf{R}$ abgeleitet werden, sofern die Beobachtungen unabhängig voneinander sind. Es gilt in diesem Fall $f = n - u + d = \tr \mathbf{R}$.
Mit Ausnahme von beobachteten Höhenunterschieden sind die Beobachtungsgleichungen für terrestrische Messungen und GNSS-Beobachtungen in der Netzausgleichung nicht-linear. Da der Ausgleichungsalgorithmus – formuliert als Gauß-Markov-Modell –, jedoch ein lineares Modell voraussetzt, sind die nicht-linearen Beobachtungsgleichungen zu linearisieren. Hierzu werden die Beobachtungsgleichungen $\phi$ am Entwicklungspunkt $\mathbf{x}_0$ unter Vernachlässigung von Termen höherer Ordnung mittels Taylorreihe entwickelt,
$$\phi(\mathbf{x})=\phi(\mathbf{x_0}+\mathbf{\hat x})=\phi(\mathbf{x_0})+\mathbf{A\hat x} + \frac{1}2\mathbf{H}[\mathbf{\hat x}\otimes\mathbf{\hat x}] + ...$$
worin $\mathbf{A}$ die Jacobi-Matrix mit den ersten partiellen Ableitungen und $\mathbf{H}$ die Hesse-Matrix mit den zweiten partiellen Ableitungen an der Entwicklungsstelle $\mathbf{x}_0$ sind und $\otimes$ den Kroneckerprodukt-Operator symbolisiert.
Diese Vorgehensweise ist zulässig, wenn die Näherungslösung der gesuchten Parameter $\mathbf{x}_0$ bereits sehr dicht an der tatsächlichen Lösung von $\mathbf{x}$ liegt und somit die Terme höherer Ordnung gegenüber dem linearen und dem konstanten Term vernachlässigt werden können.
$$\phi(\mathbf{x})=\phi(\mathbf{x_0}+\mathbf{\hat x})\approx\phi(\mathbf{x_0})+\mathbf{A\hat x}$$
Die Lösung dieser Gleichung liefert – Konvergenz vorausgesetzt – eine verbesserte Schätzung der Parameter $\mathbf{x}$ gegenüber den zuvor gewählten Näherungswerten $\mathbf{x}_0$.
$$\mathbf{x}=\mathbf{x_0}+\mathbf{\hat x}$$
Zur Bestimmung optimaler Parameter $\mathbf{x}$ lässt sich somit eine einfache Iterationsvorschrift formulieren: Die Taylorreihenentwicklung der Beobachtungsgleichung $\phi$ wird solange an der jeweils verbesserten Entwicklungsstelle $\mathbf{x}_0 := \mathbf{x}$ durchgeführt, bis der Zuschlagsvektor ein Nullvektor ist und somit
$$\mathbf{x}=\mathbf{x_0}$$
gilt.
Liegt die erste Näherung bereits sehr dicht am Optimum, so sind nur wenige Iterationsschritte nötig, um eine Lösung zu erzielen. Bei ungünstig gewählten Näherungswerten sind entsprechend mehr Iterationen notwendig. Im ungünstigsten Fall divergiert das Gleichungssystem sogar und liefert eine falsche bzw. unbrauchbare Lösung. Die Güte der Näherungskoordinaten ist somit entscheidend für das Konvergenzverhalten.